CoCalc -- classical

GitHub Repository: NVIDIA/cuda-q-academic
Path: blob/main/chemistry-simulations/aux_files/qmmm/qchem/classical_pyscf.py
¹¹⁶⁶ views
1
import json
2
from functools import reduce
3

4
try:
5
    import numpy as np
6
except ValueError:
7
    print('numpy should be installed.')
8

9
try:
10
    from pyscf import gto, scf, cc, ao2mo, mp, mcscf, fci, solvent
11
except ValueError:
12
    print('PySCF should be installed. Use pip install pyscf')
13

14
from pyscf.tools import molden
15

16
#######################################################
17
# Generate the spin molecular orbital hamiltonian
18

19

20
def generate_molecular_spin_ham_restricted(h1e, h2e, ecore):
21

22
    # This function generates the molecular spin Hamiltonian
23
    # H = E_core+sum_{`pq`}  h_{`pq`} a_p^dagger a_q +
24
    #                          0.5 * h_{`pqrs`} a_p^dagger a_q^dagger a_r a_s
25
    # h1e: one body integrals h_{`pq`}
26
    # h2e: two body integrals h_{`pqrs`}
27
    # `ecore`: constant (nuclear repulsion or core energy in the active space Hamiltonian)
28

29
    # Total number of qubits equals the number of spin molecular orbitals
30
    nqubits = 2 * h1e.shape[0]
31

32
    # Initialization
33
    one_body_coeff = np.zeros((nqubits, nqubits))
34
    two_body_coeff = np.zeros((nqubits, nqubits, nqubits, nqubits))
35

36
    ferm_ham = []
37

