diff --git a/README.md b/README.md index e7eab9b4..c9de93e0 100644 --- a/README.md +++ b/README.md @@ -34,7 +34,7 @@ processors. ##The modified version of [libDMET](https://github.com/gkclab/libdmet_preview) available at [here](https://github.com/oimeitei/libdmet_preview) is -recommended to run periodic BE using QuEmb. +recommended to run periodic BE using QuEmb. &&Wannier90 code is interfaced via [libDMET](https://github.com/gkclab/libdmet_preview) in QuEmb ### Steps @@ -43,16 +43,16 @@ recommended to run periodic BE using QuEmb. ```bash git clone https://github.com/oimeitei/quemb.git cd quemb - + 2. Install QuEmb using one of the following approaches: ```bash - pip install . + pip install . ``` or simply add `path/to/quemd` to `PYTHONPATH` ```bash export PYTHONPATH=/path/to/quemb:$PYTHONPATH ``` - + For conda (or virtual environment) installations, after creating your environment, specify the path to mol-be source as a path file, as in: ```bash echo path/to/quemb > $(python -c "from distutils.sysconfig import get_python_lib; print(get_python_lib())")/quemb.pth @@ -94,6 +94,6 @@ Alternatively, you can view the latest documentation online [here](https://quemb ## References The methods implemented in this code are described in details in the following papers: -- OR Meitei, T Van Voorhis, Periodic bootstrap embedding, [JCTC 19 3123 2023](https://doi.org/10.1021/acs.jctc.3c00069) -- OR Meitei, T Van Voorhis, Electron correlation in 2D periodic systems, [arXiv:2308.06185](https://arxiv.org/abs/2308.06185) +- OR Meitei, T Van Voorhis, Periodic bootstrap embedding, [JCTC 19 3123 2023](https://doi.org/10.1021/acs.jctc.3c00069) +- OR Meitei, T Van Voorhis, Electron correlation in 2D periodic systems, [arXiv:2308.06185](https://arxiv.org/abs/2308.06185) - HZ Ye, HK Tran, T Van Voorhis, Bootstrap embedding for large molecular systems, [JCTC 16 5035 2020](https://doi.org/10.1021/acs.jctc.0c00438) diff --git a/docs/source/index.rst b/docs/source/index.rst index 033826db..3247dabe 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -24,7 +24,7 @@ References 2. OR Meitei, T Van Voorhis, Electron correlation in 2D periodic systems, `arXiv:2308.06185 `_ 3. HZ Ye, HK Tran, T Van Voorhis, Bootstrap embedding for large molecular systems, `JCTC 16 5035 2020 `_ - + .. toctree:: :maxdepth: 1 @@ -37,4 +37,4 @@ References solvers misc - + diff --git a/docs/source/install.rst b/docs/source/install.rst index d2ad8adf..2542ec68 100644 --- a/docs/source/install.rst +++ b/docs/source/install.rst @@ -25,18 +25,18 @@ pip install ----------- :: - + pip install . Add to ``PYTHONPATH`` --------------------- Simply add ``path/to/quemb`` to ``PYTHONPATH`` :: - + export PYTHONPATH=/path/to/quemb:$PYTHONPATH Conda or virtual environment ---------------------------- For conda (or virtual environment) installations, after creating your environment, specify the path to mol-be source as a path file, as in:: - + echo path/to/quemb > $(python -c "from distutils.sysconfig import get_python_lib; print(get_python_lib())")/quemb.pth diff --git a/docs/source/kernel.rst b/docs/source/kernel.rst index 24e1a0c1..ff30b78d 100644 --- a/docs/source/kernel.rst +++ b/docs/source/kernel.rst @@ -31,5 +31,5 @@ Serial BE Solver .. toctree:: :maxdepth: 4 .. autofunction:: molbe.solver.be_func - + diff --git a/docs/source/optimize.rst b/docs/source/optimize.rst index d7fef76c..00efcc1c 100644 --- a/docs/source/optimize.rst +++ b/docs/source/optimize.rst @@ -7,7 +7,7 @@ Main BE optimization .. toctree:: :maxdepth: 4 .. autoclass:: molbe._opt.BEOPT - :members: + :members: Quasi-Newton optimization ========================= diff --git a/docs/source/pfrag.rst b/docs/source/pfrag.rst index 376154cf..941469ea 100644 --- a/docs/source/pfrag.rst +++ b/docs/source/pfrag.rst @@ -17,4 +17,4 @@ Periodic fragments :maxdepth: 4 .. automodule:: kbe.pfrag :members: - + diff --git a/docs/source/solvers.rst b/docs/source/solvers.rst index 69715a91..db50ecc3 100644 --- a/docs/source/solvers.rst +++ b/docs/source/solvers.rst @@ -14,7 +14,7 @@ Crystalline orbital localization -------------------------------- .. autofunction:: kbe.lo.localize - + Density Matching Error ====================== @@ -24,9 +24,9 @@ Interface to Quantum Chemistry Methods ====================================== .. autofunction:: molbe.solver.solve_mp2 - + .. autofunction:: molbe.solver.solve_ccsd - + .. autofunction:: molbe.helper.get_scfObj Schmidt Decomposition @@ -44,15 +44,15 @@ Periodic Schmidt decomposition Handling Hamiltonian ==================== - + .. autofunction:: molbe.helper.get_eri - + .. autofunction:: molbe.helper.get_core Build molecular HF potential ---------------------------- - + .. autofunction:: molbe.helper.get_veff Build perioidic HF potential @@ -63,9 +63,9 @@ Build perioidic HF potential Handling Energies ================= - + .. autofunction:: molbe.helper.get_frag_energy - + .. autofunction:: molbe.rdm.compute_energy_full Handling Densities diff --git a/docs/source/usage.rst b/docs/source/usage.rst index 8cafa85b..4237bee0 100644 --- a/docs/source/usage.rst +++ b/docs/source/usage.rst @@ -15,9 +15,9 @@ Simple example of BE calculation on molecular system:: # Perform pyscf calculations to get mol, mf objects # See quemb/example/molbe_h8_density_matching.py - # get mol: pyscf.gto.M - # get mf: pyscf.scf.RHF - + # get mol: pyscf.gto.M + # get mf: pyscf.scf.RHF + # Define fragments myFrag = fragpart(be_type='be2', mol=mol) @@ -26,7 +26,7 @@ Simple example of BE calculation on molecular system:: # Perform density matching in BE mybe.optimize(solver='CCSD') - + Simple example of periodic BE calculation on 1D periodic system:: @@ -47,4 +47,4 @@ Simple example of periodic BE calculation on 1D periodic system:: # Perform density matching in BE mybe.optimize(solver='CCSD') - + diff --git a/example/data/hexene.xyz b/example/data/hexene.xyz index bbe63b0e..ddbf3d7a 100644 --- a/example/data/hexene.xyz +++ b/example/data/hexene.xyz @@ -1,20 +1,20 @@ 18 -C -5.6502267899 0.7485927383 -0.0074809907 -C -4.5584842828 1.7726977952 -0.0418619714 -H -5.5515181382 0.0602177800 -0.8733001951 -H -6.6384111226 1.2516490350 -0.0591711493 -H -5.5928112720 0.1649656434 0.9355613930 -C -3.2701647911 1.4104028701 0.0085107804 -C -2.1789947571 2.4456245961 -0.0265301736 -C -0.7941361337 1.7933691827 0.0427863465 -H -2.3064879754 3.1340763205 0.8376047229 -H -2.2652945230 3.0294980525 -0.9691823088 -C 0.3121144607 2.8438276122 0.0072105803 -H -0.6626010212 1.1042272672 -0.8201573062 -H -0.7037327970 1.2086599200 0.9845037688 -H 0.2581952957 3.4296615741 -0.9351274363 -H 1.3021608456 2.3437050700 0.0587054221 -H 0.2170130582 3.5342406229 0.8723087816 -H -4.8264911382 2.8238751191 -0.1088423637 -H -3.0132059318 0.3554469837 0.0754505007 +C -5.6502267899 0.7485927383 -0.0074809907 +C -4.5584842828 1.7726977952 -0.0418619714 +H -5.5515181382 0.0602177800 -0.8733001951 +H -6.6384111226 1.2516490350 -0.0591711493 +H -5.5928112720 0.1649656434 0.9355613930 +C -3.2701647911 1.4104028701 0.0085107804 +C -2.1789947571 2.4456245961 -0.0265301736 +C -0.7941361337 1.7933691827 0.0427863465 +H -2.3064879754 3.1340763205 0.8376047229 +H -2.2652945230 3.0294980525 -0.9691823088 +C 0.3121144607 2.8438276122 0.0072105803 +H -0.6626010212 1.1042272672 -0.8201573062 +H -0.7037327970 1.2086599200 0.9845037688 +H 0.2581952957 3.4296615741 -0.9351274363 +H 1.3021608456 2.3437050700 0.0587054221 +H 0.2170130582 3.5342406229 0.8723087816 +H -4.8264911382 2.8238751191 -0.1088423637 +H -3.0132059318 0.3554469837 0.0754505007 diff --git a/example/kbe_polyacetylene.py b/example/kbe_polyacetylene.py index c77de2c0..c3214e7c 100644 --- a/example/kbe_polyacetylene.py +++ b/example/kbe_polyacetylene.py @@ -20,14 +20,14 @@ cell.a = lat cell.atom=''' -H 1.4285621630072645 0.0 -0.586173422487319 -C 0.3415633681566205 0.0 -0.5879921146011252 -H -1.4285621630072645 0.0 0.586173422487319 -C -0.3415633681566205 0.0 0.5879921146011252 -H 1.4285621630072645 0.0 1.868826577512681 -C 0.3415633681566205 0.0 1.867007885398875 -H -1.4285621630072645 0.0 3.041173422487319 -C -0.3415633681566205 0.0 3.0429921146011254 +H 1.4285621630072645 0.0 -0.586173422487319 +C 0.3415633681566205 0.0 -0.5879921146011252 +H -1.4285621630072645 0.0 0.586173422487319 +C -0.3415633681566205 0.0 0.5879921146011252 +H 1.4285621630072645 0.0 1.868826577512681 +C 0.3415633681566205 0.0 1.867007885398875 +H -1.4285621630072645 0.0 3.041173422487319 +C -0.3415633681566205 0.0 3.0429921146011254 ''' cell.unit='Angstrom' diff --git a/example/molbe_dmrg_block2.py b/example/molbe_dmrg_block2.py index 817763e4..09d6ae02 100644 --- a/example/molbe_dmrg_block2.py +++ b/example/molbe_dmrg_block2.py @@ -1,5 +1,5 @@ # A run through the `QuEmb` interface with `block2` for performing BE-DMRG. -# `block2` is a DMRG and sparse tensor network library developed by the +# `block2` is a DMRG and sparse tensor network library developed by the # Garnet-Chan group at Caltech: https://block2.readthedocs.io/en/latest/index.html import os, numpy, sys @@ -13,7 +13,7 @@ fci_ecorr, ccsd_ecorr, ccsdt_ecorr, bedmrg_ecorr = [], [], [], [] # Specify a scratch directory for fragment DMRG files: -scratch = os.getcwd() +scratch = os.getcwd() for a in seps: # Hartree-Fock serves as the starting point for all BE calculations: @@ -31,42 +31,42 @@ # Exact diagonalization (FCI) will provide the reference: mc = fci.FCI(mf) fci_ecorr.append(mc.kernel()[0] - mf.e_tot) - + # CCSD and CCSD(T) are good additional points of comparison: mycc = cc.CCSD(mf).run() et_correction = mycc.ccsd_t() ccsd_ecorr.append(mycc.e_tot - mf.e_tot) ccsdt_ecorr.append(mycc.e_tot + et_correction - mf.e_tot) - # Define BE1 fragments. Prior to DMRG the localization of MOs - # is usually necessary. While there doesn't appear to be + # Define BE1 fragments. Prior to DMRG the localization of MOs + # is usually necessary. While there doesn't appear to be # any clear advantage to using any one scheme over another, # the Pipek-Mezey scheme continues to be the most popular. With # BE-DMRG, localization takes place prior to fragmentation: fobj = fragpart(be_type='be1', mol=mol) mybe = BE( - mf, - fobj, + mf, + fobj, lo_method='pipek-mezey', # or 'lowdin', 'iao', 'boys' pop_method='lowdin' # or 'meta-lowdin', 'mulliken', 'iao', 'becke' ) - + # Next, run BE-DMRG with default parameters and maxM=100. mybe.oneshot( solver='block2', # or 'DMRG', 'DMRGSCF', 'DMRGCI' scratch=scratch, # Scratch dir for fragment DMRG maxM=100, # Max fragment bond dimension force_cleanup=True, # Remove all fragment DMRG tmpfiles - ) + ) bedmrg_ecorr.append(mybe.ebe_tot - mf.e_tot) - # Setting `force_cleanup=True` will clean the scratch directory after each + # Setting `force_cleanup=True` will clean the scratch directory after each # fragment DMRG calculation finishes. DMRG tempfiles can be quite large, so # be sure to keep an eye on them if `force_cleanup=False` (default). - + # NOTE: This does *not* delete the log files `dmrg.conf` and `dmrg.out`for each frag, # which can still be found in `/scratch/`. - + # Finally, plot the resulting potential energy curves: fig, ax = plt.subplots() @@ -81,7 +81,7 @@ # (See ../quemb/example/figures/BEDMRG_H8_PES20.png for an example.) # For larger fragments, you'll want greater control over the fragment -# DMRG calculations. Using the same setup as above for a single geometry: +# DMRG calculations. Using the same setup as above for a single geometry: mol = gto.M() mol.atom = [['H', (0.,0.,i*1.2)] for i in range(8)] mol.basis = 'sto-3g' @@ -90,15 +90,15 @@ mol.build() fobj = fragpart(be_type='be2', mol=mol) mybe = BE( - mf, - fobj, + mf, + fobj, lo_method='pipek-mezey', pop_method='lowdin' ) # We automatically construct the fragment DMRG schedules based on user keywords. The following # input, for example, yields a 60 sweep schedule which uses the two-dot algorithm from sweeps 0-49, -# and the one-dot algo from 50-60. The fragment orbitals are also reordered according the Fiedler +# and the one-dot algo from 50-60. The fragment orbitals are also reordered according the Fiedler # vector procedure, along with a few other tweaks: mybe.optimize( @@ -114,7 +114,7 @@ block_extra_keyword=['fiedler'], # Specify orbital reordering algorithm force_cleanup=True, # Remove all fragment DMRG tmpfiles only_chem=True, -) +) # Or, alternatively, we can construct a full schedule by hand: schedule={ @@ -127,16 +127,16 @@ # and pass it to the fragment solver through `schedule_kwargs`: mybe.optimize( solver='block2', - scratch=scratch, + scratch=scratch, schedule_kwargs=schedule, block_extra_keyword=['fiedler'], force_cleanup=True, only_chem=True, ) -# To make sure the calculation is proceeding as expected, make sure to check `[scratch]/dmrg.conf` +# To make sure the calculation is proceeding as expected, make sure to check `[scratch]/dmrg.conf` # and `[scratch]/dmrg.out`, which are the fragment DMRG inputs and outputs, respectively, used -# by `block2`. +# by `block2`. #NOTE: Parameters in `schedule_kwargs` will overwrite any other DMRG kwargs. #NOTE: The DMRG schedule kwargs and related syntax follows the standard notation used in block2. \ No newline at end of file diff --git a/example/molbe_h8_chemical_potential.py b/example/molbe_h8_chemical_potential.py index dc3e8818..ab148dcf 100644 --- a/example/molbe_h8_chemical_potential.py +++ b/example/molbe_h8_chemical_potential.py @@ -8,9 +8,9 @@ mol = gto.M(atom=''' H 0. 0. 0. H 0. 0. 1. -H 0. 0. 2. +H 0. 0. 2. H 0. 0. 3. -H 0. 0. 4. +H 0. 0. 4. H 0. 0. 5. H 0. 0. 6. H 0. 0. 7. @@ -32,7 +32,7 @@ # Initialize BE mybe = BE(mf, fobj) # Perform chemical potential optimization -mybe.optimize(solver='FCI', only_chem=True) +mybe.optimize(solver='FCI', only_chem=True) # Compute BE error be_ecorr = mybe.ebe_tot - mybe.ebe_hf @@ -41,8 +41,8 @@ # Define BE2 fragments fobj = fragpart(be_type='be2', mol=mol) -mybe = BE(mf, fobj) -mybe.optimize(solver='FCI', only_chem=True) +mybe = BE(mf, fobj) +mybe.optimize(solver='FCI', only_chem=True) # Compute BE error be_ecorr = mybe.ebe_tot - mybe.ebe_hf @@ -51,8 +51,8 @@ # Define BE3 fragments fobj = fragpart(be_type='be3', mol=mol) -mybe = BE(mf, fobj) -mybe.optimize(solver='FCI', only_chem=True) +mybe = BE(mf, fobj) +mybe.optimize(solver='FCI', only_chem=True) # Compute BE error be_ecorr = mybe.ebe_tot - mybe.ebe_hf diff --git a/example/molbe_h8_density_matching.py b/example/molbe_h8_density_matching.py index 2a79605b..ec529ac2 100644 --- a/example/molbe_h8_density_matching.py +++ b/example/molbe_h8_density_matching.py @@ -8,9 +8,9 @@ mol = gto.M(atom=''' H 0. 0. 0. H 0. 0. 1. -H 0. 0. 2. +H 0. 0. 2. H 0. 0. 3. -H 0. 0. 4. +H 0. 0. 4. H 0. 0. 5. H 0. 0. 6. H 0. 0. 7. @@ -41,7 +41,7 @@ # Define BE2 fragments fobj = fragpart(be_type='be2', mol=mol) -mybe = BE(mf, fobj) +mybe = BE(mf, fobj) mybe.optimize(solver='FCI') # Compute BE error @@ -51,7 +51,7 @@ # Define BE3 fragments fobj = fragpart(be_type='be3', mol=mol) -mybe = BE(mf, fobj) +mybe = BE(mf, fobj) mybe.optimize(solver='FCI') # Compute BE error diff --git a/example/molbe_oneshot_rbe_hcore.py b/example/molbe_oneshot_rbe_hcore.py index 972d995f..3ec7db48 100644 --- a/example/molbe_oneshot_rbe_hcore.py +++ b/example/molbe_oneshot_rbe_hcore.py @@ -77,7 +77,7 @@ def pyscf2lint(mol, hcore_pyscf): nproc = 1, # number of processors to parallize across ompnum = 2, be_type = 'be2', # BE type: this sets the fragment size. - frozen_core = True, # Frozen core + frozen_core = True, # Frozen core unrestricted = False, # specify restricted calculation from_chk = False, # can save the RHF as PySCF checkpoint. # Set to true if running from converged UHF chk diff --git a/example/molbe_oneshot_rbe_qmmm-fromchk.py b/example/molbe_oneshot_rbe_qmmm-fromchk.py index 2ecefa31..32c5d6cd 100644 --- a/example/molbe_oneshot_rbe_qmmm-fromchk.py +++ b/example/molbe_oneshot_rbe_qmmm-fromchk.py @@ -30,14 +30,14 @@ # using checkfile from converged RHF be_energy = be2puffin(structure, # the QM region XYZ geometry 'sto-3g', # the chosen basis set - pts_and_charges = [coords, charges], # the loaded hamiltonian + pts_and_charges = [coords, charges], # the loaded hamiltonian use_df = False, # density fitting charge = 0, # charge of QM region spin = 0, # spin of QM region nproc = 1, # number of processors to parallize across ompnum = 2, be_type = 'be2', # BE type: this sets the fragment size. - frozen_core = False, # Frozen core + frozen_core = False, # Frozen core unrestricted = False, # specify restricted calculation from_chk = True, # can save the RHF as PySCF checkpoint. # Set to true if running from converged UHF chk @@ -48,7 +48,7 @@ """ To not use or generate checkfile: from_chk = False -checkfile = None +checkfile = None To generate checkfile in be2puffin: from_chk = False diff --git a/example/molbe_oneshot_ube_qmmm.py b/example/molbe_oneshot_ube_qmmm.py index 4da7d69b..564d5df7 100644 --- a/example/molbe_oneshot_ube_qmmm.py +++ b/example/molbe_oneshot_ube_qmmm.py @@ -36,7 +36,7 @@ nproc = 1, # number of processors to parallize across ompnum = 2, # number of nodes to parallelize across be_type = 'be2', # BE type: this sets the fragment size. - frozen_core = False, # keep this to False for non-minimal basis: localization and + frozen_core = False, # keep this to False for non-minimal basis: localization and # numerical problems for ruthenium systems in non-minimal basis unrestricted = True, # specify unrestricted calculation from_chk = False, # can save the UHF as PySCF checkpoint. diff --git a/kbe/_opt.py b/kbe/_opt.py index 43452b60..1be41d91 100644 --- a/kbe/_opt.py +++ b/kbe/_opt.py @@ -61,9 +61,9 @@ def optimize(self, solver='MP2',method='QN', if only_chem: J0 = [[0.]] J0 = self.get_be_error_jacobian(jac_solver='HF') - J0 = [[J0[-1,-1]]] + J0 = [[J0[-1,-1]]] else: - J0 = self.get_be_error_jacobian(jac_solver='HF') + J0 = self.get_be_error_jacobian(jac_solver='HF') # Perform the optimization be_.optimize(method, J0=J0) diff --git a/kbe/autofrag.py b/kbe/autofrag.py index 8dafcd43..6b452c47 100644 --- a/kbe/autofrag.py +++ b/kbe/autofrag.py @@ -27,15 +27,15 @@ def add_check_k(min1, flist, sts, ksts, nk_): ksts.append(nk_) def nearestof2coord(coord1, coord2 , bond=2.6*1.88973): - + mind = 50000 lmin = () - for idx, i in enumerate(coord1): + for idx, i in enumerate(coord1): for jdx, j in enumerate(coord2): if idx == jdx: continue dist = numpy.linalg.norm(i - j) - + if dist < mind or dist-mind < 0.1: if dist <= bond: lmin = (idx, jdx) @@ -55,10 +55,10 @@ def nearestof2coord(coord1, coord2 , bond=2.6*1.88973): continue dist = numpy.linalg.norm(i - j) if dist-mind<0.1 and dist <= bond: - + lunit_.append(idx) runit_.append(jdx) - + return(lunit_, runit_) @@ -68,7 +68,7 @@ def sidefunc(cell, Idx, unit1, unit2, main_list, sub_list, coord, if ext_list == []: main_list.extend(unit2[numpy.where(unit1==Idx)[0]]) - sub_list.extend(unit2[numpy.where(unit1==Idx)[0]]) + sub_list.extend(unit2[numpy.where(unit1==Idx)[0]]) else: for sub_i in unit2[numpy.where(unit1==Idx)[0]]: if sub_i in rlist: continue @@ -81,11 +81,11 @@ def sidefunc(cell, Idx, unit1, unit2, main_list, sub_list, coord, sub_list.append(sub_i) else: main_list.append(sub_i) - sub_list.append(sub_i) - - closest = sub_list.copy() + sub_list.append(sub_i) + + closest = sub_list.copy() close_be3 = [] - + if be_type == 'be3' or be_type == 'be4': for lmin1 in unit2[numpy.where(unit1==Idx)[0]]: for jdx, j in enumerate(coord): @@ -97,7 +97,7 @@ def sidefunc(cell, Idx, unit1, unit2, main_list, sub_list, coord, main_list.append(jdx) sub_list.append(jdx) close_be3.append(jdx) - + if be_type == 'be4': for kdx, k in enumerate(coord): if kdx == jdx: @@ -127,25 +127,25 @@ def surround(cell, sidx, unit1, unit2, flist, coord, be_type, ext_list, for tmpi in sublist_: if not tmpi in rlist: for tmp_jdx, tmpj in enumerate(ext_list): - if tmpj == tmpi and klist[tmp_jdx] == NK: + if tmpj == tmpi and klist[tmp_jdx] == NK: break - else: + else: sublist.append(tmpi) - + ext_list.extend(sublist) for kdx in sublist: klist.append(NK) - + def kfrag_func(site_list, numk, nk1, uNs, Ns, nk2=None, debug=False, debug1=False, shift=False, debug2=False): if nk2 is None: nk2 = nk1 frglist = [] - for pq in site_list: + for pq in site_list: if numk > nk1 * 2: - if not uNs == Ns: + if not uNs == Ns: nk_ = numk - (nk1*2) frglist.append(uNs*nk1*2 + (nk_*uNs) + pq) else: @@ -159,15 +159,15 @@ def kfrag_func(site_list, numk, nk1, uNs, Ns, nk2=None, debug=False, else: frglist.append(uNs*numk + pq) elif numk >nk1: - if uNs == Ns: + if uNs == Ns: frglist.append(uNs*numk + pq) else: - nk_ = numk - nk1 - frglist.append(uNs*nk1 + (nk_*uNs)+pq) + nk_ = numk - nk1 + frglist.append(uNs*nk1 + (nk_*uNs)+pq) else: if not Ns==uNs: frglist.append(uNs*numk + pq) - + else: frglist.append(uNs*numk + pq) if debug: print(uNs*numk + pq, end=' ') @@ -183,7 +183,7 @@ def be_reduce(be_type, N=1): elif be_type == 'be4': return 'be3' -def autogen(mol, kpt, frozen_core=True, be_type='be2', +def autogen(mol, kpt, frozen_core=True, be_type='be2', write_geom=False, unitcell=1, gamma_2d=False, gamma_1d=False, long_bond=False, perpend_dist = 4.0, perpend_dist_tol = 1e-3, @@ -195,8 +195,8 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', Partitions a unitcell into overlapping fragments as defined in BE atom-based fragmentations. It automatically detects branched chemical chains and ring systems and partitions accordingly. - For efficiency, it only checks two atoms for connectivity (chemical bond) if they are within 3.5 Angstrom. - This value is hardcoded as normdist. Two atoms are defined as bonded if they are within 1.8 Angstrom (1.2 for Hydrogen atom). + For efficiency, it only checks two atoms for connectivity (chemical bond) if they are within 3.5 Angstrom. + This value is hardcoded as normdist. Two atoms are defined as bonded if they are within 1.8 Angstrom (1.2 for Hydrogen atom). This is also hardcoded as bond & hbond. Neighboring unitcells are used in the fragmentation, exploiting translational symmetry. Parameters @@ -243,7 +243,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', """ from pyscf import lib from .misc import sgeom - + if not float(unitcell).is_integer(): print('Fractional unitcell is not supported!') sys.exit() @@ -251,16 +251,16 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', if not nx and not ny and not nz: nx_ = unitcell if kpt[0] > 1 else 1 ny_ = unitcell if kpt[1] > 1 else 1 - nz_ = unitcell if kpt[2] > 1 else 1 + nz_ = unitcell if kpt[2] > 1 else 1 else: nx_ = unitcell if nx else 1 ny_ = unitcell if ny else 1 nz_ = unitcell if nz else 1 kmesh = [nx_, ny_, nz_] - + cell = sgeom(mol, kmesh=kmesh) cell.build() - + print(flush=True) print('No. of cells used in building fragments : {:>3}'.format(unitcell),flush=True) @@ -273,11 +273,11 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', sys.exit() else: cell = mol.copy() - + ncore, no_core_idx, core_list = get_core(cell) - + coord = cell.atom_coords() - + ang2bohr = 1.88973 normdist = 3.5 * ang2bohr bond = 1.8 * ang2bohr @@ -291,9 +291,9 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', dunit = [] munit = [] tunit = [] - + lnext = [i for i in kpt if i>1] - + if not len(lnext) == 0: nk1 = lnext[0] twoD = False @@ -312,7 +312,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', nk2 = 1 twoD = True lkpt = [2,2,1] - + if frozen_core: ncore__, no_core_idx__, core_list__ = get_core(mol) Nsite = mol.aoslice_by_atom()[-1][3]-ncore__ -1 @@ -320,47 +320,47 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', Nsite = mol.aoslice_by_atom()[-1][3]-1 # kmesh is lkpt # Building neighbouring cells - + if twoD: lnk2 = 2 else: lnk2 = 1 - - lattice_vector = cell.lattice_vectors() + + lattice_vector = cell.lattice_vectors() Ts = lib.cartesian_prod((numpy.arange(lkpt[0]), numpy.arange(lkpt[1]), numpy.arange(lkpt[2]))) - + Ls = numpy.dot(Ts, lattice_vector) - + # 1-2-(1-2)-1-2 - # * * + # * * - lcoord = Ls.reshape(-1,1,3)[lnk2]*-1 + coord + lcoord = Ls.reshape(-1,1,3)[lnk2]*-1 + coord rcoord = Ls.reshape(-1,1,3)[lnk2] + coord - + lunit, runit = nearestof2coord(coord, lcoord, bond=bond) lunit_, runit_ = nearestof2coord(coord, rcoord, bond=bond) - + if not set(lunit) == set(runit_) or not set(runit) == set(lunit_): print('Fragmentation error : wrong connection of unit cells ') sys.exit() - + if sum(i>1 for i in kpt) > 1 or gamma_2d: # only 2D is supported # neighbours up-down or right-left - + # *1-2 # | # (1-2) # | - # *1-2 + # *1-2 lcoord2 = Ls.reshape(-1,1,3)[1]*-1 + coord rcoord2 = Ls.reshape(-1,1,3)[1] + coord dunit, uunit = nearestof2coord(coord, lcoord2, bond=bond) dunit_, uunit_ = nearestof2coord(coord, rcoord2, bond=bond) - + if not set(uunit) == set(dunit_) or not set(dunit) == set(uunit_): print('Fragmentation error : wrong connection of unit cells ') sys.exit() @@ -371,24 +371,24 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', # \ / # (1-2) # / \ - # 1-2* 1-2 + # 1-2* 1-2 lcoord3 = Ls.reshape(-1,1,3)[lnk2+1]*-1 + coord rcoord3 = Ls.reshape(-1,1,3)[lnk2+1] + coord - - + + munit, tunit = nearestof2coord(coord, lcoord3, bond=bond) munit_, tunit_ = nearestof2coord(coord, rcoord3, bond=bond) - + if not set(munit) == set(tunit_) or not set(tunit) == set(munit_): print('Fragmentation error : wrong connection of unit cells ') sys.exit() # kmesh lkpt ends - + ## starts here normlist = [] for i in coord: normlist.append(numpy.linalg.norm(i)) - Frag = [] + Frag = [] pedge = [] cen = [] lsites = [] @@ -397,14 +397,14 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', dsites = [] msites = [] tsites = [] - + klsites = [] krsites = [] kusites = [] kdsites = [] kmsites = [] ktsites = [] - + lunit = numpy.asarray(lunit) runit = numpy.asarray(runit) uunit = numpy.asarray(uunit) @@ -423,8 +423,8 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', for ajdx, aj in enumerate(coord): if aidx==ajdx: continue if ai[inter_layer_axis] == aj[inter_layer_axis]: continue - dist = numpy.linalg.norm(ai-aj) - if dist > bond: + dist = numpy.linalg.norm(ai-aj) + if dist > bond: if inter_dist > dist: inter_dist = dist inter_idx = ajdx @@ -438,14 +438,14 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', inter_dist_.append(dist) inter_idx_.append(ajdx) inter_layer_dict.append([inter_idx_, inter_dist_]) - + # Assumes - minimum atom in a ring is 5 if be_type == 'be4': if twoD: print('*********************') print('USE BE4 WITH CAUTION') print('*********************') - + for idx, i in enumerate(normlist): if cell.atom_pure_symbol(idx) == 'H': continue @@ -453,27 +453,27 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', lsts = [] rsts = [] usts = [] - dsts = [] + dsts = [] msts = [] tsts = [] - + klsts = [] krsts = [] kusts = [] - kdsts = [] + kdsts = [] kmsts = [] ktsts = [] - + tmplist = normlist - i tmplist = list(tmplist) - + clist = [] cout = 0 for jdx,j in enumerate(tmplist): if not idx==jdx and not cell.atom_pure_symbol(jdx) == 'H': - if abs(j)< normdist: + if abs(j)< normdist: clist.append(jdx) - + pedg = [] flist = [] @@ -485,7 +485,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', # For systems like Graphene, BN, SiC, hexagonal 2D sheets, # BE4 can give wrong fragmentations # the following code adds in redundant atoms - + if idx in lunit: closest, close_be3 = sidefunc(cell, idx, lunit, runit, flist, lsts, coord, be_type, bond=bond) @@ -494,7 +494,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', klsts.append(nk1*(nk2-1)+1) else: klsts.append(nk1) - + for lsidx in closest: if lsidx in lunit or lsidx in munit: warn_large_fragment() if lsidx in runit: @@ -509,8 +509,8 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', if lsidx in tunit: surround(cell, lsidx, tunit, munit, flist, coord, be_type, lsts, klsts, 2, bond=bond) - - if be_type == 'be4': + + if be_type == 'be4': for lsidx in close_be3: if lsidx in lunit or lsidx in munit: warn_large_fragment() if lsidx in runit: @@ -519,15 +519,15 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', if lsidx in uunit: surround(cell, lsidx, uunit, dunit, flist, coord, 'be3', lsts, klsts, nk1*(nk2-1)+2, bond=bond) - if lsidx in dunit: + if lsidx in dunit: surround(cell, lsidx, dunit, uunit, flist, coord, 'be3', lsts, klsts, nk1*nk2, bond=bond) - if lsidx in tunit: + if lsidx in tunit: surround(cell, lsidx, tunit, munit, flist, coord, 'be3', lsts, klsts, 2, bond=bond) for i1dx, i1 in enumerate(klsts): if i1 == 1: clist_check.append(lsts[i1dx]) - + if idx in uunit: closest, close_be3 = sidefunc(cell, idx, uunit, dunit, flist, usts, coord, be_type, bond=bond) @@ -549,7 +549,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', kusts, 1, rlist=[idx]+clist_check, bond=bond) if be_type == 'be4': for usidx in close_be3: - if usidx in uunit or usidx in tunit: warn_large_fragment() + if usidx in uunit or usidx in tunit: warn_large_fragment() if usidx in lunit: surround(cell, usidx, lunit, runit, flist, coord, 'be3', usts, kusts, nk1*(nk2-1)+2, bond=bond) @@ -564,7 +564,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', kusts, 1, rlist=[idx], bond=bond) for i1dx, i1 in enumerate(kusts): if i1 == 1: clist_check.append(usts[i1dx]) - + if idx in munit: closest, close_be3 = sidefunc(cell, idx, munit, tunit, flist, msts, coord, be_type, bond=bond) @@ -618,7 +618,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', if rsidx in dunit: surround(cell, rsidx, dunit, uunit, flist, coord, be_type, rsts, krsts, nk1*2, bond=bond) - if be_type == 'be4': + if be_type == 'be4': for rsidx in close_be3: if rsidx in runit or rsidx in tunit: warn_large_fragment() if rsidx in lunit: @@ -640,7 +640,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', coord, be_type, bond=bond) for kdx in dsts: kdsts.append(nk1) - for dsidx in closest: + for dsidx in closest: if dsidx in munit or dsidx in dunit: warn_large_fragment() if dsidx in lunit: surround(cell, dsidx, lunit, runit, flist, coord, be_type, dsts, kdsts, @@ -655,7 +655,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', surround(cell, dsidx, tunit, munit, flist, coord, be_type, dsts, kdsts, nk1+1, bond=bond) if be_type == 'be4': - for dsidx in close_be3: + for dsidx in close_be3: if dsidx in munit or dsidx in dunit: warn_large_fragment() if dsidx in lunit: surround(cell, dsidx, lunit, runit, flist, coord, 'be3', dsts, kdsts, @@ -671,7 +671,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', nk1+1, bond=bond) for i1dx, i1 in enumerate(kdsts): if i1 == 1: clist_check.append(dsts[i1dx]) - + if idx in tunit: closest, close_be3 = sidefunc(cell, idx, tunit, munit, flist, tsts, coord, be_type, bond=bond) @@ -702,13 +702,13 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', nk1+1, bond=bond) for i1dx, i1 in enumerate(ktsts): if i1 == 1: clist_check.append(tsts[i1dx]) - + flist.append(idx) - cen.append(idx) - - for jdx in clist: + cen.append(idx) + + for jdx in clist: dist = numpy.linalg.norm(coord[idx] - coord[jdx]) - + if (dist <= bond) or (interlayer and dist in inter_layer_dict[idx][1] \ and jdx in inter_layer_dict[idx][0] and \ dist < perpend_dist*ang2bohr and not jdx in pedg): @@ -716,7 +716,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', flist.append(jdx) pedg.append(jdx) if dist > bond: continue - if be_type=='be3' or be_type == 'be4': + if be_type=='be3' or be_type == 'be4': if jdx in lunit: lmin1 = runit[numpy.where(lunit==jdx)[0]] if not twoD: @@ -742,16 +742,16 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', if twoD: klsts.append(nk1*(nk2-1)+1) else: - klsts.append(nk1) + klsts.append(nk1) for lsidx in lmin1: if lsidx in lunit or lsidx in munit: warn_large_fragment() if lsidx in uunit: surround(cell, lsidx, uunit, dunit, flist, coord, 'be3', lsts, klsts, nk1*(nk2-1)+2, bond=bond) - if lsidx in dunit: + if lsidx in dunit: surround(cell, lsidx, dunit, uunit, flist, coord, 'be3', lsts, klsts, nk1*nk2, bond=bond) - if lsidx in tunit: + if lsidx in tunit: surround(cell, lsidx, tunit, munit, flist, coord, 'be3', lsts, klsts, 2, bond=bond) if jdx in runit: @@ -765,7 +765,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', else: krsts.append(2) else: - add_check_k(rmin1, flist, rsts, krsts, nk1+1) + add_check_k(rmin1, flist, rsts, krsts, nk1+1) if be_type == 'be4': for kdx, k in enumerate(coord): if kdx == jdx or kdx in rmin1 or cell.atom_pure_symbol(kdx) == 'H': @@ -781,7 +781,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', else: krsts.append(2) for rsidx in rmin1: - if rsidx in runit or rsidx in tunit: warn_large_fragment() + if rsidx in runit or rsidx in tunit: warn_large_fragment() if rsidx in munit: surround(cell, rsidx, munit, tunit, flist, coord, 'be3', rsts, krsts, nk1, bond=bond) @@ -807,8 +807,8 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', flist.append(kdx) usts.append(kdx) kusts.append(2) - for usidx in umin1: - if usidx in uunit or usidx in tunit: warn_large_fragment() + for usidx in umin1: + if usidx in uunit or usidx in tunit: warn_large_fragment() if usidx in lunit: surround(cell, usidx, lunit, runit, flist, coord, 'be3', usts, kusts, nk1*(nk2-1)+2, bond=bond) @@ -817,11 +817,11 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', usts, kusts, nk1+2, bond=bond) if usidx in munit: surround(cell, usidx, munit, tunit, flist, coord, 'be3', - usts, kusts, nk1*(nk2-1)+1, bond=bond) + usts, kusts, nk1*(nk2-1)+1, bond=bond) if jdx in dunit: dmin1 = uunit[numpy.where(dunit==jdx)[0]] add_check_k(dmin1, flist, dsts, kdsts, nk1) - + if be_type == 'be4': for kdx, k in enumerate(coord): if kdx == jdx or kdx in dmin1 or \ @@ -834,7 +834,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', flist.append(kdx) dsts.append(kdx) kdsts.append(nk1) - for dsidx in dmin1: + for dsidx in dmin1: if dsidx in munit or dsidx in dunit: warn_large_fragment() if dsidx in lunit: surround(cell, dsidx, lunit, runit, flist, coord, 'be3', @@ -869,11 +869,11 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', if msidx in uunit: surround(cell, msidx, uunit, dunit, flist, coord, 'be3', msts, kmsts, nk1*(nk2-1)+1, bond=bond) - + if jdx in tunit: tmin1 = munit[numpy.where(tunit==jdx)[0]] add_check_k(tmin1, flist, tsts, ktsts, nk1+2) - + if be_type == 'be4': for kdx, k in enumerate(coord): if kdx == jdx or kdx in tmin1 or \ @@ -895,8 +895,8 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', if tsidx in dunit: surround(cell, tsidx, dunit, uunit, flist, coord, 'be3', tsts, ktsts, nk1+1, bond=bond) - - + + for kdx in clist: if not kdx == jdx: dist = numpy.linalg.norm(coord[jdx] - coord[kdx]) @@ -907,7 +907,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', flist.append(kdx) pedg.append(kdx) if be_type=='be4': - if kdx in lunit: + if kdx in lunit: lmin1 = runit[numpy.where(lunit==kdx)[0]] for zdx in lmin1: if zdx in lsts and klsts[lsts.index(zdx)] == nk1*(nk2-1)+1 \ @@ -918,7 +918,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', klsts.append(nk1*(nk2-1)+1) else: klsts.append(nk1) - if kdx in runit: + if kdx in runit: rmin1 = lunit[numpy.where(runit==kdx)[0]] for zdx in rmin1: if zdx in rsts and krsts[rsts.index(zdx)] == nk1+1 and twoD: @@ -929,7 +929,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', krsts.append(nk1+1) else: krsts.append(2) - if kdx in uunit: + if kdx in uunit: umin1 = dunit[numpy.where(uunit==kdx)[0]] for zdx in umin1: if zdx in usts and kusts[usts.index(zdx)] == 2 and twoD: @@ -960,8 +960,8 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', continue flist.append(zdx) tsts.append(zdx) - ktsts.append(nk1+2) - + ktsts.append(nk1+2) + for ldx, l in enumerate(coord): if ldx==kdx or ldx==jdx or cell.atom_pure_symbol(ldx) == 'H'or ldx in pedg: continue @@ -969,7 +969,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', if dist <= bond: flist.append(ldx) pedg.append(ldx) - + lsites.append(lsts) rsites.append(rsts) usites.append(usts) @@ -982,25 +982,25 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', kdsites.append(kdsts) kmsites.append(kmsts) ktsites.append(ktsts) - + Frag.append(flist) pedge.append(pedg) hlist = [[] for i in coord] for idx, i in enumerate(normlist): if cell.atom_pure_symbol(idx) == 'H': - + tmplist = normlist - i tmplist = list(tmplist) - + clist = [] for jdx,j in enumerate(tmplist): if not idx==jdx and not cell.atom_pure_symbol(jdx) == 'H': if abs(j)< normdist: clist.append(jdx) - + for jdx in clist: - + dist = numpy.linalg.norm(coord[idx] - coord[jdx]) if dist <= hbond: hlist[jdx].append(idx) @@ -1010,7 +1010,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', print('--------------------------',flush=True) print('Fragment | Center | Edges ',flush=True) print('--------------------------',flush=True) - + for idx,i in enumerate(Frag): print(' {:>4} | {:>5} |'.format(idx, cell.atom_pure_symbol(cen[idx])+str(cen[idx]+1)), @@ -1032,7 +1032,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', w = open('fragments.xyz','w') for idx,i in enumerate(Frag): w.write(str(len(i)+len(hlist[cen[idx]])+len(hlist[j]))+'\n') - w.write('Fragment - '+str(idx)+'\n') + w.write('Fragment - '+str(idx)+'\n') for j in hlist[cen[idx]]: w.write(' {:>3} {:>10.7f} {:>10.7f} {:>10.7f} \n'.format(cell.atom_pure_symbol(j), coord[j][0]/ang2bohr, @@ -1055,42 +1055,42 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', pao = True else: pao = False - + if pao: cell2 = cell.copy() cell2.basis = valence_basis cell2.build() - bas2list = cell2.aoslice_by_atom() + bas2list = cell2.aoslice_by_atom() nbas2 = [0 for i in range(cell.natm)] - + baslist = cell.aoslice_by_atom() sites__ = [[] for i in coord] coreshift = 0 hshift = [0 for i in coord] - + for adx in range(cell.natm): - + if not cell.atom_pure_symbol(adx) == 'H': - bas = baslist[adx] + bas = baslist[adx] start_ = bas[2] stop_ = bas[3] if pao: bas2 = bas2list[adx] nbas2[adx] += bas2[3] - bas2[2] - + if frozen_core: - start_ -= coreshift + start_ -= coreshift ncore_ = core_list[adx] stop_ -= coreshift+ncore_ - if pao: nbas2[adx] -= ncore_ + if pao: nbas2[adx] -= ncore_ coreshift += ncore_ b1list = [i for i in range(start_, stop_)] sites__[adx] = b1list else: hshift[adx] = coreshift - + hsites = [[] for i in coord] nbas2H = [0 for i in coord] for hdx, h in enumerate(hlist): @@ -1104,19 +1104,19 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', if pao: bas2H = bas2list[hidx] nbas2H[hdx] += bas2H[3] - bas2H[2] - - if frozen_core: + + if frozen_core: startH -= hshift[hidx] stopH -= hshift[hidx] - + b1list = [i for i in range(startH, stopH)] hsites[hdx].extend(b1list) - + max_site = max([j for i in sites__ for j in i]) if any(hsites): maxH = max([j for i in hsites for j in i]) max_site = max(max_site, maxH) - + fsites = [] edgsites = [] edge_idx = [] @@ -1130,15 +1130,15 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', Ns = Nsite+1 uNs = Ns*unitcell Sites = sites__.copy() - - for j_ in range(unitcell-1): - for idx, i in enumerate(Sites): + + for j_ in range(unitcell-1): + for idx, i in enumerate(Sites): if any (i>=unitcell*Ns + unitcell*Ns*j_): sites__[idx] = [ (nk1*Ns)+j+(nk1*Ns*j_ - unitcell*Ns*j_) - (unitcell*Ns) for j in i] - + for idx, i in enumerate(Frag): ftmp = [] - ftmpe = [] + ftmpe = [] indix = 0 edind = [] edg = [] @@ -1147,14 +1147,14 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', if not conmax: frglist = sites__[jdx].copy() frglist.extend(hsites[jdx]) - elif nkcon: + elif nkcon: numk = klsites[idx][jdx_]-1 frglist = kfrag_func(sites__[jdx]+hsites[jdx], numk, nk1, uNs, Ns) else: frglist = [pq-max_site-1 for pq in sites__[jdx]] frglist.extend([pq-max_site-1 for pq in hsites[jdx]]) ftmp.extend(frglist) - ls = len(sites__[jdx])+len(hsites[jdx]) + ls = len(sites__[jdx])+len(hsites[jdx]) if not pao: if not conmax: edglist = sites__[jdx].copy() @@ -1164,7 +1164,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', numk, nk1, uNs, Ns) else: edglist = [pq-max_site-1 for pq in sites__[jdx]] - edglist.extend([pq-max_site-1 for pq in hsites[jdx]]) + edglist.extend([pq-max_site-1 for pq in hsites[jdx]]) ftmpe.append(edglist) edind.append([pq for pq in range(indix,indix+ls)]) else: @@ -1179,7 +1179,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', edglist.extend([pq-max_site-1 for pq in hsites[jdx]][:nbas2H[jdx]]) ftmpe.append(edglist) - ind__ = [ indix+frglist.index(pq) for pq in edglist] + ind__ = [ indix+frglist.index(pq) for pq in edglist] edind.append(ind__) indix += ls @@ -1194,8 +1194,8 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', else: frglist = [pq-max_site-1 for pq in sites__[jdx]] frglist.extend([pq-max_site-1 for pq in hsites[jdx]]) - - ftmp.extend(frglist) + + ftmp.extend(frglist) ls = len(sites__[jdx])+len(hsites[jdx]) if not pao: if not conmax: @@ -1206,8 +1206,8 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', numk, nk1, uNs, Ns) else: edglist = [pq-max_site-1 for pq in sites__[jdx]] - edglist.extend([pq-max_site-1 for pq in hsites[jdx]]) - + edglist.extend([pq-max_site-1 for pq in hsites[jdx]]) + ftmpe.append(edglist) edind.append([pq for pq in range(indix,indix+ls)]) else: @@ -1220,12 +1220,12 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', else: edglist = [pq-max_site-1 for pq in sites__[jdx]][:nbas2[jdx]] edglist.extend([pq-max_site-1 for pq in hsites[jdx]][:nbas2H[jdx]]) - + ftmpe.append(edglist) - ind__ = [ indix+frglist.index(pq) for pq in edglist] + ind__ = [ indix+frglist.index(pq) for pq in edglist] edind.append(ind__) indix += ls - + for jdx_,jdx in enumerate(msites[idx]): edg.append(jdx) if not conmax: @@ -1235,11 +1235,11 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', numk = kmsites[idx][jdx_]-1 frglist = kfrag_func(sites__[jdx]+hsites[jdx], numk, nk1, uNs, Ns) - + else: frglist = [pq-max_site-1 for pq in sites__[jdx]] frglist.extend([pq-max_site-1 for pq in hsites[jdx]]) - + ftmp.extend(frglist) ls = len(sites__[jdx])+len(hsites[jdx]) if not pao: @@ -1251,7 +1251,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', numk, nk1, uNs, Ns) else: edglist = [pq-max_site-1 for pq in sites__[jdx]] - edglist.extend([pq-max_site-1 for pq in hsites[jdx]]) + edglist.extend([pq-max_site-1 for pq in hsites[jdx]]) ftmpe.append(edglist) edind.append([pq for pq in range(indix,indix+ls)]) else: @@ -1264,29 +1264,29 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', else: edglist = [pq-max_site-1 for pq in sites__[jdx]][:nbas2[jdx]] edglist.extend([pq-max_site-1 for pq in hsites[jdx]][:nbas2H[jdx]]) - + ftmpe.append(edglist) - ind__ = [ indix+frglist.index(pq) for pq in edglist] + ind__ = [ indix+frglist.index(pq) for pq in edglist] edind.append(ind__) - indix += ls - + indix += ls + frglist = sites__[cen[idx]].copy() frglist.extend(hsites[cen[idx]]) ftmp.extend(frglist) ls = len(sites__[cen[idx]])+len(hsites[cen[idx]]) if not pao: - centerf_idx.append([pq for pq in range(indix,indix+ls)]) + centerf_idx.append([pq for pq in range(indix,indix+ls)]) else: cntlist = sites__[cen[idx]].copy()[:nbas2[cen[idx]]] cntlist.extend(hsites[cen[idx]][:nbas2H[cen[idx]]]) - ind__ = [ indix+frglist.index(pq) for pq in cntlist] + ind__ = [ indix+frglist.index(pq) for pq in cntlist] centerf_idx.append(ind__) indix += ls - + for jdx in pedge[idx]: edg.append(jdx) - frglist = sites__[jdx].copy() + frglist = sites__[jdx].copy() frglist.extend(hsites[jdx]) ftmp.extend(frglist) ls = len(sites__[jdx]) + len(hsites[jdx]) @@ -1298,13 +1298,13 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', else: edglist = sites__[jdx][:nbas2[jdx]].copy() edglist.extend(hsites[jdx][:nbas2H[jdx]]) - + ftmpe.append(edglist) - ind__ = [ indix+frglist.index(pq) for pq in edglist] - edind.append(ind__) - indix += ls + ind__ = [ indix+frglist.index(pq) for pq in edglist] + edind.append(ind__) + indix += ls + - for jdx_,jdx in enumerate(rsites[idx]): edg.append(jdx) if not conmax: @@ -1318,9 +1318,9 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', frglist = [pq+max_site+1 for pq in sites__[jdx]] frglist.extend([pq+max_site+1 for pq in hsites[jdx]]) ftmp.extend(frglist) - + ls = len(sites__[jdx])+len(hsites[jdx]) - if not pao: + if not pao: if not conmax: edglist = sites__[jdx].copy() edglist.extend(hsites[jdx]) @@ -1329,7 +1329,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', numk, nk1, uNs, Ns) else: edglist = [pq-max_site+1 for pq in sites__[jdx]] - edglist.extend([pq-max_site+1 for pq in hsites[jdx]]) + edglist.extend([pq-max_site+1 for pq in hsites[jdx]]) ftmpe.append(edglist) edind.append([pq for pq in range(indix,indix+ls)]) else: @@ -1344,7 +1344,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', edglist.extend([pq+max_site+1 for pq in hsites[jdx]][:nbas2H[jdx]]) ftmpe.append(edglist) - ind__ = [ indix+frglist.index(pq) for pq in edglist] + ind__ = [ indix+frglist.index(pq) for pq in edglist] edind.append(ind__) indix += ls @@ -1353,16 +1353,16 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', if not conmax: frglist = sites__[jdx].copy() frglist.extend(hsites[jdx]) - elif nkcon: + elif nkcon: numk = kdsites[idx][jdx_]-1 frglist = kfrag_func(sites__[jdx]+hsites[jdx], numk, nk1, uNs, Ns) else: frglist = [pq+max_site+1 for pq in sites__[jdx]] frglist.extend([pq+max_site+1 for pq in hsites[jdx]]) - - ftmp.extend(frglist) + + ftmp.extend(frglist) ls = len(sites__[jdx])+len(hsites[jdx]) - if not pao: + if not pao: if not conmax: edglist = sites__[jdx].copy() edglist.extend(hsites[jdx]) @@ -1372,7 +1372,7 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', else: edglist = [pq-max_site+1 for pq in sites__[jdx]] edglist.extend([pq-max_site+1 for pq in hsites[jdx]]) - + ftmpe.append(edglist) edind.append([pq for pq in range(indix,indix+ls)]) else: @@ -1385,9 +1385,9 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', else: edglist = [pq+max_site+1 for pq in sites__[jdx]][:nbas2[jdx]] edglist.extend([pq+max_site+1 for pq in hsites[jdx]][:nbas2H[jdx]]) - + ftmpe.append(edglist) - ind__ = [ indix+frglist.index(pq) for pq in edglist] + ind__ = [ indix+frglist.index(pq) for pq in edglist] edind.append(ind__) indix += ls @@ -1402,19 +1402,19 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', else: frglist = [pq+max_site+1 for pq in sites__[jdx]] frglist.extend([pq+max_site+1 for pq in hsites[jdx]]) - - ftmp.extend(frglist) + + ftmp.extend(frglist) ls = len(sites__[jdx])+len(hsites[jdx]) - if not pao: + if not pao: if not conmax: edglist = sites__[jdx].copy() edglist.extend(hsites[jdx]) elif nkcon: edglist = kfrag_func(sites__[jdx]+hsites[jdx], - numk, nk1, uNs, Ns) + numk, nk1, uNs, Ns) else: edglist = [pq-max_site+1 for pq in sites__[jdx]] - edglist.extend([pq-max_site+1 for pq in hsites[jdx]]) + edglist.extend([pq-max_site+1 for pq in hsites[jdx]]) ftmpe.append(edglist) edind.append([pq for pq in range(indix,indix+ls)]) else: @@ -1423,38 +1423,38 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', edglist.extend(hsites[jdx][:nbas2H[jdx]]) elif nkcon: edglist =kfrag_func(sites__[jdx][:nbas2[jdx]]+hsites[jdx][:nbas2H[jdx]], - numk, nk1, uNs, Ns) + numk, nk1, uNs, Ns) else: edglist = [pq+max_site+1 for pq in sites__[jdx]][:nbas2[jdx]] edglist.extend([pq+max_site+1 for pq in hsites[jdx]][:nbas2H[jdx]]) - + ftmpe.append(edglist) - ind__ = [ indix+frglist.index(pq) for pq in edglist] + ind__ = [ indix+frglist.index(pq) for pq in edglist] edind.append(ind__) - indix += ls + indix += ls edge.append(edg) fsites.append(ftmp) edgsites.append(ftmpe) edge_idx.append(edind) - + center = [] for ix in edge: cen_ = [] for jx in ix: cen_.append(cen.index(jx)) center.append(cen_) - - Nfrag = len(fsites) + + Nfrag = len(fsites) ebe_weight=[] # Use IAO+PAO for computing energy for ix, i in enumerate(fsites): - + tmp_ = [i.index(pq) for pq in sites__[cen[ix]]] - tmp_.extend([i.index(pq) for pq in hsites[cen[ix]]]) + tmp_.extend([i.index(pq) for pq in hsites[cen[ix]]]) ebe_weight.append([1.0, tmp_]) - + # Center of a fragment are defined in cen[idx] # center[[idx,jdx]] defines fragments idx,jdx who's cen[idx],cen[jdx] \\ # centers are matched to the edges. @@ -1470,9 +1470,9 @@ def autogen(mol, kpt, frozen_core=True, be_type='be2', cntlist = sites__[cen[j]].copy()[:nbas2[cen[j]]] cntlist.extend(hsites[cen[j]][:nbas2H[cen[j]]]) idx.append([fsites[j].index(k) for k in cntlist]) - + center_idx.append(idx) - + return(fsites, edgsites, center, edge_idx, center_idx, centerf_idx, ebe_weight) diff --git a/kbe/chain.py b/kbe/chain.py index e77fefa0..ad481b65 100644 --- a/kbe/chain.py +++ b/kbe/chain.py @@ -19,10 +19,10 @@ def findH(mol,nh, tmphlist = []): hidx.append(idx) return hidx - + def polychain(self, mol, frozen_core=False, unitcell=1): """ - Hard coded fragmentation for polymer chains. This is not recommended for any production level calculations. + Hard coded fragmentation for polymer chains. This is not recommended for any production level calculations. """ # group H @@ -33,91 +33,91 @@ def polychain(self, mol, frozen_core=False, unitcell=1): mol2.basis = 'sto-3g' mol2.build() bas2list = mol2.aoslice_by_atom() - + sites__ = [] sites2 = [] coreshift = 0 tmphlist = [] for adx in range(mol.natm): - + if not mol.atom_pure_symbol(adx) == 'H': h1 = findH(mol, mol.atom_coord(adx),tmphlist = tmphlist ) tmphlist.extend(h1) - + bas = baslist[adx] bas2 = bas2list[adx] - + start_ = bas[2] stop_ = bas[3] start2_ = bas2[2] stop2_ = bas2[3] - + sites_ = [] sites2_ = [] for hidx in h1: - + basH = baslist[hidx] basH2 = bas2list[hidx] - + startH_ = basH[2] stopH_ = basH[3] - + startH2_ = basH2[2] - stopH2_ = basH2[3] - + stopH2_ = basH2[3] + if frozen_core: startH_ -= coreshift stopH_ -= coreshift - + startH2_ -= coreshift stopH2_ -= coreshift - + b1list = [i for i in range(startH_, stopH_)] b2list = [i for i in range(startH2_, stopH2_)] sites2_.extend([ True if i in b2list else False for i in b1list]) sites_.extend(b1list) - + if frozen_core: start_ -= coreshift start2_ -= coreshift ncore_ = self.core_list[adx] stop_ -= coreshift+ncore_ stop2_ -= coreshift+ncore_ - coreshift += ncore_ + coreshift += ncore_ b1list = [i for i in range(start_, stop_)] b2list = [i for i in range(start2_, stop2_)] - + sites_.extend(b1list) sites2_.extend([ True if i in b2list else False for i in b1list]) - + sites__.append(sites_) sites2.append(sites2_) - + if unitcell > 1: sites = [] - if isinstance(unitcell, int): + if isinstance(unitcell, int): for i in range(unitcell): - if i: + if i: ns_ = nsites[-1][-1] nsites = [] for p in sites: nsites.append([q+ns_+1 for q in p]) else: - nsites = sites__ + nsites = sites__ sites = [*sites, *nsites] else: - int_sites = int(unitcell) + int_sites = int(unitcell) frac_sites = int(len(sites__)*(unitcell-int_sites)) for i in range(int_sites): - if i: + if i: ns_ = nsites[-1][-1] nsites = [] for p in sites: nsites.append([q+ns_+1 for q in p]) else: - nsites = sites__ + nsites = sites__ sites = [*sites, *nsites] ns_ = sites[-1][-1] nsites = [] @@ -125,48 +125,48 @@ def polychain(self, mol, frozen_core=False, unitcell=1): nsites.append([q+ns_+1 for q in sites__[i]]) sites = [*sites, *nsites] elif unitcell <1: - frac_sites =int(len(sites__)*(unitcell)) + frac_sites =int(len(sites__)*(unitcell)) ns_ = sites__[-1][-1] - nsites = [] + nsites = [] for i in range(frac_sites): nsites.append([q for q in sites__[i]]) - sites = nsites + sites = nsites else: - sites = sites__ - + sites = sites__ + if self.be_type=='be2': - fs=[] + fs=[] if not self.self_match and not self.allcen: - + for i in range(len(sites)-2): self.fsites.append(sites[i]+ sites[i+1]+ sites[i+2]) fs.append([ sites[i], sites[i+1], sites[i+2]]) - + self.Nfrag = len(self.fsites) self.edge.append([fs[0][2]]) for i in fs[1:-1]: self.edge.append([i[0],i[-1]]) self.edge.append([fs[-1][0]]) - + self.center.append([1]) for i in range(self.Nfrag-2): self.center.append([i,i+2]) - self.center.append([self.Nfrag-2]) + self.center.append([self.Nfrag-2]) for ix, i in enumerate(self.fsites): tmp_ = [] elist_ = [ xx for yy in self.edge[ix] for xx in yy] - + for j in i: if not j in elist_: tmp_.append(i.index(j)) - + self.ebe_weight.append([1.0, tmp_]) - elif self.allcen: + elif self.allcen: sites_left = [-i-1 for i in sites[0]] ns = len(sites[-1]) sites_right = [i+ns for i in sites[-1]] - + self.fsites.append(sites_left+sites[0]+sites[1]) fs.append([sites_left, sites[0], sites[1]]) for i in range(len(sites)-2): @@ -177,27 +177,27 @@ def polychain(self, mol, frozen_core=False, unitcell=1): self.Nfrag = len(self.fsites) for i in fs: #[1:-1]: - self.edge.append([i[0],i[-1]]) - + self.edge.append([i[0],i[-1]]) + self.center.append([self.Nfrag-1, 1]) for i in range(self.Nfrag-2): self.center.append([i,i+2]) self.center.append([self.Nfrag-2, 0]) - + for ix, i in enumerate(self.fsites): - tmp_ = [] - + tmp_ = [] + elist_ = [xx for yy in self.edge[ix] for xx in yy] for j in i: if not j in elist_: tmp_.append(i.index(j)) - self.ebe_weight.append([1.0, tmp_]) + self.ebe_weight.append([1.0, tmp_]) else: - + for i in range(len(sites)-2): self.fsites.append(sites[i]+sites[i+1]+sites[i+2]) fs.append([ sites[i], sites[i+1], sites[i+2]]) self.Nfrag = len(self.fsites) - + #self.edge.append([fs[0][2]]) for i in fs:#[1:-1]: self.edge.append([i[0],i[-1]]) @@ -225,24 +225,24 @@ def polychain(self, mol, frozen_core=False, unitcell=1): for i in range(self.Nfrag): self.centerf_idx.append([self.fsites[i].index(j) for j in fs[i][1]]) print(self.fsites) - + elif self.be_type=='be3': - fs=[] - if not self.self_match and not self.allcen: + fs=[] + if not self.self_match and not self.allcen: # change back to i,i+1,i+2 for i in range(len(sites)-4): self.fsites.append(sites[i]+ sites[i+1]+ sites[i+2]+ sites[i+3]+ sites[i+4]) fs.append([ sites[i], sites[i+1], sites[i+2], sites[i+3], sites[i+4]]) - - self.Nfrag = len(self.fsites) + + self.Nfrag = len(self.fsites) self.edge.append([fs[0][3], fs[0][4]]) # change back 0->1 for i in fs[1:-1]: self.edge.append([i[0],i[1], i[-2], i[-1]]) self.edge.append([fs[-1][0], fs[-1][1]]) - + self.center.append([1,2]) self.center.append([0,0,2,3]) for i in range(self.Nfrag-4): @@ -250,9 +250,9 @@ def polychain(self, mol, frozen_core=False, unitcell=1): self.center.append([self.Nfrag-4, self.Nfrag-3, self.Nfrag-1, self.Nfrag-1]) self.center.append([self.Nfrag-3, self.Nfrag-2]) - + for i in range(self.Nfrag): - self.centerf_idx.append([self.fsites[i].index(j) for j in fs[i][2]]) + self.centerf_idx.append([self.fsites[i].index(j) for j in fs[i][2]]) for ix, i in enumerate(self.fsites): tmp_ = [] elist_ = [ xx for yy in self.edge[ix] for xx in yy] @@ -280,7 +280,7 @@ def polychain(self, mol, frozen_core=False, unitcell=1): self.fsites.append(sites[-3]+sites[-2]+sites[-1]+sites_right0+sites_right1) fs.append([sites[-4],sites[-3],sites[-2],sites[-1],sites_right0]) fs.append([sites[-3],sites[-2],sites[-1],sites_right0,sites_right1]) - + self.Nfrag = len(self.fsites) #self.edge.append([fs[0][3], fs[0][4]]) @@ -304,9 +304,9 @@ def polychain(self, mol, frozen_core=False, unitcell=1): for j in i: if not j in elist_: tmp_.append(i.index(j)) self.ebe_weight.append([1.0, tmp_]) - + else: - + for i in range(len(sites)-4): self.fsites.append(sites[i]+ sites[i+1]+ sites[i+2]+ sites[i+3]+ sites[i+4]) @@ -327,7 +327,7 @@ def polychain(self, mol, frozen_core=False, unitcell=1): if Nfrag > 2: self.center.append([Nfrag-4, Nfrag-3, Nfrag-1, 0]) self.center.append([Nfrag-3, Nfrag-2, 0, 1]) - + for ix, i in enumerate(self.fsites): tmp_ = [] if ix==0: @@ -338,36 +338,36 @@ def polychain(self, mol, frozen_core=False, unitcell=1): if not j in elist_: tmp_.append(i.index(j)) if ix == self.Nfrag-1: for edg_ix in range(2): - + tmp_.extend([i.index(k) for k in self.edge[ix][edg_ix+2]]) self.ebe_weight.append([1.0, tmp_]) - + # center on each fragments ?? do we add more for PBC for i in range(self.Nfrag): self.centerf_idx.append([self.fsites[i].index(j) for j in fs[i][2]]) - + print(flush=True) print(' No. of fragments : ',self.Nfrag,flush=True) - print(flush=True) - + print(flush=True) + if not self.be_type=='be1': for i in range(self.Nfrag): idx = [] for j in self.edge[i]: idx.append([self.fsites[i].index(k) for k in j]) - self.edge_idx.append(idx) + self.edge_idx.append(idx) if not self.self_match and not self.allcen: for i in range(self.Nfrag): idx = [] for j in range(len(self.center[i])): - + idx.append([self.fsites[self.center[i][j]].index(k) for k in self.edge[i][j]]) - self.center_idx.append(idx) + self.center_idx.append(idx) else: for i in range(self.Nfrag): - idx = [] - + idx = [] + for j in self.center[i]: idx__ = [] for fidx, fval in enumerate(self.fsites[j]): @@ -377,5 +377,5 @@ def polychain(self, mol, frozen_core=False, unitcell=1): self.center_idx.append(idx) - - + + diff --git a/kbe/fragment.py b/kbe/fragment.py index c6dc2f9c..8ab72eba 100644 --- a/kbe/fragment.py +++ b/kbe/fragment.py @@ -9,7 +9,7 @@ def print_mol_missing(): flush=True) print('exiting',flush=True) sys.exit() - + class fragpart: def __init__(self, natom=0, dim=1, frag_type='autogen', @@ -22,7 +22,7 @@ def __init__(self, natom=0, dim=1, frag_type='autogen', """Fragment/partitioning definition Interfaces two main fragmentation functions (autogen & polychain) in MolBE. It defines edge & - center for density matching and energy estimation. It also forms the base for IAO/PAO partitioning for + center for density matching and energy estimation. It also forms the base for IAO/PAO partitioning for a large basis set bootstrap calculation. Fragments are constructed based on atoms within a unitcell. Parameters @@ -42,7 +42,7 @@ def __init__(self, natom=0, dim=1, frag_type='autogen', valence_basis: str Name of minimal basis set for IAO scheme. 'sto-3g' suffice for most cases. valence_only: bool - If this option is set to True, all calculation will be performed in the valence basis in the IAO partitioning. + If this option is set to True, all calculation will be performed in the valence basis in the IAO partitioning. This is an experimental feature. frozen_core: bool Whether to invoke frozen core approximation. This is set to False by default @@ -77,11 +77,11 @@ def __init__(self, natom=0, dim=1, frag_type='autogen', self.valence_basis=valence_basis self.kpt = kpt self.molecule = False ### remove this - + # Check for frozen core approximation if frozen_core: self.ncore, self.no_core_idx, self.core_list = get_core(mol) - + if frag_type=='polychain': if mol is None: print_mol_missing() self.polychain(mol, frozen_core=frozen_core, unitcell=unitcell) @@ -92,21 +92,21 @@ def __init__(self, natom=0, dim=1, frag_type='autogen', flush=True) print('exiting',flush=True) sys.exit() - + fgs = autogen(mol, kpt, be_type=be_type, frozen_core=frozen_core, valence_basis=valence_basis, unitcell=unitcell, nx=nx, ny=ny, nz=nz, long_bond=long_bond, perpend_dist = perpend_dist, perpend_dist_tol = perpend_dist_tol, gamma_2d=gamma_2d, gamma_1d=gamma_1d,interlayer=interlayer) - - self.fsites, self.edge, self.center, self.edge_idx, self.center_idx, self.centerf_idx, self.ebe_weight = fgs + + self.fsites, self.edge, self.center, self.edge_idx, self.center_idx, self.centerf_idx, self.ebe_weight = fgs self.Nfrag = len(self.fsites) - + else: print('Fragmentation type = ',frag_type,' not implemented!', flush=True) print('exiting',flush=True) sys.exit() - + from .chain import polychain diff --git a/kbe/helper.py b/kbe/helper.py index fee2dbf2..6a2a88a9 100644 --- a/kbe/helper.py +++ b/kbe/helper.py @@ -25,7 +25,7 @@ def get_veff(eri_, dm, S, TA, hf_veff, return_veff0=False): Hartree-Fock effective potential for the full system. """ - + # construct rdm nk, nao, neo = TA.shape P_ = numpy.zeros((neo, neo), dtype=numpy.complex128) @@ -40,18 +40,18 @@ def get_veff(eri_, dm, S, TA, hf_veff, return_veff0=False): eri_ = numpy.asarray(eri_, dtype=numpy.double) vj, vk = scf.hf.dot_eri_dm(eri_, P_, hermi=1, with_j=True, with_k=True) Veff_ = vj - 0.5*vk - + # remove core contribution from hf_veff - + Veff0 = numpy.zeros((neo, neo), dtype=numpy.complex128) for k in range(nk): Veff0 += functools.reduce(numpy.dot, (TA[k].conj().T, hf_veff[k], TA[k])) Veff0 /= float(nk) - + Veff = Veff0 - Veff_ if return_veff0: return(Veff0, Veff) - + return Veff diff --git a/kbe/lo.py b/kbe/lo.py index 93afecf8..f51b7c3b 100644 --- a/kbe/lo.py +++ b/kbe/lo.py @@ -35,7 +35,7 @@ def localize(self, lo_method, mol=None, valence_basis='sto-3g', iao_wannier=True valence_basis: str Name of minimal basis set for IAO scheme. 'sto-3g' suffice for most cases. valence_only: bool - If this option is set to True, all calculation will be performed in the valence basis in the IAO partitioning. + If this option is set to True, all calculation will be performed in the valence basis in the IAO partitioning. This is an experimental feature. iao_wannier : bool Whether to perform Wannier localization in the IAO space @@ -54,9 +54,9 @@ def localize(self, lo_method, mol=None, valence_basis='sto-3g', iao_wannier=True elif iao_val_core: sys.exit('valence_basis='+valence_basis+' not supported for iao_val_core=True') - + if lo_method == 'lowdin': - + # Lowdin orthogonalization with k-points W = numpy.zeros_like(self.S) nk, nao, nmo = self.C.shape @@ -67,18 +67,18 @@ def localize(self, lo_method, mol=None, valence_basis='sto-3g', iao_wannier=True else: lmo_coeff = numpy.zeros_like(self.C) cinv_ = numpy.zeros((nk, nmo, nao), dtype=numpy.complex128) - + for k in range(self.nkpt): - + es_, vs_ = scipy.linalg.eigh(self.S[k]) edx = es_ > 1.e-14 - + W[k] = numpy.dot(vs_[:,edx]/numpy.sqrt(es_[edx]), vs_[:,edx].conj().T) for i in range(W[k].shape[1]): if W[k][i,i] < 0: W[:,i] *= -1 if self.frozen_core: - + pcore = numpy.eye(W[k].shape[0]) - numpy.dot(self.P_core[k], self.S[k]) C_ = numpy.dot(pcore, W[k]) @@ -88,18 +88,18 @@ def localize(self, lo_method, mol=None, valence_basis='sto-3g', iao_wannier=True Cpop = functools.reduce(numpy.dot, (C_.conj().T, self.S[k], C_)) Cpop = numpy.diag(Cpop.real) - + no_core_idx = numpy.where(Cpop > 0.7)[0] C_ = C_[:,no_core_idx] - + S_ = functools.reduce(numpy.dot, (C_.conj().T, self.S[k], C_)) - + es_, vs_ = scipy.linalg.eigh(S_) edx = es_ > 1.e-14 W_ = numpy.dot(vs_[:,edx]/numpy.sqrt(es_[edx]), vs_[:,edx].conj().T) W_nocore[k] = numpy.dot(C_,W_) - + lmo_coeff[k] = functools.reduce(numpy.dot, (W_nocore[k].conj().T,self.S[k], self.C[k][:,self.ncore:])) @@ -114,49 +114,49 @@ def localize(self, lo_method, mol=None, valence_basis='sto-3g', iao_wannier=True else: self.W = W self.lmo_coeff = lmo_coeff - self.cinv = cinv_ + self.cinv = cinv_ elif lo_method=='iao': from libdmet.lo import pywannier90 from pyscf import lo import os, h5py - + from .lo_k import get_xovlp_k, get_iao_k, get_pao_k, get_pao_native_k,symm_orth_k , remove_core_mo_k from pyscf import gto, lo if not iao_val_core or not self.frozen_core: Co = self.C[:,:,:self.Nocc].copy() S12, S2 = get_xovlp_k(self.cell, self.kpts, basis=valence_basis) - ciao_ = get_iao_k(Co, S12, self.S, S2=S2) - + ciao_ = get_iao_k(Co, S12, self.S, S2=S2) + arrange_by_atom = True # tmp - aos are not rearrange and so below is not necessary if arrange_by_atom: - + nk, nao, nlo = ciao_.shape Ciao_ = numpy.zeros((nk, nao, nlo), dtype=numpy.complex128) for k in range(self.nkpt): - aoind_by_atom = get_aoind_by_atom(self.cell) + aoind_by_atom = get_aoind_by_atom(self.cell) ctmp, iaoind_by_atom = reorder_by_atom_(ciao_[k], aoind_by_atom, self.S[k]) Ciao_[k] = ctmp else: Ciao_ = ciao_.copy() - + # get_pao_k returns canonical orthogonalized orbitals #Cpao = get_pao_k(Ciao, self.S, S12, S2, self.cell) # get_pao_native_k returns symm orthogonalized orbitals cpao_ = get_pao_native_k(Ciao_, self.S, self.cell, valence_basis, self.kpts) - + if arrange_by_atom: nk, nao, nlo = cpao_.shape Cpao_ = numpy.zeros((nk, nao, nlo), dtype=numpy.complex128) for k in range(self.nkpt): aoind_by_atom = get_aoind_by_atom(self.cell) ctmp, paoind_by_atom = reorder_by_atom_(cpao_[k], aoind_by_atom, self.S[k]) - Cpao_[k] = ctmp + Cpao_[k] = ctmp else: Cpao_ = cpao_.copy() - + nk, nao, nlo = Ciao_.shape if self.frozen_core: nk, nao, nlo = Ciao_.shape @@ -209,35 +209,35 @@ def localize(self, lo_method, mol=None, valence_basis='sto-3g', iao_wannier=True c_core_val = numpy.zeros((nk_, nao, Ciao_.shape[-1]+self.ncore), dtype=Ciao_.dtype) for k in range(nk_): c_core_val[k] = numpy.hstack((ciao_core[k], Ciao_[k])) - + arrange_by_atom = True # tmp - aos are not rearrange and so below is not necessary (iaoind_by_atom is used to stack iao|pao later) if arrange_by_atom: nk, nao, nlo = c_core_val.shape for k in range(self.nkpt): - aoind_by_atom = get_aoind_by_atom(self.cell) + aoind_by_atom = get_aoind_by_atom(self.cell) ctmp, iaoind_by_atom = reorder_by_atom_(c_core_val[k], aoind_by_atom, self.S[k]) - - cpao_ = get_pao_native_k(c_core_val, self.S, self.cell, valence_basis, self.kpts, ortho=True) + + cpao_ = get_pao_native_k(c_core_val, self.S, self.cell, valence_basis, self.kpts, ortho=True) if arrange_by_atom: nk, nao, nlo = cpao_.shape Cpao_ = numpy.zeros((nk, nao, nlo), dtype=numpy.complex128) for k in range(self.nkpt): aoind_by_atom = get_aoind_by_atom(self.cell) ctmp, paoind_by_atom = reorder_by_atom_(cpao_[k], aoind_by_atom, self.S[k]) - Cpao_[k] = ctmp - + Cpao_[k] = ctmp + Cpao = Cpao_.copy() Ciao = Ciao_.copy() if iao_wannier: - mo_energy_ = [] + mo_energy_ = [] for k in range(nk): fock_iao = reduce(numpy.dot, (Ciao_[k].conj().T, self.FOCK[k], Ciao_[k])) S_iao = reduce(numpy.dot, (Ciao_[k].conj().T, self.S[k], Ciao_[k])) - e_iao, v_iao = scipy.linalg.eigh(fock_iao, S_iao) + e_iao, v_iao = scipy.linalg.eigh(fock_iao, S_iao) mo_energy_.append(e_iao) - iaomf = KMF(self.mol, kpts = self.kpts, mo_coeff = Ciao_, mo_energy = mo_energy_) + iaomf = KMF(self.mol, kpts = self.kpts, mo_coeff = Ciao_, mo_energy = mo_energy_) num_wann = numpy.asarray(iaomf.mo_coeff).shape[2] keywords = ''' @@ -263,55 +263,55 @@ def localize(self, lo_method, mol=None, valence_basis='sto-3g', iao_wannier=True else: ovlp_ciao = uciao[k].conj().T @ self.S[k] @ Ciao[k] A_matrix[k] = ovlp_ciao - A_matrix = A_matrix.transpose(1,2,0) - + A_matrix = A_matrix.transpose(1,2,0) + w90.kernel(A_matrix=A_matrix) u_mat = numpy.array(w90.U_matrix.transpose(2,0,1), order='C', dtype=numpy.complex128) os.system('cp wannier90.wout wannier90_iao.wout') os.system('rm wannier90.*') - + nk, nao, nlo = Ciao_.shape Ciao = numpy.zeros((nk, nao, nlo), dtype=numpy.complex128) - + for k in range(self.nkpt): - Ciao[k] = numpy.dot(Ciao_[k], u_mat[k]) + Ciao[k] = numpy.dot(Ciao_[k], u_mat[k]) # Stack Ciao|Cpao Wstack = numpy.zeros((self.nkpt, Ciao.shape[1], Ciao.shape[2]+Cpao.shape[2]), - dtype=numpy.complex128) + dtype=numpy.complex128) if self.frozen_core: for k in range(self.nkpt): shift = 0 - ncore = 0 + ncore = 0 for ix in range(self.cell.natm): nc = ncore_(self.cell.atom_charge(ix)) ncore += nc - niao = len(iaoind_by_atom[ix]) + niao = len(iaoind_by_atom[ix]) iaoind_ix = [i_ - ncore for i_ in iaoind_by_atom[ix][nc:]] Wstack[k][:, shift:shift+niao-nc] = Ciao[k][:, iaoind_ix] shift += niao-nc npao = len(paoind_by_atom[ix]) - + Wstack[k][:,shift:shift+npao] = Cpao[k][:, paoind_by_atom[ix]] shift += npao else: for k in range(self.nkpt): - shift = 0 + shift = 0 for ix in range(self.cell.natm): niao = len(iaoind_by_atom[ix]) Wstack[k][:, shift:shift+niao] = Ciao[k][:, iaoind_by_atom[ix]] shift += niao npao = len(paoind_by_atom[ix]) Wstack[k][:,shift:shift+npao] = Cpao[k][:, paoind_by_atom[ix]] - shift += npao + shift += npao self.W = Wstack - + nmo = self.C.shape[2] - self.ncore nlo = self.W.shape[2] nao = self.S.shape[2] - + lmo_coeff = numpy.zeros((self.nkpt, nlo, nmo), dtype=numpy.complex128) cinv_ = numpy.zeros((self.nkpt, nlo, nao), dtype=numpy.complex128) @@ -336,9 +336,9 @@ def localize(self, lo_method, mol=None, valence_basis='sto-3g', iao_wannier=True for k in range(self.nkpt): lmo_coeff[k] = reduce(numpy.dot, (self.W[k].conj().T, self.S[k], self.C[k][:,self.ncore:])) cinv_[k] = numpy.dot(self.W[k].conj().T, self.S[k]) - + assert(numpy.allclose(lmo_coeff[k].conj().T @ lmo_coeff[k], numpy.eye(lmo_coeff[k].shape[1]))) - + self.lmo_coeff = lmo_coeff self.cinv = cinv_ @@ -356,8 +356,8 @@ def localize(self, lo_method, mol=None, valence_basis='sto-3g', iao_wannier=True if self.frozen_core: Ccore = self.C[k][:,:self.ncore] - lorb_nocore[k] = remove_core_mo_k(lorb[k], Ccore, self.S[k]) - + lorb_nocore[k] = remove_core_mo_k(lorb[k], Ccore, self.S[k]) + if not self.frozen_core: lmf = KMF(self.mol, kpts=self.kpts, mo_coeff = lorb, mo_energy = self.mo_energy) else: @@ -368,7 +368,7 @@ def localize(self, lo_method, mol=None, valence_basis='sto-3g', iao_wannier=True e__, v__ = scipy.linalg.eigh(fock_lnc, S_lnc) mo_energy_nc.append(e__) lmf = KMF(self.mol, kpts=self.kpts, mo_coeff = lorb_nocore, mo_energy = mo_energy_nc) - + num_wann = lmf.mo_coeff.shape[2] keywords = ''' num_iter = 10000 @@ -378,7 +378,7 @@ def localize(self, lo_method, mol=None, valence_basis='sto-3g', iao_wannier=True iprint = 3 kmesh_tol = 0.00001 ''' - + w90 = pywannier90.W90(lmf, self.kmesh, num_wann, other_keywords=keywords) A_matrix = numpy.zeros((self.nkpt, num_wann, num_wann), dtype=numpy.complex128) @@ -388,7 +388,7 @@ def localize(self, lo_method, mol=None, valence_basis='sto-3g', iao_wannier=True A_matrix[k] = numpy.eye(num_wann, dtype=numpy.complex128) A_matrix = A_matrix.transpose(1,2,0) - + w90.kernel(A_matrix=A_matrix) u_mat = numpy.array(w90.U_matrix.transpose(2,0,1), order='C', dtype=numpy.complex128) @@ -396,18 +396,18 @@ def localize(self, lo_method, mol=None, valence_basis='sto-3g', iao_wannier=True W = numpy.zeros((nk, nao, nlo), dtype=numpy.complex128) for k in range(nk): W[k] = numpy.dot(lmf.mo_coeff[k], u_mat[k]) - + self.W = W lmo_coeff = numpy.zeros((self.nkpt, nlo, nmo-self.ncore), dtype=numpy.complex128) cinv_ = numpy.zeros((self.nkpt, nlo, nao), dtype=numpy.complex128) - + for k in range(nk): lmo_coeff[k] = reduce(numpy.dot, (self.W[k].conj().T, self.S[k], self.C[k][:,self.ncore:])) cinv_[k] = numpy.dot(self.W[k].conj().T, self.S[k]) assert(numpy.allclose(lmo_coeff[k].conj().T @ lmo_coeff[k], numpy.eye(lmo_coeff[k].shape[1]))) self.lmo_coeff = lmo_coeff self.cinv = cinv_ - + else: print('lo_method = ',lo_method,' not implemented!',flush=True) print('exiting',flush=True) diff --git a/kbe/lo_k.py b/kbe/lo_k.py index 1e250f84..265f2a15 100644 --- a/kbe/lo_k.py +++ b/kbe/lo_k.py @@ -1,6 +1,6 @@ # Author(s): Henry Tran # Oinam Meitei -# +# import numpy, sys, scipy from functools import reduce @@ -36,7 +36,7 @@ def get_symm_orth_mat_k(A, thr=1.E-7, ovlp=None): if int(numpy.sum(e < thr)) > 0: raise ValueError("Linear dependence is detected in the column space of A: smallest eigenvalue (%.3E) is less than thr (%.3E). Please use 'cano_orth' instead." % (numpy.min(e), thr)) U = reduce(numpy.dot, (u, numpy.diag(e**-0.5), u.conj().T)) - #U = reduce(numpy.dot, (u/numpy.sqrt(e), u.conj().T)) + #U = reduce(numpy.dot, (u/numpy.sqrt(e), u.conj().T)) return U def symm_orth_k(A, thr=1.E-7, ovlp=None): @@ -113,18 +113,18 @@ def get_iao_k(Co, S12, S1, S2=None, ortho=True): Cotil = reduce(numpy.dot, (P1[k], S12[k], P2[k], S12[k].conj().T, Co[k])) ptil = numpy.dot(P1[k], S12[k]) Stil = reduce(numpy.dot, (Cotil.conj().T, S1[k], Cotil)) - + Po = numpy.dot(Co[k], Co[k].conj().T) - + Stil_inv = numpy.linalg.inv(Stil) Potil = reduce(numpy.dot, (Cotil, Stil_inv, Cotil.conj().T)) - + Ciao[k] = (numpy.eye(nao, dtype=numpy.complex128) - \ numpy.dot((Po + Potil - 2.* reduce(numpy.dot,(Po, S1[k], Potil))), S1[k])) @ ptil if ortho: Ciao[k] = symm_orth_k(Ciao[k], ovlp=S1[k]) - + rep_err = numpy.linalg.norm(Ciao[k] @ Ciao[k].conj().T @ S1[k] @ Po - Po) if rep_err > 1.E-10: raise RuntimeError @@ -152,12 +152,12 @@ def get_pao_k(Ciao, S, S12, S2): Piao = Ciao[k] @ Ciao[k].conj().T @ S[k] cpao_ = (numpy.eye(nao) - Piao)@ nonval - - numpy.o0 = cpao_.shape[-1] + + numpy.o0 = cpao_.shape[-1] Cpao.append(cano_orth(cpao_,ovlp=S[k])) numpy.o1 = Cpao[k].shape[-1] Cpao = numpy.asarray(Cpao) - + return Cpao def get_pao_native_k(Ciao, S, mol, valence_basis, kpts, ortho=True): @@ -170,9 +170,9 @@ def get_pao_native_k(Ciao, S, mol, valence_basis, kpts, ortho=True): Return: Cpao (symmetrically orthogonalized) """ - + nk, nao, niao = Ciao.shape - + # Form a mol object with the valence basis for the ao_labels mol_alt = mol.copy() mol_alt.basis = valence_basis @@ -199,8 +199,8 @@ def get_pao_native_k(Ciao, S, mol, valence_basis, kpts, ortho=True): print("# of PAO: %d --> %d" % (npao0,npao1), flush=True) print("", flush=True) else: - Cpao[k] = cpao_.copy() - + Cpao[k] = cpao_.copy() + return Cpao diff --git a/kbe/misc.py b/kbe/misc.py index 33d55b82..459ddaeb 100644 --- a/kbe/misc.py +++ b/kbe/misc.py @@ -16,11 +16,11 @@ def sgeom(cell, kmesh=None): Number of k-points in each lattice vector dimension """ from pyscf.pbc import tools - + scell = tools.super_cell(cell, kmesh) - + return scell - + def get_phase(cell, kpts, kmesh): a_vec = cell.lattice_vectors() @@ -30,10 +30,10 @@ def get_phase(cell, kpts, kmesh): Rs = numpy.dot(Ts, a_vec) tmp_ = numpy.dot(Rs, kpts.T) - + NRs = Rs.shape[0] phase = 1/numpy.sqrt(NRs) * numpy.exp(1j* numpy.dot(Rs, kpts.T)) - + return phase def get_phase1(cell, kpts, kmesh): diff --git a/kbe/pbe.py b/kbe/pbe.py index 5c82440b..bac49679 100644 --- a/kbe/pbe.py +++ b/kbe/pbe.py @@ -7,7 +7,7 @@ from pyscf import lib import h5py, os -from .misc import storePBE +from .misc import storePBE class BE: """ @@ -28,11 +28,11 @@ class BE: lo_method : str Method for orbital localization, default is 'lowdin'. """ - def __init__(self, mf, fobj, eri_file='eri_file.h5', - lo_method='lowdin',compute_hf=True, + def __init__(self, mf, fobj, eri_file='eri_file.h5', + lo_method='lowdin',compute_hf=True, restart=False, save=False, restart_file='storebe.pk', - mo_energy = None, + mo_energy = None, save_file='storebe.pk',hci_pt=False, nproc=1, ompnum=4, hci_cutoff=0.001, ci_coeff_cutoff = None, select_cutoff=None, @@ -98,7 +98,7 @@ def __init__(self, mf, fobj, eri_file='eri_file.h5', # Fragment information from fobj self.frag_type=fobj.frag_type - self.Nfrag = fobj.Nfrag + self.Nfrag = fobj.Nfrag self.fsites = fobj.fsites self.edge = fobj.edge self.center = fobj.center @@ -111,13 +111,13 @@ def __init__(self, mf, fobj, eri_file='eri_file.h5', self.mol = fobj.mol self.cell = fobj.mol self.kmesh = fobj.kpt - + unitcell_nkpt = 1 for i in self.kmesh: if i>1: unitcell_nkpt *= self.unitcell self.unitcell_nkpt = unitcell_nkpt self.ebe_hf = 0. - + nkpts_ = 1 for i in self.kmesh: if i>1: nkpts_ *= i @@ -129,22 +129,22 @@ def __init__(self, mf, fobj, eri_file='eri_file.h5', self.ci_coeff_cutoff = ci_coeff_cutoff self.select_cutoff = select_cutoff self.hci_pt=hci_pt - - if not restart: - self.mo_energy = mf.mo_energy - mf.exxdiv = None + + if not restart: + self.mo_energy = mf.mo_energy + mf.exxdiv = None self.mf = mf self.Nocc = mf.cell.nelectron//2 - self.enuc = mf.energy_nuc() + self.enuc = mf.energy_nuc() self.hcore = mf.get_hcore() - self.S = mf.get_ovlp() - self.C = numpy.array(mf.mo_coeff) + self.S = mf.get_ovlp() + self.C = numpy.array(mf.mo_coeff) self.hf_dm = mf.make_rdm1() self.hf_veff = mf.get_veff(self.cell, dm_kpts = self.hf_dm, hermi=1, kpts=self.kpts, kpts_band=None) self.hf_etot = mf.e_tot self.W = None self.lmo_coeff = None - + self.print_ini() self.Fobjs = [] self.pot = initialize_pot(self.Nfrag, self.edge_idx) @@ -168,7 +168,7 @@ def __init__(self, mf, fobj, eri_file='eri_file.h5', if cderi: self.cderi = be_var.SCRATCH+str(jobid)+'/'+cderi - + if exxdiv == 'ewald': if not restart: self.ek = self.ewald_sum(kpts=self.kpts) @@ -183,7 +183,7 @@ def __init__(self, mf, fobj, eri_file='eri_file.h5', print('exxdiv = ',exxdiv,'not implemented!',flush=True) print('Energy may diverse.',flush=True) print(flush=True) - + self.frozen_core = False if not fobj.frozen_core else True self.ncore = 0 if not restart: @@ -191,29 +191,29 @@ def __init__(self, mf, fobj, eri_file='eri_file.h5', self.C_core = None self.P_core = None self.core_veff = None - + if self.frozen_core: # Handle frozen core orbitals self.ncore = fobj.ncore self.no_core_idx = fobj.no_core_idx self.core_list = fobj.core_list - + if not restart: self.Nocc -=self.ncore - + nk, nao, nao = self.hf_dm.shape - + dm_nocore = numpy.zeros((nk, nao, nao), dtype=numpy.result_type(self.C, self.C)) C_core = numpy.zeros((nk, nao, self.ncore), dtype=self.C.dtype) P_core = numpy.zeros((nk, nao, nao), dtype=numpy.result_type(self.C, self.C)) - + for k in range(nk): dm_nocore[k]+= 2.*numpy.dot(self.C[k][:,self.ncore:self.ncore+self.Nocc], - self.C[k][:,self.ncore:self.ncore+self.Nocc].conj().T) + self.C[k][:,self.ncore:self.ncore+self.Nocc].conj().T) C_core[k] += self.C[k][:,:self.ncore] P_core[k] += numpy.dot(C_core[k], C_core[k].conj().T) - - + + self.C_core = C_core self.P_core = P_core self.hf_dm = dm_nocore @@ -228,24 +228,24 @@ def __init__(self, mf, fobj, eri_file='eri_file.h5', ecore_h1 /= float(nk) ecore_veff /= float(nk) - E_core = ecore_h1 + ecore_veff + E_core = ecore_h1 + ecore_veff if numpy.abs(E_core.imag).max() < 1.e-10: self.E_core = E_core.real else: - + print('Imaginary density in E_core ', numpy.abs(E_core.imag).max()) sys.exit() for k in range(nk): self.hf_veff[k] -= self.core_veff[k] self.hcore[k] += self.core_veff[k] - + # Needed for Wannier localization if lo_method=='wannier' or iao_wannier: self.FOCK = self.mf.get_fock(self.hcore, self.S, self.hf_veff, self.hf_dm) - + if not restart: # Localize orbitals self.localize(lo_method, mol=self.cell, valence_basis=fobj.valence_basis, @@ -262,13 +262,13 @@ def __init__(self, mf, fobj, eri_file='eri_file.h5', if not restart : self.initialize(mf._eri,compute_hf) - - + + from ._opt import optimize # this is a molbe method not BEOPT from molbe.external.optqn import get_be_error_jacobian from .lo import localize - + def print_ini(self): """ Print initialization banner for the kBE calculation. @@ -284,26 +284,26 @@ def print_ini(self): print(' PP PP BB B EE ',flush=True) print(' PP PP BBBBBBB EEEEEEE ',flush=True) print(flush=True) - - print(' PERIODIC BOOTSTRAP EMBEDDING',flush=True) + + print(' PERIODIC BOOTSTRAP EMBEDDING',flush=True) print(' BEn = ',self.be_type,flush=True) print('-----------------------------------------------------------', flush=True) print(flush=True) - + def ewald_sum(self, kpts=None): from pyscf.pbc.df.df_jk import _ewald_exxdiv_for_G0 - + dm_ = self.mf.make_rdm1() nk, nao, nao = dm_.shape - + vk_kpts = numpy.zeros(dm_.shape) * 1j _ewald_exxdiv_for_G0(self.mf.cell, self.kpts, dm_.reshape(-1, nk, nao, nao), vk_kpts.reshape(-1, nk, nao, nao), self.kpts) e_ = numpy.einsum("kij,kji->",vk_kpts,dm_)*0.25 e_ /= float(nk) - + return e_.real def initialize(self, eri_,compute_hf, restart=False): @@ -321,25 +321,25 @@ def initialize(self, eri_,compute_hf, restart=False): """ from molbe.helper import get_scfObj from multiprocessing import Pool - + import h5py, os, logging from pyscf import ao2mo from pyscf.pbc.df import fft_ao2mo from pyscf.pbc.df import df_ao2mo - from pyscf.pbc import ao2mo as pao2mo + from pyscf.pbc import ao2mo as pao2mo from libdmet.basis_transform.eri_transform import get_emb_eri_fast_gdf - + if compute_hf: E_hf = 0. EH1 = 0. ECOUL = 0. EF = 0. - + # Create a file to store ERIs if not restart: file_eri = h5py.File(self.eri_file,'w') lentmp = len(self.edge_idx) transform_parallel=False # hard set for now - for I in range(self.Nfrag): + for I in range(self.Nfrag): if lentmp: fobjs_ = Frags(self.fsites[I], I, edge=self.edge[I], eri_file=self.eri_file, @@ -353,29 +353,29 @@ def initialize(self, eri_,compute_hf, restart=False): edge_idx=[],center_idx=[],centerf_idx=[], efac=self.ebe_weight[I], unitcell=self.unitcell, unitcell_nkpt=self.unitcell_nkpt) - + fobjs_.sd(self.W, self.lmo_coeff, self.Nocc, kmesh=self.kmesh, cell=self.cell, frag_type=self.frag_type, kpts=self.kpts, h1=self.hcore) - + fobjs_.cons_h1(self.hcore) fobjs_.heff = numpy.zeros_like(fobjs_.h1) fobjs_.dm_init = fobjs_.get_nsocc(self.S, self.C, self.Nocc, ncore=self.ncore) if self.cderi is None: if not restart: - + eri = get_emb_eri_fast_gdf(self.mf.cell, self.mf.with_df, t_reversal_symm=True, symmetry=4, C_ao_emb=fobjs_.TA)[0] - + file_eri.create_dataset(fobjs_.dname, data=eri) eri = ao2mo.restore(8, eri, fobjs_.nao) fobjs_.cons_fock(self.hf_veff, self.S, self.hf_dm, eri_=eri, lmo=self.lmo_coeff, nocc=self.Nocc) else: eri=None self.Fobjs.append(fobjs_) - + # ERI & Fock parallelization for periodic calculations if self.cderi: if self.nproc == 1: @@ -418,14 +418,14 @@ def initialize(self, eri_,compute_hf, restart=False): sys.exit() self.Fobjs[frg].fock = self.Fobjs[frg].h1 + veff_.real veffs = None - + # SCF parallelized if self.nproc == 1 and not transform_parallel: for frg in range(self.Nfrag): # SCF self.Fobjs[frg].scf(fs=True, dm0 = self.Fobjs[frg].dm_init) - else: + else: nprocs = int(self.nproc/self.ompnum) pool_ = Pool(nprocs) os.system('export OMP_NUM_THREADS='+str(self.ompnum)) @@ -444,37 +444,37 @@ def initialize(self, eri_,compute_hf, restart=False): pool_.close() for frg in range(self.Nfrag): self.Fobjs[frg]._mo_coeffs = mo_coeffs[frg] - + for frg in range(self.Nfrag): self.Fobjs[frg].dm0 = numpy.dot( self.Fobjs[frg]._mo_coeffs[:,:self.Fobjs[frg].nsocc], self.Fobjs[frg]._mo_coeffs[:,:self.Fobjs[frg].nsocc].conj().T) *2. - + # energy - if compute_hf: + if compute_hf: eh1, ecoul, ef = self.Fobjs[frg].energy_hf(return_e1=True) E_hf += self.Fobjs[frg].ebe_hf - + print(flush=True) if not restart: file_eri.close() - + if compute_hf: E_hf /= self.unitcell_nkpt - hf_err = self.hf_etot-(E_hf+self.enuc+self.E_core) - + hf_err = self.hf_etot-(E_hf+self.enuc+self.E_core) + self.ebe_hf = E_hf+self.enuc+self.E_core-self.ek print('HF-in-HF error : {:>.4e} Ha'. format(hf_err), flush=True) - + if abs(hf_err)>1.e-5: print('WARNING!!! Large HF-in-HF energy error') - + couti = 0 for fobj in self.Fobjs: fobj.udim = couti couti = fobj.set_udim(couti) - + def oneshot(self, solver='MP2', nproc=1, ompnum=4, calc_frag_energy=False, clean_eri=False): """ Perform a one-shot bootstrap embedding calculation. @@ -579,9 +579,9 @@ def read_heff(self, heff_file='bepotfile.h5'): for fobj in self.Fobjs: fobj.heff = filepot.get(fobj.dname) filepot.close() - - - + + + def initialize_pot(Nfrag, edge_idx): """ Initialize the potential array for bootstrap embedding. @@ -605,7 +605,7 @@ def initialize_pot(Nfrag, edge_idx): Initialized potential array with zeros. """ pot_=[] - + if not len(edge_idx) == 0: for I in range(Nfrag): for i in edge_idx[I]: @@ -614,7 +614,7 @@ def initialize_pot(Nfrag, edge_idx): if j>k: continue pot_.append(0.0) - + pot_.append(0.) return pot_ @@ -661,5 +661,5 @@ def parallel_scf_wrapper(dname, nao, nocc, h1, dm_init, eri_file): from .helper import get_eri, get_scfObj eri = get_eri(dname, nao, eri_file=eri_file) mf_ = get_scfObj(h1, eri, nocc, dm_init) - + return mf_.mo_coeff diff --git a/kbe/pfrag.py b/kbe/pfrag.py index d648d3f4..145b90f7 100644 --- a/kbe/pfrag.py +++ b/kbe/pfrag.py @@ -19,7 +19,7 @@ def __init__(self, fsites, ifrag, edge=None, center=None, edge_idx=None, center_idx=None, efac=None, eri_file='eri_file.h5',unitcell_nkpt=1, ewald_ek=None, centerf_idx=None, unitcell=1): - """Constructor function for `Frags` class. + """Constructor function for `Frags` class. Parameters ---------- @@ -42,7 +42,7 @@ def __init__(self, fsites, ifrag, edge=None, center=None, centerf_idx : list, optional indices of the center site atoms in the fragment, by default None """ - + self.fsites = fsites self.unitcell=unitcell self.unitcell_nkpt=unitcell_nkpt @@ -53,19 +53,19 @@ def __init__(self, fsites, ifrag, edge=None, center=None, self.ifrag = ifrag self.dname = 'f'+str(ifrag) self.nao = None - self.mo_coeffs = None - self._mo_coeffs = None + self.mo_coeffs = None + self._mo_coeffs = None self.nsocc = None self._mf = None self._mc = None - + # CCSD self.t1 = None self.t2 = None - + self.heff = None self.edge = edge - self.center = center + self.center = center self.edge_idx = edge_idx self.center_idx = center_idx self.centerf_idx = centerf_idx @@ -91,7 +91,7 @@ def __init__(self, fsites, ifrag, edge=None, center=None, self.ebe_hf0 = 0. self.rdm1_lo_k = None - def sd(self, lao, lmo, nocc, + def sd(self, lao, lmo, nocc, frag_type='autogen', cell=None, kpts = None, kmesh=None, h1=None): """ @@ -112,10 +112,10 @@ def sd(self, lao, lmo, nocc, kmesh : list of int Number of k-points in each lattice vector direction """ - + from .misc import get_phase - - nk, nao, nlo = lao.shape + + nk, nao, nlo = lao.shape rdm1_lo_k = numpy.zeros((nk, nlo, nlo), dtype=numpy.result_type(lmo, lmo)) for k in range(nk): @@ -124,47 +124,47 @@ def sd(self, lao, lmo, nocc, phase = get_phase(cell, kpts, kmesh) supcell_rdm = numpy.einsum('Rk,kuv,Sk->RuSv', phase, rdm1_lo_k, phase.conj()) supcell_rdm = supcell_rdm.reshape(nk*nlo, nk*nlo) - + if numpy.abs(supcell_rdm.imag).max() < 1.e-6: supcell_rdm = supcell_rdm.real else: print('Imaginary density in Full SD', numpy.abs(supcell_rdm.imag).max()) sys.exit() - - Sites = [i+(nlo*0) for i in self.fsites] + + Sites = [i+(nlo*0) for i in self.fsites] if not frag_type == 'autogen': Sites.sort() - + TA_R = schmidt_decomp_svd(supcell_rdm, Sites) teo = TA_R.shape[-1] TA_R = TA_R.reshape(nk, nlo, teo) - + phase1 = get_phase1(cell, kpts, kmesh) - TA_k = numpy.einsum('Rim, Rk -> kim', TA_R, phase1) + TA_k = numpy.einsum('Rim, Rk -> kim', TA_R, phase1) self.TA_lo_eo = TA_k - + TA_ao_eo_k = numpy.zeros((nk, nao, teo), - dtype=numpy.result_type(lao.dtype, TA_k.dtype)) + dtype=numpy.result_type(lao.dtype, TA_k.dtype)) for k in range(nk): TA_ao_eo_k[k] = numpy.dot(lao[k], TA_k[k]) - + self.TA = TA_ao_eo_k self.nao = TA_ao_eo_k.shape[-1] - # useful for debugging -- + # useful for debugging -- rdm1_eo = numpy.zeros((teo, teo), dtype=numpy.complex128) for k in range(nk): rdm1_eo += functools.reduce(numpy.dot, (TA_k[k].conj().T, rdm1_lo_k[k], TA_k[k])) rdm1_eo /= float(nk) - + h1_eo = numpy.zeros((teo, teo), dtype=numpy.complex128) for k in range(nk): h1_eo += functools.reduce(numpy.dot, (self.TA[k].conj().T, h1[k], self.TA[k])) - h1_eo /= float(nk) + h1_eo /= float(nk) e1 = 2.0 *numpy.einsum("ij,ij->i", h1_eo[:self.nfsites], rdm1_eo[:self.nfsites]) e_h1 = 0. @@ -180,21 +180,21 @@ def cons_h1(self, h1): h1 : numpy.ndarray One-electron Hamiltonian matrix. """ - + nk, nao, teo = self.TA.shape h1_eo = numpy.zeros((teo, teo), dtype=numpy.complex128) - for k in range(nk): + for k in range(nk): h1_eo += functools.reduce(numpy.dot, (self.TA[k].conj().T, h1[k], - self.TA[k])) + self.TA[k])) h1_eo /= float(nk) - + if numpy.abs(h1_eo.imag).max() < 1.e-7: self.h1 = h1_eo.real else: print('Imaginary Hcore ', numpy.abs(h1_eo.imag).max()) - sys.exit() - + sys.exit() + def cons_fock(self, hf_veff, S, dm, eri_=None): """ Construct the Fock matrix for the fragment. @@ -210,18 +210,18 @@ def cons_fock(self, hf_veff, S, dm, eri_=None): eri_ : numpy.ndarray, optional Electron repulsion integrals, by default None. """ - + if eri_ is None: eri_ = get_eri(self.dname, self.TA.shape[1], ignore_symm=True, eri_file=self.eri_file) - + veff0, veff_ = get_veff(eri_, dm, S, self.TA, hf_veff, return_veff0=True) - if numpy.abs(veff_.imag).max() < 1.e-6: + if numpy.abs(veff_.imag).max() < 1.e-6: self.veff = veff_.real self.veff0 = veff0.real else: print('Imaginary Veff ', numpy.abs(veff_.imag).max()) sys.exit() - self.fock = self.h1 + veff_.real + self.fock = self.h1 + veff_.real def get_nsocc(self, S, C, nocc,ncore=0): """ @@ -244,7 +244,7 @@ def get_nsocc(self, S, C, nocc,ncore=0): Projected density matrix. """ import scipy.linalg - + nk, nao, neo = self.TA.shape dm_ = numpy.zeros((nk, nao, nao), dtype=numpy.result_type(C,C)) for k in range(nk): @@ -260,14 +260,14 @@ def get_nsocc(self, S, C, nocc,ncore=0): P_ = P_.real else: print('Imaginary density in get_nsocc ', numpy.abs(P_.imag).max()) - sys.exit() + sys.exit() nsocc_ = numpy.trace(P_) nsocc = int(numpy.round(nsocc_.real)/2) - + self.nsocc = nsocc return P_ - - + + def scf(self, heff=None, fs=False, eri=None, pert_h=False,pert_list=None, save_chkfile=False, dm0 = None): @@ -289,14 +289,14 @@ def scf(self, heff=None, fs=False, eri=None, if self._mf is not None: self._mf = None if self._mc is not None: self._mc = None if heff is None: heff = self.heff - + if eri is None: eri = get_eri(self.dname, self.nao, eri_file=self.eri_file) if dm0 is None: dm0 = numpy.dot( self._mo_coeffs[:,:self.nsocc], self._mo_coeffs[:,:self.nsocc].conj().T) *2. - + mf_ = get_scfObj(self.fock + heff, eri, self.nsocc, dm0 = dm0, fname = self.dname, @@ -305,16 +305,16 @@ def scf(self, heff=None, fs=False, eri=None, if pert_h: return mf_ - + if not fs: self._mf = mf_ self.mo_coeffs = mf_.mo_coeff.copy() else: self._mo_coeffs = mf_.mo_coeff.copy() - dm0 = mf_.make_rdm1() + dm0 = mf_.make_rdm1() mf_= None - + def update_heff(self,u, cout = None, return_heff=False, be_iter=None, no_chempot=False, @@ -324,7 +324,7 @@ def update_heff(self,u, cout = None, return_heff=False, Update the effective Hamiltonian for the fragment. """ import h5py - + heff_ = numpy.zeros_like(self.h1) if cout is None: @@ -332,10 +332,10 @@ def update_heff(self,u, cout = None, return_heff=False, else: cout = cout if not no_chempot: - for i,fi in enumerate(self.fsites): - if not any(i in sublist for sublist in self.edge_idx): + for i,fi in enumerate(self.fsites): + if not any(i in sublist for sublist in self.edge_idx): heff_[i,i] -= u[-1] - + if only_chem: self.heff = heff_ if return_heff: @@ -344,18 +344,18 @@ def update_heff(self,u, cout = None, return_heff=False, else: return(cout, heff_) return cout - + for idx,i in enumerate(self.edge_idx): for j in range(len(i)): for k in range(len(i)): if j>k : - continue - + continue + heff_[i[j], i[k]] = u[cout] heff_[i[k], i[j]] = u[cout] - + cout += 1 - + self.heff = heff_ if return_heff: if cout is None: @@ -376,7 +376,7 @@ def set_udim(self, cout): def energy_hf(self, rdm_hf=None, mo_coeffs = None, eri=None, return_e1=False, unrestricted = False): if mo_coeffs is None: mo_coeffs = self._mo_coeffs - + if rdm_hf is None: rdm_hf = numpy.dot(mo_coeffs[:,:self.nsocc], mo_coeffs[:,:self.nsocc].conj().T) @@ -385,10 +385,10 @@ def energy_hf(self, rdm_hf=None, mo_coeffs = None, eri=None, return_e1=False, un e1 = unrestricted*numpy.einsum("ij,ij->i", self.h1[:self.nfsites], rdm_hf[:self.nfsites]) - + ec = 0.5 * unrestricted * numpy.einsum("ij,ij->i",self.veff[:self.nfsites], rdm_hf[:self.nfsites]) - + if self.TA.ndim == 3: jmax = self.TA[0].shape[1] else: @@ -399,18 +399,18 @@ def energy_hf(self, rdm_hf=None, mo_coeffs = None, eri=None, return_e1=False, un r.close() - + e2 = numpy.zeros_like(e1) for i in range(self.nfsites): for j in range(jmax): ij = i*(i+1)//2+j if i > j else j*(j+1)//2+i Gij = (2.*rdm_hf[i,j]*rdm_hf - - numpy.outer(rdm_hf[i], rdm_hf[j]))[:jmax,:jmax] + numpy.outer(rdm_hf[i], rdm_hf[j]))[:jmax,:jmax] Gij[numpy.diag_indices(jmax)] *= 0.5 - Gij += Gij.T + Gij += Gij.T e2[i] += 0.5 * unrestricted * Gij[numpy.tril_indices(jmax)] @ eri[ij] - e_ = e1+e2+ec + e_ = e1+e2+ec etmp = 0. e1_ = 0. e2_ = 0. @@ -420,15 +420,15 @@ def energy_hf(self, rdm_hf=None, mo_coeffs = None, eri=None, return_e1=False, un e1_ += self.efac[0]*e1[i] e2_ += self.efac[0]*e2[i] ec_ += self.efac[0]*ec[i] - + self.ebe_hf = etmp if return_e1: e_h1 = 0. e_coul = 0. for i in self.efac[1]: - + e_h1 += self.efac[0]*e1[i] e_coul += self.efac[0]*(e2[i]+ec[i]) return(e_h1,e_coul, e1+e2+ec) - - return e1+e2+ec + + return e1+e2+ec diff --git a/kbe/solver.py b/kbe/solver.py index b03ed2e3..5eecce61 100644 --- a/kbe/solver.py +++ b/kbe/solver.py @@ -21,25 +21,25 @@ def schmidt_decomp_svd(rdm, Frag_sites): ------- numpy.ndarray Transformation matrix (TA) including both fragment and entangled bath orbitals. - """ + """ import scipy.linalg import functools - + thres = 1.0e-10 - Tot_sites = rdm.shape[0] - + Tot_sites = rdm.shape[0] + Fragsites = [i if i>=0 else Tot_sites+i for i in Frag_sites] - + Env_sites1 = numpy.array([i for i in range(Tot_sites) if not i in Fragsites]) nfs = len(Frag_sites) - - Denv = rdm[Env_sites1][:, Fragsites] + + Denv = rdm[Env_sites1][:, Fragsites] U, sigma, V = scipy.linalg.svd(Denv, full_matrices=False, lapack_driver='gesvd') - nbath = ( sigma >= thres).sum() + nbath = ( sigma >= thres).sum() TA = numpy.zeros((Tot_sites, nfs + nbath), dtype=numpy.complex128) TA[Fragsites, :nfs] = numpy.eye(nfs) TA[Env_sites1, nfs:] = U[:,:nbath] - + return TA diff --git a/molbe/_opt.py b/molbe/_opt.py index a1b3c0c1..8ba0c3f4 100644 --- a/molbe/_opt.py +++ b/molbe/_opt.py @@ -23,7 +23,7 @@ class BEOPT: enuc : float Nuclear component of the energy. solver : str - High-level solver in bootstrap embedding. 'MP2', 'CCSD', 'FCI' are supported. Selected CI versions, + High-level solver in bootstrap embedding. 'MP2', 'CCSD', 'FCI' are supported. Selected CI versions, 'HCI', 'SHCI', & 'SCI' are also supported. Defaults to 'MP2' only_chem : bool Whether to perform chemical potential optimization only. Refer to bootstrap embedding literatures. @@ -39,15 +39,15 @@ class BEOPT: ebe_hf : float Hartree-Fock energy. Defaults to 0.0 """ - + def __init__(self, pot, Fobjs, Nocc, enuc,solver='MP2', ecore=0., nproc=1,ompnum=4, only_chem=False, hf_veff = None, hci_pt=False,hci_cutoff=0.001, ci_coeff_cutoff = None, select_cutoff=None, max_space=500, conv_tol = 1.e-6,relax_density = False, ebe_hf =0., scratch_dir=None, **solver_kwargs): - - # Initialize class attributes + + # Initialize class attributes self.ebe_hf=ebe_hf self.hf_veff = hf_veff self.pot = pot @@ -83,7 +83,7 @@ def objfunc(self, xk): ---------- xk : list Current potentials in the BE optimization. - + Returns ------- list @@ -94,7 +94,7 @@ def objfunc(self, xk): if self.nproc == 1: err_, errvec_,ebe_ = be_func(xk, self.Fobjs, self.Nocc, self.solver, self.enuc, eeval=True, return_vec=True, hf_veff = self.hf_veff, - only_chem=self.only_chem, + only_chem=self.only_chem, hci_cutoff=self.hci_cutoff, nproc=self.ompnum, relax_density=self.relax_density, ci_coeff_cutoff = self.ci_coeff_cutoff, @@ -109,14 +109,14 @@ def objfunc(self, xk): hci_cutoff=self.hci_cutoff,relax_density=self.relax_density, ci_coeff_cutoff = self.ci_coeff_cutoff, select_cutoff = self.select_cutoff, - ecore=self.ecore, ebe_hf=self.ebe_hf, be_iter=self.iter, + ecore=self.ecore, ebe_hf=self.ebe_hf, be_iter=self.iter, scratch_dir=self.scratch_dir, **self.solver_kwargs) - + # Update error and BE energy self.err = err_ self.Ebe = ebe_ return errvec_ - + def optimize(self, method, J0 = None, trust_region=False): """Main kernel to perform BE optimization @@ -132,7 +132,7 @@ def optimize(self, method, J0 = None, trust_region=False): """ from molbe.external.optqn import FrankQN import sys - + print('-----------------------------------------------------', flush=True) print(' Starting BE optimization ', flush=True) @@ -142,13 +142,13 @@ def optimize(self, method, J0 = None, trust_region=False): print('-----------------------------------------------------', flush=True) print(flush=True) - + if method=='QN': - - print('-- In iter ',self.iter, flush=True) + + print('-- In iter ',self.iter, flush=True) # Initial step - f0 = self.objfunc(self.pot) + f0 = self.objfunc(self.pot) print('Error in density matching : {:>2.4e}'.format(self.err), flush=True) print(flush=True) @@ -157,7 +157,7 @@ def optimize(self, method, J0 = None, trust_region=False): optQN = FrankQN(self.objfunc, numpy.array(self.pot), f0, J0, max_space=self.max_space) - + if self.err < self.conv_tol: print(flush=True) print('CONVERGED w/o Optimization Steps',flush=True) @@ -178,9 +178,9 @@ def optimize(self, method, J0 = None, trust_region=False): else: print('This optimization method for BE is not supported') sys.exit() - - - + + + def optimize(self, solver='MP2',method='QN', only_chem=False, conv_tol = 1.e-6, relax_density=False, use_cumulant=True, @@ -244,9 +244,9 @@ def optimize(self, solver='MP2',method='QN', if only_chem: J0 = [[0.]] J0 = self.get_be_error_jacobian(jac_solver='HF') - J0 = [[J0[-1,-1]]] + J0 = [[J0[-1,-1]]] else: - J0 = self.get_be_error_jacobian(jac_solver='HF') + J0 = self.get_be_error_jacobian(jac_solver='HF') # Perform the optimization be_.optimize(method, J0=J0, trust_region=trust_region) @@ -256,4 +256,4 @@ def optimize(self, solver='MP2',method='QN', else: print('This optimization method for BE is not supported') sys.exit() - + diff --git a/molbe/autofrag.py b/molbe/autofrag.py index 5ce966c2..28c9f815 100644 --- a/molbe/autofrag.py +++ b/molbe/autofrag.py @@ -4,7 +4,7 @@ import numpy from .helper import get_core -def autogen(mol, frozen_core=True, be_type='be2', +def autogen(mol, frozen_core=True, be_type='be2', write_geom=False, valence_basis = None, valence_only = False, print_frags=True): @@ -13,8 +13,8 @@ def autogen(mol, frozen_core=True, be_type='be2', Partitions a molecule into overlapping fragments as defined in BE atom-based fragmentations. It automatically detects branched chemical chains and ring systems and partitions accordingly. - For efficiency, it only checks two atoms for connectivity (chemical bond) if they are within 3.5 Angstrom. - This value is hardcoded as normdist. Two atoms are defined as bonded if they are within 1.8 Angstrom (1.2 for Hydrogen atom). + For efficiency, it only checks two atoms for connectivity (chemical bond) if they are within 3.5 Angstrom. + This value is hardcoded as normdist. Two atoms are defined as bonded if they are within 1.8 Angstrom (1.2 for Hydrogen atom). This is also hardcoded as bond & hbond. Parameters @@ -70,8 +70,8 @@ def autogen(mol, frozen_core=True, be_type='be2', cell = mol.copy() cell.basis = valence_basis cell.build() - - ncore, no_core_idx, core_list = get_core(cell) + + ncore, no_core_idx, core_list = get_core(cell) coord = cell.atom_coords() ang2bohr = 1.88973 normdist = 3.5 * ang2bohr @@ -82,7 +82,7 @@ def autogen(mol, frozen_core=True, be_type='be2', normlist = [] for i in coord: normlist.append(numpy.linalg.norm(i)) - Frag = [] + Frag = [] pedge = [] cen = [] @@ -95,33 +95,33 @@ def autogen(mol, frozen_core=True, be_type='be2', open_frag = [] open_frag_cen = [] - - # Assumes that there can be only 5 member connected system + + # Assumes that there can be only 5 member connected system for idx, i in enumerate(normlist): if cell.atom_pure_symbol(idx) == 'H' and not hchain: continue - + tmplist = normlist - i tmplist = list(tmplist) - + clist = [] cout = 0 for jdx,j in enumerate(tmplist): if not idx==jdx and (not cell.atom_pure_symbol(jdx) == 'H' or hchain): - if abs(j)< normdist: + if abs(j)< normdist: clist.append(jdx) pedg = [] - flist = [] + flist = [] flist.append(idx) if not be_type == 'be1': - for jdx in clist: - dist = numpy.linalg.norm(coord[idx] - coord[jdx]) + for jdx in clist: + dist = numpy.linalg.norm(coord[idx] - coord[jdx]) if dist <= bond: flist.append(jdx) pedg.append(jdx) if be_type=='be3' or be_type == 'be4': - + for kdx in clist: if not kdx == jdx: dist = numpy.linalg.norm(coord[jdx] - coord[kdx]) @@ -129,31 +129,31 @@ def autogen(mol, frozen_core=True, be_type='be2', if not kdx in pedg: flist.append(kdx) pedg.append(kdx) - if be_type=='be4': + if be_type=='be4': for ldx, l in enumerate(coord): if ldx==kdx or ldx==jdx or\ (cell.atom_pure_symbol(ldx) == 'H' and not hchain)\ - or ldx in pedg: + or ldx in pedg: continue dist = numpy.linalg.norm(coord[kdx] - l) if dist <= bond: flist.append(ldx) - pedg.append(ldx) + pedg.append(ldx) # Update fragment and edge lists based on current partitioning for pidx, frag_ in enumerate(Frag): - if set(flist).issubset(frag_): + if set(flist).issubset(frag_): open_frag.append(pidx) open_frag_cen.append(idx) break - elif set(frag_).issubset(flist): - open_frag = [ oidx-1 if oidx > pidx else oidx for oidx in open_frag] + elif set(frag_).issubset(flist): + open_frag = [ oidx-1 if oidx > pidx else oidx for oidx in open_frag] open_frag.append(len(Frag)-1) open_frag_cen.append(cen[pidx]) del cen[pidx] del Frag[pidx] del pedge[pidx] - else: + else: Frag.append(flist) pedge.append(pedg) cen.append(idx) @@ -164,15 +164,15 @@ def autogen(mol, frozen_core=True, be_type='be2', hlist = [[] for i in coord] if not hchain: for idx, i in enumerate(normlist): - if cell.atom_pure_symbol(idx) == 'H': + if cell.atom_pure_symbol(idx) == 'H': tmplist = normlist - i - tmplist = list(tmplist) + tmplist = list(tmplist) clist = [] for jdx,j in enumerate(tmplist): if not idx==jdx and not cell.atom_pure_symbol(jdx) == 'H': if abs(j)< normdist: - clist.append(jdx) - for jdx in clist: + clist.append(jdx) + for jdx in clist: dist = numpy.linalg.norm(coord[idx] - coord[jdx]) if dist <= hbond: hlist[jdx].append(idx) @@ -184,7 +184,7 @@ def autogen(mol, frozen_core=True, be_type='be2', print('--------------------------',flush=True) print('Fragment | Center | Edges ',flush=True) print('--------------------------',flush=True) - + for idx,i in enumerate(Frag): print(' {:>4} | {:>5} |'.format(idx, cell.atom_pure_symbol(cen[idx])+str(cen[idx]+1)),end=' ', flush=True) for j in hlist[cen[idx]]: @@ -205,7 +205,7 @@ def autogen(mol, frozen_core=True, be_type='be2', w = open('fragments.xyz','w') for idx,i in enumerate(Frag): w.write(str(len(i)+len(hlist[cen[idx]])+len(hlist[j]))+'\n') - w.write('Fragment - '+str(idx)+'\n') + w.write('Fragment - '+str(idx)+'\n') for j in hlist[cen[idx]]: w.write(' {:>3} {:>10.7f} {:>10.7f} {:>10.7f} \n'.format(cell.atom_pure_symbol(j), coord[j][0]/ang2bohr, @@ -225,15 +225,15 @@ def autogen(mol, frozen_core=True, be_type='be2', # Prepare for PAO basis if requested pao = bool(valence_basis and not valence_only) - + if pao: cell2 = cell.copy() cell2.basis = valence_basis cell2.build() - bas2list = cell2.aoslice_by_atom() + bas2list = cell2.aoslice_by_atom() nbas2 = [0 for i in range(cell.natm)] - + baslist = cell.aoslice_by_atom() sites__ = [[] for i in coord] coreshift = 0 @@ -241,8 +241,8 @@ def autogen(mol, frozen_core=True, be_type='be2', # Process each atom to determine core and valence basis sets for adx in range(cell.natm): - if hchain: - bas = baslist[adx] + if hchain: + bas = baslist[adx] start_ = bas[2] stop_ = bas[3] if pao: @@ -251,27 +251,27 @@ def autogen(mol, frozen_core=True, be_type='be2', b1list = [i for i in range(start_, stop_)] sites__[adx] = b1list continue - + if not cell.atom_pure_symbol(adx) == 'H' and not hchain: - bas = baslist[adx] + bas = baslist[adx] start_ = bas[2] stop_ = bas[3] if pao: bas2 = bas2list[adx] nbas2[adx] += bas2[3] - bas2[2] - + if frozen_core: - start_ -= coreshift + start_ -= coreshift ncore_ = core_list[adx] stop_ -= coreshift+ncore_ - if pao: nbas2[adx] -= ncore_ + if pao: nbas2[adx] -= ncore_ coreshift += ncore_ b1list = [i for i in range(start_, stop_)] sites__[adx] = b1list else: hshift[adx] = coreshift - + hsites = [[] for i in coord] nbas2H = [0 for i in coord] for hdx, h in enumerate(hlist): @@ -283,14 +283,14 @@ def autogen(mol, frozen_core=True, be_type='be2', if pao: bas2H = bas2list[hidx] nbas2H[hdx] += bas2H[3] - bas2H[2] - - if frozen_core: + + if frozen_core: startH -= hshift[hidx] stopH -= hshift[hidx] - + b1list = [i for i in range(startH, stopH)] hsites[hdx].extend(b1list) - + fsites = [] edgsites = [] edge_idx = [] @@ -300,11 +300,11 @@ def autogen(mol, frozen_core=True, be_type='be2', # Create fragments and edges based on partitioning for idx, i in enumerate(Frag): ftmp = [] - ftmpe = [] + ftmpe = [] indix = 0 edind = [] - edg = [] - + edg = [] + frglist = sites__[cen[idx]].copy() frglist.extend(hsites[cen[idx]]) @@ -317,31 +317,31 @@ def autogen(mol, frozen_core=True, be_type='be2', frglist.extend(hsites[open_frag_cen[pidx__]]) ls += len(sites__[open_frag_cen[pidx__]]) +\ len(hsites[open_frag_cen[pidx__]]) - + ftmp.extend(frglist) if not pao: ls_ = len(sites__[cen[idx]])+len(hsites[cen[idx]]) - centerf_idx.append([pq for pq in range(indix,indix+ls_)]) + centerf_idx.append([pq for pq in range(indix,indix+ls_)]) else: cntlist = sites__[cen[idx]].copy()[:nbas2[cen[idx]]] cntlist.extend(hsites[cen[idx]][:nbas2H[cen[idx]]]) - ind__ = [ indix+frglist.index(pq) for pq in cntlist] + ind__ = [ indix+frglist.index(pq) for pq in cntlist] centerf_idx.append(ind__) indix += ls if not be_type=='be1': for jdx in pedge[idx]: - + if idx in open_frag: - + if jdx == open_frag_cen[open_frag.index(idx)]: continue if jdx in open_frag_cen: - continue + continue edg.append(jdx) - frglist = sites__[jdx].copy() + frglist = sites__[jdx].copy() frglist.extend(hsites[jdx]) - + ftmp.extend(frglist) ls = len(sites__[jdx]) + len(hsites[jdx]) if not pao: @@ -352,11 +352,11 @@ def autogen(mol, frozen_core=True, be_type='be2', else: edglist = sites__[jdx][:nbas2[jdx]].copy() edglist.extend(hsites[jdx][:nbas2H[jdx]]) - + ftmpe.append(edglist) - ind__ = [ indix+frglist.index(pq) for pq in edglist] - edind.append(ind__) - indix += ls + ind__ = [ indix+frglist.index(pq) for pq in edglist] + edind.append(ind__) + indix += ls edge.append(edg) edgsites.append(ftmpe) edge_idx.append(edind) @@ -373,24 +373,24 @@ def autogen(mol, frozen_core=True, be_type='be2', print(' This is more complicated than I can handle') sys.exit() center.append(cen_) - + Nfrag = len(fsites) add_centers=[[] for x in range(Nfrag)] # additional centers for mixed-basis ebe_weight=[] # Compute weights for each fragment - for ix, i in enumerate(fsites): + for ix, i in enumerate(fsites): tmp_ = [i.index(pq) for pq in sites__[cen[ix]]] - tmp_.extend([i.index(pq) for pq in hsites[cen[ix]]]) + tmp_.extend([i.index(pq) for pq in hsites[cen[ix]]]) if ix in open_frag: for pidx__,pid__ in enumerate(open_frag): if ix == pid__: add_centers[pid__].append(open_frag_cen[pidx__]) tmp_.extend([i.index(pq) for pq in sites__[open_frag_cen[pidx__]]]) - tmp_.extend([i.index(pq) for pq in hsites[open_frag_cen[pidx__]]]) + tmp_.extend([i.index(pq) for pq in hsites[open_frag_cen[pidx__]]]) ebe_weight.append([1.0, tmp_]) - + center_idx = [] if not be_type=='be1': for i in range(Nfrag): @@ -411,8 +411,8 @@ def autogen(mol, frozen_core=True, be_type='be2', idx.append([fsites[j].index(k) for k in cntlist]) jdx_continue = True break - - if jdx_continue: continue + + if jdx_continue: continue if not pao: cntlist = sites__[cen[j]].copy() cntlist.extend(hsites[cen[j]]) @@ -421,7 +421,7 @@ def autogen(mol, frozen_core=True, be_type='be2', cntlist = sites__[cen[j]].copy()[:nbas2[cen[j]]] cntlist.extend(hsites[cen[j]][:nbas2H[cen[j]]]) idx.append([fsites[j].index(k) for k in cntlist]) - + center_idx.append(idx) return(fsites, edgsites, center, edge_idx, center_idx, centerf_idx, ebe_weight, Frag, cen, hlist, add_centers) diff --git a/molbe/be_parallel.py b/molbe/be_parallel.py index 5653f8e9..dae60dc5 100644 --- a/molbe/be_parallel.py +++ b/molbe/be_parallel.py @@ -71,15 +71,15 @@ def run_solver(h1, dm0, dname, nao, nocc, nfsites, # Initialize SCF object mf_ = get_scfObj(h1, eri, nocc, dm0=dm0) rdm_return = False - + if relax_density: rdm_return = True # Select solver - if solver=='MP2': + if solver=='MP2': mc_ = solve_mp2(mf_, mo_energy=mf_.mo_energy) rdm1_tmp = mc_.make_rdm1() - + elif solver=='CCSD': if not rdm_return: t1, t2 = solve_ccsd(mf_, @@ -94,7 +94,7 @@ def run_solver(h1, dm0, dname, nao, nocc, nfsites, relax=True) elif solver=='FCI': from pyscf import fci - + mc_ = fci.FCI(mf_, mf_.mo_coeff) efci, civec = mc_.kernel() rdm1_tmp = mc_.make_rdm1(civec, mc_.norb, mc_.nelec) @@ -103,25 +103,25 @@ def run_solver(h1, dm0, dname, nao, nocc, nfsites, nao, nmo = mf_.mo_coeff.shape eri = ao2mo.kernel(mf_._eri, mf_.mo_coeff, aosym='s4', - compact=False).reshape(4*((nmo),)) + compact=False).reshape(4*((nmo),)) ci_ = hci.SCI(mf_.mol) if select_cutoff is None and ci_coeff_cutoff is None: select_cutoff = hci_cutoff ci_coeff_cutoff = hci_cutoff elif select_cutoff is None or ci_coeff_cutoff is None: sys.exit() - + ci_.select_cutoff = select_cutoff ci_.ci_coeff_cutoff = ci_coeff_cutoff - - nelec = (nocc, nocc) + + nelec = (nocc, nocc) h1_ = functools.reduce(numpy.dot, (mf_.mo_coeff.T, h1, mf_.mo_coeff)) eci, civec = ci_.kernel(h1_, eri, nmo, nelec) civec = numpy.asarray(civec) - + (rdm1a_, rdm1b_), (rdm2aa, rdm2ab, rdm2bb) = ci_.make_rdm12s(civec, nmo, nelec) rdm1_tmp = rdm1a_ + rdm1b_ - rdm2s = rdm2aa + rdm2ab + rdm2ab.transpose(2,3,0,1) + rdm2bb + rdm2s = rdm2aa + rdm2ab + rdm2ab.transpose(2,3,0,1) + rdm2bb elif solver=='SHCI': from pyscf.shciscf import shci @@ -140,7 +140,7 @@ def run_solver(h1, dm0, dname, nao, nocc, nfsites, mch.fcisolver.restart=True mch.mc1step() rdm1_tmp, rdm2s = mch.fcisolver.make_rdm12(0, nmo, nelec) - + elif solver == 'SCI': from pyscf import cornell_shci from pyscf import ao2mo, mcscf @@ -160,23 +160,23 @@ def run_solver(h1, dm0, dname, nao, nocc, nfsites, ci.config['get_2rdm_csv'] = True ci.kernel(h1, eri, nmo, nelec) rdm1_tmp, rdm2s = ci.make_rdm12(0,nmo,nelec) - + else: print('Solver not implemented',flush=True) - print('exiting',flush=True) + print('exiting',flush=True) sys.exit() # Compute RDM1 rdm1 = functools.reduce(numpy.dot, (mf_.mo_coeff, rdm1_tmp, - mf_.mo_coeff.T))*0.5 + mf_.mo_coeff.T))*0.5 if eeval: if solver =='CCSD' and not rdm_return: with_dm1 = True if use_cumulant: with_dm1 = False rdm2s = make_rdm2_urlx(t1, t2, with_dm1 = with_dm1) - + elif solver == 'MP2': rdm2s = mc_.make_rdm2() elif solver == 'FCI': @@ -199,9 +199,9 @@ def run_solver(h1, dm0, dname, nao, nocc, nfsites, if return_rdm_ao: return(e_f, mf_.mo_coeff, rdm1, rdm2s, rdm1_tmp) - + return (e_f, mf_.mo_coeff, rdm1, rdm2s) - + def run_solver_u(fobj_a, fobj_b, solver, enuc, hf_veff, frag_energy=True, relax_density=False, frozen=False, @@ -212,7 +212,7 @@ def run_solver_u(fobj_a, fobj_b, solver, enuc, hf_veff, Parameters ---------- fobj_a : - Alpha spin molbe.pfrag.Frags object + Alpha spin molbe.pfrag.Frags object fobj_b : Beta spin molbe.pfrag.Frags object solver : str @@ -228,7 +228,7 @@ def run_solver_u(fobj_a, fobj_b, solver, enuc, hf_veff, frozen : bool, optional If True, uses frozen core, defaults to False eri_file : str, optional - Filename for the electron repulsion integrals. Default is 'eri_file.h5'. + Filename for the electron repulsion integrals. Default is 'eri_file.h5'. use_cumulant : bool, optional If True, uses the cumulant approximation for RDM2. Default is True. ereturn : bool, optional @@ -236,7 +236,7 @@ def run_solver_u(fobj_a, fobj_b, solver, enuc, hf_veff, Returns ------- - + As implemented, only returns the UCCSD fragment energy """ print("obj type", type(fobj_a)) @@ -285,7 +285,7 @@ def run_solver_u(fobj_a, fobj_b, solver, enuc, hf_veff, else: raise NotImplementedError("RDM Return not Implemented") - fobj_a.__rdm2 = rdm2s[0].copy() + fobj_a.__rdm2 = rdm2s[0].copy() fobj_b.__rdm2 = rdm2s[1].copy() # Calculate energy on a per-fragment basis @@ -313,24 +313,24 @@ def run_solver_u(fobj_a, fobj_b, solver, enuc, hf_veff, else: return NotImplementedError("Energy only calculated on a per-fragment basis") - + def be_func_parallel(pot, Fobjs, Nocc, solver, enuc, hf_veff=None, nproc=1, ompnum=4, only_chem=False,relax_density=False,use_cumulant=True, - eeval=False, ereturn=False, frag_energy=False, + eeval=False, ereturn=False, frag_energy=False, hci_cutoff=0.001, ci_coeff_cutoff = None, select_cutoff=None, return_vec=False, ecore=0., ebe_hf=0., be_iter=None, writeh1=False): """ Embarrassingly Parallel High-Level Computation - Performs high-level bootstrap embedding (BE) computation for each fragment. Computes 1-RDMs - and 2-RDMs for each fragment. It also computes error vectors in BE density match. For selected + Performs high-level bootstrap embedding (BE) computation for each fragment. Computes 1-RDMs + and 2-RDMs for each fragment. It also computes error vectors in BE density match. For selected CI solvers, this function exposes thresholds used in selected CI calculations (hci_cutoff, ci_coeff_cutoff, select_cutoff). Parameters ---------- pot : list of float - Potentials (local & global) that are added to the 1-electron Hamiltonian component. + Potentials (local & global) that are added to the 1-electron Hamiltonian component. The last element in the list is the chemical potential. Fobjs : list of MolBE.fragpart Fragment definitions. @@ -368,7 +368,7 @@ def be_func_parallel(pot, Fobjs, Nocc, solver, enuc, hf_veff=None, Returns ------- float or tuple - Depending on the parameters, returns the error norm or a tuple containing the error norm, + Depending on the parameters, returns the error norm or a tuple containing the error norm, error vector, and the computed energy. """ from multiprocessing import Pool @@ -386,16 +386,16 @@ def be_func_parallel(pot, Fobjs, Nocc, solver, enuc, hf_veff=None, nprocs = int(nproc/ompnum) # Update the effective Hamiltonian with potentials - if not pot is None: + if not pot is None: for fobj in Fobjs: fobj.update_heff(pot, only_chem=only_chem) - + pool_ = Pool(nprocs) - results = [] + results = [] rdms = [] # Run solver in parallel for each fragment - for nf in range(nfrag): + for nf in range(nfrag): h1 = Fobjs[nf].fock + Fobjs[nf].heff dm0 = Fobjs[nf].dm0.copy() dname = Fobjs[nf].dname @@ -412,13 +412,13 @@ def be_func_parallel(pot, Fobjs, Nocc, solver, enuc, hf_veff=None, solver,Fobjs[nf].eri_file, veff0, hci_cutoff, ci_coeff_cutoff,select_cutoff, ompnum, writeh1, True, True, use_cumulant, relax_density, frag_energy]) - + results.append(result) # Collect results [rdms.append(result.get()) for result in results] pool_.close() - + if frag_energy: # Compute and return fragment energy # rdms are the returned energies, not density matrices! @@ -426,50 +426,50 @@ def be_func_parallel(pot, Fobjs, Nocc, solver, enuc, hf_veff=None, e_2 = 0. e_c = 0. for i in range(len(rdms)): - e_1 += rdms[i][0] + e_1 += rdms[i][0] e_2 += rdms[i][1] e_c += rdms[i][2] return (e_1+e_2+e_c, (e_1, e_2, e_c)) - + # Compute total energy e_1 = 0. e_2 = 0. e_c = 0. for idx, fobj in enumerate(Fobjs): - e_1 += rdms[idx][0][0] + e_1 += rdms[idx][0][0] e_2 += rdms[idx][0][1] - e_c += rdms[idx][0][2] + e_c += rdms[idx][0][2] fobj.mo_coeffs = rdms[idx][1] fobj._rdm1 = rdms[idx][2] fobj.__rdm2 = rdms[idx][3] fobj.__rdm1 = rdms[idx][4] - - del rdms + + del rdms ernorm, ervec = solve_error(Fobjs,Nocc, only_chem=only_chem) - + if return_vec: return (ernorm, ervec, [e_1+e_2+e_c, [e_1, e_2, e_c]]) - + if eeval: print('Error in density matching : {:>2.4e}'.format(ernorm), flush=True) - + return ernorm def be_func_parallel_u(pot, Fobjs, solver, enuc, hf_veff=None, nproc=1, ompnum=4, relax_density=False,use_cumulant=True, - eeval=False, ereturn=False, frag_energy=False, + eeval=False, ereturn=False, frag_energy=False, ecore=0., ebe_hf=0., frozen=False): """ Embarrassingly Parallel High-Level Computation - Performs high-level unrestricted bootstrap embedding (UBE) computation for each fragment. Computes 1-RDMs + Performs high-level unrestricted bootstrap embedding (UBE) computation for each fragment. Computes 1-RDMs and 2-RDMs for each fragment to return the energy. As such, this currently is equipped for one-shot U-CCSD BE. Parameters ---------- pot : list of float - Potentials (local & global) that are added to the 1-electron Hamiltonian component. + Potentials (local & global) that are added to the 1-electron Hamiltonian component. The last element in the list is the chemical potential. Should always be 0, as this is still a one-shot only implementation Fobjs : list of tuples of MolBE.fragpart @@ -500,7 +500,7 @@ def be_func_parallel_u(pot, Fobjs, solver, enuc, hf_veff=None, Returns ------- float - Returns the computed energy + Returns the computed energy """ from multiprocessing import Pool import os @@ -538,7 +538,7 @@ def be_func_parallel_u(pot, Fobjs, solver, enuc, hf_veff=None, e_2 = 0. e_c = 0. for i in range(len(energy_list)): - e_1 += energy_list[i][0] + e_1 += energy_list[i][0] e_2 += energy_list[i][1] e_c += energy_list[i][2] return (e_1+e_2+e_c, (e_1, e_2, e_c)) diff --git a/molbe/eri_onthefly.py b/molbe/eri_onthefly.py index 687f76c7..b9b2fd36 100644 --- a/molbe/eri_onthefly.py +++ b/molbe/eri_onthefly.py @@ -40,7 +40,7 @@ def calculate_pqL(aux_range): ints = getints3c(mf.mol._add_suffix('int3c2e'), atm, bas, env, shls_slice, 1, 's1', ao_loc, cintopt, out=None) if be_var.PRINT_LEVEL > 10: print('Finish calculating (μν|P) for range', aux_range, flush=True) return ints - + def block_step_size(nfrag, naux, nao): """ Internal function to calculate the block step size for the 3-center integrals calculation diff --git a/molbe/external/ccsd_rdm.py b/molbe/external/ccsd_rdm.py index 6377e7a5..a0305b4f 100644 --- a/molbe/external/ccsd_rdm.py +++ b/molbe/external/ccsd_rdm.py @@ -33,19 +33,19 @@ def make_rdm2_urlx(t1, t2, with_dm1=True): if with_dm1: dm1 = make_rdm1_ccsd_t1(t1) dm1[numpy.diag_indices(nocc)] -= 2 - + for i in range(nocc): dm2[i,i,:,:] += dm1 * 2 dm2[:,:,i,i] += dm1 * 2 dm2[:,i,i,:] -= dm1 dm2[i,:,:,i] -= dm1.T - + for i in range(nocc): for j in range(nocc): dm2[i,i,j,j] += 4 dm2[i,j,j,i] -= 2 - return dm2 + return dm2 def make_rdm1_uccsd(ucc, relax=False): from pyscf.cc.uccsd_rdm import make_rdm1 diff --git a/molbe/external/lo_helper.py b/molbe/external/lo_helper.py index 6f8faa04..e701f0b4 100644 --- a/molbe/external/lo_helper.py +++ b/molbe/external/lo_helper.py @@ -35,11 +35,11 @@ def get_aoind_by_atom(mol, atomind_by_motif=None): for ia in range(natom)] # if motif info is provided, group lo by motif if atomind_by_motif is None: - + aoind_by_atom = [list(range(*aoshift_by_atom[ia:ia+2])) for ia in range(natom)] else: - + nmotif = len(atomind_by_motif) assert( set([ia for im in range(nmotif) @@ -50,19 +50,19 @@ def get_aoind_by_atom(mol, atomind_by_motif=None): for ia in atomind_by_motif[im]: aoind_by_atom[im] += list( range(*aoshift_by_atom[ia:ia+2])) - + return aoind_by_atom def reorder_by_atom_(Clo, aoind_by_atom, S, thr=0.5): import numpy as np - + natom = len(aoind_by_atom) nlo = Clo.shape[1] X = get_symm_mat_pow(S, 0.5) - + Clo_soao = X @ Clo - + loind_reorder = [] loind_by_atom = [None] * natom loshift = 0 diff --git a/molbe/external/optqn.py b/molbe/external/optqn.py index 6342ca02..56227170 100644 --- a/molbe/external/optqn.py +++ b/molbe/external/optqn.py @@ -20,24 +20,24 @@ def line_search_LF(func, xold, fold, dx, iter_): xk = xold + dx lcout = 0 - - fk = func(xk) + + fk = func(xk) lcout += 1 - + norm_dx = numpy.linalg.norm(dx) norm_fk = numpy.linalg.norm(fk) norm_fold = numpy.linalg.norm(fold) alp = 1. - + if norm_fk > rho*norm_fold - sigma2*norm_dx**2.: while norm_fk > (1.+eta)*norm_fold - sigma1*alp**2.*norm_dx**2.: - + alp *= beta xk = xold + alp * dx - - + + fk = func(xk) - + lcout += 1 norm_fk = numpy.linalg.norm(fk) if lcout == 20: @@ -48,7 +48,7 @@ def line_search_LF(func, xold, fold, dx, iter_): return alp, xk, fk def trustRegion(func, xold, fold, Binv, c = 0.5): - """Perform Trust Region Optimization. See "A Broyden Trust Region Quasi-Newton Method + """Perform Trust Region Optimization. See "A Broyden Trust Region Quasi-Newton Method for Nonlinear Equations" (https://www.iaeng.org/IJCS/issues_v46/issue_3/IJCS_46_3_09.pdf) Algorithm 1 for more information @@ -123,17 +123,17 @@ class FrankQN: Performs quasi newton optimization. Interfaces many functionalities of the frankestein code originaly written by Hong-Zhou Ye - + """ - + def __init__(self, func, x0, f0, J0, trust=0.5, max_space=500): - + self.x0 = x0 - self.n = x0.size + self.n = x0.size self.f0 = f0 self.func = func - + self.B0 = numpy.linalg.pinv(J0) self.iter_ = 0 @@ -143,19 +143,19 @@ def __init__(self, func, x0, f0, J0, trust=0.5, max_space=500): self.xold = None # old errvec self.fnew = None # new jacobian? self.fold = None # old jacobian? - self.max_subspace = max_space + self.max_subspace = max_space self.dxs = numpy.empty([self.max_subspace,self.n]) self.fs = numpy.empty([self.max_subspace,self.n]) self.us = numpy.empty([self.max_subspace,self.n]) # u_m = B_m @ f_m self.vs = numpy.empty([self.max_subspace,self.n]) # v_m = B_0 @ f_{m+1} self.B = None self.trust = trust - + def next_step(self, trust_region=False): - if self.iter_ == 0: - self.xnew = self.x0 - self.fnew = self.func(self.xnew) if self.f0 is None else self.f0 + if self.iter_ == 0: + self.xnew = self.x0 + self.fnew = self.func(self.xnew) if self.f0 is None else self.f0 self.fs[0] = self.fnew.copy() self.us[0] = numpy.dot(self.B0, self.fnew) self.Binv = self.B0.copy() @@ -164,26 +164,26 @@ def next_step(self, trust_region=False): if not self.iter_ == 0: dx_i = self.xnew - self.xold df_i = self.fnew - self.fold - + self.xold = self.xnew.copy() self.fold = self.fnew.copy() - + if not self.iter_ ==0: - + tmp__ = numpy.outer(dx_i - self.Binv@df_i, dx_i@self.Binv)/\ (dx_i@self.Binv@df_i) - self.Binv += tmp__ + self.Binv += tmp__ us_tmp = self.Binv@self.fnew if trust_region: self.xnew, self.fnew = trustRegion(self.func, self.xold, self.fold, self.Binv, c = self.trust) else: self.us[self.iter_] = self.get_Bnfn(self.iter_) - + alp, self.xnew, self.fnew = line_search_LF( self.func, self.xold, self.fold, -self.us[self.iter_], self.iter_) - - # udpate vs, dxs, and fs + + # udpate vs, dxs, and fs self.vs[self.iter_] = numpy.dot(self.B0, self.fnew) self.dxs[self.iter_] = self.xnew - self.xold self.fs[self.iter_+1] = self.fnew.copy() @@ -204,11 +204,11 @@ def get_Bnfn(self, n): for j in range(n-i+1): a = vs[j] b = vs[n-i] - un_ - + vps[j] = a + (dxn_@a)/(dxn_@b) * (dxn_ - b) - + vs = vps - + return vs[0] def get_be_error_jacobian(self,jac_solver='HF'): @@ -226,22 +226,22 @@ def get_be_error_jacobian(self,jac_solver='HF'): res_func = ccsdres_func elif jac_solver=='HF': res_func = hfres_func - + Ncout = [None] * self.Nfrag for A in range(self.Nfrag): Jes[A], Jcs[A], xes[A], xcs[A], ys[A], alphas[A], Ncout[A] = \ get_atbe_Jblock_frag(self.Fobjs[A], res_func) - + alpha = sum(alphas) - + # build Jacobian ''' ignore! F0-M1 F1-M2M2 F2-M3M3 F3-M4 M1 M2 M2 M3 M3 M4 - M1 E0 C1-1 + M1 E0 C1-1 M2 C0-0 E1 E1 - M2 E1 E1 C2-2 + M2 E1 E1 C2-2 M3 C1-1 E2 E2 M3 E2 E2 C3-3 M4 C2-2 E3 @@ -249,7 +249,7 @@ def get_be_error_jacobian(self,jac_solver='HF'): N_ = sum(Ncout) J = numpy.zeros((N_+1, N_+1)) cout = 0 - + for findx, fobj in enumerate(self.Fobjs): J[cout:Ncout[findx]+cout, cout:Ncout[findx]+cout] = Jes[findx] J[cout:Ncout[findx]+cout, N_:] = numpy.array(xes[findx]).reshape(-1,1) @@ -265,9 +265,9 @@ def get_be_error_jacobian(self,jac_solver='HF'): J[cout+coutc_: cout+coutc, N_:] += numpy.array( xcs[fobj.center[cindx]]).reshape(-1,1) coutc_ = coutc - cout += Ncout[findx] + cout += Ncout[findx] J[N_:,N_:] = alpha - + return J def get_atbe_Jblock_frag(fobj, res_func): @@ -275,13 +275,13 @@ def get_atbe_Jblock_frag(fobj, res_func): vpots = get_vpots_frag(fobj.nao, fobj.edge_idx, fobj.fsites) - eri_ = get_eri(fobj.dname, fobj.nao, eri_file=fobj.eri_file) + eri_ = get_eri(fobj.dname, fobj.nao, eri_file=fobj.eri_file) dm0 = numpy.dot( fobj._mo_coeffs[:,:fobj.nsocc], - fobj._mo_coeffs[:,:fobj.nsocc].T) *2. + fobj._mo_coeffs[:,:fobj.nsocc].T) *2. mf_ = get_scfObj(fobj.fock+fobj.heff, eri_, fobj.nsocc, dm0=dm0) - - dPs, dP_mu = res_func(mf_, vpots, eri_, fobj.nsocc) - + + dPs, dP_mu = res_func(mf_, vpots, eri_, fobj.nsocc) + Je = [] Jc = [] y=[] @@ -290,64 +290,64 @@ def get_atbe_Jblock_frag(fobj, res_func): cout = 0 for edge in fobj.edge_idx: - + for j_ in range(len(edge)): for k_ in range(len(edge)): if j_>k_: continue ## response w.r.t matching pot - # edges + # edges tmpje_ = [] - + for edge_ in fobj.edge_idx: lene = len(edge_) - + for j__ in range(lene): for k__ in range(lene): if j__>k__: continue - - tmpje_.append(dPs[cout][edge_[j__],edge_[k__]]) + + tmpje_.append(dPs[cout][edge_[j__],edge_[k__]]) y_ = 0. for fidx, fval in enumerate(fobj.fsites): - + if not any(fidx in sublist for sublist in fobj.edge_idx): y_ += dPs[cout][fidx, fidx] - + y.append(y_) - - tmpjc_ = [] + + tmpjc_ = [] # center on the same fragment - #for cen in fobj.efac[1]: + #for cen in fobj.efac[1]: for j__ in fobj.centerf_idx: for k__ in fobj.centerf_idx: if j__>k__: - continue + continue tmpjc_.append(-dPs[cout][j__,k__]) - + Je.append(tmpje_) - + Jc.append(tmpjc_) - + ## response w.r.t. chem pot # edge - xe.append(dP_mu[edge[j_],edge[k_]]) - cout += 1 + xe.append(dP_mu[edge[j_],edge[k_]]) + cout += 1 Je = numpy.array(Je).T Jc = numpy.array(Jc).T - + alpha = 0. - for fidx, fval in enumerate(fobj.fsites): + for fidx, fval in enumerate(fobj.fsites): if not any(fidx in sublist for sublist in fobj.edge_idx): - alpha += dP_mu[fidx, fidx] - + alpha += dP_mu[fidx, fidx] + for j__ in fobj.centerf_idx: for k__ in fobj.centerf_idx: if j__>k__: continue - xc.append(-dP_mu[j__,k__]) + xc.append(-dP_mu[j__,k__]) return Je, Jc, xe, xc, y, alpha, cout @@ -366,19 +366,19 @@ def get_be_error_jacobian_selffrag(self,jac_solver='HF'): res_func = ccsdres_func elif jac_solver=='HF': res_func = hfres_func - + Jes, Jcs, xes, xcs, ys, alphas, Ncout = \ get_atbe_Jblock_frag(self.Fobjs[0], res_func) - - N_ = Ncout + + N_ = Ncout J = numpy.zeros((N_+1, N_+1)) - + J[:Ncout, :Ncout] = Jes J[:Ncout, N_:] = numpy.array(xes).reshape(-1,1) J[N_:, :Ncout] = ys J[:Ncout, N_:] += numpy.array([*xcs, *xcs]).reshape(-1,1) - J[N_:,N_:] = alphas - + J[N_:,N_:] = alphas + return J @@ -392,23 +392,23 @@ def hfres_func(mf, vpots, eri, nsocc): us = cphf_kernel_batch(C, moe, eri, no, vpots) dPs = [get_rhf_dP_from_u(C, no, us[I]) for I in range(len(vpots)-1)] dP_mu = get_rhf_dP_from_u(C, no, us[-1]) - + return dPs, dP_mu def mp2res_func(mf, vpots, eri, nsocc): from molbe.external.cpmp2_utils import get_dPmp2_batch_r - + C = mf.mo_coeff moe = mf.mo_energy eri = mf._eri no = nsocc - + dPs_an = get_dPmp2_batch_r(C, moe, eri, no, vpots, aorep=True) dPs_an = numpy.array([dp_*0.5 for dp_ in dPs_an]) dP_mu = dPs_an[-1] - + return dPs_an[:-1], dP_mu - + def ccsdres_func(mf, vpots, eri, nsocc): from molbe.external.jac_utils import get_dPccsdurlx_batch_u @@ -427,7 +427,7 @@ def ccsdres_func(mf, vpots, eri, nsocc): def get_vpots_frag(nao, edge_idx, fsites): vpots = [] - + for edge_ in edge_idx: lene = len(edge_) for j__ in range(lene): @@ -443,13 +443,13 @@ def get_vpots_frag(nao, edge_idx, fsites): # outer edges not included tmppot = numpy.zeros((nao,nao)) for fidx, fval in enumerate(fsites): - + if not any(fidx in sublist for sublist in edge_idx): tmppot[fidx, fidx] = -1 - + vpots.append(tmppot) return vpots - + diff --git a/molbe/external/uccsd_eri.py b/molbe/external/uccsd_eri.py index dbe0278a..9c40ed5d 100644 --- a/molbe/external/uccsd_eri.py +++ b/molbe/external/uccsd_eri.py @@ -4,7 +4,7 @@ def make_eris_incore(mycc, Vss, Vos, mo_coeff=None, ao2mofn=None, frozen=False): vhf = frank_get_veff(mycc, mycc._scf.make_rdm1(mycc.mo_coeff, mycc.mo_occ), Vss, Vos) fockao = frank_get_fock(mycc, vhf, frozen) - mycc._scf.get_veff = lambda *args, **kwargs: vhf + mycc._scf.get_veff = lambda *args, **kwargs: vhf mycc._scf.get_fock = lambda *args, **kwargs: fockao return _make_eris_incore(mycc, mo_coeff=mo_coeff, ao2mofn=ao2mofn) #, frozen) diff --git a/molbe/external/unrestricted_utils.py b/molbe/external/unrestricted_utils.py index 1251ed11..5a44e66a 100644 --- a/molbe/external/unrestricted_utils.py +++ b/molbe/external/unrestricted_utils.py @@ -32,7 +32,7 @@ def make_uhf_obj(fobj_a, fobj_b, frozen=False): full_uhf.TA = [fobj_a.TA, fobj_b.TA] - # Build core components + # Build core components if frozen: full_uhf.gcores_raw = [fobj_a.TA.T @ (fobj_a.hf_veff-fobj_a.core_veff) @ fobj_a.TA, fobj_b.TA.T @ (fobj_b.hf_veff-fobj_b.core_veff) @ fobj_b.TA] diff --git a/molbe/fragment.py b/molbe/fragment.py index 5bc8433e..cd2ca22b 100644 --- a/molbe/fragment.py +++ b/molbe/fragment.py @@ -8,7 +8,7 @@ class fragpart: """Fragment/partitioning definition Interfaces two main fragmentation functions (autogen & chain) in MolBE. It defines edge & - center for density matching and energy estimation. It also forms the base for IAO/PAO partitioning for + center for density matching and energy estimation. It also forms the base for IAO/PAO partitioning for a large basis set bootstrap calculation. Parameters @@ -28,7 +28,7 @@ class fragpart: valence_basis: str Name of minimal basis set for IAO scheme. 'sto-3g' suffice for most cases. valence_only: bool - If this option is set to True, all calculation will be performed in the valence basis in the IAO partitioning. + If this option is set to True, all calculation will be performed in the valence basis in the IAO partitioning. This is an experimental feature. frozen_core: bool Whether to invoke frozen core approximation. This is set to False by default @@ -43,7 +43,7 @@ def __init__(self, frag_type='autogen', valence_basis=None,valence_only=False, print_frags=True, write_geom=False, be_type='be2', mol=None, frozen_core=False): - + # Initialize class attributes self.mol = mol self.frag_type = frag_type @@ -56,7 +56,7 @@ def __init__(self, frag_type='autogen', self.center_idx = [] self.centerf_idx = [] self.be_type = be_type - self.frozen_core = frozen_core + self.frozen_core = frozen_core self.valence_basis = valence_basis self.valence_only = valence_only @@ -80,26 +80,26 @@ def __init__(self, frag_type='autogen', flush=True) print('exiting',flush=True) sys.exit() - self.chain(mol, frozen_core=frozen_core,closed=closed) + self.chain(mol, frozen_core=frozen_core,closed=closed) elif frag_type=='autogen': if mol is None: print('Provide pyscf gto.M object in fragpart() and restart!', flush=True) print('exiting',flush=True) sys.exit() - + fgs = autogen(mol, be_type=be_type, frozen_core=frozen_core,write_geom=write_geom, valence_basis=valence_basis, valence_only=valence_only, print_frags=print_frags) - + self.fsites, self.edge, self.center, self.edge_idx, self.center_idx, self.centerf_idx, self.ebe_weight, self.Frag_atom, self.center_atom, self.hlist_atom, self.add_center_atom = fgs self.Nfrag = len(self.fsites) - + else: print('Fragmentation type = ',frag_type,' not implemented!', flush=True) print('exiting',flush=True) sys.exit() - + from .lchain import chain def hchain_simple(self): """Hard coded fragmentation feature @@ -110,7 +110,7 @@ def hchain_simple(self): self.fsites.append([i]) self.edge.append([]) self.Nfrag = len(self.fsites) - + elif self.be_type=='be2': for i in range(self.natom-2): self.fsites.append([i, i+1, i+2]) @@ -121,47 +121,47 @@ def hchain_simple(self): for i in self.fsites[1:-1]: self.edge.append([[i[0]],[i[-1]]]) self.edge.append([[self.fsites[-1][0]]]) - + self.center.append([1]) for i in range(self.Nfrag-2): self.center.append([i,i+2]) self.center.append([self.Nfrag-2]) - + elif self.be_type=='be3': for i in range(self.natom-4): self.fsites.append([i, i+1, i+2, i+3, i+4]) self.centerf_idx.append([2]) self.Nfrag = len(self.fsites) - + self.edge.append([[3],[4]]) for i in self.fsites[1:-1]: self.edge.append([[i[0]],[i[1]],[i[-2]],[i[-1]]]) self.edge.append([[self.fsites[-1][0]],[self.fsites[-1][1]]]) - + self.center.append([1,2]) self.center.append([0,0,2,3]) for i in range(self.Nfrag-4): self.center.append([i,i+1, i+3,i+4]) - + self.center.append([self.Nfrag-4,self.Nfrag-3, self.Nfrag-1,self.Nfrag-1]) - self.center.append([self.Nfrag-3,self.Nfrag-2]) - + self.center.append([self.Nfrag-3,self.Nfrag-2]) + for ix, i in enumerate(self.fsites): tmp_ = [] elist_ = [ xx for yy in self.edge[ix] for xx in yy] for j in i: if not j in elist_: tmp_.append(i.index(j)) - self.ebe_weight.append([1.0, tmp_]) - + self.ebe_weight.append([1.0, tmp_]) + if not self.be_type=='be1': for i in range(self.Nfrag): idx = [] for j in self.edge[i]: - + idx.append([self.fsites[i].index(k) for k in j]) self.edge_idx.append(idx) - + for i in range(self.Nfrag): idx = [] for j in range(len(self.center[i])): diff --git a/molbe/helper.py b/molbe/helper.py index 5b4bb00a..e6be4208 100644 --- a/molbe/helper.py +++ b/molbe/helper.py @@ -29,14 +29,14 @@ def get_veff(eri_, dm, S, TA, hf_veff): numpy.ndarray Effective HF potential in the embedding basis. """ - + import functools from pyscf import scf # Transform the density matrix ST = numpy.dot(S, TA) P_ = functools.reduce(numpy.dot,(ST.T, dm, ST)) - + # Ensure the transformed density matrix and ERI are real and double-precision P_ = numpy.asarray(P_.real, dtype=numpy.double) eri_ = numpy.asarray(eri_, dtype=numpy.double) @@ -45,7 +45,7 @@ def get_veff(eri_, dm, S, TA, hf_veff): vj, vk = scf.hf.dot_eri_dm(eri_, P_, hermi=1, with_j=True, with_k=True) Veff_ = vj - 0.5 * vk Veff = functools.reduce(numpy.dot,(TA.T, hf_veff, TA)) - Veff_ - + return Veff # create pyscf pbc scf object @@ -79,7 +79,7 @@ def get_scfObj(h1, Eri, nocc, dm0=None, enuc=0., """ # from 40-customizing_hamiltonian.py in pyscf examples from pyscf import gto, scf - + nao = h1.shape[0] # Initialize a dummy molecule with the required number of electrons @@ -87,22 +87,22 @@ def get_scfObj(h1, Eri, nocc, dm0=None, enuc=0., mol = gto.M() mol.nelectron = nocc * 2 mol.incore_anyway = True - + # Initialize an RHF object mf_ = scf.RHF(mol) mf_.get_hcore = lambda *args:h1 mf_.get_ovlp = lambda *args: S mf_._eri = Eri mf_.incore_anyway = True - mf_.max_cycle=50 + mf_.max_cycle=50 mf_.verbose=0 - + # Run the SCF calculation if dm0 is None: mf_.kernel() else: mf_.kernel(dm0=dm0) - + # Check if the SCF calculation converged if not mf_.converged: print(flush=True) @@ -110,7 +110,7 @@ def get_scfObj(h1, Eri, nocc, dm0=None, enuc=0., print(flush=True) mf_.verbose=0 mf_.level_shift=0.2 - mf_.diis_space=25 + mf_.diis_space=25 if dm0 is None: mf_.kernel() else: @@ -123,7 +123,7 @@ def get_scfObj(h1, Eri, nocc, dm0=None, enuc=0., print(flush=True) print('SCF Converged!',flush=True) print(flush=True) - + return mf_ @@ -160,7 +160,7 @@ def get_eri(i_frag, Nao, symm = 8, ignore_symm = False, eri_file='eri_file.h5'): # Optionally restore the symmetry of the ERI if not ignore_symm: - + # Set the number of threads for the library to 1 lib.num_threads(1) eri__ = ao2mo.restore(symm, eri__, Nao) @@ -170,7 +170,7 @@ def get_eri(i_frag, Nao, symm = 8, ignore_symm = False, eri_file='eri_file.h5'): return eri__ def ncore_(z): - + if 1<= z<=2: nc = 0 elif 2<=z<=5: @@ -190,7 +190,7 @@ def ncore_(z): flush=True) print('exiting',flush=True) sys.exit() - + return nc @@ -257,10 +257,10 @@ def get_frag_energy(mo_coeffs, nsocc, nfsites, efac, TA, h1, hf_veff, rdm1, rdm2 List containing the energy contributions: [e1_tmp, e2_tmp, ec_tmp]. """ - # Rotate the RDM1 into the MO basis + # Rotate the RDM1 into the MO basis rdm1s_rot = mo_coeffs @ rdm1 @ mo_coeffs.T * 0.5 - # Construct the Hartree-Fock 1-RDM + # Construct the Hartree-Fock 1-RDM hf_1rdm = numpy.dot(mo_coeffs[:,:nsocc], mo_coeffs[:,:nsocc].conj().T) @@ -296,9 +296,9 @@ def get_frag_energy(mo_coeffs, nsocc, nfsites, efac, TA, h1, hf_veff, rdm1, rdm2 for i in range(nfsites): for j in range(jmax): ij = i*(i+1)//2+j if i > j else j*(j+1)//2+i - Gij = rdm2s[i,j,:jmax,:jmax].copy() + Gij = rdm2s[i,j,:jmax,:jmax].copy() Gij[numpy.diag_indices(jmax)] *= 0.5 - Gij += Gij.T + Gij += Gij.T e2[i] += Gij[numpy.tril_indices(jmax)] @ eri[ij] # Sum the energy contributions @@ -313,13 +313,13 @@ def get_frag_energy(mo_coeffs, nsocc, nfsites, efac, TA, h1, hf_veff, rdm1, rdm2 # Calculate the total energy contribution for the specified fragment indices for i in efac[1]: etmp += efac[0]*e_[i] - e1_tmp += efac[0]*e1[i] - e2_tmp += efac[0]*e2[i] + e1_tmp += efac[0]*e1[i] + e2_tmp += efac[0]*e2[i] ec_tmp += efac[0]*ec[i] - + return [e1_tmp,e2_tmp,ec_tmp] -def get_frag_energy_u(mo_coeffs, nsocc, nfsites, efac, TA, h1, hf_veff, rdm1, rdm2s, dname, +def get_frag_energy_u(mo_coeffs, nsocc, nfsites, efac, TA, h1, hf_veff, rdm1, rdm2s, dname, eri_file='eri_file.h5', gcores=None, frozen=False, veff0=None): """ Compute the fragment energy for unrestricted calculations @@ -351,7 +351,7 @@ def get_frag_energy_u(mo_coeffs, nsocc, nfsites, efac, TA, h1, hf_veff, rdm1, rd Dataset name in the HDF5 file. eri_file : str, optional Filename of the HDF5 file containing the electron repulsion integrals. Defaults to 'eri_file.h5'. - gcores : + gcores : frozen : bool, optional Indicate frozen core. Default is False @@ -385,7 +385,7 @@ def get_frag_energy_u(mo_coeffs, nsocc, nfsites, efac, TA, h1, hf_veff, rdm1, rd # Calculate the one-electron and effective potential energy contributions e1 = [numpy.einsum("ij,ij->i",h1[s][:nfsites[s]], delta_rdm1[s][:nfsites[s]]) for s in [0,1]] ec = [numpy.einsum("ij,ij->i",veff0[s][:nfsites[s]], delta_rdm1[s][:nfsites[s]]) for s in [0,1]] - + jmax = [TA[0].shape[1],TA[1].shape[1]] # Load ERIs from the HDF5 file diff --git a/molbe/lchain.py b/molbe/lchain.py index b67a4a74..df52d0c0 100644 --- a/molbe/lchain.py +++ b/molbe/lchain.py @@ -9,7 +9,7 @@ def chain(self,mol, frozen_core=False, closed=False): for adx, bas in enumerate(mol.aoslice_by_atom()): start_ = bas[2] stop_ = bas[3] - + if not frozen_core: sites.append([i for i in range(start_, stop_)]) else: @@ -27,7 +27,7 @@ def chain(self,mol, frozen_core=False, closed=False): sys.exit() Ns = mol.aoslice_by_atom()[-1][3] - + if self.be_type=='be2': fs=[] if closed: @@ -36,7 +36,7 @@ def chain(self,mol, frozen_core=False, closed=False): #sites_left = [i for i in sites[-1]] ns = len(sites[-1]) sites_right = [i+ns for i in sites[-1]] - + if closed: self.fsites.append(sites_left+sites[0]+sites[1]) fs.append([sites_left, sites[0], sites[1]]) @@ -46,7 +46,7 @@ def chain(self,mol, frozen_core=False, closed=False): if closed: self.fsites.append(sites[-2]+sites[-1]+sites_right) fs.append([sites[-2],sites[-1],sites_right]) - + self.Nfrag = len(self.fsites) if closed: @@ -74,24 +74,24 @@ def chain(self,mol, frozen_core=False, closed=False): if closed: for ix, i in enumerate(self.fsites): - tmp_ = [] + tmp_ = [] elist_ = [xx for yy in self.edge[ix] for xx in yy] for j in i: if not j in elist_: tmp_.append(i.index(j)) self.ebe_weight.append([1.0, tmp_]) - + for i in range(self.Nfrag): self.centerf_idx.append([self.fsites[i].index(j) for j in fs[i][1]]) - + if self.be_type=='be3' and not closed: fs=[] - + for i in range(mol.natm-4): self.fsites.append(sites[i]+ sites[i+1]+ sites[i+2]+ sites[i+3]+ sites[i+4]) fs.append([ sites[i], sites[i+1], sites[i+2], sites[i+3], sites[i+4]]) - + self.Nfrag = len(self.fsites) self.edge.append([fs[0][3], fs[0][4]]) @@ -105,11 +105,11 @@ def chain(self,mol, frozen_core=False, closed=False): self.center.append([0,0,2,3]) for i in range(self.Nfrag-4): self.center.append([i, i+1, i+3, i+4]) - + self.center.append([self.Nfrag-4, self.Nfrag-3, self.Nfrag-1, self.Nfrag-1]) self.center.append([self.Nfrag-3, self.Nfrag-2]) - + for i in range(self.Nfrag): self.centerf_idx.append([self.fsites[i].index(j) for j in fs[i][2]]) @@ -134,7 +134,7 @@ def chain(self,mol, frozen_core=False, closed=False): self.fsites.append(sites[-3]+sites[-2]+sites[-1]+sites_right0+sites_right1) fs.append([sites[-4],sites[-3],sites[-2],sites[-1],sites_right0]) fs.append([sites[-3],sites[-2],sites[-1],sites_right0,sites_right1]) - + self.Nfrag = len(self.fsites) #self.edge.append([fs[0][3], fs[0][4]]) @@ -160,7 +160,7 @@ def chain(self,mol, frozen_core=False, closed=False): self.ebe_weight.append([1.0, tmp_]) for i in range(self.Nfrag): self.centerf_idx.append([self.fsites[i].index(j) for j in fs[i][2]]) - + if self.be_type=='be4' and not closed: fs=[] for i in range(mol.natm-6): @@ -181,40 +181,40 @@ def chain(self,mol, frozen_core=False, closed=False): self.center.append([0,0,1,3,4,5]) for i in range(self.Nfrag-6): self.center.append([i, i+1,i+2, i+4, i+5,i+6]) - + self.center.append([self.Nfrag-6, self.Nfrag-5, self.Nfrag-4, self.Nfrag-2, self.Nfrag-1, self.Nfrag-1]) self.center.append([self.Nfrag-5, self.Nfrag-4, self.Nfrag-3, self.Nfrag-1, self.Nfrag-1, self.Nfrag-1]) self.center.append([self.Nfrag-4, self.Nfrag-3, self.Nfrag-2]) - + for i in range(self.Nfrag): - self.centerf_idx.append([self.fsites[i].index(j) for j in fs[i][3]]) + self.centerf_idx.append([self.fsites[i].index(j) for j in fs[i][3]]) if self.be_type=='be4' and closed: print('Will add this soon!') sys.exit() - if not closed: + if not closed: for ix, i in enumerate(self.fsites): tmp_ = [] elist_ = [ xx for yy in self.edge[ix] for xx in yy] for j in i: if not j in elist_: tmp_.append(i.index(j)) self.ebe_weight.append([1.0, tmp_]) - + if not self.be_type=='be1': for i in range(self.Nfrag): idx = [] for j in self.edge[i]: - + idx.append([self.fsites[i].index(k) for k in j]) self.edge_idx.append(idx) if closed: for i in range(self.Nfrag): - idx = [] - + idx = [] + for j in self.center[i]: idx__ = [] for fidx, fval in enumerate(self.fsites[j]): @@ -226,9 +226,9 @@ def chain(self,mol, frozen_core=False, closed=False): for i in range(self.Nfrag): idx = [] for j in range(len(self.center[i])): - + idx.append([self.fsites[self.center[i][j]].index(k) for k in self.edge[i][j]]) - self.center_idx.append(idx) - + self.center_idx.append(idx) + diff --git a/molbe/lo.py b/molbe/lo.py index 01e5301c..2ab34cff 100644 --- a/molbe/lo.py +++ b/molbe/lo.py @@ -37,7 +37,7 @@ def get_symm_orth_mat(A, thr=1.E-6, ovlp=None): def symm_orth(A, thr=1.E-6, ovlp=None): """ Symmetrically orthogonalize columns of A - """ + """ U = get_symm_orth_mat(A, thr, ovlp) return A @ U @@ -45,23 +45,23 @@ def symm_orth(A, thr=1.E-6, ovlp=None): def remove_core_mo(Clo, Ccore, S, thr=0.5): assert(numpy.allclose(Clo.T@S@Clo,numpy.eye(Clo.shape[1]))) assert(numpy.allclose(Ccore.T@S@Ccore,numpy.eye(Ccore.shape[1]))) - + n,nlo = Clo.shape ncore = Ccore.shape[1] Pcore = Ccore@Ccore.T @ S - Clo1 = (numpy.eye(n) - Pcore) @ Clo + Clo1 = (numpy.eye(n) - Pcore) @ Clo pop = numpy.diag(Clo1.T @ S @ Clo1) idx_keep = numpy.where(pop>thr)[0] - assert(len(idx_keep) == nlo-ncore) + assert(len(idx_keep) == nlo-ncore) Clo2 = symm_orth(Clo1[:,idx_keep], ovlp=S) return Clo2 - + def reorder_lo(C, S, idao_by_atom, atom_by_motif, motifname, ncore_by_motif, thresh=0.5, verbose=3): - """ Reorder localized orbitals - + """ Reorder localized orbitals + This function reorders the IAOs and PAOs so that the IAOs and PAOs for each atom are grouped together. """ @@ -70,7 +70,7 @@ def reorder_lo(C, S, idao_by_atom, atom_by_motif, motifname, for idao in idao_by_atom: pop = numpy.sum(pop_by_ao[idao], axis = 0) - + def get_xovlp(mol, basis='sto-3g'): """ Gets set of valence orbitals based on smaller (should be minimal) basis @@ -108,7 +108,7 @@ def get_iao(Co, S12, S1, S2 = None): S2 = S12.T @ numpy.linalg.inv(S1) @ S12 P1 = numpy.linalg.inv(S1) P2 = numpy.linalg.inv(S2) - + # depolarized occ mo Cotil = P1 @ S12 @ P2 @ S12.T @ Co @@ -121,8 +121,8 @@ def get_iao(Co, S12, S1, S2 = None): Ciao = (numpy.eye(n) - (Po + Potil - 2 * Po @ S1 @ Potil) @ S1) @ ptil Ciao = symm_orth(Ciao, ovlp = S1) - - # check span + + # check span rep_err = numpy.linalg.norm(Ciao @ Ciao.T @ S1 @ Po - Po) if rep_err > 1.E-10: raise RuntimeError @@ -143,16 +143,16 @@ def get_pao(Ciao, S, S12, S2, mol): s12 = numpy.linalg.inv(S) @ S12 nonval = numpy.eye(n) - s12 @ s12.T # set of orbitals minus valence (orth in working basis) - + Piao = Ciao @ Ciao.T @ S # projector into IAOs Cpao = (numpy.eye(n) - Piao) @ nonval # project out IAOs from non-valence basis - + # begin canonical orthogonalization to get rid of redundant orbitals numpy.o0 = Cpao.shape[1] Cpao = cano_orth(Cpao, ovlp=S) numpy.o1 = Cpao.shape[1] - + return Cpao def get_pao_native(Ciao, S, mol, valence_basis): @@ -166,7 +166,7 @@ def get_pao_native(Ciao, S, mol, valence_basis): Cpao (symmetrically orthogonalized) """ n = Ciao.shape[0] - + # Form a mol object with the valence basis for the ao_labels mol_alt = mol.copy() mol_alt.basis = valence_basis @@ -201,12 +201,12 @@ def get_loc(mol, C, method, pop_method=None, init_guess=None): from pyscf.lo import Boys as Localizer else: raise NotImplementedError('Localization scheme not understood') - + mlo = Localizer(mol, C) if pop_method is not None: mlo.pop_method=str(pop_method) - mlo.init_guess = init_guess + mlo.init_guess = init_guess C_ = mlo.kernel() return C_ @@ -228,7 +228,7 @@ def localize(self, lo_method, mol=None, valence_basis='sto-3g', valence_basis: str Name of minimal basis set for IAO scheme. 'sto-3g' suffice for most cases. valence_only: bool - If this option is set to True, all calculation will be performed in the valence basis in the IAO partitioning. + If this option is set to True, all calculation will be performed in the valence basis in the IAO partitioning. This is an experimental feature. """ from numpy.linalg import eigh @@ -236,11 +236,11 @@ def localize(self, lo_method, mol=None, valence_basis='sto-3g', from pyscf.lo import orth import scipy.linalg,functools from .helper import ncore_ - + if lo_method == 'lowdin': - + es_, vs_ = eigh(self.S) - edx = es_ > 1.e-15 + edx = es_ > 1.e-15 self.W = numpy.dot(vs_[:,edx]/numpy.sqrt(es_[edx]), vs_[:,edx].T) if self.frozen_core: if self.unrestricted: @@ -275,7 +275,7 @@ def localize(self, lo_method, mol=None, valence_basis='sto-3g', s_ = numpy.diag(1.0/s_) W_ = functools.reduce(numpy.dot, (vs_, s_, vs_.T)) - self.W = numpy.dot(C_, W_) + self.W = numpy.dot(C_, W_) if self.unrestricted: if self.frozen_core: @@ -294,18 +294,18 @@ def localize(self, lo_method, mol=None, valence_basis='sto-3g', (self.W.T, self.S, self.C[:,self.ncore:])) else: self.lmo_coeff = functools.reduce(numpy.dot, - (self.W.T, self.S, self.C)) + (self.W.T, self.S, self.C)) elif lo_method in ['pipek-mezey','pipek', 'PM']: - + es_, vs_ = eigh(self.S) - edx = es_ > 1.e-15 + edx = es_ > 1.e-15 self.W = numpy.dot(vs_[:,edx]/numpy.sqrt(es_[edx]), vs_[:,edx].T) - + es_, vs_ = eigh(self.S) - edx = es_ > 1.e-15 + edx = es_ > 1.e-15 W_ = numpy.dot(vs_[:,edx]/numpy.sqrt(es_[edx]), vs_[:,edx].T) - if self.frozen_core: + if self.frozen_core: P_core = numpy.eye(W_.shape[0]) - numpy.dot(self.P_core, self.S) C_ = numpy.dot(P_core, W_) Cpop = functools.reduce(numpy.dot, @@ -319,72 +319,72 @@ def localize(self, lo_method, mol=None, valence_basis='sto-3g', s_ = numpy.diag(1.0/s_) W_ = functools.reduce(numpy.dot, (vs_, s_, vs_.T)) - W_ = numpy.dot(C_, W_) + W_ = numpy.dot(C_, W_) self.W = get_loc(self.mol, W_, 'PM', pop_method=pop_method, init_guess=init_guess) - + if not self.frozen_core: self.lmo_coeff = self.W.T @ self.S @ self.C - else: + else: self.lmo_coeff = self.W.T @ self.S @ self.C[:,self.ncore:] - elif lo_method=='iao': + elif lo_method=='iao': from pyscf import lo import os, h5py - + loc_type = 'SO' val_basis = 'sto-3g' - + # Occupied mo_coeff (with core) Co = self.C[:,:self.Nocc] # Get necessary overlaps, second arg is IAO basis S12, S2 = get_xovlp(self.mol, basis=val_basis) # Use these to get IAOs Ciao = get_iao(Co, S12, self.S, S2 = S2) - + if not valence_only: # Now get PAOs if loc_type.upper() != 'SO': Cpao = get_pao(Ciao, self.S, S12, S2, self.mol) elif loc_type.upper() == 'SO': Cpao = get_pao_native(Ciao, self.S, self.mol, valence_basis=val_basis) - + # rearrange by atom aoind_by_atom = get_aoind_by_atom(self.mol) Ciao, iaoind_by_atom = reorder_by_atom_(Ciao, aoind_by_atom, self.S) if not valence_only: Cpao, paoind_by_atom = reorder_by_atom_(Cpao, aoind_by_atom, self.S) - + if self.frozen_core: # Remove core MOs Cc = self.C[:,:self.ncore] # Assumes core are first Ciao = remove_core_mo(Ciao, Cc, self.S) - + # Localize orbitals beyond symm orth if loc_type.upper() != 'SO': Ciao = get_loc(self.mol, Ciao, loc_type) if not valence_only: Cpao = get_loc(self.mol, Cpao, loc_type) - + shift = 0 ncore = 0 - if not valence_only: + if not valence_only: Wstack = numpy.zeros((Ciao.shape[0], Ciao.shape[1]+Cpao.shape[1])) #-self.ncore)) else: Wstack = numpy.zeros((Ciao.shape[0], Ciao.shape[1])) - - if self.frozen_core: + + if self.frozen_core: for ix in range(self.mol.natm): nc = ncore_(self.mol.atom_charge(ix)) ncore += nc niao = len(iaoind_by_atom[ix]) iaoind_ix = [ i_ - ncore for i_ in iaoind_by_atom[ix][nc:]] - Wstack[:, shift:shift+niao-nc] = Ciao[:, iaoind_ix] + Wstack[:, shift:shift+niao-nc] = Ciao[:, iaoind_ix] shift += niao-nc if not valence_only: npao = len(paoind_by_atom[ix]) - Wstack[:,shift:shift+npao] = Cpao[:, paoind_by_atom[ix]] + Wstack[:,shift:shift+npao] = Cpao[:, paoind_by_atom[ix]] shift += npao else: if not hstack: @@ -397,9 +397,9 @@ def localize(self, lo_method, mol=None, valence_basis='sto-3g', Wstack[:,shift:shift+npao] = Cpao[:, paoind_by_atom[ix]] shift += npao else: - Wstack = numpy.hstack((Ciao, Cpao)) - if not nosave: - self.W = Wstack + Wstack = numpy.hstack((Ciao, Cpao)) + if not nosave: + self.W = Wstack assert(numpy.allclose(self.W.T @ self.S @ self.W, numpy.eye(self.W.shape[1]))) else: assert(numpy.allclose(Wstack.T @ self.S @ Wstack, numpy.eye(Wstack.shape[1]))) @@ -408,7 +408,7 @@ def localize(self, lo_method, mol=None, valence_basis='sto-3g', nlo = self.W.shape[1] if not valence_only: - if nmo > nlo: + if nmo > nlo: Co_nocore = self.C[:,self.ncore:self.Nocc] Cv = self.C[:,self.Nocc:] # Ensure that the LOs span the occupied space @@ -422,15 +422,15 @@ def localize(self, lo_method, mol=None, valence_basis='sto-3g', self.lmo_coeff = self.W.T @ self.S @ C_ else: self.lmo_coeff = self.W.T @ self.S @ self.C[:,self.ncore:] - else: + else: self.lmo_coeff = self.W.T @ self.S @ self.C[:,self.ncore:] - + elif lo_method == 'boys': from pyscf.lo import Boys es_, vs_ = eigh(self.S) - edx = es_ > 1.e-15 + edx = es_ > 1.e-15 W_ = numpy.dot(vs_[:,edx]/numpy.sqrt(es_[edx]), vs_[:,edx].T) - if self.frozen_core: + if self.frozen_core: P_core = numpy.eye(W_.shape[0]) - numpy.dot(self.P_core, self.S) C_ = numpy.dot(P_core, W_) Cpop = functools.reduce(numpy.dot, @@ -444,15 +444,15 @@ def localize(self, lo_method, mol=None, valence_basis='sto-3g', s_ = numpy.diag(1.0/s_) W_ = functools.reduce(numpy.dot, (vs_, s_, vs_.T)) - W_ = numpy.dot(C_, W_) - + W_ = numpy.dot(C_, W_) + self.W = get_loc(self.mol, W_, 'BOYS') - + if not self.frozen_core: self.lmo_coeff = self.W.T @ self.S @ self.C - else: - self.lmo_coeff = self.W.T @ self.S @ self.C[:,self.ncore:] - + else: + self.lmo_coeff = self.W.T @ self.S @ self.C[:,self.ncore:] + else: print('lo_method = ',lo_method,' not implemented!',flush=True) print('exiting',flush=True) diff --git a/molbe/mbe.py b/molbe/mbe.py index 95aa2608..d9221a8e 100644 --- a/molbe/mbe.py +++ b/molbe/mbe.py @@ -10,7 +10,7 @@ class storeBE: def __init__(self, Nocc, hf_veff, hcore, S, C, hf_dm, hf_etot, W, lmo_coeff, - enuc, + enuc, E_core, C_core, P_core, core_veff, mo_energy): self.Nocc = Nocc self.hf_veff = hf_veff @@ -48,11 +48,11 @@ class BE: Method for orbital localization, default is 'lowdin'. """ - def __init__(self, mf, fobj, eri_file='eri_file.h5', - lo_method='lowdin', pop_method=None, compute_hf=True, + def __init__(self, mf, fobj, eri_file='eri_file.h5', + lo_method='lowdin', pop_method=None, compute_hf=True, restart=False, save=False, restart_file='storebe.pk', - mo_energy = None, + mo_energy = None, save_file='storebe.pk',hci_pt=False, nproc=1, ompnum=4, scratch_dir=None, hci_cutoff=0.001, ci_coeff_cutoff = None, select_cutoff=None, @@ -91,7 +91,7 @@ def __init__(self, mf, fobj, eri_file='eri_file.h5', auxbasis : str, optional Auxiliary basis for density fitting, by default None (uses default auxiliary basis defined in PySCF). """ - + if restart: # Load previous calculation data from restart file with open(restart_file, 'rb') as rfile: @@ -112,7 +112,7 @@ def __init__(self, mf, fobj, eri_file='eri_file.h5', self.P_core = store_.P_core self.core_veff = store_.core_veff self.mo_energy = store_.mo_energy - + self.unrestricted = False self.nproc = nproc self.ompnum = ompnum @@ -121,7 +121,7 @@ def __init__(self, mf, fobj, eri_file='eri_file.h5', # Fragment information from fobj self.frag_type=fobj.frag_type - self.Nfrag = fobj.Nfrag + self.Nfrag = fobj.Nfrag self.fsites = fobj.fsites self.edge = fobj.edge self.center = fobj.center @@ -131,27 +131,27 @@ def __init__(self, mf, fobj, eri_file='eri_file.h5', self.ebe_weight = fobj.ebe_weight self.be_type = fobj.be_type self.mol = fobj.mol - + self.ebe_hf = 0. self.ebe_tot = 0. - + # HCI parameters self.hci_cutoff = hci_cutoff self.ci_coeff_cutoff = ci_coeff_cutoff self.select_cutoff = select_cutoff self.hci_pt=hci_pt - - self.mf = mf - if not restart: + + self.mf = mf + if not restart: self.mo_energy = mf.mo_energy - + self.mf = mf - self.Nocc = mf.mol.nelectron//2 + self.Nocc = mf.mol.nelectron//2 self.enuc = mf.energy_nuc() - + self.hcore = mf.get_hcore() self.S = mf.get_ovlp() - self.C = numpy.array(mf.mo_coeff) + self.C = numpy.array(mf.mo_coeff) self.hf_dm = mf.make_rdm1() self.hf_veff = mf.get_veff() self.hf_etot = mf.e_tot @@ -164,7 +164,7 @@ def __init__(self, mf, fobj, eri_file='eri_file.h5', self.pot = initialize_pot(self.Nfrag, self.edge_idx) self.eri_file = eri_file self.scratch_dir = scratch_dir - + # Set scratch directory jobid='' if be_var.CREATE_SCRATCH_DIR: @@ -172,7 +172,7 @@ def __init__(self, mf, fobj, eri_file='eri_file.h5', jobid = str(os.environ['SLURM_JOB_ID']) except: jobid = '' - if not be_var.SCRATCH=='': + if not be_var.SCRATCH=='': self.scratch_dir = be_var.SCRATCH+str(jobid) os.system('mkdir -p '+self.scratch_dir) else: @@ -181,7 +181,7 @@ def __init__(self, mf, fobj, eri_file='eri_file.h5', self.eri_file = be_var.SCRATCH+eri_file else: self.eri_file = self.scratch_dir+'/'+eri_file - + self.frozen_core = False if not fobj.frozen_core else True self.ncore = 0 if not restart: @@ -189,57 +189,57 @@ def __init__(self, mf, fobj, eri_file='eri_file.h5', self.C_core = None self.P_core = None self.core_veff = None - + if self.frozen_core: # Handle frozen core orbitals self.ncore = fobj.ncore self.no_core_idx = fobj.no_core_idx self.core_list = fobj.core_list - + if not restart: - self.Nocc -=self.ncore + self.Nocc -=self.ncore self.hf_dm = 2.*numpy.dot(self.C[:,self.ncore:self.ncore+self.Nocc], self.C[:,self.ncore:self.ncore+self.Nocc].T) self.C_core = self.C[:,:self.ncore] self.P_core = numpy.dot(self.C_core, self.C_core.T) self.core_veff = mf.get_veff(dm = self.P_core*2.) - self.E_core = numpy.einsum('ji,ji->',2.*self.hcore+self.core_veff, self.P_core) + self.E_core = numpy.einsum('ji,ji->',2.*self.hcore+self.core_veff, self.P_core) self.hf_veff -= self.core_veff self.hcore += self.core_veff - + if not restart: # Localize orbitals self.localize(lo_method, pop_method=pop_method, mol=self.mol, valence_basis=fobj.valence_basis, valence_only=fobj.valence_only) - + if fobj.valence_only and lo_method=='iao': self.Ciao_pao = self.localize(lo_method, pop_method=pop_method, mol=self.mol, valence_basis=fobj.valence_basis, hstack=True, valence_only=False, nosave=True) - + if save: # Save intermediate results for restart store_ = storeBE(self.Nocc, self.hf_veff, self.hcore, self.S, self.C, self.hf_dm, self.hf_etot, - self.W, self.lmo_coeff, self.enuc, + self.W, self.lmo_coeff, self.enuc, self.E_core, self.C_core, self.P_core, self.core_veff, self.mo_energy) with open(save_file, 'wb') as rfile: pickle.dump(store_, rfile, pickle.HIGHEST_PROTOCOL) rfile.close() - - + + if not restart : # Initialize fragments and perform initial calculations self.initialize(mf._eri,compute_hf) - else: + else: self.initialize(None,compute_hf, restart=True) - - + + from ._opt import optimize from molbe.external.optqn import get_be_error_jacobian from .lo import localize from .rdm import rdm1_fullbasis, compute_energy_full - + def print_ini(self): """ Print initialization banner for the MOLBE calculation. @@ -254,14 +254,14 @@ def print_ini(self): print(' M M OO OO LL BB B EE ',flush=True) print(' M M OO OO LL BB B EE ',flush=True) print(' M M OOOO LLLLLL BBBBBBB EEEEEEE ',flush=True) - + print(flush=True) - print(' MOLECULAR BOOTSTRAP EMBEDDING',flush=True) + print(' MOLECULAR BOOTSTRAP EMBEDDING',flush=True) print(' BEn = ',self.be_type,flush=True) print('-----------------------------------------------------------', flush=True) print(flush=True) - + def initialize(self, eri_,compute_hf, restart=False): """ @@ -276,19 +276,19 @@ def initialize(self, eri_,compute_hf, restart=False): restart : bool, optional Whether to restart from a previous calculation, by default False. """ - from .helper import get_scfObj + from .helper import get_scfObj import h5py from pyscf import ao2mo from multiprocessing import Pool - + if compute_hf: E_hf = 0. - + # Create a file to store ERIs if not restart: file_eri = h5py.File(self.eri_file,'w') lentmp = len(self.edge_idx) for I in range(self.Nfrag): - + if lentmp: fobjs_ = Frags(self.fsites[I], I, edge=self.edge[I], eri_file=self.eri_file, @@ -301,9 +301,9 @@ def initialize(self, eri_,compute_hf, restart=False): edge_idx=[],center_idx=[],centerf_idx=[], efac=self.ebe_weight[I]) fobjs_.sd(self.W, self.lmo_coeff, self.Nocc) - + self.Fobjs.append(fobjs_) - + if not restart: # Transform ERIs for each fragment and store in the file # ERI Transform Decision Tree @@ -333,50 +333,50 @@ def initialize(self, eri_,compute_hf, restart=False): return NotImplementedError else: eri=None - + for fobjs_ in self.Fobjs: # Process each fragment eri = numpy.array(file_eri.get(fobjs_.dname)) dm_init = fobjs_.get_nsocc(self.S, self.C, self.Nocc, ncore=self.ncore) - + fobjs_.cons_h1(self.hcore) - + if not restart: eri = ao2mo.restore(8, eri, fobjs_.nao) - + fobjs_.cons_fock(self.hf_veff, self.S, self.hf_dm, eri_=eri) - + fobjs_.heff = numpy.zeros_like(fobjs_.h1) fobjs_.scf(fs=True, eri=eri) - + fobjs_.dm0 = numpy.dot( fobjs_._mo_coeffs[:,:fobjs_.nsocc], fobjs_._mo_coeffs[:,:fobjs_.nsocc].conj().T) *2. - + if compute_hf: - + eh1, ecoul, ef = fobjs_.energy_hf(return_e1=True) E_hf += fobjs_.ebe_hf if not restart: file_eri.close() - + if compute_hf: - + self.ebe_hf = E_hf+self.enuc+self.E_core hf_err = self.hf_etot - self.ebe_hf print('HF-in-HF error : {:>.4e} Ha'. format(hf_err), flush=True) if abs(hf_err)>1.e-5: print('WARNING!!! Large HF-in-HF energy error') - + print(flush=True) - + couti = 0 for fobj in self.Fobjs: fobj.udim = couti couti = fobj.set_udim(couti) - - def oneshot(self, solver='MP2', nproc=1, ompnum=4, calc_frag_energy=False, clean_eri=False, + + def oneshot(self, solver='MP2', nproc=1, ompnum=4, calc_frag_energy=False, clean_eri=False, scratch_dir=None, **solver_kwargs): """ Perform a one-shot bootstrap embedding calculation. @@ -396,10 +396,10 @@ def oneshot(self, solver='MP2', nproc=1, ompnum=4, calc_frag_energy=False, clean """ from .solver import be_func from .be_parallel import be_func_parallel - + self.scratch_dir = scratch_dir self.solver_kwargs = solver_kwargs - + print("Calculating Energy by Fragment? ", calc_frag_energy) if nproc == 1: rets = be_func(None, self.Fobjs, self.Nocc, solver, self.enuc, hf_veff=self.hf_veff, @@ -484,9 +484,9 @@ def read_heff(self, heff_file='bepotfile.h5'): for fobj in self.Fobjs: fobj.heff = filepot.get(fobj.dname) filepot.close() - - - + + + def initialize_pot(Nfrag, edge_idx): """ Initialize the potential array for bootstrap embedding. @@ -510,7 +510,7 @@ def initialize_pot(Nfrag, edge_idx): Initialized potential array with zeros. """ pot_=[] - + if not len(edge_idx) == 0: for I in range(Nfrag): for i in edge_idx[I]: @@ -519,7 +519,7 @@ def initialize_pot(Nfrag, edge_idx): if j>k: continue pot_.append(0.) - + pot_.append(0.) return pot_ diff --git a/molbe/misc.py b/molbe/misc.py index 6d35aa63..c52b632f 100644 --- a/molbe/misc.py +++ b/molbe/misc.py @@ -287,7 +287,7 @@ def be2puffin( Unrestricted vs restricted HF and CCSD, by default False from_chk: bool, optional Run calculation from converged RHF/UHF checkpoint. By default False - checkfile: string, optional + checkfile: string, optional if not None: - if from_chk: specify the checkfile to run the embedding calculation - if not from_chk: specify where to save the checkfile @@ -372,7 +372,7 @@ def be2puffin( if not hcore is None: mf.get_hcore = lambda *args: hcore_pyscf if not jk is None: mf.get_jk = lambda *args: jk_pyscf - if checkfile: + if checkfile: print("Saving checkfile to:", checkfile) mf.chkfile = checkfile time_pre_mf = time.time() @@ -427,13 +427,13 @@ def be2puffin( def print_energy(ecorr, e_V_Kapprox, e_F_dg, e_hf): - + # Print energy results print('-----------------------------------------------------', flush=True) print(' BE ENERGIES with cumulant-based expression', flush=True) - + print('-----------------------------------------------------', flush=True) print(' E_BE = E_HF + Tr(F del g) + Tr(V K_approx)', flush=True) @@ -444,5 +444,5 @@ def print_energy(ecorr, e_V_Kapprox, e_F_dg, e_hf): print(' Ecorr BE : {:>14.8f} Ha'.format(ecorr), flush=True) print('-----------------------------------------------------', flush=True) - + print(flush=True) diff --git a/molbe/pfrag.py b/molbe/pfrag.py index 62151f1d..e9c45540 100644 --- a/molbe/pfrag.py +++ b/molbe/pfrag.py @@ -15,13 +15,13 @@ class Frags: This class contains various functionalities required for managing and manipulating fragments for BE calculations. """ - + def __init__(self, fsites, ifrag, edge=None, center=None, edge_idx=None, center_idx=None, efac=None, eri_file='eri_file.h5', centerf_idx=None, unrestricted=False): - """Constructor function for `Frags` class. + """Constructor function for `Frags` class. Parameters ---------- @@ -46,7 +46,7 @@ def __init__(self, fsites, ifrag, edge=None, center=None, unrestricted : bool, optional unrestricted calculation, by default False """ - + self.fsites = fsites self.nfsites = len(fsites) self.TA = None @@ -58,19 +58,19 @@ def __init__(self, fsites, ifrag, edge=None, center=None, else: self.dname = 'f'+str(ifrag) self.nao = None - self.mo_coeffs = None - self._mo_coeffs = None + self.mo_coeffs = None + self._mo_coeffs = None self.nsocc = None self._mf = None self._mc = None - + # CCSD self.t1 = None self.t2 = None - + self.heff = None self.edge = edge - self.center = center + self.center = center self.edge_idx = edge_idx self.center_idx = center_idx self.centerf_idx = centerf_idx @@ -105,14 +105,14 @@ def sd(self, lao, lmo, nocc, norb=None, return_orb_count=False): nocc : int Number of occupied orbitals. norb : int, optional - Specify number of bath orbitals. + Specify number of bath orbitals. Used for UBE, where different number of alpha and beta orbitals Default is None, allowing orbitals to be chosen by threshold return_orb_count : bool, optional Retrun the number of orbitals in each space, for UBE use/ Default is False """ - + if return_orb_count: TA, n_f, n_b = schmidt_decomposition(lmo, nocc, self.fsites, norb=norb, return_orb_count=return_orb_count) else: @@ -133,11 +133,11 @@ def cons_h1(self, h1): h1 : numpy.ndarray One-electron Hamiltonian matrix. """ - + h1_tmp = functools.reduce(numpy.dot, (self.TA.T, h1, self.TA)) self.h1 = h1_tmp - + def cons_fock(self, hf_veff, S, dm, eri_=None): """ Construct the Fock matrix for the fragment. @@ -156,11 +156,11 @@ def cons_fock(self, hf_veff, S, dm, eri_=None): if eri_ is None: eri_ = get_eri(self.dname, self.TA.shape[1], ignore_symm=True, eri_file=self.eri_file) - + veff_ = get_veff(eri_, dm, S, self.TA, hf_veff) - self.veff = veff_.real + self.veff = veff_.real self.fock = self.h1 + veff_.real - + def get_nsocc(self, S, C, nocc,ncore=0): """ @@ -183,23 +183,23 @@ def get_nsocc(self, S, C, nocc,ncore=0): Projected density matrix. """ import scipy.linalg - + C_ = functools.reduce(numpy.dot,(self.TA.T, S, C[:,ncore:ncore+nocc])) P_ = numpy.dot(C_, C_.T) nsocc_ = numpy.trace(P_) nsocc = int(numpy.round(nsocc_)) - + try: mo_coeffs = scipy.linalg.svd(C_)[0] except: - + mo_coeffs = scipy.linalg.eigh(C_)[1][:,-nsocc:] - + self._mo_coeffs = mo_coeffs - self.nsocc = nsocc + self.nsocc = nsocc return P_ - - + + def scf(self, heff=None, fs=False, eri=None, dm0 = None, unrestricted=False, spin_ind=None): """ Perform self-consistent field (SCF) calculation for the fragment. @@ -219,7 +219,7 @@ def scf(self, heff=None, fs=False, eri=None, dm0 = None, unrestricted=False, spi spin_ind : int, optional Alpha (0) or beta (1) spin for unrestricted calculation, by default None """ - + import copy if self._mf is not None: self._mf = None if self._mc is not None: self._mc = None @@ -228,14 +228,14 @@ def scf(self, heff=None, fs=False, eri=None, dm0 = None, unrestricted=False, spi if eri is None: if unrestricted: dname = self.dname[spin_ind] - else: + else: dname = self.dname eri = get_eri(dname, self.nao, eri_file=self.eri_file) if dm0 is None: dm0 = numpy.dot( self._mo_coeffs[:,:self.nsocc], self._mo_coeffs[:,:self.nsocc].conj().T) *2. - + mf_ = get_scfObj(self.fock + heff, eri, self.nsocc, dm0 = dm0) if not fs: self._mf = mf_ @@ -243,21 +243,21 @@ def scf(self, heff=None, fs=False, eri=None, dm0 = None, unrestricted=False, spi else: self._mo_coeffs = mf_.mo_coeff.copy() mf_= None - + def update_heff(self,u, cout = None, return_heff=False, only_chem=False): """ Update the effective Hamiltonian for the fragment. """ - + heff_ = numpy.zeros_like(self.h1) if cout is None: cout = self.udim else: cout = cout - - for i,fi in enumerate(self.fsites): - if not any(i in sublist for sublist in self.edge_idx): + + for i,fi in enumerate(self.fsites): + if not any(i in sublist for sublist in self.edge_idx): heff_[i,i] -= u[-1] if only_chem: @@ -273,14 +273,14 @@ def update_heff(self,u, cout = None, return_heff=False, only_chem=False): for j in range(len(i)): for k in range(len(i)): if j>k :#or j==k: - continue - + continue + heff_[i[j], i[k]] = u[cout] heff_[i[k], i[j]] = u[cout] - + cout += 1 - - + + self.heff = heff_ if return_heff: if cout is None: @@ -297,29 +297,29 @@ def set_udim(self, cout): continue cout += 1 return cout - + def energy(self,rdm2s, eri=None, print_fragE=False): ## This function uses old energy expression and will be removed rdm2s = numpy.einsum("ijkl,pi,qj,rk,sl->pqrs", 0.5*rdm2s, - *([self.mo_coeffs]*4),optimize=True) - + *([self.mo_coeffs]*4),optimize=True) + e1 = 2.*numpy.einsum("ij,ij->i",self.h1[:self.nfsites], self._rdm1[:self.nfsites]) ec = numpy.einsum("ij,ij->i",self.veff[:self.nfsites],self._rdm1[:self.nfsites]) - + if self.TA.ndim == 3: jmax = self.TA[0].shape[1] else: jmax = self.TA.shape[1] - + if eri is None: r = h5py.File(self.eri_file,'r') eri = r[self.dname][()] r.close() - - + + e2 = numpy.zeros_like(e1) for i in range(self.nfsites): for j in range(jmax): @@ -327,10 +327,10 @@ def energy(self,rdm2s, eri=None, print_fragE=False): Gij = rdm2s[i,j,:jmax,:jmax].copy() Gij[numpy.diag_indices(jmax)] *= 0.5 Gij += Gij.T - + e2[i] += Gij[numpy.tril_indices(jmax)] @ eri[ij] - - e_ = e1+e2+ec + + e_ = e1+e2+ec etmp = 0. e1_ = 0. ec_ = 0. @@ -340,10 +340,10 @@ def energy(self,rdm2s, eri=None, print_fragE=False): e1_ += self.efac[0] * e1[i] ec_ += self.efac[0] * ec[i] e2_ += self.efac[0] * e2[i] - + print('BE Energy Frag-{:>3} {:>12.7f} {:>12.7f} {:>12.7f}; Total : {:>12.7f}'. format(self.dname, e1_, ec_, e2_, etmp)) - + self.ebe = etmp return (e1+e2+ec) @@ -376,35 +376,35 @@ def energy_hf(self, rdm_hf=None, mo_coeffs = None, eri=None, return_e1=False, un r.close() - + e2 = numpy.zeros_like(e1) for i in range(self.nfsites): for j in range(jmax): ij = i*(i+1)//2+j if i > j else j*(j+1)//2+i Gij = (2.*rdm_hf[i,j]*rdm_hf - - numpy.outer(rdm_hf[i], rdm_hf[j]))[:jmax,:jmax] + numpy.outer(rdm_hf[i], rdm_hf[j]))[:jmax,:jmax] Gij[numpy.diag_indices(jmax)] *= 0.5 - Gij += Gij.T + Gij += Gij.T if unrestricted: #unrestricted ERI file has 3 spin components: a, b, ab e2[i] += 0.5 * unrestricted_fac * Gij[numpy.tril_indices(jmax)] @ eri[spin_ind][ij] else: e2[i] += 0.5 * unrestricted_fac * Gij[numpy.tril_indices(jmax)] @ eri[ij] - e_ = e1+e2+ec + e_ = e1+e2+ec etmp = 0. for i in self.efac[1]: etmp += self.efac[0]*e_[i] - + self.ebe_hf = etmp - + if return_e1: e_h1 = 0. e_coul = 0. for i in self.efac[1]: - + e_h1 += self.efac[0]*e1[i] e_coul += self.efac[0]*(e2[i]+ec[i]) return(e_h1,e_coul, e1+e2+ec) - + return e1+e2+ec - + diff --git a/molbe/rdm.py b/molbe/rdm.py index b0906f5b..72902668 100644 --- a/molbe/rdm.py +++ b/molbe/rdm.py @@ -7,7 +7,7 @@ def rdm1_fullbasis(self, return_ao=True, only_rdm1=False, return_RDM2=True, print_energy=False): """ Compute the one-particle and two-particle reduced density matrices (RDM1 and RDM2). - + Parameters: ----------- return_ao : bool, optional @@ -22,7 +22,7 @@ def rdm1_fullbasis(self, return_ao=True, only_rdm1=False, Whether to return the two-particle RDM (RDM2). Default is True. print_energy : bool, optional Whether to print the energy contributions. Default is False. - + Returns: -------- rdm1AO : numpy.ndarray @@ -41,13 +41,13 @@ def rdm1_fullbasis(self, return_ao=True, only_rdm1=False, from pyscf import scf, ao2mo # Copy the molecular orbital coefficients - C_mo = self.C.copy() + C_mo = self.C.copy() nao, nmo = C_mo.shape # Initialize density matrices for atomic orbitals (AO) rdm1AO = numpy.zeros((nao, nao)) rdm2AO = numpy.zeros((nao, nao, nao, nao)) - + for fobjs in self.Fobjs: if return_RDM2: # Adjust the one-particle reduced density matrix (RDM1) @@ -58,18 +58,18 @@ def rdm1_fullbasis(self, return_ao=True, only_rdm1=False, dm_nc = numpy.einsum('ij,kl->ijkl', drdm1, drdm1, dtype=numpy.float64, optimize=True) - \ 0.5*numpy.einsum('ij,kl->iklj', drdm1, drdm1, dtype=numpy.float64,optimize=True) fobjs.__rdm2 -= dm_nc - - # Generate the projection matrix + + # Generate the projection matrix cind = [ fobjs.fsites[i] for i in fobjs.efac[1]] Pc_ = fobjs.TA.T @ self.S @ self.W[:, cind] @ self.W[:, cind].T @ self.S @ fobjs.TA if not only_rdm2: # Compute RDM1 in the localized orbital (LO) basis and transform to AO basis - rdm1_eo = fobjs.mo_coeffs @ fobjs.__rdm1 @ fobjs.mo_coeffs.T + rdm1_eo = fobjs.mo_coeffs @ fobjs.__rdm1 @ fobjs.mo_coeffs.T rdm1_center = Pc_ @ rdm1_eo rdm1_ao = fobjs.TA @ rdm1_center @ fobjs.TA.T rdm1AO += rdm1_ao - + if not only_rdm1: # Transform RDM2 to AO basis rdm2s = numpy.einsum("ijkl,pi,qj,rk,sl->pqrs", fobjs.__rdm2, @@ -100,7 +100,7 @@ def rdm1_fullbasis(self, return_ao=True, only_rdm1=False, rdm2LO = numpy.einsum("ijkl,pi,qj,rk,sl->pqrs", rdm2AO, CloT_S, CloT_S, CloT_S, CloT_S, optimize=True) - + if not only_rdm2: # Symmetrize RDM1 rdm1AO = (rdm1AO + rdm1AO.T)/2. @@ -118,23 +118,23 @@ def rdm1_fullbasis(self, return_ao=True, only_rdm1=False, Eh1 = numpy.einsum('ij,ij', self.hcore, rdm1AO, optimize=True) eri = ao2mo.restore(1,self.mf._eri, self.mf.mo_coeff.shape[1]) E2 = 0.5*numpy.einsum('pqrs,pqrs', eri,rdm2AO, optimize=True) - print(flush=True) + print(flush=True) print('-----------------------------------------------------', flush=True) print(' BE ENERGIES with cumulant-based expression', flush=True) - + print('-----------------------------------------------------', flush=True) - + print(' 1-elec E : {:>15.8f} Ha'.format(Eh1), flush=True) print(' 2-elec E : {:>15.8f} Ha'.format(E2), flush=True) E_tot = Eh1+E2+self.E_core + self.enuc print(' E_BE : {:>15.8f} Ha'.format(E_tot), flush=True) print(' Ecorr BE : {:>15.8f} Ha'.format((E_tot)-self.ebe_hf), flush=True) print('-----------------------------------------------------', - flush=True) + flush=True) print(flush=True) - + if only_rdm1: if return_ao: return rdm1AO @@ -170,25 +170,25 @@ def compute_energy_full(self, approx_cumulant=False, use_full_rdm=False, return_ Returns ------- tuple of numpy.ndarray or None - If `return_rdm` is True, returns a tuple containing the one-particle and two-particle + If `return_rdm` is True, returns a tuple containing the one-particle and two-particle reduced density matrices (RDM1 and RDM2). Otherwise, returns None. Notes ----- - This function computes the total energy in the full basis, with options to use - approximate or true cumulants, and to return the reduced density matrices (RDMs). The + This function computes the total energy in the full basis, with options to use + approximate or true cumulants, and to return the reduced density matrices (RDMs). The energy components are printed as part of the function's output. """ - + from pyscf import scf, ao2mo # Compute the one-particle reduced density matrix (RDM1) and the cumulant (Kumul) in the full basis rdm1f, Kumul, rdm1_lo, rdm2_lo = self.rdm1_fullbasis(return_lo=True, return_RDM2=False) - + if not approx_cumulant: # Compute the true two-particle reduced density matrix (RDM2) if not using approximate cumulant Kumul_T = self.rdm1_fullbasis(only_rdm2=True) - + if return_rdm: # Construct the full RDM2 from RDM1 RDM2_full = numpy.einsum('ij,kl->ijkl', rdm1f, rdm1f, dtype=numpy.float64, optimize=True) - \ @@ -227,15 +227,15 @@ def compute_energy_full(self, approx_cumulant=False, use_full_rdm=False, return_ if not approx_cumulant: # Compute the true two-electron energy if not using approximate cumulant EKumul_T = numpy.einsum('pqrs,pqrs', eri,Kumul_T, optimize=True) - + if use_full_rdm and return_rdm: # Compute the full two-electron energy using the full RDM2 E2 = numpy.einsum('pqrs,pqrs', eri,RDM2_full, optimize=True) # Compute the approximate BE total energy - EKapprox = self.ebe_hf + Eh1_dg + Eveff_dg + EKumul/2. + EKapprox = self.ebe_hf + Eh1_dg + Eveff_dg + EKumul/2. self.ebe_tot = EKapprox - + if not approx_cumulant: # Compute the true BE total energy if not using approximate cumulant EKtrue = Eh1 + EVeff/2. + EKumul_T/2. + self.enuc + self.E_core @@ -245,7 +245,7 @@ def compute_energy_full(self, approx_cumulant=False, use_full_rdm=False, return_ print('-----------------------------------------------------', flush=True) print(' BE ENERGIES with cumulant-based expression', flush=True) - + print('-----------------------------------------------------', flush=True) print(' E_BE = E_HF + Tr(F del g) + Tr(V K_approx)', flush=True) @@ -254,7 +254,7 @@ def compute_energy_full(self, approx_cumulant=False, use_full_rdm=False, return_ print(' Tr(V K_aprrox) : {:>14.8f} Ha'.format(EKumul/2.), flush=True) print(' E_BE : {:>14.8f} Ha'.format(EKapprox), flush=True) print(' Ecorr BE : {:>14.8f} Ha'.format(EKapprox-self.ebe_hf), flush=True) - + if not approx_cumulant: print(flush=True) print(' E_BE = Tr(F[g] g) + Tr(V K_true)', flush=True) @@ -262,7 +262,7 @@ def compute_energy_full(self, approx_cumulant=False, use_full_rdm=False, return_ print(' Tr(Veff[g] g) : {:>14.8f} Ha'.format(EVeff/2.), flush=True) print(' Tr(V K_true) : {:>14.8f} Ha'.format(EKumul_T/2.), flush=True) print(' E_BE : {:>14.8f} Ha'.format(EKtrue), flush=True) - if use_full_rdm and return_rdm: + if use_full_rdm and return_rdm: print(' E(g+G) : {:>14.8f} Ha'.format(Eh1 + 0.5*E2 + self.E_core + self.enuc), flush=True) print(' Ecorr BE : {:>14.8f} Ha'.format(EKtrue-self.ebe_hf), flush=True) @@ -270,9 +270,9 @@ def compute_energy_full(self, approx_cumulant=False, use_full_rdm=False, return_ print(' True - approx : {:>14.4e} Ha'.format(EKtrue-EKapprox)) print('-----------------------------------------------------', flush=True) - + print(flush=True) - + # Return the RDMs if requested if return_rdm: return(rdm1f, RDM2_full) - + diff --git a/molbe/solver.py b/molbe/solver.py index 5c0e3cec..15714554 100644 --- a/molbe/solver.py +++ b/molbe/solver.py @@ -12,7 +12,7 @@ def be_func(pot, Fobjs, Nocc, solver, enuc, hf_veff=None, only_chem = False, nproc=4, hci_pt=False, hci_cutoff=0.001, ci_coeff_cutoff = None, select_cutoff=None, eeval=False, ereturn=False, frag_energy=False, relax_density=False, - return_vec=False, ecore=0., ebe_hf=0., be_iter=None, use_cumulant=True, + return_vec=False, ecore=0., ebe_hf=0., be_iter=None, use_cumulant=True, scratch_dir=None, **solver_kwargs): """ Perform bootstrap embedding calculations for each fragment. @@ -81,14 +81,14 @@ def be_func(pot, Fobjs, Nocc, solver, enuc, hf_veff=None, only_chem=only_chem) else: heff_ = None - + # Compute the one-electron Hamiltonian h1_ = fobj.fock + fobj.heff # Perform SCF calculation fobj.scf() # Solve using the specified solver - if solver=='MP2': + if solver=='MP2': fobj._mc = solve_mp2(fobj._mf, mo_energy=fobj._mf.mo_energy) elif solver=='CCSD': if rdm_return: @@ -102,47 +102,47 @@ def be_func(pot, Fobjs, Nocc, solver, enuc, hf_veff=None, mo_energy=fobj._mf.mo_energy, rdm_return=False) rdm1_tmp = make_rdm1_ccsd_t1(fobj.t1) - + elif solver=='FCI': from .helper import get_eri import scipy - + mc = fci.FCI(fobj._mf, fobj._mf.mo_coeff) - efci, civec = mc.kernel() + efci, civec = mc.kernel() rdm1_tmp = mc.make_rdm1(civec, mc.norb, mc.nelec) elif solver=='HCI': from pyscf import hci from .helper import get_eri # pilot pyscf.hci only in old versions - + nao, nmo = fobj._mf.mo_coeff.shape eri = ao2mo.kernel(fobj._mf._eri, fobj._mf.mo_coeff, aosym='s4', compact=False).reshape(4*((nmo),)) - + ci_ = hci.SCI(fobj._mf.mol) if select_cutoff is None and ci_coeff_cutoff is None: select_cutoff = hci_cutoff ci_coeff_cutoff = hci_cutoff elif select_cutoff is None or ci_coeff_cutoff is None: sys.exit() - + ci_.select_cutoff = select_cutoff ci_.ci_coeff_cutoff = ci_coeff_cutoff - - nelec = (fobj.nsocc, fobj.nsocc) + + nelec = (fobj.nsocc, fobj.nsocc) h1_ = fobj.fock+fobj.heff h1_ = functools.reduce(numpy.dot, (fobj._mf.mo_coeff.T, h1_, fobj._mf.mo_coeff)) eci, civec = ci_.kernel(h1_, eri, nmo, nelec) civec = numpy.asarray(civec) - - + + (rdm1a_, rdm1b_), (rdm2aa, rdm2ab, rdm2bb) = ci_.make_rdm12s(civec, nmo, nelec) rdm1_tmp = rdm1a_ + rdm1b_ rdm2s = rdm2aa + rdm2ab + rdm2ab.transpose(2,3,0,1) + rdm2bb - + elif solver=='SHCI': from pyscf.shciscf import shci @@ -155,13 +155,13 @@ def be_func(pot, Fobjs, Nocc, solver, enuc, hf_veff=None, if not os.path.isdir(tmp): os.system('mkdir -p '+tmp) nao, nmo = fobj._mf.mo_coeff.shape - + nelec = (fobj.nsocc, fobj.nsocc) mch = shci.SHCISCF(fobj._mf, nmo, nelec, orbpath=fobj.dname) mch.fcisolver.mpiprefix = 'mpirun -np '+str(nproc) if hci_pt: mch.fcisolver.stochastic = False - mch.fcisolver.epsilon2 = hci_cutoff + mch.fcisolver.epsilon2 = hci_cutoff else: mch.fcisolver.stochastic = True # this is for PT and doesnt add PT to rdm mch.fcisolver.nPTiter = 0 @@ -171,7 +171,7 @@ def be_func(pot, Fobjs, Nocc, solver, enuc, hf_veff=None, mch.fcisolver.scratchDirectory = scratch_dir mch.mc1step() rdm1_tmp, rdm2s = mch.fcisolver.make_rdm12(0, nmo, nelec) - + elif solver == 'SCI': from pyscf import cornell_shci from pyscf import ao2mo, mcscf @@ -191,9 +191,9 @@ def be_func(pot, Fobjs, Nocc, solver, enuc, hf_veff=None, ci.config['get_2rdm_csv'] = True ci.kernel(h1, eri, nmo, nelec) rdm1_tmp, rdm2s = ci.make_rdm12(0,nmo,nelec) - + elif solver in ['block2', 'DMRG','DMRGCI','DMRGSCF']: - + solver_kwargs_ = solver_kwargs.copy() if scratch_dir is None and be_var.CREATE_SCRATCH_DIR: tmp = os.path.join(be_var.SCRATCH, str(os.getpid()), str(fobj.dname)) @@ -201,15 +201,15 @@ def be_func(pot, Fobjs, Nocc, solver, enuc, hf_veff=None, tmp = os.path.join(scratch_dir, str(os.getpid()), str(fobj.dname)) if not os.path.isdir(tmp): os.system('mkdir -p '+tmp) - + try: - rdm1_tmp, rdm2s = solve_block2(fobj._mf, - fobj.nsocc, + rdm1_tmp, rdm2s = solve_block2(fobj._mf, + fobj.nsocc, frag_scratch = tmp, **solver_kwargs_) except Exception as inst: print(f"Fragment DMRG solver failed with Exception: {type(inst)}\n", inst, flush=True) - + finally: if solver_kwargs_.pop('force_cleanup', False): os.system('rm -r '+ os.path.join(tmp,'F.*')) @@ -253,22 +253,22 @@ def be_func(pot, Fobjs, Nocc, solver, enuc, hf_veff=None, rdm2s -= nc fobj.__rdm2 = rdm2s.copy() if frag_energy or eeval: - # Find the energy of a given fragment, with the cumulant definition. + # Find the energy of a given fragment, with the cumulant definition. # Return [e1, e2, ec] as e_f and add to the running total_e. e_f = get_frag_energy(fobj.mo_coeffs, fobj.nsocc, fobj.nfsites, fobj.efac, fobj.TA, fobj.h1, hf_veff, rdm1_tmp, rdm2s, fobj.dname, eri_file=fobj.eri_file, veff0=fobj.veff0) total_e = [sum(x) for x in zip(total_e, e_f)] fobj.energy_hf() - + if frag_energy or eeval: Ecorr = sum(total_e) - + if frag_energy: return (Ecorr, total_e) - + ernorm, ervec = solve_error(Fobjs,Nocc, only_chem=only_chem) - + if return_vec: return (ernorm, ervec, [Ecorr, total_e]) @@ -277,9 +277,9 @@ def be_func(pot, Fobjs, Nocc, solver, enuc, hf_veff=None, return ernorm - + def be_func_u(pot, Fobjs, solver, enuc, hf_veff=None, - eeval=False, ereturn=False, frag_energy=True, + eeval=False, ereturn=False, frag_energy=True, relax_density=False, ecore=0., ebe_hf=0., scratch_dir=None, use_cumulant=True, frozen=False): """ @@ -428,7 +428,7 @@ def solve_error(Fobjs, Nocc, only_chem=False): Error vector. """ import math - + err_edge = [] err_chempot = 0. @@ -439,11 +439,11 @@ def solve_error(Fobjs, Nocc, only_chem=False): err_chempot += fobj._rdm1[i,i] err_chempot /= Fobjs[0].unitcell_nkpt err = err_chempot - Nocc - + return abs(err), numpy.asarray([err]) # Compute edge and chemical potential errors - for fobj in Fobjs: + for fobj in Fobjs: #match rdm-edge for edge in fobj.edge_idx: @@ -452,36 +452,36 @@ def solve_error(Fobjs, Nocc, only_chem=False): if j_>k_: continue err_edge.append(fobj._rdm1[edge[j_], edge[k_]]) - #chem potential + #chem potential for i in fobj.efac[1]: err_chempot += fobj._rdm1[i,i] - + err_chempot /= Fobjs[0].unitcell_nkpt err_edge.append(err_chempot) # far-end edges are included as err_chempot # Compute center errors err_cen = [] - for findx, fobj in enumerate(Fobjs): + for findx, fobj in enumerate(Fobjs): # Match RDM for centers - for cindx, cens in enumerate(fobj.center_idx): + for cindx, cens in enumerate(fobj.center_idx): lenc = len(cens) for j_ in range(lenc): for k_ in range(lenc): if j_>k_: - continue + continue err_cen.append(Fobjs[fobj.center[cindx]]._rdm1[cens[j_], - cens[k_]]) - + cens[k_]]) + err_cen.append(Nocc) err_edge = numpy.array(err_edge) - err_cen = numpy.array(err_cen) + err_cen = numpy.array(err_cen) # Compute the error vector err_vec = err_edge - err_cen # Compute the norm of the error vector norm_ = numpy.mean(err_vec * err_vec)**0.5 - + return norm_, err_vec def solve_mp2(mf, frozen=None, mo_coeff=None, mo_occ=None, mo_energy=None): @@ -627,7 +627,7 @@ def solve_ccsd(mf, frozen=None, mo_coeff=None,relax=False, use_cumulant=False, w def solve_block2(mf:object, nocc:int, frag_scratch:str = None, **solver_kwargs): """ DMRG fragment solver using the pyscf.dmrgscf wrapper. - + Parameters ---------- mf: pyscf.scf.hf.RHF @@ -650,19 +650,19 @@ def solve_block2(mf:object, nocc:int, frag_scratch:str = None, **solver_kwargs): Sweep index at which to transition to one-dot DMRG algorithm. All sweeps prior to this will use the two-dot algorithm. block_extra_keyword: list(str), optional Other keywords to be passed to block2. See: https://block2.readthedocs.io/en/latest/user/keywords.html - + Returns ------- rdm1: numpy.ndarray 1-Particle reduced density matrix for fragment. rdm2: numpy.ndarray 2-Particle reduced density matrix for fragment. - + Other Parameters ---------------- schedule_kwargs: dict, optional Dictionary containing DMRG scheduling parameters to be passed to block2. - + e.g. The default schedule used here would be equivalent to the following: schedule_kwargs = { 'scheduleSweeps': [0, 10, 20, 30, 40, 50], @@ -670,7 +670,7 @@ def solve_block2(mf:object, nocc:int, frag_scratch:str = None, **solver_kwargs): 'scheduleTols': [1e-5,1e-5, 1e-6, 1e-6, 1e-8, 1e-8], 'scheduleNoises': [0.01, 0.01, 0.001, 0.001, 1e-4, 0.0], } - + Raises ------ @@ -692,48 +692,48 @@ def solve_block2(mf:object, nocc:int, frag_scratch:str = None, **solver_kwargs): root = solver_kwargs.pop("root", 0) block_extra_keyword = solver_kwargs.pop("block_extra_keyword", ['fiedler']) schedule_kwargs = solver_kwargs.pop("schedule_kwargs", {}) - + if norb <= 2: block_extra_keyword = ['noreorder'] #Other reordering algorithms explode if the network is too small. - + if lo_method is None: orbs = mf.mo_coeff elif isinstance(lo_method, str): raise NotImplementedError("Localization within the fragment+bath subspace is currently not supported.") - + mc = mcscf.CASCI(mf, norb, nelec) mc.fcisolver = dmrgscf.DMRGCI(mf.mol) ###Sweep scheduling - mc.fcisolver.scheduleSweeps = schedule_kwargs.pop("scheduleSweeps", - [(1*max_iter)//6, - (2*max_iter)//6, - (3*max_iter)//6, - (4*max_iter)//6, - (5*max_iter)//6, + mc.fcisolver.scheduleSweeps = schedule_kwargs.pop("scheduleSweeps", + [(1*max_iter)//6, + (2*max_iter)//6, + (3*max_iter)//6, + (4*max_iter)//6, + (5*max_iter)//6, max_iter] ) - mc.fcisolver.scheduleMaxMs = schedule_kwargs.pop("scheduleMaxMs", - [startM if (startMiklj',hf_dm, hf_dm) + \ numpy.einsum('ij,kl->iklj',hf_dm, del_rdm1) + \ numpy.einsum('ij,kl->iklj',del_rdm1, hf_dm))*0.5 - + rdm2 -= nc - + return rdm1, rdm2 def solve_uccsd(mf, eris_inp, frozen=None, mo_coeff=None, relax=False, @@ -780,7 +780,7 @@ def solve_uccsd(mf, eris_inp, frozen=None, mo_coeff=None, relax=False, ---------- mf : pyscf.scf.hf.UHF Mean-field object from PySCF. Constructed with make_uhf_obj - eris_inp : + eris_inp : Custom fragment ERIs object frozen : list or int, optional List of frozen orbitals or number of frozen core orbitals. Defaults to None. @@ -897,7 +897,7 @@ def schmidt_decomposition(mo_coeff, nocc, Frag_sites, cinv = None, rdm=None, no Specifies number of bath orbitals. Used for UBE to make alpha and beta spaces the same size. Defaults to None return_orb_count : bool, optional - Return more information about the number of orbitals. Used in UBE. + Return more information about the number of orbitals. Used in UBE. Defaults to False Returns @@ -917,12 +917,12 @@ def schmidt_decomposition(mo_coeff, nocc, Frag_sites, cinv = None, rdm=None, no # Compute the reduced density matrix (RDM) if not provided if not mo_coeff is None: - C = mo_coeff[:,:nocc] + C = mo_coeff[:,:nocc] if rdm is None: Dhf = numpy.dot(C, C.T) if not cinv is None: Dhf = functools.reduce(numpy.dot, - (cinv, Dhf, cinv.conj().T)) + (cinv, Dhf, cinv.conj().T)) else: Dhf = rdm @@ -957,16 +957,16 @@ def schmidt_decomposition(mo_coeff, nocc, Frag_sites, cinv = None, rdm=None, no Bidx.append(i) else: for i in range(len(Eval)): - if thres < numpy.abs(Eval[i]) < 1.0 - thres: + if thres < numpy.abs(Eval[i]) < 1.0 - thres: Bidx.append(i) # Initialize the transformation matrix (TA) TA = numpy.zeros([Tot_sites, len(Frag_sites) + len(Bidx)]) TA[Frag_sites, :len(Frag_sites)] = numpy.eye(len(Frag_sites)) # Fragment part TA[Env_sites1,len(Frag_sites):] = Evec[:,Bidx] # Environment part - + if return_orb_count: - # return TA, norbs_frag, norbs_bath + # return TA, norbs_frag, norbs_bath return TA, Frag_sites1.shape[0], len(Bidx) else: return TA diff --git a/molbe/ube.py b/molbe/ube.py index 67ed5dcc..9d4547e7 100644 --- a/molbe/ube.py +++ b/molbe/ube.py @@ -141,7 +141,7 @@ def __init__( jobid = str(os.environ['SLURM_JOB_ID']) except: jobid = '' - if not be_var.SCRATCH=='': + if not be_var.SCRATCH=='': self.scratch_dir = be_var.SCRATCH+str(jobid) os.system('mkdir -p '+self.scratch_dir) else: @@ -317,7 +317,7 @@ def initialize(self, eri_, compute_hf): print("| Fragment | Nocc | Fragment Orbs | Bath Orbs | Schmidt Space |", flush=True) print("____________________________________________________________________", flush=True) for I in range(self.Nfrag): - print('| {:>2} | ({:>3},{:>3}) | ({:>3},{:>3}) | ({:>3},{:>3}) | ({:>3},{:>3}) |'.format(I, all_noccs[I][0],all_noccs[I][1], + print('| {:>2} | ({:>3},{:>3}) | ({:>3},{:>3}) | ({:>3},{:>3}) | ({:>3},{:>3}) |'.format(I, all_noccs[I][0],all_noccs[I][1], orb_count_a[I][0], orb_count_b[I][0], orb_count_a[I][1], orb_count_b[I][1], orb_count_a[I][0]+orb_count_a[I][1], @@ -373,7 +373,7 @@ def oneshot(self, solver="UCCSD", nproc=1, ompnum=4, calc_frag_energy=False, cle relax_density=False, frag_energy=calc_frag_energy, frozen=self.frozen_core, - nproc=nproc, + nproc=nproc, ompnum=ompnum) @@ -386,7 +386,7 @@ def oneshot(self, solver="UCCSD", nproc=1, ompnum=4, calc_frag_energy=False, cle self.ebe_tot = E + self.uhf_full_e print("Total Energy : {:>12.8f} Ha".format((self.ebe_tot), flush=True)) print("Corr Energy : {:>12.8f} Ha".format((E), flush=True)) - + if clean_eri == True: try: os.remove(self.eri_file) diff --git a/tests/chempot_molBE_test.py b/tests/chempot_molBE_test.py index df375592..447189a0 100644 --- a/tests/chempot_molBE_test.py +++ b/tests/chempot_molBE_test.py @@ -20,7 +20,7 @@ def test_h8_sto3g_ben(self): mol.build() self.molecular_restricted_test(mol, 'be2', 'H8 (BE2)', 'hchain_simple', -4.30628355, only_chem = True) self.molecular_restricted_test(mol, 'be3', 'H8 (BE3)', 'hchain_simple', -4.30649890, only_chem = True) - + @unittest.skipIf(os.getenv("GITHUB_ACTIONS") == "true", "Skip expensive tests on Github Actions") def test_octane_sto3g_ben(self): # Octane, STO-3G diff --git a/tests/dm_molBE_test.py b/tests/dm_molBE_test.py index 053f12a7..4387a280 100644 --- a/tests/dm_molBE_test.py +++ b/tests/dm_molBE_test.py @@ -18,7 +18,7 @@ def test_h8_sto3g_ben_trustRegion(self): mol.charge = 0.; mol.spin = 0. mol.build() self.molecular_QN_test(mol, 'be2', 'H8 (BE2)', 'hchain_simple', only_chem = False) - + def molecular_QN_test(self, mol, be_type, test_name, frag_type, delta = 1e-6, only_chem = True): mf = scf.RHF(mol); mf.max_cycle = 100; mf.kernel() fobj = fragpart(frag_type=frag_type, be_type=be_type, mol = mol) diff --git a/tests/dmrg_molBE_test.py b/tests/dmrg_molBE_test.py index 6bd2a8f1..2e2bdf12 100644 --- a/tests/dmrg_molBE_test.py +++ b/tests/dmrg_molBE_test.py @@ -1,5 +1,5 @@ """ -This script tests the QuEmb to block2 interface for performing ground state BE-DMRG. +This script tests the QuEmb to block2 interface for performing ground state BE-DMRG. Author(s): Shaun Weatherly """ @@ -14,7 +14,7 @@ class TestBE_DMRG(unittest.TestCase): from pyscf import dmrgscf except ImportError: dmrgscf=None - + @unittest.skipIf(dmrgscf is None, "Module `pyscf.dmrgscf` not imported correctly.") def test_h8_sto3g_pipek(self): mol = gto.M() @@ -30,7 +30,7 @@ def molecular_DMRG_test(self, mol, be_type, maxM, test_name, frag_type, target, mf = scf.RHF(mol); mf.kernel() fobj = fragpart(frag_type=frag_type, be_type=be_type, mol=mol) mybe = BE(mf, fobj, lo_method='pipek', pop_method='lowdin') - mybe.oneshot(solver='block2', + mybe.oneshot(solver='block2', scratch_dir=str(tmp), maxM=int(maxM), maxIter=30, diff --git a/tests/hf-in-hf_BE_test.py b/tests/hf-in-hf_BE_test.py index 09f4dc82..0fa5fdeb 100644 --- a/tests/hf-in-hf_BE_test.py +++ b/tests/hf-in-hf_BE_test.py @@ -19,7 +19,7 @@ def test_h8_sto3g_ben(self): self.molecular_restricted_test(mol, 'be1', 'H8 (BE1)') self.molecular_restricted_test(mol, 'be2', 'H8 (BE2)') self.molecular_restricted_test(mol, 'be3', 'H8 (BE3)') - + def test_h8_ccpvdz_ben(self): # Linear Equidistant (r=1Å) H8 Chain, cc-pVDZ mol = gto.M() @@ -41,7 +41,7 @@ def test_octane_sto3g_ben(self): self.molecular_restricted_test(mol, 'be1', 'Octane (BE1)') self.molecular_restricted_test(mol, 'be2', 'Octane (BE2)') self.molecular_restricted_test(mol, 'be3', 'Octane (BE3)') - + def molecular_restricted_test(self, mol, be_type, test_name, delta = 1e-5): mf = scf.RHF(mol); mf.kernel() fobj = fragpart(frag_type='autogen', be_type=be_type, mol = mol) diff --git a/tests/xyz/hexene.xyz b/tests/xyz/hexene.xyz index bbe63b0e..ddbf3d7a 100644 --- a/tests/xyz/hexene.xyz +++ b/tests/xyz/hexene.xyz @@ -1,20 +1,20 @@ 18 -C -5.6502267899 0.7485927383 -0.0074809907 -C -4.5584842828 1.7726977952 -0.0418619714 -H -5.5515181382 0.0602177800 -0.8733001951 -H -6.6384111226 1.2516490350 -0.0591711493 -H -5.5928112720 0.1649656434 0.9355613930 -C -3.2701647911 1.4104028701 0.0085107804 -C -2.1789947571 2.4456245961 -0.0265301736 -C -0.7941361337 1.7933691827 0.0427863465 -H -2.3064879754 3.1340763205 0.8376047229 -H -2.2652945230 3.0294980525 -0.9691823088 -C 0.3121144607 2.8438276122 0.0072105803 -H -0.6626010212 1.1042272672 -0.8201573062 -H -0.7037327970 1.2086599200 0.9845037688 -H 0.2581952957 3.4296615741 -0.9351274363 -H 1.3021608456 2.3437050700 0.0587054221 -H 0.2170130582 3.5342406229 0.8723087816 -H -4.8264911382 2.8238751191 -0.1088423637 -H -3.0132059318 0.3554469837 0.0754505007 +C -5.6502267899 0.7485927383 -0.0074809907 +C -4.5584842828 1.7726977952 -0.0418619714 +H -5.5515181382 0.0602177800 -0.8733001951 +H -6.6384111226 1.2516490350 -0.0591711493 +H -5.5928112720 0.1649656434 0.9355613930 +C -3.2701647911 1.4104028701 0.0085107804 +C -2.1789947571 2.4456245961 -0.0265301736 +C -0.7941361337 1.7933691827 0.0427863465 +H -2.3064879754 3.1340763205 0.8376047229 +H -2.2652945230 3.0294980525 -0.9691823088 +C 0.3121144607 2.8438276122 0.0072105803 +H -0.6626010212 1.1042272672 -0.8201573062 +H -0.7037327970 1.2086599200 0.9845037688 +H 0.2581952957 3.4296615741 -0.9351274363 +H 1.3021608456 2.3437050700 0.0587054221 +H 0.2170130582 3.5342406229 0.8723087816 +H -4.8264911382 2.8238751191 -0.1088423637 +H -3.0132059318 0.3554469837 0.0754505007