@@ -595,16 +595,141 @@ def orthonormalize_basis_orbitals(self):
595595 # backward compatibility
596596 return self .use_native_orbitals ()
597597
598- def use_native_orbitals (self , inplace = False ):
598+ def use_native_orbitals (self , inplace = False , core : list = None , * args , ** kwargs ):
599599 """
600+ Parameters
601+ ----------
602+ inplace:
603+ update current molecule or return a new instance
604+ core:
605+ list of core orbital indices (optional) — orbitals will be frozen and treated as doubly occupied. The orbitals correspond to
606+ the currently used orbitals of the molecule (default is usually canonical HF), see mol.print_basis_info() if unsure. Providing core
607+ orbitals is optional; the default is inherited from the active space set in self.integral_manager. If core is provided, the
608+ corresponding active native orbitals will be chosen based on their overlap with the core orbitals.
609+ active(in kwargs):
610+ list the active orbital indices (optional, in kwargs) - on the Native Orbital schema. Default: All orbitals, if core (see above) is provided,
611+ then the default is to automatically select the active orbitals based on their overlap with the provided core orbitals (selectint the N-|core|
612+ orbitals that have smallest overlap with coree).
613+ As an example, Assume the input geometry was H, He, H. active=[0,1,2] is selecting the (orthonormalized) atomic 1s (left H), 1s (He), 1s (right H).
614+ If core=[0] and active is not set, then active=[0,2] will be selected automatically (as the 1s He atomic orbital will have the largest overlap
615+ with the lowest energy HF orbital).
600616 Returns
601617 -------
602618 New molecule in the native (orthonormalized) basis given
603619 e.g. for standard basis sets the orbitals are orthonormalized Gaussian Basis Functions
604620 """
605- if not self .integral_manager .active_space_is_trivial ():
606- warnings .warn ("orthonormalize_basis_orbitals: active space is set and might lead to inconsistent behaviour" , TequilaWarning )
607-
621+ c = copy .deepcopy (self .integral_manager .orbital_coefficients )
622+ s = self .integral_manager .overlap_integrals
623+ d = self .integral_manager .get_orthonormalized_orbital_coefficients ()
624+ def inner (a , b ,s ):
625+ return numpy .sum (numpy .multiply (numpy .outer (a ,b ),s ))
626+ def orthogonalize (c ,d ,s ):
627+ '''
628+ :return: orthogonalized orbitals with core HF orbitals and active Orthongonalized Native orbitals.
629+ '''
630+ ### Computing Core-Active overlap Matrix
631+ # sbar_{ki} = \langle \phi_k | \varphi_i \rangle = \sum_{m,n} c_{nk}d_{mi}\langle \chi_n | \chi_m \rangle
632+ # c_{nk} = HF coeffs, d_{mi} = nat orb coef s_{mn} = Atomic Overlap Matrix
633+ # k \in active orbs, i \in core orbs, m,n \in basis coeffs
634+ # sbar = np.einsum('nk,mi,nm->ki', c, d, s) #only works if active == to_active
635+ c = c .T
636+ d = d .T
637+ sbar = numpy .zeros (shape = s .shape )
638+ for k in active :
639+ for i in core :
640+ sbar [i ][to_active [k ]] = inner (c [i ],d [k ],s )
641+ ### Projecting out Core orbitals from the Native ones
642+ # dbar_{ji} = d_{ji} - \sum_k sbar_{ki}c_{jk}
643+ # k \in active, i \in core, j in basis coeffs
644+ dbar = numpy .zeros (shape = s .shape )
645+
646+ for j in active :
647+ dbar [to_active [j ]] = d [j ]
648+ for i in core :
649+ temp = sbar [i ][to_active [j ]] * c [i ]
650+ dbar [to_active [j ]] -= temp
651+ ### Projected-out Nat Orbs Normalization
652+ for i in to_active .values ():
653+ norm = numpy .sqrt (numpy .sum (numpy .multiply (numpy .outer (dbar [i ], dbar [i ]), s .T )))
654+ if not numpy .isclose (norm , 0 ):
655+ dbar [i ] = dbar [i ] / norm
656+ ### Reintroducing the New Coeffs on the HF coeff matrix
657+ for j in to_active .values ():
658+ c [j ] = dbar [j ]
659+ ### Compute new orbital overlap matrix:
660+ sprima = numpy .eye (len (c ))
661+ for idx , i in enumerate (to_active .values ()):
662+ for j in [* to_active .values ()][idx :]:
663+ sprima [i ][j ] = inner (c [i ],c [j ],s )
664+ sprima [j ][i ] = sprima [i ][j ]
665+ ### Symmetric orthonormalization
666+ lam_s , l_s = numpy .linalg .eigh (sprima )
667+ lam_s = lam_s * numpy .eye (len (lam_s ))
668+ lam_sqrt_inv = numpy .sqrt (numpy .linalg .inv (lam_s ))
669+ symm_orthog = numpy .dot (l_s , numpy .dot (lam_sqrt_inv , l_s .T ))
670+ return symm_orthog .dot (c ).T
671+ def get_active (core ):
672+ ov = numpy .zeros (shape = (len (self .integral_manager .orbitals )))
673+ for i in core :
674+ for j in range (len (d )):
675+ ov [j ] += numpy .abs (inner (c .T [i ], d .T [j ],s ))
676+ act = []
677+ for i in range (len (self .integral_manager .orbitals ) - len (core )):
678+ idx = numpy .argmin (ov )
679+ act .append (idx )
680+ ov [idx ] = 1 * len (core )
681+ act .sort ()
682+ return act
683+
684+ def get_core (active ):
685+ ov = numpy .zeros (shape = (len (self .integral_manager .orbitals )))
686+ for i in active :
687+ for j in range (len (d )):
688+ ov [j ] += numpy .abs (inner (d .T [i ], c .T [j ], s ))
689+ co = []
690+ for i in range (len (self .integral_manager .orbitals ) - len (active )):
691+ idx = numpy .argmin (ov )
692+ co .append (idx )
693+ ov [idx ] = 1 * len (active )
694+ co .sort ()
695+ return co
696+
697+ active = None
698+ if not self .integral_manager .active_space_is_trivial () and core is None :
699+ core = [i .idx_total for i in self .integral_manager .orbitals if i .idx is None ]
700+ if 'active' in kwargs :
701+ active = kwargs ['active' ]
702+ kwargs .pop ('active' )
703+ if core is None :
704+ core = get_core (active )
705+ else :
706+ if active is None :
707+ if core is None :
708+ core = []
709+ active = [i for i in range (len (self .integral_manager .orbitals ))]
710+ else :
711+ if isinstance (core ,int ):
712+ core = [core ]
713+ active = get_active (core )
714+ assert len (active ) + len (core ) == len (self .integral_manager .orbitals )
715+ to_active = [i for i in range (len (self .integral_manager .orbitals )) if i not in core ]
716+ to_active = {active [i ]: to_active [i ] for i in range (len (active ))}
717+ if len (core ):
718+ coeff = orthogonalize (c ,d ,s )
719+ if inplace :
720+ self .integral_manager = self .initialize_integral_manager (one_body_integrals = self .integral_manager .one_body_integrals ,
721+ two_body_integrals = self .integral_manager .two_body_integrals ,constant_term = self .integral_manager .constant_term ,
722+ active_orbitals = [* to_active .values ()],reference_orbitals = [i .idx_total for i in self .integral_manager .reference_orbitals ]
723+ , frozen_orbitals = core , orbital_coefficients = coeff , overlap_integrals = s )
724+ return self
725+ else :
726+ integral_manager = self .initialize_integral_manager (one_body_integrals = self .integral_manager .one_body_integrals ,
727+ two_body_integrals = self .integral_manager .two_body_integrals ,constant_term = self .integral_manager .constant_term
728+ , active_orbitals = [* to_active .values ()],reference_orbitals = [i .idx_total for i in self .integral_manager .reference_orbitals ]
729+ , frozen_orbitals = core , orbital_coefficients = coeff , overlap_integrals = s )
730+ parameters = copy .deepcopy (self .parameters )
731+ result = QuantumChemistryBase (parameters = parameters , integral_manager = integral_manager ,transformation = self .transformation ,active_orbitals = [* to_active .values ()])
732+ return result
608733 # can not be an instance of a specific backend (otherwise we get inconsistencies with classical methods in the backend)
609734 if inplace :
610735 self .integral_manager .transform_to_native_orbitals ()
0 commit comments