Source code for libra_py.workflows.nbra.mapping

#*********************************************************************************
#* Copyright (C) 2018-2019 Brendan A. Smith, Wei Li, Alexey V. Akimov
#*
#* This file is distributed under the terms of the GNU General Public License
#* as published by the Free Software Foundation, either version 2 of
#* the License, or (at your option) any later version.
#* See the file LICENSE in the root directory of this distribution
#* or <http://www.gnu.org/licenses/>.
#*
#*********************************************************************************/
"""

.. module:: mapping
   :platform: Unix, Windows
   :synopsis: This module implements the methods to map single-electron properties (e.g. KS basis)
       to many-electron ones (e.g. Slater Determinat basis)

.. moduleauthor:: Brendan Smith, Wei Li, and Alexey V. Akimov

"""

import os
import sys
import math

# Fisrt, we add the location of the library to test to the PYTHON path
if sys.platform=="cygwin":
    from cyglibra_core import *
elif sys.platform=="linux" or sys.platform=="linux2":
    from liblibra_core import *

#from libra_py import *


[docs]def sd2indx(inp,nbasis, do_sort=False): """ This function maps a list of integers defining the occupied spin-orbitals in a SD onto the global indices of the orbitals that are being occupied. Args: inp ( list of ints ): indices of the occupied spin-orbitals. Indexing starts with 1, not 0! Positive number ```n``` correspond to an alpha electron occupying the orbital (whether it is an alpha-spatial or beta-spatial component) ```n``` Negative number ```-n``` correspond to a beta electron occupying the orbital (whether it is an alpha-spatial or beta-spatial component) ```n``` nbasis ( int ): the total number of orbitals orbitals (both alpha- and beta- spatial components ) in the selected active space [currently not really used!] do_sort ( Boolean ): the flag to tell whether the indices should be sorted in a particular order: - True - for new scheme (needed for the SAC) [default] - use with False for Pyxaid mapping! Returns: list of ints: the indices of the orbitals occupied in this SD (different convention) Example: Assume out active space consists of 4 orbitals: |1^a>, |1^b>, |2^a>, |2^b>. Here, the superscript indicates the spin-component of the given orbital. In the step2 they are grouped as |1^a>, |2^a>, |1^b>, |2^b> And indexed (overall) as: 1 2 3 4 (as used in the input) 0 1 2 3 (as used in the output) Note: this indexing convention assumes starting from 1, not from 0! So the ```nbasis``` should be = 4 Input SD = [1, -3] means that the orbital |1^a> is occupied by an alpha electron |a> and the orbital |1^b> is occupied by a beta electron |b>, so the SD is: SD = 1/sqrt(2) * | |1^a>|a> |1^b>|b> | This function will return the global indices of the occupied orbitals in the set of active orbitals. That is the conversion will be: [1, -3] -> [0, 2] Note that the output indices start from 0 """ sz = len(inp) spat = [0] * sz for i in range(0,sz): # alpha if inp[i] > 0: spat[i] = inp[i] - 1 # beta else: spat[i] = (abs(inp[i])) - 1 #+ nbasis/2 - 1 # Rearrange in ascending order: this is needed for # a consistency of the final results among different orbitals # Warning! But the reordering messes up the mapping procedure, so it # is better not to have it. Note sure what effect it might have on spin-adaptation # though out = list(spat) if do_sort==True: out = sorted(spat) return out
[docs]def energy_arb(SD, e): """Computes the energy of a Slater Determinat (SD) state Args: SD ( list of ints ): indices of the orbitals occupied in this SD SeeAlso: ```inp``` in the ```sd2indx(inp,nbasis)``` function e ( MATRIX(nbasis, nbasis) ): energies of all 1-electron levels (alpha- and beta- components) Returns: double: the sum of the orbital energies for the occupied spin-orbitals - This is the approximation of the SD energy """ nbasis = e.num_of_rows sd = sd2indx(SD,nbasis) res = 0.0+0.0j for i in sd: res = res + e.get(i,i) return res
[docs]def energy_mat_arb(SD, e, dE): """Computes a matrix of the SD energies Args: SD (list of lists of ints ): a list of SD determinants, such that: SD[iSD] is a list of integers defining which orbitals are occupied in SD with index ```iSD``` and how SeeAlso: ```inp``` in the ```sd2indx(inp,nbasis)``` function e ( MATRIX(nbasis, nbasis) ): energies of all 1-electron levels (alpha- and beta- components) dE ( list of doubles ): energy corrections added to each SD Returns: CMATRIX(N,N): the matrix of energies in the SD basis. Here, N = len(SD) - the number of SDs. """ n = len(SD) E = CMATRIX(n,n) E0 = 0.0 #energy_arb(SD[0], e) + dE[0]*(1.0+0.0j) for i in range(0,n): E.set(i,i, energy_arb(SD[i], e) + dE[i]*(1.0+0.0j) - E0 ) return E
[docs]def orbs2spinorbs(s): """ This function converts the matrices in the orbital basis (e.g. old PYXAID style) to the spin-orbital basis. Essentially, it makes a block matrix of a double dimension: ( s 0 ) s --> ( 0 s ) This is meant to be used for backward compatibility with PYXIAD-generated data Args: s ( CMATRIX(N,N) ): whatever matrix in the orbital basis, N - is the number of doubly-occupied orbitals Returns: CMATRIX(2N, 2N): whatever matrix in the spin-orbital basis. It has a doubled dimensionality """ sz = s.num_of_cols zero = CMATRIX(sz, sz) act_sp1 = list(range(0, sz)) act_sp2 = list(range(sz, 2*sz)) S = CMATRIX(2*sz, 2*sz) push_submatrix(S, s, act_sp1, act_sp1) push_submatrix(S, zero, act_sp1, act_sp2) push_submatrix(S, zero, act_sp2, act_sp1) push_submatrix(S, s, act_sp2, act_sp2) return S
[docs]def ovlp_arb(SD1, SD2, S, use_minimal=False): """Compute the overlap of two generic SDs: <SD1|SD2> Args: SD1 ( lists of ints ): first SD, such that: SeeAlso: ```inp``` in the ```sd2indx(inp,nbasis)``` function SD2 ( lists of ints ): second SD, such that: SeeAlso: ```inp``` in the ```sd2indx(inp,nbasis)``` function S ( CMATRIX(N,N) ): is the matrix in the space of 1-el spin-orbitals (either spin-diabatic or spin-adiabatic or both) , N - is the number of 1-el orbitals. use_minimal ( Boolean ): If True, use the minimal subset of Kohn-Sham orbitals needed to describe the overlap of the SDs Returns: complex: the overlap of the two determinants <SD1|SD2> """ nbasis = S.num_of_rows # Converting the SDs provided by the user into the internal format to be read by Libra sd1_tmp = sd2indx(SD1,nbasis) sd2_tmp = sd2indx(SD2,nbasis) # Now, we will either keep the full SDs (this may cause artifacts when computing overlaps) # Or will use a smaller subset of these. sd1, sd2 = [], [] if use_minimal == True: #print("sd1_tmp = ", sd1_tmp) #print("sd2_tmp = ", sd2_tmp) for i in range( len(sd1_tmp) ): if sd1_tmp[i] not in sd2_tmp: sd1.append( sd1_tmp[i] ) for i in range( len(sd1_tmp) ): if sd2_tmp[i] not in sd1_tmp: sd2.append( sd2_tmp[i] ) #print("sd1 = ", sd1) #print("sd2 = ", sd2) # if both sd1 and sd2 are empty set them with sd1_tmp and sd2_tmp if len(sd1) == 0 and len(sd2) == 0: sd1 = sd1_tmp sd2 = sd2_tmp else: sd1 = sd1_tmp sd2 = sd2_tmp s = CMATRIX(len(sd1),len(sd2)) # Forming the overlap of the SDs for i in range(0,len(sd1)): for j in range(0,len(sd2)): # The overlap is non-zero only if the orbitals # are occupied with the same-spin electrons. if (SD1[i] * SD2[j]) > 0: s.set(i,j,S.get(sd1[i],sd2[j])) else: s.set(i,j,0.0,0.0) # Checking if the matrix is square #print("\nChecking if s is a square matrix ...") if s.num_of_rows != s.num_of_cols: print("\nWARNING: the matarix of Kohn-Sham orbitial overlaps is not a square matrix") print("\nExiting now ..") print("sd1 = ", sd1) print("sd2 = ", sd2) print("s.num_of_rows = ", s.num_of_rows) print("s.num_of_cols = ", s.num_of_cols) sys.exit(0) return det(s)
[docs]def ovlp_mat_arb(SD1, SD2, S, use_minimal=True): """Compute a matrix of overlaps in the SD basis Args: SD1 ( list of lists of N ints ): a list of N SD determinants, such that: SD1[iSD] is a list of integers defining which orbitals are occupied in SD with index ```iSD``` and how SeeAlso: ```inp``` in the ```sd2indx(inp,nbasis)``` function SD2 ( list of lists of M ints ): a list of M SD determinants, such that: SD2[iSD] is a list of integers defining which orbitals are occupied in SD with index ```iSD``` and how SeeAlso: ```inp``` in the ```sd2indx(inp,nbasis)``` function S ( CMATRIX(K,K) ): is the matrix in the space of 1-el spin-orbitals (either spin-diabatic or spin-adiabatic or both) , K - is the number of 1-el orbitals. use_minimal ( Boolean ): If True, use the minimal subset of Kohn-Sham orbitals needed to describe the overlap of the SDs. This is passed to the funciton that computes the overlaps Returns: CMATRIX(N, M): overlap matrix composed of elements <SD1(i)|SD2(j)> """ N, M = len(SD1), len(SD2) res = CMATRIX(N,M) for n in range(0,N): for m in range(0,M): res.set(n,m, ovlp_arb(SD1[n], SD2[m], S, use_minimal)) return res