# Copyright 2018-2024 Xanadu Quantum Technologies Inc.# Licensed under the Apache License, Version 2.0 (the "License");# you may not use this file except in compliance with the License.# You may obtain a copy of the License at# http://www.apache.org/licenses/LICENSE-2.0# Unless required by applicable law or agreed to in writing, software# distributed under the License is distributed on an "AS IS" BASIS,# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.# See the License for the specific language governing permissions and# limitations under the License."""This module contains functions to perform VSCF calculation."""importnumpyasnp# pylint: disable=too-many-arguments, too-many-positional-argumentsdef_build_fock(mode,active_ham_terms,active_terms,modals,h_mat,mode_rot):r"""Builds the fock matrix for each vibrational mode. Args: mode (int): index of the active mode for which the Fock matrix is being calculated active_ham_terms (dict): dictionary containing the active Hamiltonian terms active_terms (list): list containing the active Hamiltonian terms for the given mode modals (list[int]): list containing the maximum number of modals to consider for each mode h_mat (array[array[float]]): the Hamiltonian matrix mode_rot (TensorLike[float]): rotation matrix for the given vibrational mode Returns: array[array[float]]: fock matrix for the given mode """fock_matrix=np.zeros((modals[mode],modals[mode]))forterminactive_terms:ham_term=active_ham_terms[term]h_val,modal_indices,excitation_indices=ham_term[0],list(ham_term[1]),list(ham_term[2])modal_idx=modal_indices.index(mode)i=excitation_indices[modal_idx]j=excitation_indices[modal_idx+len(modal_indices)]ifi<modals[mode]andj<modals[mode]:h_val=h_val*np.prod([h_mat[term,modal]formodalinmodal_indicesifmodal!=mode])fock_matrix[i,j]+=h_valfock_matrix=mode_rot.T@fock_matrix@mode_rotreturnfock_matrixdef_update_h(h_mat,mode,active_ham_terms,mode_rot,active_terms,modals):r"""Updates value of Hamiltonian matrix for a mode after other modes have been rotated. Args: h_mat (array[array[float]]): the Hamiltonian matrix mode (int): index of the active mode for which the Hamiltonian matrix is being updated active_ham_terms (dict): dictionary containing the active Hamiltonian terms mode_rot (array[array[float]]): rotation matrix for a given mode active_terms (list): list containing active Hamiltonian terms for the given mode modals (list[int]): list containing the maximum number of modals to consider for each mode Returns: array[array[float]]: the updated Hamiltonian matrix """u_coeffs=mode_rot[:,0]fortinactive_terms:ham_term=active_ham_terms[t]modal_indices,excitation_indices=ham_term[1],ham_term[2]modal_idx=modal_indices.index(mode)i=excitation_indices[modal_idx]j=excitation_indices[modal_idx+len(modal_indices)]ifi<modals[mode]andj<modals[mode]:h_mat[t,mode]=u_coeffs[i]*u_coeffs[j]returnh_matdef_find_active_terms(h_integrals,modals,cutoff):r"""Identifies the active terms in the Hamiltonian, following the equations 20-22 in `J. Chem. Theory Comput. 2010, 6, 235–248 <https://pubs.acs.org/doi/10.1021/ct9004454>`_. The equation suggests that if mode m is not contained in a Hamiltonian term, it evaluates to zero. Args: h_integrals (list[TensorLike[float]]): list of n-mode expansion of Hamiltonian integrals modals (list[int]): list containing the maximum number of modals to consider for each mode cutoff (float): threshold value for including matrix elements into operator Returns: tuple: A tuple containing the following: - dict: dictionary containing the active Hamiltonian terms - dict: dictionary containing the active Hamiltonian terms for all vibrational modes - int: total number of active terms in the Hamiltonian """active_ham_terms={}active_num=0nmodes=np.shape(h_integrals[0])[0]forn,ham_ninenumerate(h_integrals):foridx,h_valinnp.ndenumerate(ham_n):ifnp.abs(h_val)<cutoff:continuemodal_indices=idx[:n+1]ifall(modal_indices[i]>modal_indices[i+1]foriinrange(n)):excitation_indices=idx[n+1:]exc_in_modals=Trueform_idx,minenumerate(modal_indices):i=excitation_indices[m_idx]j=excitation_indices[m_idx+len(modal_indices)]ifi>=modals[m]orj>=modals[m]:exc_in_modals=Falsebreakifexc_in_modals:active_ham_terms[active_num]=[h_val,modal_indices,excitation_indices]active_num+=1active_mode_terms={m:[tfortinrange(active_num)ifminactive_ham_terms[t][1]]forminrange(nmodes)}returnactive_ham_terms,active_mode_terms,active_numdef_fock_energy(h_mat,active_ham_terms,active_mode_terms,modals,mode_rots):r"""Calculates vibrational energy. Args: h_mat (array[array[float]]): the Hamiltonian matrix active_ham_terms (dict): dictionary containing the active Hamiltonian terms active_mode_terms (dict): list containing the active Hamiltonian terms for all vibrational modes modals (list[int]): list containing the maximum number of modals to consider for each mode mode_rots (List[TensorLike[float]]): list of rotation matrices for each vibrational mode Returns: float: vibrational energy """nmodes=h_mat.shape[1]e0s=np.zeros(nmodes)formodeinrange(nmodes):fock_mat=_build_fock(mode,active_ham_terms,active_mode_terms[mode],modals,h_mat,mode_rots[mode])e0s[mode]=fock_mat[0,0]returnnp.sum(e0s)def_vscf(h_integrals,modals,cutoff,tol=1e-8,max_iters=10000):r"""Performs the VSCF calculation. Args: h_integrals (list(TensorLike[float])): list containing Hamiltonian integral matrices modals (list(int)): list containing the maximum number of modals to consider for each mode cutoff (float): threshold value for including matrix elements into operator. tol (float): convergence tolerance for vscf calculation max_iters (int): maximum number of iterations for vscf to converge Returns: tuple: A tuple of the following: - float: vscf energy - list[TensorLike[float]]: list of rotation matrices for all vibrational modes """nmodes=np.shape(h_integrals[0])[0]active_ham_terms,active_mode_terms,active_num=_find_active_terms(h_integrals,modals,cutoff)mode_rots=[np.eye(modals[mode])formodeinrange(nmodes)]h_mat=np.zeros((active_num,nmodes))formodeinrange(nmodes):h_mat=_update_h(h_mat,mode,active_ham_terms,mode_rots[mode],active_mode_terms[mode],modals)e0=_fock_energy(h_mat,active_ham_terms,active_mode_terms,modals,mode_rots)enew=e0+2*tolforcurr_iterinrange(max_iters):ifcurr_iter!=0:e0=enewformodeinrange(nmodes):fock=_build_fock(mode,active_ham_terms,active_mode_terms[mode],modals,h_mat,mode_rots[mode])_,eigvec=np.linalg.eigh(fock)mode_rots[mode]=mode_rots[mode]@eigvech_mat=_update_h(h_mat,mode,active_ham_terms,mode_rots[mode],active_mode_terms[mode],modals)fock=np.transpose(eigvec)@fock@eigvecenew=_fock_energy(h_mat,active_ham_terms,active_mode_terms,modals,mode_rots)ifnp.abs(enew-e0)<=tol:breakreturnenew,mode_rotsdef_rotate_one_body(h1,nmodes,mode_rots,modals):r"""Rotates one body integrals. Args: h1 (TensorLike[float]): one-body integrals nmodes (int): number of vibrational modes mode_rots (list[TensorLike[float]]): list of rotation matrices for all vibrational modes modals (list[int]): list containing the maximum number of modals to consider for each mode Returns: TensorLike[float]: rotated one-body integrals """imax=np.max(modals)h1_rot=np.zeros((nmodes,imax,imax))forminrange(nmodes):h1_rot[m,:modals[m],:modals[m]]=np.einsum("ij,ia,jb->ab",h1[m,:,:],mode_rots[m][:,:modals[m]],mode_rots[m][:,:modals[m]])returnh1_rotdef_rotate_two_body(h2,nmodes,mode_rots,modals):r"""Rotates two body integrals. Args: h2 (TensorLike[float]): two-body integrals nmodes (int): number of vibrational modes mode_rots (list[TensorLike[float]]): list of rotation matrices for all vibrational modes modals (list[int]): list containing the maximum number of modals to consider for each mode Returns: TensorLike[float]: rotated two-body integrals """imax=np.max(modals)U_mats=[mode_rots[m]forminrange(nmodes)]h2_rot=np.zeros((nmodes,nmodes,imax,imax,imax,imax))form1inrange(nmodes):form2inrange(nmodes):h2_rot[m1,m2,:modals[m1],:modals[m2],:modals[m1],:modals[m2]]=np.einsum("ijkl,ia,jb,kc,ld->abcd",h2[m1,m2,:,:,:,:],U_mats[m1][:,:modals[m1]],U_mats[m2][:,:modals[m2]],U_mats[m1][:,:modals[m1]],U_mats[m2][:,:modals[m2]],)returnh2_rotdef_rotate_three_body(h3,nmodes,mode_rots,modals):r"""Rotates three body integrals. Args: h3 (TensorLike[float]): three-body integrals nmodes (int): number of vibrational modes mode_rots (list[TensorLike[float]]): list of rotation matrices for all vibrational modes modals (list[int]): list containing the maximum number of modals to consider for each mode Returns: TensorLike[float]: rotated three-body integrals """imax=np.max(modals)h3_rot=np.zeros((nmodes,nmodes,nmodes,imax,imax,imax,imax,imax,imax))form1inrange(nmodes):form2inrange(nmodes):form3inrange(nmodes):h3_rot[m1,m2,m3,:modals[m1],:modals[m2],:modals[m3],:modals[m1],:modals[m2],:modals[m3],]=np.einsum("ijklmn,ia,jb,kc,ld,me,nf->abcdef",h3[m1,m2,m3,:,:,:,:,:,:],mode_rots[m1][:,:modals[m1]],mode_rots[m2][:,:modals[m2]],mode_rots[m3][:,:modals[m3]],mode_rots[m1][:,:modals[m1]],mode_rots[m2][:,:modals[m2]],mode_rots[m3][:,:modals[m3]],)returnh3_rotdef_rotate_dipole(d_integrals,mode_rots,modals):r"""Generates VSCF rotated dipole. Args: d_integrals (list[TensorLike[float]]): list of n-mode expansion of dipole integrals mode_rots (list[TensorLike[float]]): list of rotation matrices for all vibrational modes modals (list[int]): list containing the maximum number of modals to consider for each mode Returns: tuple(TensorLike[float]): a tuple of rotated dipole integrals """n=len(d_integrals)nmodes=np.shape(d_integrals[0])[0]imax=np.max(modals)d1_rot=np.zeros((3,nmodes,imax,imax))foralphainrange(3):d1_rot[alpha,::]=_rotate_one_body(d_integrals[0][alpha,::],nmodes,mode_rots,modals=modals)dip_data=[d1_rot]ifn>1:d2_rot=np.zeros((3,nmodes,nmodes,imax,imax,imax,imax))foralphainrange(3):d2_rot[alpha,::]=_rotate_two_body(d_integrals[1][alpha,::],nmodes,mode_rots,modals=modals)dip_data=[d1_rot,d2_rot]ifn>2:d3_rot=np.zeros((3,nmodes,nmodes,nmodes,imax,imax,imax,imax,imax,imax))foralphainrange(3):d3_rot[alpha,::]=_rotate_three_body(d_integrals[2][alpha,::],nmodes,mode_rots,modals=modals)dip_data=[d1_rot,d2_rot,d3_rot]returndip_datadef_rotate_hamiltonian(h_integrals,mode_rots,modals):r"""Generates VSCF rotated Hamiltonian. Args: h_integrals (list[TensorLike[float]]): list of n-mode expansion of Hamiltonian integrals mode_rots (list[TensorLike[float]]): list of rotation matrices for all vibrational modes modals (list[int]): list containing the maximum number of modals to consider for each mode Returns: tuple(TensorLike[float]): tuple of rotated Hamiltonian integrals """n=len(h_integrals)nmodes=np.shape(h_integrals[0])[0]h1_rot=_rotate_one_body(h_integrals[0],nmodes,mode_rots,modals)h_data=[h1_rot]ifn>1:h2_rot=_rotate_two_body(h_integrals[1],nmodes,mode_rots,modals)h_data=[h1_rot,h2_rot]ifn>2:h3_rot=_rotate_three_body(h_integrals[2],nmodes,mode_rots,modals)h_data=[h1_rot,h2_rot,h3_rot]returnh_data
[docs]defvscf_integrals(h_integrals,d_integrals=None,modals=None,cutoff=None,cutoff_ratio=1e-6):r"""Generates vibrational self-consistent field rotated integrals. This function converts the Christiansen vibrational Hamiltonian integrals obtained in the harmonic oscillator basis to integrals in the vibrational self-consistent field (VSCF) basis. The implementation is based on the method described in `J. Chem. Theory Comput. 2010, 6, 235–248 <https://pubs.acs.org/doi/10.1021/ct9004454>`_. Args: h_integrals (list[TensorLike[float]]): List of Hamiltonian integrals for up to 3 coupled vibrational modes. Look at the Usage Details for more information. d_integrals (list[TensorLike[float]]): List of dipole integrals for up to 3 coupled vibrational modes. Look at the Usage Details for more information. modals (list[int]): list containing the maximum number of modals to consider for each vibrational mode. Default value is the maximum number of modals. cutoff (float): threshold value for including matrix elements into operator cutoff_ratio (float): ratio for discarding elements with respect to biggest element in the integrals. Default value is ``1e-6``. Returns: tuple: a tuple containing: - list[TensorLike[float]]: List of Hamiltonian integrals in VSCF basis for up to 3 coupled vibrational modes. - list[TensorLike[float]]: List of dipole integrals in VSCF basis for up to 3 coupled vibrational modes. ``None`` is returned if ``d_integrals`` is ``None``. **Example** >>> h1 = np.array([[[0.00968289, 0.00233724, 0.0007408, 0.00199125], ... [0.00233724, 0.02958449, 0.00675431, 0.0021936], ... [0.0007408, 0.00675431, 0.0506012, 0.01280986], ... [0.00199125, 0.0021936, 0.01280986, 0.07282307]]]) >>> qml.qchem.vscf_integrals(h_integrals=[h1], modals=[4,4,4]) ([array([[[ 9.36124041e-03, -4.20128342e-19, 3.25260652e-19, 1.08420217e-18], [-9.21571847e-19, 2.77803512e-02, -3.46944695e-18, 5.63785130e-18], [-3.25260652e-19, -8.67361738e-19, 4.63297357e-02, -1.04083409e-17], [ 1.30104261e-18, 5.20417043e-18, -1.38777878e-17, 7.92203227e-02]]])], None) .. details:: :title: Usage Details The ``h_integral`` tensor must have one of these dimensions: - 1-mode coupled integrals: `(n, m)` - 2-mode coupled integrals: `(n, n, m, m, m, m)` - 3-mode coupled integrals: `(n, n, n, m, m, m, m, m, m)` where ``n`` is the number of vibrational modes in the molecule and ``m`` represents the number of modals. The ``d_integral`` tensor must have one of these dimensions: - 1-mode coupled integrals: `(3, n, m)` - 2-mode coupled integrals: `(3, n, n, m, m, m, m)` - 3-mode coupled integrals: `(3, n, n, n, m, m, m, m, m, m)` where ``n`` is the number of vibrational modes in the molecule, ``m`` represents the number of modals and the first axis represents the ``x, y, z`` component of the dipole. Default is ``None``. """iflen(h_integrals)>3:raiseValueError(f"Building n-mode Hamiltonian is not implemented for n equal to {len(h_integrals)}.")ifd_integralsisnotNone:iflen(d_integrals)>3:raiseValueError(f"Building n-mode dipole is not implemented for n equal to {len(d_integrals)}.")nmodes=np.shape(h_integrals[0])[0]imax=np.shape(h_integrals[0])[1]max_modals=nmodes*[imax]ifmodalsisNone:modals=max_modalselse:ifnp.max(modals)>imax:raiseValueError("Number of maximum modals cannot be greater than the modals for unrotated integrals.")imax=np.max(modals)ifcutoffisNone:max_val=np.max([np.max(np.abs(H))forHinh_integrals])cutoff=max_val*cutoff_ratio_,mode_rots=_vscf(h_integrals,modals=max_modals,cutoff=cutoff)h_data=_rotate_hamiltonian(h_integrals,mode_rots,modals)ifd_integralsisnotNone:dip_data=_rotate_dipole(d_integrals,mode_rots,modals)returnh_data,dip_datareturnh_data,None