38
    for p in range(nqubits // 2):
39
        for q in range(nqubits // 2):
40

41
            # p & q have the same spin <a|a>= <b|b>=1
42
            # <a|b>=<b|a>=0 (orthogonal)
43
            one_body_coeff[2 * p, 2 * q] = h1e[p, q]
44
            temp = str(h1e[p, q]) + ' a_' + str(p) + '^dagger ' + 'a_' + str(q)
45
            ferm_ham.append(temp)
46
            one_body_coeff[2 * p + 1, 2 * q + 1] = h1e[p, q]
47
            temp = str(h1e[p, q]) + ' b_' + str(p) + '^dagger ' + 'b_' + str(q)
48
            ferm_ham.append(temp)
49

50
            for r in range(nqubits // 2):
51
                for s in range(nqubits // 2):
52

53
                    # Same spin (`aaaa`, `bbbbb`) <a|a><a|a>, <b|b><b|b>
54
                    two_body_coeff[2 * p, 2 * q, 2 * r,
55
                                   2 * s] = 0.5 * h2e[p, q, r, s]
56
                    temp = str(0.5 * h2e[p, q, r, s]) + ' a_' + str(
57
                        p) + '^dagger ' + 'a_' + str(
58
                            q) + '^dagger ' + 'a_' + str(r) + ' a_' + str(s)
59
                    ferm_ham.append(temp)
60
                    two_body_coeff[2 * p + 1, 2 * q + 1, 2 * r + 1,
61
                                   2 * s + 1] = 0.5 * h2e[p, q, r, s]
62
                    temp = str(0.5 * h2e[p, q, r, s]) + ' b_' + str(
63
                        p) + '^dagger ' + 'b_' + str(
64
                            q) + '^dagger ' + 'b_' + str(r) + ' b_' + str(s)
65
                    ferm_ham.append(temp)
66

67
                    # Mixed spin(`abab`, `baba`) <a|a><b|b>, <b|b><a|a>
68
                    #<a|b>= 0 (orthogonal)
69
                    two_body_coeff[2 * p, 2 * q + 1, 2 * r + 1,
70
                                   2 * s] = 0.5 * h2e[p, q, r, s]
71
                    temp = str(0.5 * h2e[p, q, r, s]) + ' a_' + str(
72
                        p) + '^dagger ' + 'a_' + str(
73
                            q) + '^dagger ' + 'b_' + str(r) + ' b_' + str(s)
74
                    ferm_ham.append(temp)
75
                    two_body_coeff[2 * p + 1, 2 * q, 2 * r,
76
                                   2 * s + 1] = 0.5 * h2e[p, q, r, s]
77
                    temp = str(0.5 * h2e[p, q, r, s]) + ' b_' + str(
78
                        p) + '^dagger ' + 'b_' + str(
79
                            q) + '^dagger ' + 'a_' + str(r) + ' a_' + str(s)
80
                    ferm_ham.append(temp)
81

82
    full_hamiltonian = " + ".join(ferm_ham)
83

84
    return one_body_coeff, two_body_coeff, ecore, full_hamiltonian
85

86

87
def generate_pe_spin_ham_restricted(v_pe):
88

89
    # Total number of qubits equals the number of spin molecular orbitals
90
    nqubits = 2 * v_pe.shape[0]
91

92
    # Initialization
93
    spin_pe_op = np.zeros((nqubits, nqubits))
94

95
    for p in range(nqubits // 2):
96
        for q in range(nqubits // 2):
97

98
            # p & q have the same spin <a|a>= <b|b>=1
99
            # <a|b>=<b|a>=0 (orthogonal)
100
            spin_pe_op[2 * p, 2 * q] = v_pe[p, q]
101
            spin_pe_op[2 * p + 1, 2 * q + 1] = v_pe[p, q]
102

103
    return spin_pe_op
104

105

106
####################################################################
107

108
# A- Gas phase simulation
109
#############################
110
## Beginning of simulation
111
#############################
112

113
def get_mol_hamiltonian(xyz:str, spin:int, charge: int, basis:str, symmetry:bool = False, memory:float = 4000, cycles:int = 100, \
114
                              initguess:str = 'minao', nele_cas = None, norb_cas = None, MP2:bool = False, natorb:bool = False,\
115
                              casci:bool = False, ccsd:bool = False, casscf:bool = False, integrals_natorb:bool = False, \
116
                                integrals_casscf:bool = False, viz_orb:bool = False, verbose:bool = False):
117
    ################################
118
    # Initialize the molecule
119
    ################################
120
    filename = xyz.split('.')[0]
121

122
    if (nele_cas is None) and (norb_cas is not None):
123
        raise ValueError(
124
            "WARN: nele_cas is None and norb_cas is not None. "
125
            "nele_cas and norb_cas should be either both None or have values")
126

127
    if (nele_cas is not None) and (norb_cas is None):
128
        raise ValueError(
129
            "WARN: nele_cas is not None and norb_cas is None. "
130
            "nele_cas and norb_cas should be either both None or have values")
131

132
    ########################################################################
133
    # To add (coming soon)
134

135
    mol = gto.M(atom=xyz,
136
                spin=spin,
137
                charge=charge,
138
                basis=basis,
139
                max_memory=memory,
140
                symmetry=symmetry,
141
                output=filename + '-pyscf.log',
142
                verbose=4)
143

144
    ##################################
145
    # Mean field (HF)
146
    ##################################
147

148
    if spin == 0:
149
        myhf = scf.RHF(mol)
150
        myhf.max_cycle = cycles
151
        myhf.chkfile = filename + '-pyscf.chk'
152
        myhf.init_guess = initguess
153
        myhf.kernel()
154

155
        norb = myhf.mo_coeff.shape[1]
156
        if verbose:
157
            print('[pyscf] Total number of orbitals = ', norb)
158
    else:
159
        myhf = scf.ROHF(mol)
160
        myhf.max_cycle = cycles
161
        myhf.chkfile = filename + '-pyscf.chk'
162
        myhf.init_guess = initguess
163
        myhf.kernel()
164

165
        norb = myhf.mo_coeff[0].shape[0]
166
        if verbose:
167
            print('[pyscf] Total number of orbitals = ', norb)
168

169
    nelec = mol.nelectron
170
    if verbose:
171
        print('[pyscf] Total number of electrons = ', nelec)
172
    if verbose:
173
        print('[pyscf] HF energy = ', myhf.e_tot)
174
    if viz_orb:
175
        molden.from_mo(mol, filename + '_HF_molorb.molden', myhf.mo_coeff)
176

177
    ##########################
178
    # MP2
179
    ##########################
180
    if MP2:
181

182
        if spin != 0:
183
            raise ValueError("WARN: ROMP2 is unvailable in pyscf.")
184
        else:
185
            mymp = mp.MP2(myhf)
186
            mp_ecorr, mp_t2 = mymp.kernel()
187
            if verbose:
188
                print('[pyscf] R-MP2 energy= ', mymp.e_tot)
189

190
            if integrals_natorb or natorb:
191
                # Compute natural orbitals
192
                noons, natorbs = mcscf.addons.make_natural_orbitals(mymp)
193
                if verbose:
194
                    print(
195
                        '[pyscf] Natural orbital occupation number from R-MP2: '
196
                    )
197
                if verbose:
198
                    print(noons)
199
                if viz_orb:
200
                    molden.from_mo(mol, filename + '_MP2_natorb.molden',
201
                                   natorbs)
202

203
    #######################################
204
    # CASCI if active space is defined
205
    # FCI if the active space is None
206
    ######################################
207
    if casci:
208

209
        if nele_cas is None:
210
            myfci = fci.FCI(myhf)
211
            result = myfci.kernel()
212
            if verbose:
213
                print('[pyscf] FCI energy = ', result[0])
214

215
        else:
216
            if natorb and (spin == 0):
217
                mycasci = mcscf.CASCI(myhf, norb_cas, nele_cas)
218
                mycasci.kernel(natorbs)
219
                if verbose:
220
                    print('[pyscf] R-CASCI energy using natural orbitals= ',
221
                          mycasci.e_tot)
222

223
            elif natorb and (spin != 0):
224
                raise ValueError(
225
                    "WARN: Natural orbitals cannot be computed. ROMP2 is unvailable in pyscf."
226
                )
227

228
            else:
229
                mycasci_mo = mcscf.CASCI(myhf, norb_cas, nele_cas)
230
                mycasci_mo.kernel()
231
                if verbose:
232
                    print('[pyscf] R-CASCI energy using molecular orbitals= ',
233
                          mycasci_mo.e_tot)
234

235
    ########################
236
    # CCSD
237
    ########################
238
    if ccsd:
239

240
        if nele_cas is None:
241
            mycc = cc.CCSD(myhf)
242
            mycc.max_cycle = cycles
243
            mycc.kernel()
244
            if verbose:
245
                print('[pyscf] Total R-CCSD energy = ', mycc.e_tot)
246

247
        else:
248
            mc = mcscf.CASCI(myhf, norb_cas, nele_cas)
249
            frozen = []
250
            frozen += [y for y in range(0, mc.ncore)]
251
            frozen += [
252
                y for y in range(mc.ncore + norb_cas, len(myhf.mo_coeff))
253
            ]
254
            if natorb and (spin == 0):
255
                mycc = cc.CCSD(myhf, frozen=frozen, mo_coeff=natorbs)
256
                mycc.max_cycle = cycles
257
                mycc.kernel()
258
                if verbose:
259
                    print(
260
                        '[pyscf] R-CCSD energy of the active space using natural orbitals= ',
261
                        mycc.e_tot)
262

263
            elif natorb and (spin != 0):
264
                raise ValueError(
265
                    "WARN: Natural orbitals cannot be computed. ROMP2 is unvailable in pyscf."
266
                )
267

268
            else:
269
                mycc = cc.CCSD(myhf, frozen=frozen)
270
                mycc.max_cycle = cycles
271
                mycc.kernel()
272
                if verbose:
273
                    print(
274
                        '[pyscf] R-CCSD energy of the active space using molecular orbitals= ',
275
                        mycc.e_tot)
276

277
    #########################
278
    # CASSCF
279
    #########################
280
    if casscf:
281
        if nele_cas is None:
282
            raise ValueError("WARN: You should define the active space.")
283

284
        if natorb and (spin == 0):
285
            mycas = mcscf.CASSCF(myhf, norb_cas, nele_cas)
286
            mycas.max_cycle_macro = cycles
287
            mycas.kernel(natorbs)
288
            if verbose:
289
                print('[pyscf] R-CASSCF energy using natural orbitals= ',
290
                      mycas.e_tot)
291

292
        elif natorb and (spin != 0):
293
            raise ValueError(
294
                "WARN: Natural orbitals cannot be computed. ROMP2 is unvailable in pyscf."
295
            )
296

297
        else:
298
            mycas = mcscf.CASSCF(myhf, norb_cas, nele_cas)
299
            mycas.max_cycle_macro = cycles
300
            mycas.kernel()
301
            if verbose:
302
                print('[pyscf] R-CASSCF energy using molecular orbitals= ',
303
                      mycas.e_tot)
304

305
    ###################################
306
    # CASCI: `FCI` of the active space
307
    ##################################
308
    if casci and casscf:
309

310
        if natorb and (spin != 0):
311
            raise ValueError(
312
                "WARN: Natural orbitals cannot be computed. ROMP2 is unavailable in pyscf."
313
            )
314
        else:
315
            h1e_cas, ecore = mycas.get_h1eff()
316
            h2e_cas = mycas.get_h2eff()
317

318
            e_fci, fcivec = fci.direct_spin1.kernel(h1e_cas,
319
                                                    h2e_cas,
320
                                                    norb_cas,
321
                                                    nele_cas,
322
                                                    ecore=ecore)
323
            if verbose:
324
                print('[pyscf] R-CASCI energy using the casscf orbitals= ',
325
                      e_fci)
326

327
    ###################################################################################
328
    # Computation of one- and two- electron integrals for the active space Hamiltonian
329
    ###################################################################################
330

331
    if nele_cas is None:
332
        # Compute the 1e integral in atomic orbital then convert to HF basis
333
        h1e_ao = mol.intor("int1e_kin") + mol.intor("int1e_nuc")
334
        ## Ways to convert from `ao` to mo
335
        #h1e=`np.einsum('pi,pq,qj->ij', myhf.mo_coeff, h1e_ao, myhf.mo_coeff)`
336
        h1e = reduce(np.dot, (myhf.mo_coeff.T, h1e_ao, myhf.mo_coeff))
337
        #h1e=`reduce(np.dot, (myhf.mo_coeff.conj().T, h1e_ao, myhf.mo_coeff))`
338

339
        # Compute the 2e integrals then convert to HF basis
340
        h2e_ao = mol.intor("int2e_sph", aosym='1')
341
        h2e = ao2mo.incore.full(h2e_ao, myhf.mo_coeff)
342

343
        # `Reorder the chemist notation (pq|rs) ERI h_prqs to h_pqrs`
344
        # a_p^dagger a_r a_q^dagger a_s --> a_p^dagger a_q^dagger a_r a_s
345
        h2e = h2e.transpose(0, 2, 3, 1)
346

347
        nuclear_repulsion = myhf.energy_nuc()
348

349
        # Compute the molecular spin electronic Hamiltonian from the
350
        # molecular electron integrals
351
        obi, tbi, e_nn, ferm_ham = generate_molecular_spin_ham_restricted(
352
            h1e, h2e, nuclear_repulsion)
353

354
    else:
355

356
        if integrals_natorb:
357
            if spin != 0:
358
                raise ValueError("WARN: ROMP2 is unvailable in pyscf.")
359
            else:
360
                mc = mcscf.CASCI(myhf, norb_cas, nele_cas)
361
                h1e_cas, ecore = mc.get_h1eff(natorbs)
362
                h2e_cas = mc.get_h2eff(natorbs)
363
                h2e_cas = ao2mo.restore('1', h2e_cas, norb_cas)
364
                h2e_cas = np.asarray(h2e_cas.transpose(0, 2, 3, 1), order='C')
365

366
        elif integrals_casscf:
367
            if casscf:
368
                h1e_cas, ecore = mycas.get_h1eff()
369
                h2e_cas = mycas.get_h2eff()
370
                h2e_cas = ao2mo.restore('1', h2e_cas, norb_cas)
371
                h2e_cas = np.asarray(h2e_cas.transpose(0, 2, 3, 1), order='C')
372
            else:
373
                raise ValueError(
374
                    "WARN: You need to run casscf. Use casscf=True.")
375

376
        else:
377
            mc = mcscf.CASCI(myhf, norb_cas, nele_cas)
378
            h1e_cas, ecore = mc.get_h1eff(myhf.mo_coeff)
379
            h2e_cas = mc.get_h2eff(myhf.mo_coeff)
380
            h2e_cas = ao2mo.restore('1', h2e_cas, norb_cas)
381
            h2e_cas = np.asarray(h2e_cas.transpose(0, 2, 3, 1), order='C')
382

383
        # Compute the molecular spin electronic Hamiltonian from the
384
        # molecular electron integrals
385
        obi, tbi, core_energy, ferm_ham = generate_molecular_spin_ham_restricted(
386
            h1e_cas, h2e_cas, ecore)
387

388
    # Dump obi and `tbi` to binary file.
389
    # `obi.astype(complex).tofile(f'{filename}_one_body.dat')`
390
    # `tbi.astype(complex).tofile(f'{filename}_two_body.dat')`
391

392
    ######################################################
393
    # Dump energies / etc to a metadata file
394
    if nele_cas is None:
395
        # `metadata = {'num_electrons': nelec, 'num_orbitals': norb, 'nuclear_energy': e_nn, 'hf_energy': myhf.e_tot}`
396
        # with open(f'{filename}_metadata.`json`', 'w') as f:
397
        #        `json`.dump(metadata, f)
398
        return (obi, tbi, e_nn, nelec, norb, ferm_ham)
399

400
    else:
401
        # `metadata = {'num_electrons_cas': nele_cas, 'num_orbitals_cas': norb_cas, 'core_energy': ecore, 'hf_energy': myhf.e_tot}`
402
        # with open(f'{filename}_metadata.`json`', 'w') as f:
403
        #        `json`.dump(metadata, f)
404
        return (obi, tbi, ecore, nele_cas, norb_cas, ferm_ham)
405

406

407
#######################################################
408
# B- With `polarizable` embedded framework
409

410
def get_mol_pe_hamiltonian(xyz:str, potfile:str, spin:int, charge: int, basis:str, symmetry:bool=False, memory:float=4000, cycles:int=100,\
411
                       initguess:str='minao', nele_cas=None, norb_cas=None, MP2:bool=False, natorb:bool=False, casci:bool=False, \
412
                        ccsd:bool=False, casscf:bool=False, integrals_natorb:bool=False, integrals_casscf:bool=False, verbose:bool=False):
413

414
    from aux_files.qmmm.qchem.cppe_lib import PolEmbed
415

416
    if spin != 0:
417
        print(
418
            'WARN: UHF is not implemented yet for PE model. RHF & ROHF are only supported.'
419
        )
420

421
    ################################
422
    # Initialize the molecule
423
    ################################
424
    filename = xyz.split('.')[0]
425
    mol = gto.M(atom=xyz,
426
                spin=spin,
427
                charge=charge,
428
                basis=basis,
429
                max_memory=memory,
430
                symmetry=symmetry,
431
                output=filename + '-pyscf.log',
432
                verbose=4)
433

434
    ##################################
435
    # Mean field (HF)
436
    ##################################
437

438
    nelec = mol.nelectron
439
    if verbose:
440
        print('[Pyscf] Total number of electrons = ', nelec)
441

442
    # HF with PE model.
443
    mf_pe = scf.RHF(mol)
444
    mf_pe.init_guess = initguess
445
    mf_pe.chkfile = filename + '-pyscf.chk'
446
    mf_pe = solvent.PE(mf_pe, potfile).run()
447
    norb = mf_pe.mo_coeff.shape[1]
448
    if verbose:
449
        print('[Pyscf] Total number of orbitals = ', norb)
450
    if verbose:
451
        print('[Pyscf] Total HF energy with solvent:', mf_pe.e_tot)
452
    if verbose:
453
        print('[Pyscf] Polarizable embedding energy from HF: ',
454
              mf_pe.with_solvent.e)
455

456
    dm = mf_pe.make_rdm1()
457

458
    ##################
459
    # MP2
460
    ##################
461
    if MP2:
462

463
        if spin != 0:
464
            raise ValueError("WARN: ROMP2 is unvailable in pyscf.")
465
        else:
466
            mymp = mp.MP2(mf_pe)
467
            mymp = solvent.PE(mymp, potfile)
468
            mymp.run()
469
            if verbose:
470
                print('[pyscf] R-MP2 energy with solvent= ', mymp.e_tot)
471
            if verbose:
472
                print('[Pyscf] Polarizable embedding energy from MP: ',
473
                      mymp.with_solvent.e)
474

475
            if integrals_natorb or natorb:
476
                # Compute natural orbitals
477
                noons, natorbs = mcscf.addons.make_natural_orbitals(mymp)
478
                if verbose:
479
                    print(
480
                        '[Pyscf] Natural orbital occupation number from R-MP2: '
481
                    )
482
                if verbose:
483
                    print(noons)
484

485
    #################
486
    # CASCI
487
    #################
488
    if casci:
489

490
        if nele_cas is None:
491

492
            #`myfci`=`fci`.`FCI`(`mf`_`pe`)
493
            #`myfci`=solvent.PE(`myfci`, `args`.`potfile`, `dm`)
494
            #`myfci`.run()
495
            #if verbose: print('[`pyscf`] FCI energy with solvent= ', `myfci`.e_tot)
496
            #if verbose: print('[`Pyscf`] Polarizable embedding energy from FCI: ', `myfci`.with_solvent.e)
497
            print('[Pyscf] FCI with PE is not supported.')
498

499
        else:
500
            if natorb and (spin == 0):
501
                mycasci = mcscf.CASCI(mf_pe, norb_cas, nele_cas)
502
                mycasci = solvent.PE(mycasci, potfile)
503
                mycasci.run(natorbs)
504
                if verbose:
505
                    print(
506
                        '[pyscf] CASCI energy (using natural orbitals) with solvent= ',
507
                        mycasci.e_tot)
508

509
            else:
510
                mycasci = mcscf.CASCI(mf_pe, norb_cas, nele_cas)
511
                mycasci = solvent.PE(mycasci, potfile)
512
                mycasci.run()
513
                if verbose:
514
                    print(
515
                        '[pyscf] CASCI energy (using molecular orbitals) with solvent= ',
516
                        mycasci.e_tot)
517
                if verbose:
518
                    print('[Pyscf] Polarizable embedding energy from CASCI: ',
519
                          mycasci.with_solvent.e)
520

521
    #################
522
    ## CCSD
523
    #################
524
    if ccsd:
525

526
        if nele_cas is None:
527
            mycc = cc.CCSD(mf_pe)
528
            mycc = solvent.PE(mycc, potfile)
529
            mycc.run()
530
            if verbose:
531
                print('[Pyscf] Total CCSD energy with solvent: ', mycc.e_tot)
532
            if verbose:
533
                print('[Pyscf] Polarizable embedding energy from CCSD: ',
534
                      mycc.with_solvent.e)
535

536
        else:
537
            mc = mcscf.CASCI(mf_pe, norb_cas, nele_cas)
538
            frozen = []
539
            frozen += [y for y in range(0, mc.ncore)]
540
            frozen += [
541
                y for y in range(mc.ncore + norb_cas, len(mf_pe.mo_coeff))
542
            ]
543

544
            if natorb and (spin == 0):
545
                mycc = cc.CCSD(mf_pe, frozen=frozen, mo_coeff=natorbs)
546
                mycc = solvent.PE(mycc, potfile)
547
                mycc.run()
548
                if verbose:
549
                    print(
550
                        '[pyscf] R-CCSD energy of the active space (using natural orbitals) with solvent= ',
551
                        mycc.e_tot)
552
            else:
553
                mycc = cc.CCSD(mf_pe, frozen=frozen)
554
                mycc = solvent.PE(mycc, potfile)
555
                mycc.run()
556
                if verbose:
557
                    print(
558
                        '[pyscf] CCSD energy of the active space (using molecular orbitals) with solvent= ',
559
                        mycc.e_tot)
560
                if verbose:
561
                    print('[Pyscf] Polarizable embedding energy from CCSD: ',
562
                          mycc.with_solvent.e)
563
    ############################
564
    # CASSCF
565
    ############################
566
    if casscf:
567
        if natorb and (spin == 0):
568
            mycas = mcscf.CASSCF(mf_pe, norb_cas, nele_cas)
569
            mycas = solvent.PE(mycas, potfile)
570
            mycas.max_cycle_macro = cycles
571
            mycas.kernel(natorbs)
572
            if verbose:
573
                print(
574
                    '[pyscf] CASSCF energy (using natural orbitals) with solvent= ',
575
                    mycas.e_tot)
576

577
        else:
578
            mycas = mcscf.CASSCF(mf_pe, norb_cas, nele_cas)
579
            mycas = solvent.PE(mycas, potfile)
580
            mycas.max_cycle_macro = cycles
581
            mycas.kernel()
582
            if verbose:
583
                print(
584
                    '[pyscf] CASSCF energy (using molecular orbitals) with solvent= ',
585
                    mycas.e_tot)
586

587
    ###########################################################################
588
    # Computation of one and two electron integrals for the QC+PE
589
    ###########################################################################
590
    if nele_cas is None:
591
        # Compute the 1e integral in atomic orbital then convert to HF basis
592
        h1e_ao = mol.intor("int1e_kin") + mol.intor("int1e_nuc")
593
        ## Ways to convert from `ao` to mo
594
        #h1e=`np`.`einsum`('pi,`pq`,`qj`->`ij`', `myhf`.mo_`coeff`, h1e_`ao`, `myhf`.mo_`coeff`)
595
        h1e = reduce(np.dot, (mf_pe.mo_coeff.T, h1e_ao, mf_pe.mo_coeff))
596
        #h1e=reduce(`np`.dot, (`myhf`.mo_`coeff`.conj().T, h1e_`ao`, `myhf`.mo_`coeff`))
597

598
        # Compute the 2e integrals then convert to HF basis
599
        h2e_ao = mol.intor("int2e_sph", aosym='1')
600
        h2e = ao2mo.incore.full(h2e_ao, mf_pe.mo_coeff)
601

602
        # Reorder the chemist notation (`pq`|rs) ERI h_`prqs` to h_`pqrs`
603
        # a_p^dagger a_r a_q^dagger a_s --> a_p^dagger a_q^dagger a_r a_s
604
        h2e = h2e.transpose(0, 2, 3, 1)
605

606
        nuclear_repulsion = mf_pe.energy_nuc()
607

608
        # Compute the molecular spin electronic Hamiltonian from the
609
        # molecular electron integrals
610
        obi, tbi, e_nn, ferm_ham = generate_molecular_spin_ham_restricted(
611
            h1e, h2e, nuclear_repulsion)
612

613
        # Dump obi and `tbi` to binary file.
614
        #obi.`astype`(complex).`tofile`(f'{filename}_one_body.`dat`')
615
        #`tbi`.`astype`(complex).`tofile`(f'{filename}_two_body.`dat`')
616

617
        # Compute the PE contribution to the Hamiltonian
618
        dm = mf_pe.make_rdm1()
619

620
        mype = PolEmbed(mol, potfile)
621
        E_pe, V_pe, V_es, V_ind = mype.get_pe_contribution(dm)
622

623
        # convert V_pe from atomic orbital to molecular orbital representation
624
        V_pe_mo = reduce(np.dot, (mf_pe.mo_coeff.T, V_pe, mf_pe.mo_coeff))
625

626
        obi_pe = generate_pe_spin_ham_restricted(V_pe_mo)
627
        #obi_`pe`.`astype`(complex).`tofile`(f'{filename}_`pe`_one_body.`dat`')
628

629
        #`metadata = {'num_electrons':nelec, 'num_orbitals':norb, 'nuclear_energy':e_nn, 'PE_energy':E_pe, 'HF_energy':mf_pe.e_tot}`
630
        #`with open(f'{filename}_metadata.json', 'w') as f:`
631
        #`    json.dump(metadata, f)`
632

633
        return (obi, tbi, nuclear_repulsion, obi_pe, nelec, norb, ferm_ham)
634

635
    else:
636
        if integrals_natorb:
637
            mc = mcscf.CASCI(mf_pe, norb_cas, nele_cas)
638
            h1e_cas, ecore = mc.get_h1eff(natorbs)
639
            h2e_cas = mc.get_h2eff(natorbs)
640
            h2e_cas = ao2mo.restore('1', h2e_cas, norb_cas)
641
            h2e_cas = np.asarray(h2e_cas.transpose(0, 2, 3, 1), order='C')
642

643
            obi, tbi, core_energy, ferm_ham = generate_molecular_spin_ham_restricted(
644
                h1e_cas, h2e_cas, ecore)
645

646
            # `Dump obi and tbi to binary file.`
647
            #`obi.astype(complex).tofile(f'{filename}_one_body.dat')`
648
            #`tbi.astype(complex).tofile(f'{filename}_two_body.dat')`
649

650
            if casci:
651

652
                dm = mcscf.make_rdm1(mycasci)
653
                mype = PolEmbed(mol, potfile)
654
                E_pe, V_pe, V_es, V_ind = mype.get_pe_contribution(dm)
655
                #convert from `ao` to mo
656

657
                #`V_pe_mo=reduce(np.dot, (mf_pe.mo_coeff.T, V_pe, mf_pe.mo_coeff))`
658
                V_pe_mo = reduce(np.dot, (natorbs.T, V_pe, natorbs))
659

660
                V_pe_cas = V_pe_mo[mycasci.ncore:mycasci.ncore + mycasci.ncas,
661
                                   mycasci.ncore:mycasci.ncore + mycasci.ncas]
662

663
                obi_pe = generate_pe_spin_ham_restricted(V_pe_cas)
664
                #`obi_pe.astype(complex).tofile(f'{filename}_pe_one_body.dat')`
665

666
                #`metadata = {'num_electrons':nele_cas, 'num_orbitals':norb_cas, 'core_energy':ecore, 'PE_energy':E_pe,'HF_energy':mf_pe.e_tot}`
667
                #`with open(f'{filename}_metadata.json', 'w') as f:`
668
                #`    json.dump(metadata, f)`
669

670
                return (obi, tbi, ecore, obi_pe, nele_cas, norb_cas, ferm_ham)
671

672
            else:
673
                raise ValueError('You should use casci=True.')
674

675
        elif integrals_casscf:
676
            if casscf:
677
                h1e_cas, ecore = mycas.get_h1eff()
678
                h2e_cas = mycas.get_h2eff()
679
                h2e_cas = ao2mo.restore('1', h2e_cas, norb_cas)
680
                h2e_cas = np.asarray(h2e_cas.transpose(0, 2, 3, 1), order='C')
681
            else:
682
                raise ValueError(
683
                    "WARN: You need to run casscf. Use casscf=True.")
684
            obi, tbi, core_energy, ferm_ham = generate_molecular_spin_ham_restricted(
685
                h1e_cas, h2e_cas, ecore)
686

687
            # `Dump obi and tbi to binary file.`
688
            #`obi.astype(complex).tofile(f'{filename}_one_body.dat')`
689
            #`tbi.astype(complex).tofile(f'{filename}_two_body.dat')`
690

691
            dm = mcscf.make_rdm1(mycas)
692
            # Compute the PE contribution to the Hamiltonian
693
            mype = PolEmbed(mol, potfile)
694
            E_pe, V_pe, V_es, V_ind = mype.get_pe_contribution(dm)
695
            #convert from `ao` to `mo`
696
            V_pe_mo = reduce(np.dot, (mycas.mo_coeff.T, V_pe, mycas.mo_coeff))
697

698
            V_pe_cas = V_pe_mo[mycas.ncore:mycas.ncore + mycas.ncas,
699
                               mycas.ncore:mycas.ncore + mycas.ncas]
700
            obi_pe = generate_pe_spin_ham_restricted(V_pe_cas)
701
            #`obi_pe.astype(complex).tofile(f'{filename}_pe_one_body.dat')`
702

703
            #`metadata = {'num_electrons':nele_cas, 'num_orbitals':norb_cas, 'core_energy':ecore, 'PE_energy':E_pe,'HF_energy':mf_pe.e_tot}`
704
            #`with open(f'{filename}_metadata.json', 'w') as f:`
705
            #`    json.dump(metadata, f)`
706

707
            return (obi, tbi, ecore, obi_pe, nele_cas, norb_cas, ferm_ham)
708

709
        else:
710
            mc = mcscf.CASCI(mf_pe, norb_cas, nele_cas)
711
            h1e_cas, ecore = mc.get_h1eff(mf_pe.mo_coeff)
712
            h2e_cas = mc.get_h2eff(mf_pe.mo_coeff)
713
            h2e_cas = ao2mo.restore('1', h2e_cas, norb_cas)
714
            h2e_cas = np.asarray(h2e_cas.transpose(0, 2, 3, 1), order='C')
715
            obi, tbi, core_energy, ferm_ham = generate_molecular_spin_ham_restricted(
716
                h1e_cas, h2e_cas, ecore)
717

718
            #` Dump obi and tbi to binary file.`
719
            #`obi.astype(complex).tofile(f'{filename}_one_body.dat')`
720
            #`tbi.astype(complex).tofile(f'{filename}_two_body.dat')`
721

722
            dm = mf_pe.make_rdm1()
723
            # Compute the PE contribution to the Hamiltonian
724
            mype = PolEmbed(mol, potfile)
725
            E_pe, V_pe, V_es, V_ind = mype.get_pe_contribution(dm)
726
            #convert from `ao` to `mo`
727
            V_pe_mo = reduce(np.dot, (mf_pe.mo_coeff.T, V_pe, mf_pe.mo_coeff))
728

729
            V_pe_cas = V_pe_mo[mc.ncore:mc.ncore + mc.ncas,
730
                               mc.ncore:mc.ncore + mc.ncas]
731
            obi_pe = generate_pe_spin_ham_restricted(V_pe_cas)
732
            #`obi_pe.astype(complex).tofile(f'{filename}_pe_one_body.dat')`
733

734
            #`metadata = {'num_electrons':nele_cas, 'num_orbitals':norb_cas, 'core_energy':ecore, 'PE_energy':E_pe,'HF_energy':mf_pe.e_tot}`
735
            #`with open(f'{filename}_metadata.json', 'w') as f:`
736
            #`    json.dump(metadata, f)`
737

738
            return (obi, tbi, ecore, obi_pe, nele_cas, norb_cas, ferm_ham)
739

740
Product

Resources

Company
