CCL:G: Start post-Hartree Fock calculation from Localized orbitals in
PySCF
- From: "Prasanta Bandyopadhyay"
<mailprasanta88||gmail.com>
- Subject: CCL:G: Start post-Hartree Fock calculation from Localized
orbitals in PySCF
- Date: Fri, 8 Dec 2023 19:05:49 -0500
Sent to CCL by: "Prasanta Bandyopadhyay" [mailprasanta88=gmail.com]
Hello everyone. This is my first message to the community.
It is known that the Localized Molecular Orbitals can do CCSD calculations in
less time. After a canonical HF calculation, I am trying to do Boys Localization
and feed the localized orbitals to start the CCSD calculation in PySCF. I did
not find any such references or examples. So, my input may also be erroneous. If
my inputs are correct, then there may be some problems with CCSD convergence.
I have run a sample code in ipython notebook and sharing the input/output.
Please help.
Python 3.10.11 | packaged by conda-forge | (main, May 10 2023, 18:58:44) [GCC
11.3.0]
Type 'copyright', 'credits' or 'license' for more information
IPython 8.15.0 -- An enhanced Interactive Python. Type '?' for help.
...: atom = [
...: [ 'H' , 0.000000 , 0.000000 , 0.100000],
...: [ 'H1' , 0.000000 , 0.000000 , 0.860000],
...: [ 'H' , 2.000000 , 0.000000 , 0.100000],
...: [ 'H1' , 2.000000 , 0.000000 , 0.860000]],
...: basis = {'H': gto.parse('''
...: H S
...: 0.9000000000D+00 0.1000000000D+01
...: H S
...: 0.3000000000D+00 0.1000000000D+01
...: '''),'H1': gto.parse('''
...: H S
...: 0.9900000000D+00 0.1000000000D+01
...: H S
...: 0.3900000000D+00 0.1000000000D+01
...: ''')
...: },
...: cart=True
...: )
...: mf=scf.RHF(mol)
...: mf.kernel()
...: # Canonical CCSD calculation
...: ccsd_can = cc.CCSD(mf)
...: ccsd_can.kernel()
System: uname_result(system='Linux', node='prasanta',
release='6.2.0-37-generic', version='#38~22.04.1-Ubuntu SMP PREEMPT_DYNAMIC Thu
Nov 2 18:01:13 UTC 2', machine='x86_64') Threads 6
Python 3.10.11 | packaged by conda-forge | (main, May 10 2023, 18:58:44) [GCC
11.3.0]
numpy 1.24.3 scipy 1.10.1
Date: Fri Dec 8 23:45:32 2023
PySCF version 2.3.0
PySCF path /home/pro/psi4conda/lib/python3.10/site-packages/pyscf
[CONFIG] conf_file None
[INPUT] verbose = 4
[INPUT] num. atoms = 4
[INPUT] num. electrons = 4
[INPUT] charge = 0
[INPUT] spin (= nelec alpha-beta = 2S) = 0
[INPUT] symmetry False subgroup None
[INPUT] Mole.unit = angstrom
[INPUT] Cartesian GTO integrals (6d 10f)
[INPUT] Symbol X Y Z unit
X Y Z unit Magmom
[INPUT] 1 H 0.000000000000 0.000000000000 0.100000000000 AA
0.000000000000 0.000000000000 0.188972612457 Bohr 0.0
[INPUT] 2 H1 0.000000000000 0.000000000000 0.860000000000 AA
0.000000000000 0.000000000000 1.625164467126 Bohr 0.0
[INPUT] 3 H 2.000000000000 0.000000000000 0.100000000000 AA
3.779452249130 0.000000000000 0.188972612457 Bohr 0.0
[INPUT] 4 H1 2.000000000000 0.000000000000 0.860000000000 AA
3.779452249130 0.000000000000 1.625164467126 Bohr 0.0
nuclear repulsion = 2.41641498659376
number of shells = 8
number of NR pGTOs = 8
number of NR cGTOs = 8
basis = {'H': [[0, [0.9, 1.0]], [0, [0.3, 1.0]]], 'H1': [[0, [0.99, 1.0]], [0,
[0.39, 1.0]]]}
ecp = {}
CPU time: 0.76
******** <class 'pyscf.scf.hf.RHF'> ********
method = RHF
initial guess = minao
damping factor = 0
level_shift factor = 0
DIIS = <class 'pyscf.scf.diis.CDIIS'>
diis_start_cycle = 1
diis_space = 8
SCF conv_tol = 1e-09
SCF conv_tol_grad = None
SCF max_cycles = 50
direct_scf = True
direct_scf_tol = 1e-13
chkfile to save SCF result = /home/pro/tmpho51b2t8
max_memory 4000 MB (current use 133 MB)
Set gradient conv threshold to 3.16228e-05
Initial guess from minao.
init E= -1.39106822755651
HOMO = -0.433617041428962 LUMO = 0.282308074774963
cycle= 1 E= -2.06295523982167 delta_E= -0.672 |g|= 0.0727 |ddm|= 0.799
HOMO = -0.473556403806346 LUMO = 0.412321896221053
cycle= 2 E= -2.06454743844608 delta_E= -0.00159 |g|= 0.0118 |ddm|= 0.0747
HOMO = -0.469582406635784 LUMO = 0.416563705787143
cycle= 3 E= -2.06459105895304 delta_E= -4.36e-05 |g|= 0.000111 |ddm|= 0.0144
HOMO = -0.469553180529867 LUMO = 0.416580776242043
cycle= 4 E= -2.06459106222831 delta_E= -3.28e-09 |g|= 6.29e-06 |ddm|=
0.000152
HOMO = -0.469554412248582 LUMO = 0.4165811859113
cycle= 5 E= -2.06459106224563 delta_E= -1.73e-11 |g|= 3.17e-07 |ddm|=
1.39e-05
HOMO = -0.469554206106277 LUMO = 0.416581169256107
Extra cycle E= -2.06459106224566 delta_E= -3.73e-14 |g|= 5.96e-08 |ddm|=
3.61e-07
converged SCF energy = -2.06459106224566
******** <class 'pyscf.cc.ccsd.CCSD'> ********
CC2 = 0
CCSD nocc = 2, nmo = 8
max_cycle = 50
direct = 0
conv_tol = 1e-07
conv_tol_normt = 1e-05
diis_space = 6
diis_start_cycle = 0
diis_start_energy_diff = 1e+09
max_memory 4000 MB (current use 141 MB)
Init t2, MP2 energy = -2.09735912778898 E_corr(MP2) -0.0327680655433131
Init E_corr(CCSD) = -0.0327680655433144
cycle = 1 E_corr(CCSD) = -0.0433288723433739 dE = -0.0105608068 norm(t1,t2) =
0.0338082
cycle = 2 E_corr(CCSD) = -0.0469749314685062 dE = -0.00364605913 norm(t1,t2)
= 0.0129309
cycle = 3 E_corr(CCSD) = -0.0491424777892665 dE = -0.00216754632 norm(t1,t2)
= 0.00506568
cycle = 4 E_corr(CCSD) = -0.0491155494004375 dE = 2.69283888e-05 norm(t1,t2)
= 0.000214851
cycle = 5 E_corr(CCSD) = -0.0491166716803206 dE = -1.12227988e-06 norm(t1,t2)
= 2.73751e-05
cycle = 6 E_corr(CCSD) = -0.0491164004888309 dE = 2.7119149e-07 norm(t1,t2) =
7.07441e-06
cycle = 7 E_corr(CCSD) = -0.0491163959830784 dE = 4.50575245e-09 norm(t1,t2)
= 2.03493e-06
CCSD converged
E(CCSD) = -2.113707458228743 E_corr = -0.04911639598307844
Out[1]:
(-0.04911639598307844,
array([[ 9.66395727e-04, -7.15861468e-17, 3.73379105e-03,
2.44065462e-18, -1.57904882e-03, -1.32952145e-16],
[ 2.00856366e-16, -1.07204164e-03, 2.38841448e-16,
3.31332992e-03, 8.65144493e-17, -1.87163682e-03]]),
array([[[[-5.54907257e-02, -6.98304215e-17, 4.96365459e-03,
6.15454933e-17, 1.62904246e-02, 2.35860458e-16],
[-6.98304215e-17, -4.32600566e-02, -1.46096206e-17,
-6.18073805e-03, 2.82455911e-16, -1.38761644e-02],
[ 4.96365459e-03, -1.46096206e-17, -1.76383333e-02,
3.47288702e-17, -7.08184505e-04, -3.53184896e-18],
[ 6.15454933e-17, -6.18073805e-03, 3.47288702e-17,
-1.34470086e-02, 2.46462720e-17, -7.67761587e-04],
[ 1.62904246e-02, 2.82455911e-16, -7.08184505e-04,
2.46462720e-17, -1.08204284e-02, -1.05483455e-17],
[ 2.35860458e-16, -1.38761644e-02, -3.53184896e-18,
-7.67761587e-04, -1.05483455e-17, -9.06019104e-03]],
[[-2.71717382e-16, 5.62101570e-02, 8.10011333e-17,
7.13198951e-03, -2.29155940e-16, 1.58849468e-02],
[ 4.83892195e-02, 2.06293400e-16, -5.44502289e-03,
-2.97797424e-17, -1.81692514e-02, -2.42204538e-16],
[ 5.75513408e-17, -5.73018600e-03, 4.25602864e-17,
-1.63495118e-02, -1.40076323e-18, 5.15636384e-04],
[ 6.22957127e-03, -3.83416194e-17, -1.65269939e-02,
-6.88327791e-18, -1.97040906e-03, -2.58176270e-17],
[-2.00659956e-16, -1.75761121e-02, 4.95767349e-18,
-2.03951127e-03, 3.50784025e-16, -1.06212681e-02],
[ 1.48299958e-02, -2.38310658e-16, 4.72373240e-04,
-2.78193986e-17, -1.06327588e-02, -3.40485550e-16]]],
[[[-2.71717382e-16, 4.83892195e-02, 5.75513408e-17,
6.22957127e-03, -2.00659956e-16, 1.48299958e-02],
[ 5.62101570e-02, 2.06293400e-16, -5.73018600e-03,
-3.83416194e-17, -1.75761121e-02, -2.38310658e-16],
[ 8.10011333e-17, -5.44502289e-03, 4.25602864e-17,
-1.65269939e-02, 4.95767349e-18, 4.72373240e-04],
[ 7.13198951e-03, -2.97797424e-17, -1.63495118e-02,
-6.88327791e-18, -2.03951127e-03, -2.78193986e-17],
[-2.29155940e-16, -1.81692514e-02, -1.40076323e-18,
-1.97040906e-03, 3.50784025e-16, -1.06327588e-02],
[ 1.58849468e-02, -2.42204538e-16, 5.15636384e-04,
-2.58176270e-17, -1.06212681e-02, -3.40485550e-16]],
[[-5.46358905e-02, 1.30287310e-16, 5.13339295e-03,
8.49872253e-17, 1.95471507e-02, 3.66679032e-16],
[ 1.30287310e-16, -5.74430970e-02, -5.56956115e-17,
-7.49341692e-03, 2.84447125e-16, -1.68101818e-02],
[ 5.13339295e-03, -5.56956115e-17, -2.14792869e-02,
-1.04757120e-17, -6.44811716e-04, -1.83577455e-18],
[ 8.49872253e-17, -7.49341692e-03, -1.04757120e-17,
-1.65366923e-02, 1.87452778e-17, -7.38305129e-04],
[ 1.95471507e-02, 2.84447125e-16, -6.44811716e-04,
1.87452778e-17, -1.26396619e-02, -4.00584700e-17],
[ 3.66679032e-16, -1.68101818e-02, -1.83577455e-18,
-7.38305129e-04, -4.00584700e-17, -1.06101681e-02]]]]))
In [2]:
In [2]: localization_method = 'boys' # You can choose 'boys', 'iao', 'pipek', e
...: tc.
...: mo_coeff_localized = lo.orth.orth_ao(mol, localization_method)
...:
...: mo_occ_localized = mf.mo_occ
/home/pro/psi4conda/lib/python3.10/site-packages/pyscf/dft/libxc.py:772:
UserWarning: Since PySCF-2.3, B3LYP (and B3P86) are changed to the VWN-RPA
variant, the same to the B3LYP functional in Gaussian and ORCA (issue 1480). To
restore the VWN5 definition, you can put the setting "B3LYP_WITH_VWN5 =
True" in pyscf_conf.py
warnings.warn('Since PySCF-2.3, B3LYP (and B3P86) are changed to the VWN-RPA
variant, '
In [3]: ccsd_localized = cc.CCSD(mf)
...: ccsd_localized.mo_coeff = mo_coeff_localized
...: mo_occ_localized = mf.mo_occ
...: ccsd_localized.mo_occ = mo_occ_localized
...: ccsd_localized.kernel()
******** <class 'pyscf.cc.ccsd.CCSD'> ********
CC2 = 0
CCSD nocc = 2, nmo = 8
max_cycle = 50
direct = 0
conv_tol = 1e-07
conv_tol_normt = 1e-05
diis_space = 6
diis_start_cycle = 0
diis_start_energy_diff = 1e+09
max_memory 4000 MB (current use 150 MB)
Init t2, MP2 energy = 3.27530080154194 E_corr(MP2) -0.00569411069633339
Init E_corr(CCSD) = 0.216714827164018
cycle = 1 E_corr(CCSD) = -3.04461060576688 dE = -3.26132543 norm(t1,t2) =
36.1895
cycle = 2 E_corr(CCSD) = -19.1686545365852 dE = -16.1240439 norm(t1,t2) =
2491.19
cycle = 3 E_corr(CCSD) = 19.8047923107763 dE = 38.9734468 norm(t1,t2) =
2.68056e+06
cycle = 4 E_corr(CCSD) = 25.879418175191 dE = 6.07462586 norm(t1,t2) =
1.40707e+06
cycle = 5 E_corr(CCSD) = 27.2766254838092 dE = 1.39720731 norm(t1,t2) =
329057
cycle = 6 E_corr(CCSD) = -2.01780064492512 dE = -29.2944261 norm(t1,t2) =
364817
cycle = 7 E_corr(CCSD) = -2.97753280814283 dE = -0.959732163 norm(t1,t2) =
47972.5
cycle = 8 E_corr(CCSD) = -1817.22134315427 dE = -1814.24381 norm(t1,t2) =
109750
cycle = 9 E_corr(CCSD) = 0 dE = 1817.22134 norm(t1,t2) = 6.03432e+10
cycle = 10 E_corr(CCSD) = 0 dE = 0 norm(t1,t2) = 2.2255
cycle = 11 E_corr(CCSD) = 0 dE = 0 norm(t1,t2) = 2.2255
cycle = 12 E_corr(CCSD) = 0 dE = 0 norm(t1,t2) = 2.2255
cycle = 13 E_corr(CCSD) = 0 dE = 0 norm(t1,t2) = 2.2255
cycle = 14 E_corr(CCSD) = 0 dE = 0 norm(t1,t2) = 2.2255
cycle = 15 E_corr(CCSD) = 0.216714827164027 dE = 0.216714827 norm(t1,t2) =
2.2255
cycle = 16 E_corr(CCSD) = 0.108084643542402 dE = -0.108630184 norm(t1,t2) =
36.1895
cycle = 17 E_corr(CCSD) = -0.170149619083507 dE = -0.278234263 norm(t1,t2) =
35.9229
WARN: diis singular, eigh(h) [-3.72002754e-01 -1.08414679e-13 -2.85868422e-14
4.94132623e+00
2.12412279e+01 4.71628680e+01 2.69791993e+03]
---------------------------------------------------------------------------
LinAlgError Traceback (most recent call last)
Cell In[3], line 5
3 mo_occ_localized = mf.mo_occ
4 ccsd_localized.mo_occ = mo_occ_localized
----> 5 ccsd_localized.kernel()
File ~/psi4conda/lib/python3.10/site-packages/pyscf/cc/ccsd.py:1076, in
CCSD.kernel(self, t1, t2, eris)
1075 def kernel(self, t1=None, t2=None, eris=None):
-> 1076 return self.ccsd(t1, t2, eris)
File ~/psi4conda/lib/python3.10/site-packages/pyscf/cc/ccsd.py:1091, in
CCSD.ccsd(self, t1, t2, eris)
1087 if eris is None:
1088 eris = self.ao2mo(self.mo_coeff)
1090 self.converged, self.e_corr, self.t1, self.t2 = \
-> 1091 kernel(self, eris, t1, t2, max_cycle=self.max_cycle,
1092 tol=self.conv_tol, tolnormt=self.conv_tol_normt,
1093 verbose=self.verbose, callback=self.callback)
1094 self._finalize()
1095 return self.e_corr, self.t1, self.t2
File ~/psi4conda/lib/python3.10/site-packages/pyscf/cc/ccsd.py:83, in
kernel(mycc, eris, t1, t2, max_cycle, tol, tolnormt, verbose, callback)
81 t1, t2 = t1new, t2new
82 t1new = t2new = None
---> 83 t1, t2 = mycc.run_diis(t1, t2, istep, normt, eccsd-eold, adiis)
84 eold, eccsd = eccsd, mycc.energy(t1, t2, eris)
85 log.info('cycle = %d E_corr(CCSD) = %.15g dE = %.9g norm(t1,t2) =
%.6g',
86 istep+1, eccsd, eccsd - eold, normt)
File ~/psi4conda/lib/python3.10/site-packages/pyscf/cc/ccsd.py:1243, in
CCSD.run_diis(self, t1, t2, istep, normt, de, adiis)
1239 if (adiis and
1240 istep >= self.diis_start_cycle and
1241 abs(de) < self.diis_start_energy_diff):
1242 vec = self.amplitudes_to_vector(t1, t2)
-> 1243 t1, t2 = self.vector_to_amplitudes(adiis.update(vec))
1244 logger.debug1(self, 'DIIS for step %d', istep)
1245 return t1, t2
File ~/psi4conda/lib/python3.10/site-packages/pyscf/lib/diis.py:237, in
DIIS.update(self, x, xerr)
235 else:
236 self._xprev = None # release memory first
--> 237 self._xprev = xnew = self.extrapolate(nd)
239 self._store('xprev', xnew)
240 if 'xprev' not in self._buffer: # not incore
File ~/psi4conda/lib/python3.10/site-packages/pyscf/lib/diis.py:264, in
DIIS.extrapolate(self, nd)
262 except numpy.linalg.linalg.LinAlgError as e:
263 logger.warn(self, ' diis singular, eigh(h) %s', w)
--> 264 raise e
265 logger.debug1(self, 'diis-c %s', c)
267 xnew = None
File ~/psi4conda/lib/python3.10/site-packages/pyscf/lib/diis.py:261, in
DIIS.extrapolate(self, nd)
259 else:
260 try:
--> 261 c = numpy.linalg.solve(h, g)
262 except numpy.linalg.linalg.LinAlgError as e:
263 logger.warn(self, ' diis singular, eigh(h) %s', w)
File <__array_function__ internals>:200, in solve(*args, **kwargs)
File ~/psi4conda/lib/python3.10/site-packages/numpy/linalg/linalg.py:386, in
solve(a, b)
384 signature = 'DD->D' if isComplexType(t) else 'dd->d'
385 extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
--> 386 r = gufunc(a, b, signature=signature, extobj=extobj)
388 return wrap(r.astype(result_t, copy=False))
File ~/psi4conda/lib/python3.10/site-packages/numpy/linalg/linalg.py:89, in
_raise_linalgerror_singular(err, flag)
88 def _raise_linalgerror_singular(err, flag):
---> 89 raise LinAlgError("Singular matrix")
LinAlgError: Singular matrix