diff --git a/src/mqc_algebra.F03 b/src/mqc_algebra.F03 index dbfde698..c9b9863e 100644 --- a/src/mqc_algebra.F03 +++ b/src/mqc_algebra.F03 @@ -89,6 +89,8 @@ Module MQC_Algebra Procedure, Public::cval => MQC_Scalar_Get_Intrinsic_Complex !> \brief Take the absolute value of the MQC Scalar Procedure, Public::abs => MQC_Scalar_Get_ABS_Value +!> \brief Take the phase of the MQC Scalar + Procedure, Public::phase => MQC_Scalar_Get_Phase !> \brief Take the square root of the MQC Scalar Procedure, Public::sqrt => MQC_Scalar_Sqrt !> \brief Take the exponential of the MQC Scalar @@ -97,8 +99,6 @@ Module MQC_Algebra Procedure, Public::log => MQC_Scalar_Log !> \brief Return a random value to the MQC Scalar Procedure, Public::random => MQC_Scalar_Get_Random_Value -!> \brief Take the phase of the MQC Scalar - Procedure, Public::phase => MQC_Scalar_Get_Phase End Type MQC_Scalar ! ! Vectors... @@ -132,6 +132,12 @@ Module MQC_Algebra Generic, Public::at => at_sca, at_int !> \brief Returns the subvector between specified element of the MQC Vector Procedure, Public::vat => MQC_Vector_Vector_At +!> \brief Returns the VecI element of the MQC Vector + Procedure, Public::getArrayI => MQC_Vector_Integer_Get_Array +!> \brief Returns the VecR element of the MQC Vector + Procedure, Public::getArrayR => MQC_Vector_Real_Get_Array +!> \brief Returns the VecC element of the MQC Vector + Procedure, Public::getArrayC => MQC_Vector_Complex_Get_Array Procedure, Private::put_sca => MQC_Vector_Scalar_Put_Scalar Procedure, Private::put_int => MQC_Vector_Scalar_Put_Int !> \brief Updates the specified element of the MQC_Vector with the specified value @@ -653,12 +659,12 @@ Module MQC_Algebra module procedure MQC_R4Tensor_Put_Int end interface ! -!> \brief Puts a value in the R4tensor at the specied index +!> \brief Returns the identity matrix of the specified size interface identity module procedure MQC_Matrix_Identity_Func end interface ! -!> \brief Puts a value in the R4tensor at the specied index +!> \brief Initializes a matrix of the specified size interface init module procedure MQC_Matrix_Initialize_Func end interface @@ -2816,8 +2822,7 @@ Function MQC_Scalar_Get_Phase(Scalar,thresh_inp) Result(Output) else thresh=1.0E-8 end if - - if (abs(Scalar).lt.thresh) then + if (abs(Scalar).lt.thresh) then output = 0.0 else output = acos(Scalar%rval()/abs(Scalar)) @@ -8193,6 +8198,140 @@ Function MQC_Vector_Vector_At(Vec,I,J) Result(Vector) Vector%Column = Vec%Column End Function + +! +! +! PROCEDURE MQC_Vector_Integer_Get_Array +! +!> \brief MQC_Vector_Integer_Get_Array is a function that returns the VecI variable +!> of a MQC vector +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_Vector_Integer_Get_Array is a function that returns the VecI variable of a MQC +!> vector. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] Vec +!> \verbatim +!> Vec is Class(MQC_Vector) +!> The MQC_Vector to extract the VecI variable. +!> \endverbatim +!> +! +! Authors: +! ======== +!> \author M. M. F. Moraes +!> \date 2025 +! + Function MQC_Vector_Integer_Get_Array(Vec) Result(ArrayOut) +! + Implicit None + Class(MQC_Vector)::Vec + Integer(kind=int64),Dimension(:),Allocatable::ArrayOut + + if (Vec%Data_Type.eq.'Integer') then + ArrayOut = Vec%VecI + else + Call MQC_Error_A('Vector is not an Integer Data_Type in MQC_Vector_Integer_Get_Array', 6, & + 'Vec%Data_Type', Vec%Data_Type ) + EndIf + End Function MQC_Vector_Integer_Get_Array +! +! +! PROCEDURE MQC_Vector_Real_Get_Array +! +!> \brief MQC_Vector_Real_Get_Array is a function that returns the VecR variable +!> of a MQC vector +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_Vector_Real_Get_Array is a function that returns the VecR variable of a MQC +!> vector. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] Vec +!> \verbatim +!> Vec is Class(MQC_Vector) +!> The MQC_Vector to extract the VecR variable. +!> \endverbatim +!> +! +! Authors: +! ======== +!> \author M. M. F. Moraes +!> \date 2025 +! + Function MQC_Vector_Real_Get_Array(Vec) Result(ArrayOut) +! + Implicit None + Class(MQC_Vector)::Vec + Real(kind=real64),Dimension(:),Allocatable::ArrayOut + + if (Vec%Data_Type.eq.'Real') then + ArrayOut = Vec%VecR + else + Call MQC_Error_A('Vector is not a Real Data_Type in MQC_Vector_Real_Get_Array', 6, & + 'Vec%Data_Type', Vec%Data_Type ) + EndIf + End Function MQC_Vector_Real_Get_Array +! +! +! PROCEDURE MQC_Vector_Complex_Get_Array +! +!> \brief MQC_Vector_Complex_Get_Array is a function that returns the VecC variable +!> of a MQC vector +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_Vector_Complex_Get_Array is a function that returns the VecC variable of a MQC +!> vector. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] Vec +!> \verbatim +!> Vec is Class(MQC_Vector) +!> The MQC_Vector to extract the VecC variable. +!> \endverbatim +!> +! +! Authors: +! ======== +!> \author M. M. F. Moraes +!> \date 2025 +! + Function MQC_Vector_Complex_Get_Array(Vec) Result(ArrayOut) +! + Implicit None + Class(MQC_Vector)::Vec + Complex(kind=real64),Dimension(:),Allocatable::ArrayOut + + if (Vec%Data_Type.eq.'Complex') then + ArrayOut = Vec%VecC + else + Call MQC_Error_A('Vector is not a Complex Data_Type in MQC_Vector_Complex_Get_Array', 6, & + 'Vec%Data_Type', Vec%Data_Type ) + EndIf + End Function MQC_Vector_Complex_Get_Array + ! ! ! PROCEDURE MQC_Set_Vector2IntegerArray @@ -15721,6 +15860,7 @@ Subroutine MQC_Matrix_SymmMatrix_Put_Complex(mat,symmMatrixIn,label) end subroutine MQC_Matrix_SymmMatrix_Put_Complex ! ! +! ! PROCEDURE MQC_Matrix_Matrix_Put ! !> \brief MQC_Matrix_Matrix_Put is a subroutine that writes a submatrix to the @@ -23166,6 +23306,11 @@ Function MQC_MatrixMatrixDotProduct(MA,MB) Result(MC) Type(MQC_Matrix)::MC Type(MQC_Matrix)::MAreal,MBreal,Tmp1,Tmp2 integer(kind=int64)::i,j + + If (MA%NCol.eq.0.or.MB%NRow.eq.0) then + Call MQC_Allocate_Matrix(0,0,MC,'Real','StorFull') + return + end if If (MA%NCol /= MB%NRow) call MQC_Error_I('The two matrices are not conformable for multiplication', 6, & 'MA%NCol', MA%NCol, & @@ -24132,7 +24277,7 @@ Function MQC_VectormatrixDotProduct(VA,MB) Result(VC) If (VA%Column .or. VA%Length /= MB%NRow) then Call MQC_Error_i('Matrix and vector are not conformable for multiplication', 6, & 'VA%Length', VA%Length, & - 'MB%NRow', MB%NRow ) + 'MB%NRow', MB%NRow) EndIf If (MQC_Matrix_HaveReal(MB) .and. MQC_Vector_HaveReal(VA)) then @@ -26838,71 +26983,66 @@ Subroutine MQC_Matrix_SVD(A,EVals,EUVecs,EVVecs) storFlag = A%Storage 1050 Format( A, I10 ) ! - if(storFlag.eq.'StorDiag') return + if(storFlag.eq.'StorDiag') call mqc_matrix_diag2Full(A) if(storFlag.eq.'StorSymm'.or.storFlag.eq.'StorAsym'.or.& storFlag.eq.'StorHerm'.or.storFlag.eq.'StorAhrm') & call mqc_matrix_symm2Full(A) - if(storFlag.eq.'StorFull') then - if(typeFlag.eq.'Real'.or.typeFlag.eq.'Integer') then - if(typeFlag.eq.'Integer') call mqc_matrix_copy_int2Real(A) - Allocate(A_Temp(NDim1,NDim2),A_SVals(Min(NDim1,NDim2)),A_UVecs(NDim1,NDim1), & - A_VVecs(NDim2,NDim2),Work(1)) - A_Temp = A - Call DGESVD('A','A',NDim1,NDim2,A_Temp,NDim1,A_SVals,A_UVecs,NDim1,A_VVecs, & - NDim2,Work,Work_Length,IError) - Work_Length = Int(Work(1)) - DeAllocate(Work) - Allocate(Work(Work_Length)) - Call DGESVD('A','A',NDim1,NDim2,A_Temp,NDim1,A_SVals,A_UVecs,NDim1,A_VVecs, & - NDim2,Work,Work_Length,IError) - If(IError.ne.0) Write(*,1050)' SVD FAILED: IError = ',IError - DeAllocate(Work) - Work_Length = -1 - elseIf(typeFlag.eq.'Complex') then - Allocate(AC_Temp(NDim1,NDim2),A_SVals(Min(NDim1,NDim2)),AC_UVecs(NDim1,NDim1), & - AC_VVecs(NDim2,NDim2),WorkC(1),Work(5*min(NDim1,NDim2))) - AC_Temp = A - Call ZGESVD('A','A',NDim1,NDim2,AC_Temp,NDim1,A_SVals,AC_UVecs,NDim1,AC_VVecs, & - NDim2,WorkC,Work_Length,Work,IError) - Work_Length = Int(WorkC(1)) - DeAllocate(WorkC) - Allocate(WorkC(Work_Length)) - Call ZGESVD('A','A',NDim1,NDim2,AC_Temp,NDim1,A_SVals,AC_UVecs,NDim1,AC_VVecs, & - NDim2,WorkC,Work_Length,Work,IError) - If(IError.ne.0) Write(*,1050)' SVD FAILED: IError = ',IError - DeAllocate(WorkC,Work) - Work_Length = -1 - else - call mqc_error_a('typeFlag not recognised in mqc_matrix_svd', 6, & - 'typeFlag', typeFlag ) - endIf -! - if(present(EVals)) EVals = A_SVals - if(present(EUVecs)) then - if(allocated(A_UVecs)) EUVecs = A_UVecs - if(allocated(AC_UVecs)) EUVecs = AC_UVecs - endIf - if(present(EVVecs)) then - if(allocated(A_VVecs)) EVVecs = A_VVecs - if(allocated(AC_VVecs)) EVVecs = AC_VVecs - endIf - if(allocated(A_SVals)) deallocate(A_SVals) - if(allocated(A_UVecs)) deallocate(A_UVecs) - if(allocated(AC_UVecs)) deallocate(AC_UVecs) - if(allocated(A_VVecs)) deallocate(A_VVecs) - if(allocated(AC_VVecs)) deallocate(AC_VVecs) - if(allocated(A_Temp)) deallocate(A_Temp) - if(allocated(AC_Temp)) deallocate(AC_Temp) -! + if(typeFlag.eq.'Real'.or.typeFlag.eq.'Integer') then + if(typeFlag.eq.'Integer') call mqc_matrix_copy_int2Real(A) + Allocate(A_Temp(NDim1,NDim2),A_SVals(Min(NDim1,NDim2)),A_UVecs(NDim1,NDim1), & + A_VVecs(NDim2,NDim2),Work(1)) + A_Temp = A + Call DGESVD('A','A',NDim1,NDim2,A_Temp,NDim1,A_SVals,A_UVecs,NDim1,A_VVecs, & + NDim2,Work,Work_Length,IError) + Work_Length = Int(Work(1)) + DeAllocate(Work) + Allocate(Work(Work_Length)) + Call DGESVD('A','A',NDim1,NDim2,A_Temp,NDim1,A_SVals,A_UVecs,NDim1,A_VVecs, & + NDim2,Work,Work_Length,IError) + If(IError.ne.0) Write(*,1050)' SVD FAILED: IError = ',IError + DeAllocate(Work) + Work_Length = -1 + elseIf(typeFlag.eq.'Complex') then + Allocate(AC_Temp(NDim1,NDim2),A_SVals(Min(NDim1,NDim2)),AC_UVecs(NDim1,NDim1), & + AC_VVecs(NDim2,NDim2),WorkC(1),Work(5*min(NDim1,NDim2))) + AC_Temp = A + Call ZGESVD('A','A',NDim1,NDim2,AC_Temp,NDim1,A_SVals,AC_UVecs,NDim1,AC_VVecs, & + NDim2,WorkC,Work_Length,Work,IError) + Work_Length = Int(WorkC(1)) + DeAllocate(WorkC) + Allocate(WorkC(Work_Length)) + Call ZGESVD('A','A',NDim1,NDim2,AC_Temp,NDim1,A_SVals,AC_UVecs,NDim1,AC_VVecs, & + NDim2,WorkC,Work_Length,Work,IError) + If(IError.ne.0) Write(*,1050)' SVD FAILED: IError = ',IError + DeAllocate(WorkC,Work) + Work_Length = -1 else - call mqc_error_a('storFlag not recognised in mqc_matrix_svd', 6, & - 'storFlag', storFlag ) - endIf + call mqc_error_a('typeFlag not recognised in mqc_matrix_svd', 6, & + 'typeFlag', typeFlag ) + endIf +! + if(present(EVals)) EVals = A_SVals + if(present(EUVecs)) then + if(allocated(A_UVecs)) EUVecs = A_UVecs + if(allocated(AC_UVecs)) EUVecs = AC_UVecs + endIf + if(present(EVVecs)) then + if(allocated(A_VVecs)) EVVecs = A_VVecs + if(allocated(AC_VVecs)) EVVecs = AC_VVecs + endIf + if(allocated(A_SVals)) deallocate(A_SVals) + if(allocated(A_UVecs)) deallocate(A_UVecs) + if(allocated(AC_UVecs)) deallocate(AC_UVecs) + if(allocated(A_VVecs)) deallocate(A_VVecs) + if(allocated(AC_VVecs)) deallocate(AC_VVecs) + if(allocated(A_Temp)) deallocate(A_Temp) + if(allocated(AC_Temp)) deallocate(AC_Temp) ! if(typeFlag.eq.'Integer') call mqc_matrix_copy_real2Int(A) if(storFlag.eq.'StorSymm'.or.storFlag.eq.'StorAsym'.or.& storFlag.eq.'StorHerm'.or.storFlag.eq.'StorAhrm') & call mqc_matrix_full2Symm(A) + if(storFlag.eq.'StorDiag') call mqc_matrix_full2diag(A) ! Return End Subroutine MQC_Matrix_SVD @@ -27324,6 +27464,7 @@ Function MQC_Givens_Matrix(m_size,angle,p,q) !> Matrix2 is type(mqc_matrix) !> The second matrix in the outer product. !> \endverbatim +!> ! ! Authors: ! ======== @@ -27342,9 +27483,9 @@ Function MQC_Matrix_Matrix_Outer(matrix1,matrix2) result(outer) N = size(matrix1,2) O = size(matrix2,1) P = size(matrix2,2) - + call outer%init(M,N,O,P) - + do i = 1,M do j = 1,N do k = 1,O @@ -28051,8 +28192,8 @@ subroutine mqc_matrix_power(A,P) call mqc_error_I('P type not defined in MQC_Vector_Power',6) end select ! -! call A%eigensys(eigenvals=sVals,reigenvecs=uMatrix) - call A%diag(evals=sVals,evecs=uMatrix) + call A%eigensys(eigenvals=sVals,reigenvecs=uMatrix) +! call A%diag(evals=sVals,evecs=uMatrix) call sVals%power(scalar) call sVals%diag(A) A = matmul(matmul(uMatrix,A),dagger(uMatrix)) diff --git a/src/mqc_est.F03 b/src/mqc_est.F03 index 69ed509f..fde1be32 100644 --- a/src/mqc_est.F03 +++ b/src/mqc_est.F03 @@ -98,6 +98,8 @@ Module MQC_EST Procedure, Public::getLabel => mqc_integral_array_name !> \brief Set the label of the MQC_SCF_Integral object Procedure, Public::addLabel => mqc_integral_add_name +!> \brief Return a new MQC_SCF_Integral object with distinct Array_type + Procedure, Public::changeType => mqc_integral_change_type !> \brief Return a spin block of MQC_SCF_Integral object as an MQC matrix Procedure, Public::getBlock => mqc_integral_output_block !> \brief Set an MQC_SCF_Integral object as the identity matrix @@ -125,11 +127,14 @@ Module MQC_EST Procedure, Public::powerf => mqc_scf_integral_power_func !> \brief Psuedo inverse of an MQC_SCF_Integral object Procedure, Public::pinv => mqc_scf_integral_pinv +!> \brief Updates an element of an MQC_SCF_Integral object Procedure, Public::put => mqc_integral_put +!> \brief Returns an element of an MQC_SCF_Integral object Procedure, Public::at => mqc_integral_at Procedure, Public::setEList => mqc_integral_set_energy_list Procedure, Public::getEList => mqc_integral_get_energy_list Procedure, Public::deleteEList => mqc_integral_delete_energy_list + Procedure, Public::ipad => mqc_integral_pad !> \brief Returns the sum of elements in the MQC_Matrix Procedure, Public::esum => MQC_Integral_Sum_Elements Procedure, Public::pesum => MQC_Integral_Partial_Sum_Elements @@ -151,8 +156,6 @@ Module MQC_EST Procedure, Public::crossouter2ERI => mqc_integral_integral_crossouter_product !> \brief Returns the diagonal elements of an MQC_SCF_Integral object Procedure, Public::diagonal => mqc_scf_integral_diagonal_elements -!> \brief Pads the blocks of an MQC_SCF_Integral object - Procedure, Public::ipad => mqc_integral_pad End Type ! !> \brief EST intermediate diagonal matrix variable @@ -168,6 +171,8 @@ Module MQC_EST Procedure, Private::hasBeta => mqc_eigenvalues_has_beta Procedure, Private::type => mqc_eigenvalues_array_type Procedure, Private::blockSize => mqc_eigenvalues_dimension +!> \brief Initilizes the MQC_SCF_Eigenvalues object + Procedure, Public::init => mqc_eigenvalues_initialize !> \brief Return the label of the MQC_SCF_Eigenvalues object Procedure, Public::getLabel => mqc_eigenvalues_array_name !> \brief Set the label of the MQC_SCF_Eigenvalues object @@ -176,6 +181,12 @@ Module MQC_EST Procedure, Public::getBlock => mqc_eigenvalues_output_block !> \brief Returns values of an MQC_SCF_Eigenvalues object raised to a power Procedure, Public::power => mqc_scf_eigenvalues_power +!> \brief Alters the value of a specified element in an MQC_SCF_Eigenvalues +!> object + Procedure, Public::put => mqc_eigenvalues_put +!> \brief Appends the specified value to the end of the MQC_SCF_Eigenvalues +!> object + Procedure, Public::push => mqc_eigenvalues_push !> \brief Returns the value of a specified element in an MQC_SCF_Eigenvalues !> object Procedure, Public::at => mqc_eigenvalues_at @@ -189,6 +200,8 @@ Module MQC_EST Procedure, Public::minloc => MQC_Eigenvalues_MinLoc !> \brief Returns the location of the maximum value in MQC_SCF_Eigenvalues Procedure, Public::maxloc => MQC_Eigenvalues_MaxLoc +!> \brief Returns the absolute values of elements in MQC_SCF_Eigenvalues + Procedure, Public::abs => MQC_Eigenvalues_Abs End Type ! !> \brief EST intermediates of single determinant wavefunctions @@ -384,6 +397,16 @@ Module MQC_EST Module Procedure MQC_Eigenvalues_MinLoc End Interface ! +!> \brief Returns the absolute values of elements in MQC_SCF_Eigenvalues + Interface abs + Module Procedure MQC_Eigenvalues_Abs + End Interface +! +!> \brief Returns the identity MQC_SCF_Integral of the specified size + interface intiden + module procedure MQC_Integral_Identity_Func + end interface +! ! !---------------------------------------------------------------- ! | @@ -426,6 +449,11 @@ Module MQC_EST Module Procedure MQC_Integral_Scalar_Multiply End Interface ! +!> \brief Computes the element-wise product of two objects + Interface Operator (.ewp.) + Module Procedure MQC_Element_Integral_Multiply + End Interface +! ! !---------------------------------------------------------------- ! | @@ -1762,6 +1790,106 @@ subroutine mqc_integral_add_name(integral,arrayName) end subroutine mqc_integral_add_name ! ! +! +! PROCEDURE MQC_Integral_Change_Type +! +!> \brief MQC_Integral_Change_Type is a function that copy a MQC_Integral with a +!> less restrict typing. +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_Integral_Change_Type is a function that copy a MQC_Integral with a +!> less restrict typing: 'space' -> 'spin' or 'general; 'spin' -> 'general' +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] Integral +!> \verbatim +!> Integral is Class(MQC_SCF_Integral) +!> The name of the MQC_SCF_Integral variable. +!> \endverbatim +!> +!> \param[in] newType +!> \verbatim +!> newType is Character(Len=*) +!> The new type of the MQC_SCF_Integral: +!> 'space', 'spin' or 'general' +!> \endverbatim +!> +!> \param[out] IntegralOut +!> \verbatim +!> Integral is Type(MQC_SCF_Integral) +!> IntegralOut has the same orbitals as Integral, with newType type. +!> \endverbatim +!> +! +! Authors: +! ======== +!> \author M. M. F. Moraes +!> \date 2025 +! + function mqc_integral_change_type(integral,newType) result(integralOut) +! + implicit none + class(mqc_scf_integral)::integral + type(mqc_scf_integral)::integralOut + character(Len=*)::newType + character(Len=64)::newArrayType,oldArrayType + integer(kind=int64)::nDimAlpha1=0,nDimBeta1=0,nDimTotal1=0,nDimAlpha2=0,nDimBeta2=0,nDimTotal2=0,zeroI + real(kind=real64)::zeroR + complex(kind=real64)::zeroC + parameter(zeroI=0,zeroR=0.0,zeroC=(0.0,0.0)) + type(mqc_matrix)::aa,bb,ab,ba +! + call string_change_case(newType,'L',newArrayType) + oldArrayType = integral%Array_Type + if (oldArrayType.eq.newArrayType) then + integralOut = integral + return + end if + + nDimAlpha1 = integral%blockSize('Alpha',1) + nDimBeta1 = integral%blockSize('Beta',1) + nDimTotal1 = nDimAlpha1 + nDimBeta1 + nDimAlpha2 = integral%blockSize('Alpha',2) + nDimBeta2 = integral%blockSize('Beta',2) + nDimTotal2 = nDimAlpha2 + nDimBeta2 + + select case (newArrayType) + case ('space') + call mqc_error_a('New type is more restrictive than the original one in mqc_integral_change_type',& + 6,'Original Type:',oldArrayType,'New Type:',newArrayType) + case ('spin') + if (oldArrayType.eq.'general') & + call mqc_error_a('New type is more restrictive than the original one in mqc_integral_change_type',& + 6,'Original Type:',oldArrayType,'New Type:',newArrayType) + aa = integral%getBlock('alpha') + call mqc_integral_allocate(integralOut,integral%Array_Name,'spin',aa,aa) + case ('general') + if (oldArrayType.eq.'space') then + aa = integral%getBlock('alpha') + bb = aa + else + aa = integral%getBlock('alpha') + bb = integral%getBlock('beta') + end if + call ab%init(nDimBeta1,nDimAlpha2,zeroR) + call ba%init(nDimAlpha1,nDimBeta2,zeroR) + call mqc_integral_allocate(integralOut,integral%Array_Name,'general',aa,bb,ab,ba) + case default + call mqc_error_A('Unrecognised integral type in mqc_integral_change_type', 6, & + 'New Type:', newArrayType ) + end select + +! + end function mqc_integral_change_type +! +! ! PROCEDURE MQC_eigenvalues_Add_Name ! !> \brief MQC_Eigenvalues_Add_Name is a function that sets the label of an MQC @@ -1985,6 +2113,94 @@ function mqc_eigenvalues_dimension(eigenvalues,label) result(dimBlock) end select ! end function mqc_eigenvalues_dimension + +! +! +! PROCEDURE MQC_Eigenvalues_Initialize +! +!> \brief MQC_Eigenvalues_Initialize is a subroutine that initializes a MQC_SCF_Eigenvalues +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_Eigenvalues_Initialize is a subroutine that initializes a MQC_SCF_Eigenvalues. Default +!> element values are 0.0 in a 'space' type object, or otherwise vector can be initialized with optional +!> argument. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in,out] eigenvalues +!> \verbatim +!> Vector is Class(MQC_Eigenvalues_Initialize) +!> The MQC_Eigenvalues_Initialize to intialize. +!> \endverbatim +!> +!> \param[in] length +!> \verbatim +!> Length is Integer(kind=int64) +!> The length to initialize eigenvalues object +!> \endverbatim +!> +!> +!> \param[in] scalar +!> \verbatim +!> Scalar is Class(*),Optional +!> Value to set each element of Vector. If not present, +!> the value is set to 0.0. Can be of type integer, real, +!> complex or MQC_Scalar. +!> \endverbatim +! +! Authors: +! ======== +!> \author M. M. F. Moraes +!> \date 2025 +! + subroutine mqc_eigenvalues_initialize(eigenvalues,length,scalar,lengthBeta) +! +! Variable Declarations. +! + Implicit None + Real(kind=real64)::Zero=0.0d0 + Character(Len=64)::eigenvaluesType + Integer(kind=int64)::lenAlpha,lenBeta + Integer(kind=int64),Intent(In)::length + type(mqc_vector)::vecAlpha,vecBeta + Class(*),Optional::scalar + Integer(kind=int64),Optional::lengthBeta + Class(mqc_scf_eigenvalues),Intent(InOut)::eigenvalues + + If(Present(lengthBeta)) then + eigenvaluesType='spin' + lenAlpha = length + lenBeta = lengthBeta + else + eigenvaluesType='space' + lenAlpha = length + endif + eigenvalues%Array_Type = eigenvaluesType + + If(Present(Scalar)) then + call vecAlpha%init(lenAlpha,Scalar) + if(eigenvaluesType.eq.'spin') then + call vecBeta%init(lenBeta,Scalar) + endif + else + call vecAlpha%init(lenAlpha) + if(eigenvaluesType.eq.'spin') then + call vecBeta%init(lenBeta) + endif + endIf + eigenvalues%Alpha = vecAlpha + if(eigenvaluesType.eq.'spin') eigenvalues%Beta = vecBeta + + end subroutine mqc_eigenvalues_initialize + + + ! ! ! PROCEDURE MQC_ERIs_Dimension @@ -2554,6 +2770,8 @@ subroutine mqc_integral_allocate(integral,arrayName,arrayType,alpha, & if(myArrayType.eq.'general') then if(present(alphaBeta).and.allocated(alphaBeta).and.size(alphaBeta,1).ne.0.and.size(alphaBeta,2).ne.0) then integral%alphaBeta = alphaBeta + if(integral%nBRows.eq.0) integral%nBRows = size(integral%alphaBeta,1) + if(integral%nACols.eq.0) integral%nACols = size(integral%alphaBeta,2) if(size(integral%alphaBeta,1).ne.integral%nBRows.or.& size(integral%alphaBeta,2).ne.integral%nACols) & call mqc_error_i('Alpha-beta array size mismatch in & @@ -2564,6 +2782,8 @@ subroutine mqc_integral_allocate(integral,arrayName,arrayType,alpha, & endIf if(present(betaAlpha).and.allocated(betaAlpha).and.size(betaAlpha,1).ne.0.and.size(betaAlpha,2).ne.0) then integral%betaAlpha = betaAlpha + if(integral%nARows.eq.0) integral%nARows = size(integral%betaAlpha,1) + if(integral%nBCols.eq.0) integral%nBCols = size(integral%betaAlpha,2) if(size(integral%betaAlpha,1).ne.integral%nARows.or.& size(integral%betaAlpha,2).ne.integral%nBCols) & call mqc_error_i('Beta-Alpha array size mismatch in & @@ -3157,6 +3377,93 @@ subroutine mqc_integral_identity(integral,nAlpha,nBeta,label,nAlpha2,nBeta2) end subroutine mqc_integral_identity ! ! +! PROCEDURE MQC_Integral_Identity_Func +! +!> \brief MQC_Integral_Identity_Func returns an MQC integral type variable +!> as the identity matrix +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_Integral_Identity returns an MQC integral type variable of the as the +!> identity matrix. Optional argument label allows a description to be carried +!> around with the variable. Optional arguments nAlpha2 and nBeta2 allow the +!> columns of spin blocks to be sized differently to the rows. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] NAlpha +!> \verbatim +!> NAlpha is integer +!> The dimension of the alpha spin block. +!> \endverbatim +!> +!> \param[in] NBeta +!> \verbatim +!> NBeta is integer +!> The dimension of the beta spin block. +!> \endverbatim +!> +!> \param[in] Label +!> \verbatim +!> Label is Character(len=*),Optional +!> Description of the MQC_SCF_Integral contents. +!> \endverbatim +!> +!> \param[in] NAlpha2 +!> \verbatim +!> NAlpha2 is Integer,Optional +!> The column number of the alpha spin block, if a +!> rectangular spin block is desired. +!> \endverbatim +!> +!> \param[in] NBeta2 +!> \verbatim +!> NBeta2 is Integer,Optional +!> The column number of the beta spin block, if a +!> rectangular spin block is desired. +!> \endverbatim +! +! Authors: +! ======== +!> \author L. M. Thompson +!> \date 2025 +! + function mqc_integral_identity_func(nAlpha,nBeta,label,nAlpha2,nBeta2) result(integral) +! + implicit none + type(mqc_scf_integral)::integral + integer,intent(in)::nAlpha,nBeta + integer,intent(in),optional::nAlpha2,nBeta2 + character(Len=*),optional,intent(in)::label + Character(Len=64)::myLabel + integer::my_nAlpha2,my_nBeta2 +! + if(present(label)) then + call string_change_case(label,'l',myLabel) + else + myLabel = '' + endIf + if(.not.present(nAlpha2)) then + my_nAlpha2 = nAlpha + else + my_nAlpha2 = nAlpha2 + endIf + if(.not.present(nBeta2)) then + my_nBeta2 = nBeta + else + my_nBeta2 = nBeta2 + endIf +! + call integral%identity(nAlpha,nBeta,myLabel,my_nAlpha2,my_nBeta2) +! + end function mqc_integral_identity_func +! +! ! PROCEDURE MQC_Integral_Initialize ! !> \brief MQC_Integral_Initialize is used to initialize an MQC integral type @@ -3235,7 +3542,7 @@ subroutine mqc_integral_initialize(integral,nAlpha,nBeta,scalar,label,nAlpha2,nB character(Len=*),optional,intent(in)::label class(*),optional,intent(in)::Scalar integer,intent(in),optional::nAlpha2,nBeta2 - type(mqc_matrix)::tmpMatrixAlpha,tmpMatrixBeta + type(mqc_matrix)::tmpMatrixAlpha,tmpMatrixBeta,tmpMatrixAlphaBeta,tmpMatrixBetaAlpha Character(Len=64)::myLabel integer::my_nAlpha2,my_nBeta2 ! @@ -3260,7 +3567,14 @@ subroutine mqc_integral_initialize(integral,nAlpha,nBeta,scalar,label,nAlpha2,nB call mqc_integral_allocate(integral,myLabel,'space',tmpMatrixAlpha) else call tmpMatrixBeta%init(nBeta,my_nBeta2,scalar) - call mqc_integral_allocate(integral,myLabel,'spin',tmpMatrixAlpha,tmpMatrixBeta) + if(my_nBeta2.eq.0.or.nBeta.eq.0.or.my_nAlpha2.eq.0.or.nAlpha.eq.0) then + call tmpMatrixAlphaBeta%init(nAlpha,my_nBeta2,scalar) + call tmpMatrixBetaAlpha%init(nBeta,my_nAlpha2,scalar) + call mqc_integral_allocate(integral,myLabel,'general',tmpMatrixAlpha,tmpMatrixBeta,& + tmpMatrixAlphaBeta,tmpMatrixBetaAlpha) + else + call mqc_integral_allocate(integral,myLabel,'spin',tmpMatrixAlpha,tmpMatrixBeta) + endIf endIf else call tmpMatrixAlpha%init(nAlpha,my_nAlpha2) @@ -3268,7 +3582,14 @@ subroutine mqc_integral_initialize(integral,nAlpha,nBeta,scalar,label,nAlpha2,nB call mqc_integral_allocate(integral,myLabel,'space',tmpMatrixAlpha) else call tmpMatrixBeta%init(nBeta,my_nBeta2) - call mqc_integral_allocate(integral,myLabel,'spin',tmpMatrixAlpha,tmpMatrixBeta) + if(my_nBeta2.eq.0.or.nBeta.eq.0.or.my_nAlpha2.eq.0.or.nAlpha.eq.0) then + call tmpMatrixAlphaBeta%init(nAlpha,my_nBeta2) + call tmpMatrixBetaAlpha%init(nBeta,my_nAlpha2) + call mqc_integral_allocate(integral,myLabel,'general',tmpMatrixAlpha,tmpMatrixBeta,& + tmpMatrixAlphaBeta,tmpMatrixBetaAlpha) + else + call mqc_integral_allocate(integral,myLabel,'spin',tmpMatrixAlpha,tmpMatrixBeta) + endIf endIf endIf ! @@ -4172,7 +4493,7 @@ subroutine mqc_integral_combine_orbitals(integral,occupieds,virtuals,axis,arrayN call alpha%mput(tmp1,[0],[1,size(tmp1,2)]) call alpha%mput(tmp2,[0],[size(tmp1,2)+1,size(tmp1,2)+size(tmp2,2)]) endIf - call mqc_integral_allocate(integral,arrayName,'space',alpha,nDimA1=nTA1,nDimA2=nTA1,& + call mqc_integral_allocate(integral,arrayName,'space',alpha,nDimA1=nTA1,nDimA2=nTA2,& nDimB1=nTB1,nDimB2=nTB2) elseIf(occupieds%type().eq.'space'.and.virtuals%type().eq.'spin') then tmp1 = occupieds%getBlock('alpha') @@ -4193,7 +4514,7 @@ subroutine mqc_integral_combine_orbitals(integral,occupieds,virtuals,axis,arrayN call beta%mput(tmp1,[0],[1,size(tmp1,2)]) call beta%mput(tmp3,[0],[size(tmp1,2)+1,size(tmp1,2)+size(tmp3,2)]) endIf - call mqc_integral_allocate(integral,arrayName,'spin',alpha,beta,nDimA1=nTA1,nDimA2=nTA1,& + call mqc_integral_allocate(integral,arrayName,'spin',alpha,beta,nDimA1=nTA1,nDimA2=nTA2,& nDimB1=nTB1,nDimB2=nTB2) elseIf(occupieds%type().eq.'spin'.and.virtuals%type().eq.'space') then tmp1 = occupieds%getBlock('alpha') @@ -4214,7 +4535,7 @@ subroutine mqc_integral_combine_orbitals(integral,occupieds,virtuals,axis,arrayN call beta%mput(tmp2,[0],[1,size(tmp2,2)]) call beta%mput(tmp3,[0],[size(tmp2,2)+1,size(tmp2,2)+size(tmp3,2)]) endIf - call mqc_integral_allocate(integral,arrayName,'spin',alpha,beta,nDimA1=nTA1,nDimA2=nTA1,& + call mqc_integral_allocate(integral,arrayName,'spin',alpha,beta,nDimA1=nTA1,nDimA2=nTA2,& nDimB1=nTB1,nDimB2=nTB2) elseIf(occupieds%type().eq.'spin'.and.virtuals%type().eq.'spin') then tmp1 = occupieds%getBlock('alpha') @@ -4236,7 +4557,7 @@ subroutine mqc_integral_combine_orbitals(integral,occupieds,virtuals,axis,arrayN call beta%mput(tmp2,[0],[1,size(tmp2,2)]) call beta%mput(tmp4,[0],[size(tmp2,2)+1,size(tmp2,2)+size(tmp4,2)]) endIf - call mqc_integral_allocate(integral,arrayName,'spin',alpha,beta,nDimA1=nTA1,nDimA2=nTA1,& + call mqc_integral_allocate(integral,arrayName,'spin',alpha,beta,nDimA1=nTA1,nDimA2=nTA2,& nDimB1=nTB1,nDimB2=nTB2) elseIf(occupieds%type().eq.'space'.and.virtuals%type().eq.'general') then tmp1 = occupieds%getBlock('alpha') @@ -4268,7 +4589,7 @@ subroutine mqc_integral_combine_orbitals(integral,occupieds,virtuals,axis,arrayN call betaAlpha%mput(tmp5,[0],[size(tmp1,2)+1,size(tmp1,2)+size(tmp5,2)]) endIf call mqc_integral_allocate(integral,arrayName,'general',alpha,beta,alphaBeta,betaAlpha,& - nDimA1=nTA1,nDimA2=nTA1,nDimB1=nTB1,nDimB2=nTB2) + nDimA1=nTA1,nDimA2=nTA2,nDimB1=nTB1,nDimB2=nTB2) elseIf(occupieds%type().eq.'general'.and.virtuals%type().eq.'space') then tmp1 = occupieds%getBlock('alpha') tmp2 = occupieds%getBlock('beta') @@ -4299,7 +4620,7 @@ subroutine mqc_integral_combine_orbitals(integral,occupieds,virtuals,axis,arrayN call betaAlpha%mput(tmp4,[0],[1,size(tmp4,2)]) endIf call mqc_integral_allocate(integral,arrayName,'general',alpha,beta,alphaBeta,betaAlpha,& - nDimA1=nTA1,nDimA2=nTA1,nDimB1=nTB1,nDimB2=nTB2) + nDimA1=nTA1,nDimA2=nTA2,nDimB1=nTB1,nDimB2=nTB2) elseIf(occupieds%type().eq.'spin'.and.virtuals%type().eq.'general') then tmp1 = occupieds%getBlock('alpha') tmp2 = occupieds%getBlock('beta') @@ -4331,7 +4652,7 @@ subroutine mqc_integral_combine_orbitals(integral,occupieds,virtuals,axis,arrayN call betaAlpha%mput(tmp6,[0],[size(tmp2,2)+1,size(tmp2,2)+size(tmp6,2)]) endIf call mqc_integral_allocate(integral,arrayName,'general',alpha,beta,alphaBeta,betaAlpha,& - nDimA1=nTA1,nDimA2=nTA1,nDimB1=nTB1,nDimB2=nTB2) + nDimA1=nTA1,nDimA2=nTA2,nDimB1=nTB1,nDimB2=nTB2) elseIf(occupieds%type().eq.'general'.and.virtuals%type().eq.'spin') then tmp1 = occupieds%getBlock('alpha') tmp2 = occupieds%getBlock('beta') @@ -4363,7 +4684,7 @@ subroutine mqc_integral_combine_orbitals(integral,occupieds,virtuals,axis,arrayN call betaAlpha%mput(tmp4,[0],[1,size(tmp4,2)]) endIf call mqc_integral_allocate(integral,arrayName,'general',alpha,beta,alphaBeta,betaAlpha,& - nDimA1=nTA1,nDimA2=nTA1,nDimB1=nTB1,nDimB2=nTB2) + nDimA1=nTA1,nDimA2=nTA2,nDimB1=nTB1,nDimB2=nTB2) elseIf(occupieds%type().eq.'general'.and.virtuals%type().eq.'general') then tmp1 = occupieds%getBlock('alpha') tmp2 = occupieds%getBlock('beta') @@ -4401,7 +4722,7 @@ subroutine mqc_integral_combine_orbitals(integral,occupieds,virtuals,axis,arrayN call betaAlpha%mput(tmp8,[0],[size(tmp4,2)+1,size(tmp4,2)+size(tmp8,2)]) endIf call mqc_integral_allocate(integral,arrayName,'general',alpha,beta,alphaBeta,betaAlpha,& - nDimA1=nTA1,nDimA2=nTA1,nDimB1=nTB1,nDimB2=nTB2) + nDimA1=nTA1,nDimA2=nTA2,nDimB1=nTB1,nDimB2=nTB2) endIf ! end subroutine mqc_integral_combine_orbitals @@ -4726,15 +5047,16 @@ function mqc_eigenvalues_output_block(eigenvalues,blockName) result(vectorOut) if (eigenvalues%type().eq.'space') then if (eigenvalues%hasAlpha()) then if(mqc_vector_haveInteger(eigenvalues%alpha)) then - call vectorOut%init(nDimAlpha,zeroI) + call vectorOut%init(2*nDimAlpha,zeroI) elseif(mqc_vector_haveReal(eigenvalues%alpha)) then - call vectorOut%init(nDimAlpha,zeroR) + call vectorOut%init(2*nDimAlpha,zeroR) elseIf(mqc_vector_haveComplex(eigenvalues%alpha)) then - call vectorOut%init(nDimAlpha,zeroC) + call vectorOut%init(2*nDimAlpha,zeroC) else call mqc_error('unrecognized integral data type in mqc_integral_output_block') endIf call vectorOut%vput(eigenvalues%alpha,1) + call vectorOut%vput(eigenvalues%alpha,nDimAlpha+1) else call mqc_error_L('block does not exist in mqc_eigenvalues_output_block', 6, & 'eigenvalues%hasAlpha()', eigenvalues%hasAlpha() ) @@ -6265,50 +6587,217 @@ function mqc_matrix_integral_multiply(matrixA,integralB,label) result(integralOu end function mqc_matrix_integral_multiply ! ! -! PROCEDURE MQC_Integral_Sum +! PROCEDURE MQC_Integral_Pad ! -!> \brief MQC_Integral_Sum is used to sum two MQC integral type variables +!> \brief MQC_Integral_Pad is used to increase the dimension +!> of a MQC integral variable ! !> \par Purpose: ! ============= !> !> \verbatim !> -!> MQC_Integral_Sum is used to sum two MQC integral type variables. +!> MQC_Matrix_Pad is used to multiply an MQC matrix type variable +!> with a MQC integral and returns an MQC integral. !> !> \endverbatim ! ! Arguments: ! ========== -!> \param[in] IntegralA +!> \param[in] Integral !> \verbatim -!> IntegralA is type(mqc_matrix) -!> The first MQC integral variable to sum. +!> Integral is type(mqc_scf_integral) +!> Initial MQC integral variable. !> \endverbatim !> -!> \param[in] IntegralB +!> \param[in] Axis1 !> \verbatim -!> IntegralB is type(mqc_scf_integral) -!> The second MQC integral variable to sum. +!> Axis1 is integer(kind=int64) +!> Size of padded data. +!> If Axis2 is present, Axis1 is used only for rows. +!> Otherwise Axis1 is used both for rows and columns. +!> \endverbatim +!> +!> \param[in] Axis2 +!> \verbatim +!> Axis2 is integer(kind=int64),optional +!> Size of padded columns. +!> \endverbatim +!> +!> \param[in] BlockName +!> \verbatim +!> BlockName is Character(Len=*),optional,default='full' +!> Define which set of spin-orbitals should be modiffied. +!> In the general case, the off-diagonal blocks (alpha-beta +!> and beta-alpha) are modified based on the alpha-alpha and +!> beta-beta dimensions. +!> +!> 1. BlockName = 'full' outputs the full MQC_TwoERIs variable. +!> 2. BlockName = 'alpha', 'alpha-alpha' or 'aaaa' outputs the alpha-alpha spin block. +!> 3. BlockName = 'beta', 'beta-beta' or 'bbbb' outputs the beta-beta spin block. !> \endverbatim ! ! Authors: ! ======== -!> \author L. M. Thompson, M. M. F. Moraes -!> \date 2018, 2024 +!> \author M. M. F. Moraes +!> \date 2025 ! - function mqc_integral_sum(integralA,integralB) result(integralOut) + function mqc_integral_pad(integral,axis1,axis2,blockname) result(integralOut) ! implicit none - type(mqc_scf_integral),intent(in)::integralA,integralB + class(mqc_scf_integral),intent(in)::integral + integer(kind=int64),intent(in)::axis1 + integer(kind=int64),optional,intent(in)::axis2 + Character(Len=*),optional,intent(in)::blockname + integer(kind=int64)::ax1,ax2 type(mqc_scf_integral)::integralOut + type(mqc_matrix)::tmpMatrixAlpha,tmpMatrixBeta,tmpMatrixAlphaBeta,tmpMatrixBetaAlpha - integer(kind=int64)::nA1,nA2,nB1,nB2 + Character(Len=64)::myLabel ! - if(size(integralA,'alpha',1).ne.size(integralB,'alpha',1).or.size(integralA,'alpha',2).ne.& - size(integralA,'alpha',2).or.size(integralA,'beta',1).ne.size(integralA,'beta',1).or.& - size(integralA,'beta',2).ne.size(integralA,'beta',2)) call mqc_error_i('Integral blocks are & - ¬ conformable in mqc_integral_sum',6,'size(integralA,alpha,1)-size(integralB,alpha,1)',& + if(present(blockname)) then + call string_change_case(blockname,'l',myLabel) + else + myLabel = 'full' + endIf + + ax1 = axis1 + if(present(axis2)) then + ax2 = axis2 + else + ax2 = axis1 + endIf + + select case(integral%array_type) + case('space') + select case(myLabel) + case('full') + tmpMatrixAlpha = integral%getBlock('alpha') + tmpMatrixAlpha = tmpMatrixAlpha%mpad(ax1,ax2) + call mqc_integral_allocate(integralOut,integral%getLabel(),'space',tmpMatrixAlpha) + case('alpha') + tmpMatrixBeta = integral%getBlock('alpha') + tmpMatrixAlpha = tmpMatrixBeta%mpad(ax1,ax2) + call mqc_integral_allocate(integralOut,integral%getLabel(),'spin',tmpMatrixAlpha,tmpMatrixBeta) + case('beta') + tmpMatrixAlpha = integral%getBlock('alpha') + tmpMatrixBeta = tmpMatrixAlpha%mpad(ax1,ax2) + call mqc_integral_allocate(integralOut,integral%getLabel(),'spin',tmpMatrixAlpha,tmpMatrixBeta) + case default + call mqc_error_A('Unknown label in mqc_integral_pad', 6, & + 'label', myLabel ) + end select + case('spin') + select case(myLabel) + case('full') + tmpMatrixAlpha = integral%getBlock('alpha') + tmpMatrixAlpha = tmpMatrixAlpha%mpad(ax1,ax2) + tmpMatrixBeta = integral%getBlock('beta') + tmpMatrixBeta = tmpMatrixBeta%mpad(ax1,ax2) + call mqc_integral_allocate(integralOut,integral%getLabel(),'spin',tmpMatrixAlpha,tmpMatrixBeta) + case('alpha') + tmpMatrixAlpha = integral%getBlock('alpha') + tmpMatrixAlpha = tmpMatrixAlpha%mpad(ax1,ax2) + tmpMatrixBeta = integral%getBlock('beta') + call mqc_integral_allocate(integralOut,integral%getLabel(),'spin',tmpMatrixAlpha,tmpMatrixBeta) + case('beta') + tmpMatrixAlpha = integral%getBlock('alpha') + tmpMatrixBeta = integral%getBlock('beta') + tmpMatrixBeta = tmpMatrixBeta%mpad(ax1,ax2) + call mqc_integral_allocate(integralOut,integral%getLabel(),'spin',tmpMatrixAlpha,tmpMatrixBeta) + case default + call mqc_error_A('Unknown label in mqc_integral_pad', 6, & + 'label', myLabel ) + end select + case('general') + select case(myLabel) + case('full') + tmpMatrixAlpha = integral%getBlock('alpha') + tmpMatrixAlpha = tmpMatrixAlpha%mpad(ax1,ax2) + tmpMatrixBeta = integral%getBlock('beta') + tmpMatrixBeta = tmpMatrixBeta%mpad(ax1,ax2) + tmpMatrixAlphaBeta = integral%getBlock('alphabeta') + tmpMatrixAlphaBeta = tmpMatrixAlphaBeta%mpad(ax1,ax2) + tmpMatrixBetaAlpha = integral%getBlock('betaalpha') + tmpMatrixBetaAlpha = tmpMatrixBetaAlpha%mpad(ax1,ax2) + call mqc_integral_allocate(integralOut,integral%getLabel(),'general',tmpMatrixAlpha, & + tmpMatrixBeta,tmpMatrixAlphaBeta,tmpMatrixBetaAlpha) + case('alpha') + tmpMatrixAlpha = integral%getBlock('alpha') + tmpMatrixAlpha = tmpMatrixAlpha%mpad(ax1,ax2) + tmpMatrixBeta = integral%getBlock('beta') + tmpMatrixAlphaBeta = integral%getBlock('alphabeta') + tmpMatrixAlphaBeta = tmpMatrixAlphaBeta%mpad(integral%nBRows,ax2) + tmpMatrixBetaAlpha = integral%getBlock('betaalpha') + tmpMatrixBetaAlpha = tmpMatrixBetaAlpha%mpad(ax1,integral%nBCols) + call mqc_integral_allocate(integralOut,integral%getLabel(),'general',tmpMatrixAlpha, & + tmpMatrixBeta,tmpMatrixAlphaBeta,tmpMatrixBetaAlpha) + case('beta') + tmpMatrixAlpha = integral%getBlock('alpha') + tmpMatrixBeta = integral%getBlock('beta') + tmpMatrixBeta = tmpMatrixBeta%mpad(ax1,ax2) + tmpMatrixAlphaBeta = integral%getBlock('alphabeta') + tmpMatrixAlphaBeta = tmpMatrixAlphaBeta%mpad(ax1,integral%nACols) + tmpMatrixBetaAlpha = integral%getBlock('betaalpha') + tmpMatrixBetaAlpha = tmpMatrixBetaAlpha%mpad(integral%nARows,ax2) + call mqc_integral_allocate(integralOut,integral%getLabel(),'general',tmpMatrixAlpha, & + tmpMatrixBeta,tmpMatrixAlphaBeta,tmpMatrixBetaAlpha) + case default + call mqc_error_A('Unknown label in mqc_integral_pad', 6, & + 'label', myLabel ) + end select + case default + call mqc_error_A('Unknown integral type in mqc_integral_pad', 6, & + 'integral%array_type', integral%array_type ) + end select +! + end function mqc_integral_pad +! +! +! PROCEDURE MQC_Integral_Sum +! +!> \brief MQC_Integral_Sum is used to sum two MQC integral type variables +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_Integral_Sum is used to sum two MQC integral type variables. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] IntegralA +!> \verbatim +!> IntegralA is type(mqc_matrix) +!> The first MQC integral variable to sum. +!> \endverbatim +!> +!> \param[in] IntegralB +!> \verbatim +!> IntegralB is type(mqc_scf_integral) +!> The second MQC integral variable to sum. +!> \endverbatim +! +! Authors: +! ======== +!> \author L. M. Thompson, M. M. F. Moraes +!> \date 2018, 2024 +! + function mqc_integral_sum(integralA,integralB) result(integralOut) +! + implicit none + type(mqc_scf_integral),intent(in)::integralA,integralB + type(mqc_scf_integral)::integralOut + type(mqc_matrix)::tmpMatrixAlpha,tmpMatrixBeta,tmpMatrixAlphaBeta,tmpMatrixBetaAlpha + integer(kind=int64)::nA1,nA2,nB1,nB2 +! + if(size(integralA,'alpha',1).ne.size(integralB,'alpha',1).or.size(integralA,'alpha',2).ne.& + size(integralA,'alpha',2).or.size(integralA,'beta',1).ne.size(integralA,'beta',1).or.& + size(integralA,'beta',2).ne.size(integralA,'beta',2)) call mqc_error_i('Integral blocks are & + ¬ conformable in mqc_integral_sum',6,'size(integralA,alpha,1)-size(integralB,alpha,1)',& size(integralA,'alpha',1)-size(integralB,'alpha',1),'size(integralA,alpha,2)-size(integralB,alpha,2)',& size(integralA,'alpha',2)-size(integralB,'alpha',2),'size(integralA,beta,1)-size(integralB,beta,1)',& size(integralA,'beta',1)-size(integralB,'beta',1),'size(integralA,beta,2)-size(integralB,beta,2)',& @@ -6524,243 +7013,76 @@ function mqc_integral_sum(integralA,integralB) result(integralOut) end function mqc_integral_sum ! ! -! PROCEDURE MQC_Integral_Pad +! PROCEDURE MQC_ERI_Sum ! -!> \brief MQC_Integral_Pad is used to increase the dimension -!> of a MQC integral variable +!> \brief MQC_ERI_Sum is used to sum two MQC twoERIs type +!> variables ! !> \par Purpose: ! ============= !> !> \verbatim !> -!> MQC_Matrix_Pad is used to multiply an MQC matrix type variable -!> with a MQC integral and returns an MQC integral. +!> MQC_ERI_Sum is used to sum two MQC twoERIs type variables. !> !> \endverbatim ! ! Arguments: ! ========== -!> \param[in] Integral -!> \verbatim -!> Integral is type(mqc_scf_integral) -!> Initial MQC integral variable. -!> \endverbatim -!> -!> \param[in] Axis1 -!> \verbatim -!> Axis1 is integer(kind=int64) -!> Size of padded data. -!> If Axis2 is present, Axis1 is used only for rows. -!> Otherwise Axis1 is used both for rows and columns. -!> \endverbatim -!> -!> \param[in] Axis2 +!> \param[in] ERIA !> \verbatim -!> Axis2 is integer(kind=int64),optional -!> Size of padded columns. +!> ERIA is type(mqc_twoERIs) +!> The MQC twoERIs variable to sum. !> \endverbatim -!> -!> \param[in] BlockName +!> +!> \param[in] ERIB !> \verbatim -!> BlockName is Character(Len=*),optional,default='full' -!> Define which set of spin-orbitals should be modiffied. -!> In the general case, the off-diagonal blocks (alpha-beta -!> and beta-alpha) are modified based on the alpha-alpha and -!> beta-beta dimensions. -!> -!> 1. BlockName = 'full' outputs the full MQC_TwoERIs variable. -!> 2. BlockName = 'alpha', 'alpha-alpha' or 'aaaa' outputs the alpha-alpha spin block. -!> 3. BlockName = 'beta', 'beta-beta' or 'bbbb' outputs the beta-beta spin block. +!> ERIB is type(mqc_twoERIs) +!> The MQC twoERIs variable to sum. !> \endverbatim ! ! Authors: ! ======== !> \author M. M. F. Moraes -!> \date 2025 +!> \date 2024 ! - function mqc_integral_pad(integral,axis1,axis2,blockname) result(integralOut) + function mqc_eri_sum(ERIA,ERIB) result(ERIOut) ! implicit none - class(mqc_scf_integral),intent(in)::integral - integer(kind=int64),intent(in)::axis1 - integer(kind=int64),optional,intent(in)::axis2 - Character(Len=*),optional,intent(in)::blockname - integer(kind=int64)::ax1,ax2 - type(mqc_scf_integral)::integralOut - - type(mqc_matrix)::tmpMatrixAlpha,tmpMatrixBeta,tmpMatrixAlphaBeta,tmpMatrixBetaAlpha - Character(Len=64)::myLabel + type(mqc_twoERIs),intent(in)::ERIA,ERIB + type(mqc_twoERIs)::ERIOut + type(mqc_r4tensor)::aaaa,bbbb,aabb,bbaa,abab,abba,baab,baba,& + aaab,aaba,abaa,baaa,bbba,bbab,babb,abbb ! - if(present(blockname)) then - call string_change_case(blockname,'l',myLabel) - else - myLabel = 'full' - endIf - - ax1 = axis1 - if(present(axis2)) then - ax2 = axis2 - else - ax2 = axis1 - endIf - - select case(integral%array_type) + select case(ERIA%integraltype) case('space') - select case(myLabel) - case('full') - tmpMatrixAlpha = integral%getBlock('alpha') - tmpMatrixAlpha = tmpMatrixAlpha%mpad(ax1,ax2) - call mqc_integral_allocate(integralOut,integral%getLabel(),'space',tmpMatrixAlpha) - case('alpha') - tmpMatrixBeta = integral%getBlock('alpha') - tmpMatrixAlpha = tmpMatrixBeta%mpad(ax1,ax2) - call mqc_integral_allocate(integralOut,integral%getLabel(),'spin',tmpMatrixAlpha,tmpMatrixBeta) - case('beta') - tmpMatrixAlpha = integral%getBlock('alpha') - tmpMatrixBeta = tmpMatrixAlpha%mpad(ax1,ax2) - call mqc_integral_allocate(integralOut,integral%getLabel(),'spin',tmpMatrixAlpha,tmpMatrixBeta) + select case (ERIB%integraltype) + case('space') + aaaa = ERIA%getBlock('aaaa')+ERIB%getBlock('aaaa') + call mqc_twoeris_allocate(ERIOut,'full','space',aaaa) + case('spin') + aaaa = ERIA%getBlock('aaaa')+ERIB%getBlock('aaaa') + bbbb = ERIA%getBlock('aaaa')+ERIB%getBlock('bbbb') + aabb = ERIA%getBlock('aaaa')+ERIB%getBlock('aabb') + bbaa = ERIA%getBlock('aaaa')+ERIB%getBlock('bbaa') + call mqc_twoeris_allocate(ERIOut,'full','spin',aaaa,bbbb,& + aabb,bbaa) + case('general') + aaaa = ERIA%getBlock('aaaa')+ERIB%getBlock('aaaa') + bbbb = ERIA%getBlock('aaaa')+ERIB%getBlock('bbbb') + aabb = ERIA%getBlock('aaaa')+ERIB%getBlock('aabb') + bbaa = ERIA%getBlock('aaaa')+ERIB%getBlock('bbaa') + call mqc_twoeris_allocate(ERIOut,'full','general',aaaa,bbbb,& + aabb,bbaa,& + ERIB%getBlock('abab'), ERIB%getBlock('abba'),& + ERIB%getBlock('baab'), ERIB%getBlock('baba'),& + ERIB%getBlock('aaab'), ERIB%getBlock('aaba'),& + ERIB%getBlock('abaa'), ERIB%getBlock('baaa'),& + ERIB%getBlock('bbba'), ERIB%getBlock('bbab'),& + ERIB%getBlock('babb'), ERIB%getBlock('abbb')) case default - call mqc_error_A('Unknown label in mqc_integral_pad', 6, & - 'label', myLabel ) - end select - case('spin') - select case(myLabel) - case('full') - tmpMatrixAlpha = integral%getBlock('alpha') - tmpMatrixAlpha = tmpMatrixAlpha%mpad(ax1,ax2) - tmpMatrixBeta = integral%getBlock('beta') - tmpMatrixBeta = tmpMatrixBeta%mpad(ax1,ax2) - call mqc_integral_allocate(integralOut,integral%getLabel(),'spin',tmpMatrixAlpha,tmpMatrixBeta) - case('alpha') - tmpMatrixAlpha = integral%getBlock('alpha') - tmpMatrixAlpha = tmpMatrixAlpha%mpad(ax1,ax2) - tmpMatrixBeta = integral%getBlock('beta') - call mqc_integral_allocate(integralOut,integral%getLabel(),'spin',tmpMatrixAlpha,tmpMatrixBeta) - case('beta') - tmpMatrixAlpha = integral%getBlock('alpha') - tmpMatrixBeta = integral%getBlock('beta') - tmpMatrixBeta = tmpMatrixBeta%mpad(ax1,ax2) - call mqc_integral_allocate(integralOut,integral%getLabel(),'spin',tmpMatrixAlpha,tmpMatrixBeta) - case default - call mqc_error_A('Unknown label in mqc_integral_pad', 6, & - 'label', myLabel ) - end select - case('general') - select case(myLabel) - case('full') - tmpMatrixAlpha = integral%getBlock('alpha') - tmpMatrixAlpha = tmpMatrixAlpha%mpad(ax1,ax2) - tmpMatrixBeta = integral%getBlock('beta') - tmpMatrixBeta = tmpMatrixBeta%mpad(ax1,ax2) - tmpMatrixAlphaBeta = integral%getBlock('alphabeta') - tmpMatrixAlphaBeta = tmpMatrixAlphaBeta%mpad(ax1,ax2) - tmpMatrixBetaAlpha = integral%getBlock('betaalpha') - tmpMatrixBetaAlpha = tmpMatrixBetaAlpha%mpad(ax1,ax2) - call mqc_integral_allocate(integralOut,integral%getLabel(),'general',tmpMatrixAlpha, & - tmpMatrixBeta,tmpMatrixAlphaBeta,tmpMatrixBetaAlpha) - case('alpha') - tmpMatrixAlpha = integral%getBlock('alpha') - tmpMatrixAlpha = tmpMatrixAlpha%mpad(ax1,ax2) - tmpMatrixBeta = integral%getBlock('beta') - tmpMatrixAlphaBeta = integral%getBlock('alphabeta') - tmpMatrixAlphaBeta = tmpMatrixAlphaBeta%mpad(integral%nBRows,ax2) - tmpMatrixBetaAlpha = integral%getBlock('betaalpha') - tmpMatrixBetaAlpha = tmpMatrixBetaAlpha%mpad(ax1,integral%nBCols) - call mqc_integral_allocate(integralOut,integral%getLabel(),'general',tmpMatrixAlpha, & - tmpMatrixBeta,tmpMatrixAlphaBeta,tmpMatrixBetaAlpha) - case('beta') - tmpMatrixAlpha = integral%getBlock('alpha') - tmpMatrixBeta = integral%getBlock('beta') - tmpMatrixBeta = tmpMatrixBeta%mpad(ax1,ax2) - tmpMatrixAlphaBeta = integral%getBlock('alphabeta') - tmpMatrixAlphaBeta = tmpMatrixAlphaBeta%mpad(ax1,integral%nACols) - tmpMatrixBetaAlpha = integral%getBlock('betaalpha') - tmpMatrixBetaAlpha = tmpMatrixBetaAlpha%mpad(integral%nARows,ax2) - call mqc_integral_allocate(integralOut,integral%getLabel(),'general',tmpMatrixAlpha, & - tmpMatrixBeta,tmpMatrixAlphaBeta,tmpMatrixBetaAlpha) - case default - call mqc_error_A('Unknown label in mqc_integral_pad', 6, & - 'label', myLabel ) - end select - case default - call mqc_error_A('Unknown integral type in mqc_integral_pad', 6, & - 'integral%array_type', integral%array_type ) - end select -! - end function mqc_integral_pad -! -! -! PROCEDURE MQC_ERI_Sum -! -!> \brief MQC_ERI_Sum is used to sum two MQC twoERIs type -!> variables -! -!> \par Purpose: -! ============= -!> -!> \verbatim -!> -!> MQC_ERI_Sum is used to sum two MQC twoERIs type variables. -!> -!> \endverbatim -! -! Arguments: -! ========== -!> \param[in] ERIA -!> \verbatim -!> ERIA is type(mqc_twoERIs) -!> The MQC twoERIs variable to sum. -!> \endverbatim -!> -!> \param[in] ERIB -!> \verbatim -!> ERIB is type(mqc_twoERIs) -!> The MQC twoERIs variable to sum. -!> \endverbatim -! -! Authors: -! ======== -!> \author M. M. F. Moraes -!> \date 2024 -! - function mqc_eri_sum(ERIA,ERIB) result(ERIOut) -! - implicit none - type(mqc_twoERIs),intent(in)::ERIA,ERIB - type(mqc_twoERIs)::ERIOut - type(mqc_r4tensor)::aaaa,bbbb,aabb,bbaa,abab,abba,baab,baba,& - aaab,aaba,abaa,baaa,bbba,bbab,babb,abbb -! - select case(ERIA%integraltype) - case('space') - select case (ERIB%integraltype) - case('space') - aaaa = ERIA%getBlock('aaaa')+ERIB%getBlock('aaaa') - call mqc_twoeris_allocate(ERIOut,'full','space',aaaa) - case('spin') - aaaa = ERIA%getBlock('aaaa')+ERIB%getBlock('aaaa') - bbbb = ERIA%getBlock('aaaa')+ERIB%getBlock('bbbb') - aabb = ERIA%getBlock('aaaa')+ERIB%getBlock('aabb') - bbaa = ERIA%getBlock('aaaa')+ERIB%getBlock('bbaa') - call mqc_twoeris_allocate(ERIOut,'full','spin',aaaa,bbbb,& - aabb,bbaa) - case('general') - aaaa = ERIA%getBlock('aaaa')+ERIB%getBlock('aaaa') - bbbb = ERIA%getBlock('aaaa')+ERIB%getBlock('bbbb') - aabb = ERIA%getBlock('aaaa')+ERIB%getBlock('aabb') - bbaa = ERIA%getBlock('aaaa')+ERIB%getBlock('bbaa') - call mqc_twoeris_allocate(ERIOut,'full','general',aaaa,bbbb,& - aabb,bbaa,& - ERIB%getBlock('abab'), ERIB%getBlock('abba'),& - ERIB%getBlock('baab'), ERIB%getBlock('baba'),& - ERIB%getBlock('aaab'), ERIB%getBlock('aaba'),& - ERIB%getBlock('abaa'), ERIB%getBlock('baaa'),& - ERIB%getBlock('bbba'), ERIB%getBlock('bbab'),& - ERIB%getBlock('babb'), ERIB%getBlock('abbb')) - case default - call mqc_error_A('Unknown twoERI type in mqc_eri_difference', 6, & - 'ERIB%integraltype', ERIB%integraltype ) + call mqc_error_A('Unknown twoERI type in mqc_eri_difference', 6, & + 'ERIB%integraltype', ERIB%integraltype ) end select case('spin') select case (ERIB%integraltype) @@ -7451,9 +7773,21 @@ function mqc_integral_integral_multiply(integralA,integralB,label) result(integr nBB1 = size(integralB,'beta',1) nBB2 = size(integralB,'beta',2) - if(nAA2.ne.nBA1.or.nAB2.ne.nBB1) call mqc_error_i('Integral dimensions are & + if (nAA2.eq.0.or.nBA1.eq.0) then + if (nAB2.ne.0.and.nBB1.ne.0) then + if (nAB2.ne.nBB1) call mqc_error_i('Integral dimensions are & + ¬ conformable in mqc_integral_integral_multiply',6,'nAA2',nAA2,'nBA1',nBA1,& + 'nAB2',nAB2,'nBB1',nBB1) + end if + elseif (nAB2.eq.0.or.nBB1.eq.0) then + if (nAA2.ne.nBA1) call mqc_error_i('Integral dimensions are & + ¬ conformable in mqc_integral_integral_multiply',6,'nAA2',nAA2,'nBA1',nBA1,& + 'nAB2',nAB2,'nBB1',nBB1) + elseif(nAA2.ne.nBA1.or.nAB2.ne.nBB1) then + call mqc_error_i('Integral dimensions are & ¬ conformable in mqc_integral_integral_multiply',6,'nAA2',nAA2,'nBA1',nBA1,& 'nAB2',nAB2,'nBB1',nBB1) + end if select case (integralA%array_type) case('space') @@ -11699,11 +12033,11 @@ subroutine mqc_scf_integral_generalized_eigensystem(integralA,integralB,eVals,rE else tmpMat1 = integralA%getblock('alpha') tmpMat2 = integralA%getblock('beta') - call tmpMat1%eigensys(eigenvals=tmpVec1,reigenvecs=tmpMat3,leigenvecs=tmpMat4) - call tmpMat2%eigensys(eigenvals=tmpVec2,reigenvecs=tmpMat5,leigenvecs=tmpMat6) + call tmpMat1%eigensys(eigenvals=tmpVec1,reigenvecs=tmpMat3,leigenvecs=tmpMat5) + call tmpMat2%eigensys(eigenvals=tmpVec2,reigenvecs=tmpMat4,leigenvecs=tmpMat6) If(present(eVals)) call mqc_eigenvalues_allocate(eVals,'eigenvalues','spin',tmpVec1,tmpVec2) - If(present(rEVecs)) call mqc_integral_allocate(rEVecs,'eigenvectors','space',tmpMat3,tmpMat4) - If(present(lEVecs)) call mqc_integral_allocate(lEVecs,'eigenvectors','space',tmpMat5,tmpMat6) + If(present(rEVecs)) call mqc_integral_allocate(rEVecs,'eigenvectors','spin',tmpMat3,tmpMat4) + If(present(lEVecs)) call mqc_integral_allocate(lEVecs,'eigenvectors','spin',tmpMat5,tmpMat6) endIf elseIf(integralA%type().eq.'general') then nDimAlpha1 = integralA%blockSize('alpha',1) @@ -14227,18 +14561,30 @@ subroutine mqc_integral_put(integral,element,i,j,spinBlockIn) if(present(spinBlockIn)) then call string_change_case(spinBlockIn,'L',spinBlock) if(integral%array_type.eq.'space') then - if(spinBlock.eq.'alpha'.or.spinBlock.eq.'beta') then + if(spinBlock.eq.'alpha') then + integral = integral%changeType('spin') if(integral%hasAlpha()) call integral%alpha%put(element,myI,myJ) - elseIf(spinBlock.ne.'alphabeta'.and.spinBlock.ne.'betaalpha') then - call mqc_error_A('SpinBlock type not valid in MQC_Integral_Put',6,'SpinBlock',SpinBlock) + else if(spinBlock.eq.'beta') then + integral = integral%changeType('spin') + call integral%beta%put(element,myI,myJ) + elseIf(spinBlock.ne.'alphabeta') then + integral = integral%changeType('general') + call integral%alphabeta%put(element,myI,myJ) + elseIf(spinBlock.ne.'betaalpha') then + integral = integral%changeType('general') + call integral%betaalpha%put(element,myI,myJ) endIf elseIf(integral%array_type.eq.'spin') then if(spinBlock.eq.'alpha') then if(integral%hasAlpha()) call integral%alpha%put(element,myI,myJ) elseif(spinblock.eq.'beta') then if(integral%hasBeta()) call integral%beta%put(element,myI,myJ) - elseIf(spinBlock.ne.'alphabeta'.and.spinBlock.ne.'betaalpha') then - call mqc_error_A('SpinBlock type not valid in MQC_Integral_Put',6,'SpinBlock',SpinBlock) + elseIf(spinBlock.ne.'alphabeta') then + integral = integral%changeType('general') + call integral%alphabeta%put(element,myI,myJ) + elseIf(spinBlock.ne.'betaalpha') then + integral = integral%changeType('general') + call integral%betaalpha%put(element,myI,myJ) endIf elseIf(integral%array_type.eq.'general') then if(spinBlock.eq.'alpha') then @@ -14612,84 +14958,328 @@ function mqc_eigenvalues_at(eigenvalues,i,spinBlockIn) result(element) end function mqc_eigenvalues_at ! ! -! PROCEDURE MQC_SCF_Transformation_Matrix +! PROCEDURE MQC_Eigenvalues_Put ! -!> \brief MQC_SCF_Transformation_Matrix is a subroutine that returns the -!> atomic orbital to orthogonal atomic orbital transformation matrix and the -!> number of linearly independent atomic orbital basis functions +!> \brief MQC_Eigenvalues_Put is a function that updates the value of an element +!> of a MQC integral variable !> !> \par Purpose: ! ============= !> !> \verbatim !> -!> MQC_SCF_Transformation_Matrix is a subroutine that returns the atomic orbital -!> to orthogonal atomic orbital transformation matrix and the number of linearly -!> independent atomic orbital basis functions. +!> MQC_Eigenvalues_Put is a function that updates the value of (I)th element of a +!> MQC eigenvalues variable as an MQC scalar. If I is negative, the (N-I+1)th +!> index is selected. Optional argument SpinBlock specifies the spin block of MQC +!> eigenvalues. If I is greater than the dimension of the alpha spin block, the +!> remainder is used as the index of the beta spin block. The following options +!> are available: +!> +!> 1. SpinBlock = 'alpha' uses the alpha spin block. +!> 2. SpinBlock = 'beta' uses the beta spin block. !> !> \endverbatim ! ! Arguments: ! ========== -!> \param[in] Overlap +!> \param[in] eigenvalues !> \verbatim -!> Overlap is type(mqc_scf_integral) -!> The atomic orbital overlap matrix. +!> Integral is Class(MQC_SCF_Eigenvalues) +!> The MQC eigenvalues to return the value of the (I)th !> element. !> \endverbatim !> -!> \param[out] Transform_Matrix +!> \param[in] Element !> \verbatim -!> Transform_Matrix is type(mqc_scf_integral) -!> The AO -> OAO transformation matrix. +!> Element is Type(mqc_scalar) +!> The value to put in the integral. !> \endverbatim !> -!> \param[out] NBasUse +!> \param[in] I !> \verbatim -!> NBasUse is integer(kind=int64),optional -!> The number of linearly independent AO basis functions. +!> I is Integer(kind=int64) +!> The first dimension of the element in MQC eigenvalues +!> If I>0 row count is from first index +!> If I<0 row count is from last index. +!> \endverbatim +!> +!> \param[in,out] SpinBlockIn +!> \verbatim +!> SpinBlock is character(len=64),optional +!> = 'alpha': alpha spin block +!> = 'beta': beta spin block !> \endverbatim ! ! Authors: ! ======== -!> \author L. M. Thompson -!> \date 2018 +!> \author M. M. F. Moraes +!> \date 2025 ! - subroutine mqc_scf_transformation_matrix(overlap,transform_matrix,nBasUse) + subroutine mqc_eigenvalues_put(eigenvalues,element,i,spinBlockIn) ! implicit none - type(mqc_scf_integral),intent(in)::overlap - type(mqc_scf_integral),intent(out)::transform_matrix - integer(kind=int64),optional,intent(out)::nBasUse - type(mqc_matrix)::transform_matrix_tmp,SVecs - type(mqc_vector)::SVals - type(mqc_scalar)::negHalf - integer(kind=int64)::i - real(kind=real64)::thresh=1.0d-4 -! - transform_matrix_tmp = overlap%getBlock('alpha') - call transform_matrix_tmp%diag(SVals,SVecs) - if(present(nBasUse)) then - nBasUse = 0 - do i = 1, SVals%size() - if(MQC_Scalar_Get_Intrinsic_Real(SVals%at(i)).gt.thresh) nBasUse = nBasUse+1 - endDo + class(mqc_scf_eigenvalues)::eigenvalues + class(*),intent(in)::i + integer(kind=int64)::myI + type(mqc_scalar)::element + character(len=*),optional,intent(in)::spinBlockIn + character(len=64)::spinBlock + real(kind=real64)::zero=0.0d0 + + select type (i) + type is (integer) + myI = i + type is (mqc_scalar) + myI = i + class default + call mqc_error('index type not defined in MQC_twoERIs_put') + end select + + if(myI.gt.eigenvalues%blockSize('alpha')+eigenvalues%blockSize('beta')) & + call mqc_error('Requested I dimension out of bounds in mqc_eigenvalues_put') + + if(present(spinBlockIn)) then + call string_change_case(spinBlockIn,'L',spinBlock) + if(eigenvalues%array_type.eq.'space') then + if(spinBlock.eq.'alpha'.or.spinBlock.eq.'beta') then + if(eigenvalues%hasAlpha()) call eigenvalues%alpha%put(element,myI) + else + call mqc_error_A('SpinBlock type not valid in MQC_Integral_Put',6,'SpinBlock',SpinBlock) + endIf + elseIf(eigenvalues%array_type.eq.'spin') then + if(spinBlock.eq.'alpha') then + if(eigenvalues%hasAlpha()) call eigenvalues%alpha%put(element,myI) + elseif(spinblock.eq.'beta') then + if(eigenvalues%hasBeta()) call eigenvalues%beta%put(element,myI) + else + call mqc_error_A('SpinBlock type not valid in MQC_Integral_Put',6,'SpinBlock',SpinBlock) + endIf + else + call mqc_error_A('Array type not valid in MQC_Eigenvalues_Put',6,'Array_Type',eigenvalues%array_type) + endIf + else + if(eigenvalues%array_type.eq.'space') then + if(myI.le.eigenvalues%blockSize('alpha')) then + if(eigenvalues%hasAlpha()) call eigenvalues%alpha%put(element,myI) + elseIf(myI.gt.eigenvalues%blockSize('alpha')) then + myI = myI - eigenvalues%blockSize('alpha') + if(eigenvalues%hasAlpha()) call eigenvalues%alpha%put(element,myI) + endIf + elseIf(eigenvalues%array_type.eq.'spin') then + if(myI.le.eigenvalues%blockSize('alpha')) then + if(eigenvalues%hasAlpha()) call eigenvalues%alpha%put(element,myI) + elseIf(myI.gt.eigenvalues%blockSize('alpha')) then + myI = myI - eigenvalues%blockSize('alpha') + if(eigenvalues%hasBeta()) call eigenvalues%beta%put(element,myI) + endIf + else + call mqc_error_A('Array type not valid in MQC_Eigenvalues_Put',6,'Array_Type',eigenvalues%array_type) + endIf endIf -#ifdef PGI - negHalf = -(dfloat(1)/dfloat(2)) -#else - negHalf = -(float(1)/float(2)) -#endif - call mqc_vector_power(SVals,negHalf) - call MQC_Matrix_DiagMatrix_Put_Vector(SVals,transform_matrix_tmp) - transform_matrix_tmp = matmul(SVecs,transform_matrix_tmp) - call mqc_integral_allocate(transform_matrix,'transformation matrix','space',transform_matrix_tmp) -! - return - end subroutine mqc_scf_transformation_matrix + end subroutine mqc_eigenvalues_put + ! +! PROCEDURE MQC_Eigenvalues_Push ! -! PROCEDURE mqc_scf_integral_diagonal_elements +!> \brief MQC_Eigenvalues_Put is a function that adds a value to the end +!> of a MQC integral variable +!> +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_Eigenvalues_Push is a function that adds a value at the end of a MQC vector +!> of a MQC eigenvalues variable as an MQC scalar. Optional argument SpinBlock +!> specifies the spin block of MQC eigenvalues. The following options are available: +!> +!> 1. SpinBlock = 'alpha' the value is added to the alpha spin block. +!> 2. SpinBlock = 'beta' the value is added to the beta spin block. +!> 3. SpinBlock = 'both' the value is added to the both alpha and beta spin blocks. +!> +!> If MQC_Eigenvalues Type is 'space' +!> SpinBlock equal to 'alpha' or 'beta' will convert it to 'spin' +!> SpinBlock equal to 'both' or absent will keep the 'space' type +!> eseIf MQC_Eigenvalues Type is 'spin' +!> SpinBlock equal to 'alpha' or absent will append to the alpha block +!> SpinBlock equal to 'beta' will append to the beta block +!> SpinBlock equal to 'both' will append to the alpha and beta blocks +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] eigenvalues +!> \verbatim +!> Integral is Class(MQC_SCF_Eigenvalues) +!> The MQC eigenvalues to return the value of the (I)th +!> element. +!> \endverbatim +!> +!> \param[in] Element +!> \verbatim +!> Element is Class(*) +!> The value to put in the integral. +!> \endverbatim +!> +!> \param[in,out] SpinBlockIn +!> \verbatim +!> SpinBlock is character(len=64),optional +!> = 'alpha': alpha spin block +!> = 'beta': beta spin block +!> = 'both': alpha and beta spin blocks +!> \endverbatim +! +! Authors: +! ======== +!> \author M. M. F. Moraes +!> \date 2025 +! + subroutine mqc_eigenvalues_push(eigenvalues,element,spinBlockIn) +! + implicit none + class(mqc_scf_eigenvalues)::eigenvalues + type(mqc_scalar)::my_element + class(*),intent(in)::element + character(len=*),optional,intent(in)::spinBlockIn + character(len=64)::spinBlock + real(kind=real64)::zero=0.0d0 + + select type (element) + type is (integer) + my_element = element + type is (real) + my_element = element + type is (complex) + my_element = element + type is (mqc_scalar) + my_element = element + class default + call mqc_error('Scalar type not defined in mqc_vector_push',6) + end select + + + + if(present(spinBlockIn)) then + call string_change_case(spinBlockIn,'L',spinBlock) + if(eigenvalues%array_type.eq.'space') then + if(spinBlock.eq.'alpha') then + eigenvalues%beta = eigenvalues%alpha + eigenvalues%array_type = 'spin' + call eigenvalues%alpha%push(my_element) + else if (spinBlock.eq.'beta') then + eigenvalues%beta = eigenvalues%alpha + eigenvalues%array_type = 'spin' + call eigenvalues%beta%push(my_element) + else if (spinBlock.eq.'both') then + call eigenvalues%alpha%push(my_element) + else + call mqc_error_A('SpinBlock type not valid in MQC_Integral_Put',6,'SpinBlock',SpinBlock) + endIf + elseIf(eigenvalues%array_type.eq.'spin') then + if(spinBlock.eq.'alpha') then + call eigenvalues%alpha%push(my_element) + elseif(spinblock.eq.'beta') then + call eigenvalues%beta%push(my_element) + else if (spinBlock.eq.'both') then + call eigenvalues%alpha%push(my_element) + call eigenvalues%beta%push(my_element) + else + call mqc_error_A('SpinBlock type not valid in MQC_Integral_Put',6,'SpinBlock',SpinBlock) + endIf + else + call mqc_error_A('Array type not valid in MQC_Eigenvalues_Put',6,'Array_Type',eigenvalues%array_type) + endIf + else + if(eigenvalues%array_type.eq.'space'.or.eigenvalues%array_type.eq.'spin') then + call eigenvalues%alpha%push(my_element) + else + call mqc_error_A('Array type not valid in MQC_Eigenvalues_Put',6,'Array_Type',eigenvalues%array_type) + endIf + endIf + end subroutine mqc_eigenvalues_push + +! +! +! +! PROCEDURE MQC_SCF_Transformation_Matrix +! +!> \brief MQC_SCF_Transformation_Matrix is a subroutine that returns the +!> atomic orbital to orthogonal atomic orbital transformation matrix and the +!> number of linearly independent atomic orbital basis functions +!> +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_SCF_Transformation_Matrix is a subroutine that returns the atomic orbital +!> to orthogonal atomic orbital transformation matrix and the number of linearly +!> independent atomic orbital basis functions. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] Overlap +!> \verbatim +!> Overlap is type(mqc_scf_integral) +!> The atomic orbital overlap matrix. +!> element. +!> \endverbatim +!> +!> \param[out] Transform_Matrix +!> \verbatim +!> Transform_Matrix is type(mqc_scf_integral) +!> The AO -> OAO transformation matrix. +!> \endverbatim +!> +!> \param[out] NBasUse +!> \verbatim +!> NBasUse is integer(kind=int64),optional +!> The number of linearly independent AO basis functions. +!> \endverbatim +! +! Authors: +! ======== +!> \author L. M. Thompson +!> \date 2018 +! + subroutine mqc_scf_transformation_matrix(overlap,transform_matrix,nBasUse) +! + implicit none + type(mqc_scf_integral),intent(in)::overlap + type(mqc_scf_integral),intent(out)::transform_matrix + integer(kind=int64),optional,intent(out)::nBasUse + type(mqc_matrix)::transform_matrix_tmp,SVecs + type(mqc_vector)::SVals + type(mqc_scalar)::negHalf + integer(kind=int64)::i + real(kind=real64)::thresh=1.0d-4 +! + transform_matrix_tmp = overlap%getBlock('alpha') + call transform_matrix_tmp%diag(SVals,SVecs) + if(present(nBasUse)) then + nBasUse = 0 + do i = 1, SVals%size() + if(MQC_Scalar_Get_Intrinsic_Real(SVals%at(i)).gt.thresh) nBasUse = nBasUse+1 + endDo + endIf +#ifdef PGI + negHalf = -(dfloat(1)/dfloat(2)) +#else + negHalf = -(float(1)/float(2)) +#endif + call mqc_vector_power(SVals,negHalf) + call MQC_Matrix_DiagMatrix_Put_Vector(SVals,transform_matrix_tmp) + transform_matrix_tmp = matmul(SVecs,transform_matrix_tmp) + call mqc_integral_allocate(transform_matrix,'transformation matrix','space',transform_matrix_tmp) +! + return + end subroutine mqc_scf_transformation_matrix +! +! +! PROCEDURE mqc_scf_integral_diagonal_elements ! !> \brief MQC_SCF_Integral_Diagonal_Elements is a function that returns the !> diagonal of an MQC SCF integral as an MQC SCF eigenvalue object @@ -15054,6 +15644,47 @@ function mqc_eigenvalues_maxloc(eigenvalues,label) result(output) end function mqc_eigenvalues_maxloc ! ! +! PROCEDURE MQC_Eigenvalues_Abs +! +!> \brief MQC_Eigenvalues_Abs is a function that returns the absolute +!> values of all elements of an MQC_Eigenvalues variable +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_Eigenvalues_Abs is a function that returns the absolute values of +!> all elements of an MQC_Eigenvalues variable +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in,out] Eigenvalues +!> \verbatim +!> Eigenvalues is Class(MQC_Eigenvalues) +!> The name of the MQC_Eigenvalues variable. +!> \endverbatim +! +! Authors: +! ======== +!> \author L. M. Thompson +!> \date 2025 +! + function mqc_eigenvalues_abs(eigenvalues) result(output) +! + implicit none + class(mqc_scf_eigenvalues),intent(in)::eigenvalues + type(mqc_scf_eigenvalues)::output +! + output = eigenvalues + if(output%hasalpha()) output%alpha = abs(output%alpha) + if(output%hasbeta()) output%beta = abs(output%beta) +! + end function mqc_eigenvalues_abs +! +! !===================================================================== ! ! POST-SCF ROUTINES @@ -16064,7 +16695,7 @@ function mqc_detString_to_occArray(string) result(arrayOut) type(mqc_vector)::tmpArray do i = 0, string%nBits-1 - if(btest(string,i)) call tmpArray%push(i) + if(btest(string,i)) call tmpArray%push(i+1) endDo allocate(arrayOut(size(tmpArray))) arrayOut = tmpArray @@ -16369,7 +17000,7 @@ Function Slater_Condon(IOut,IPrint,Alpha_String_1,Beta_String_1, & ISgn = ISgn-2 If(Det(1).eq.Det(2).and.Det(3).eq.Det(4)) then -! Write(IOut,*) 'Entering case 1' + ! Write(IOut,*) 'Entering case 1' ! If(Spin(1).gt.Spin(2)) ISgn = ISgn+1 If(Spin(3).gt.Spin(4)) ISgn = ISgn+1 @@ -21337,8 +21968,7 @@ End Function S2_Mat_Elem !> Slater Condon rules. The spatial molecular orbital basis !> overlap (i.e. do not multiply by spin overlap which !> always leads to zero alpha-beta blocks) must be provided -!> in MO_Core_Ham argument if this option is used) must be provided -!> in MO_Core_Ham argument if this option is used.. +!> in MO_Core_Ham argument if this option is used. !> \endverbatim !> ! Authors: @@ -21655,37 +22285,54 @@ Subroutine MQC_Build_CI_Hamiltonian(IOut,IPrint,CI_Hamiltonian,Determinants,Dets End Subroutine MQC_Build_CI_Hamiltonian ! -! -! !===================================================================== ! -! PROCEDURE MQC_BUILD_CI_OPERATOR +! PROCEDURE MQC_BUILD_INDEX_DETERMINANT_STRINGS ! -!> \brief MQC_BUILD_CI_OPERATOR is a subroutine that builds an -!> operator matrix given MO integrals and determinant strings +!> \brief MQC_BUILD_INDEX_DETERMINANT_STRINGS is a subroutine that +!> a set of alpha and beta determinant strings, possibly different for +!> bra and ket, along with infomration about substitution levels etc. +!> It returns the alpha and beta determinant string for each index in the +!> operator matrix. !> !> \par Purpose: ! ============= !> !> \verbatim !> -!> MQC_BUILD_CI_OPERATOR is a subroutine that builds an operator matrix -!> given MO integrals and determinant strings. It can be used for any one, -!> two, or one and two particle operator. The routine can also be used to -!> build the S^2 operator matrix by passing the orbital overlap matrix and -!> setting the appropriate input flag. In contrast to the similar routine -!> MQC_BUILD_CI_HAMILTONIAN, MQC_BUILD_CI_OPERATOR can accept multiple -!> determinant strings corresponsing to substring blocks. Due to the -!> requirement to compute the determinant strings corresponding to a given -!> index, which has no analytic form and must be determined by searching -!> through all strings of a given substitution level, the routine is less -!> computationally efficient for a given orbital string length than -!> MQC_BUILD_CI_HAMILTONIAN. +!> MQC_BUILD_INDEX_DETERMINANT_STRINGS is a subroutine that a set of alpha +!> and beta determinant strings, possibly different for bra and ket, along +!> with infomration about substitution levels etc. It returns the alpha and +!> beta determinant string for each index in the operator matrix. !> !> \endverbatim ! ! Arguments: ! ========== +!> \param[out] L_Alpha_Strings +!> \verbatim +!> L_Alpha_Strings is type(MQC_bits),dimension(:),allocatable +!> The left alpha string array. +!> \endverbatim +!> +!> \param[out] R_Alpha_Strings +!> \verbatim +!> R_Alpha_Strings is type(MQC_bits),dimension(:),allocatable +!> The right alpha string array. +!> \endverbatim +!> +!> \param[out] L_Beta_Strings +!> \verbatim +!> L_Beta_Strings is type(MQC_bits),dimension(:),allocatable +!> The left beta string array. +!> \endverbatim +!> +!> \param[out] R_Beta_Strings +!> \verbatim +!> R_Beta_Strings is type(MQC_bits),dimension(:),allocatable +!> The right alpha string array. +!> \endverbatim +!> !> \param[in] IOut !> \verbatim !> IOut is Integer(kind=int64) @@ -21698,12 +22345,6 @@ End Subroutine MQC_Build_CI_Hamiltonian !> The print level for the subroutine. !> \endverbatim !> -!> \param[out] Operator_Matrix -!> \verbatim -!> Operator_Matrix is Type(MQC_Matrix) -!> The CI operator matrix returned. -!> \endverbatim -!> !> \param[in] Determinants !> \verbatim !> Determinants is Type(MQC_Determinant),Dimension(:) @@ -21724,18 +22365,6 @@ End Subroutine MQC_Build_CI_Hamiltonian !> should be provided for each substring set. !> \endverbatim !> -!> \param[in] MO_Core_Ham -!> \verbatim -!> MO_Core_Ham is Type(MQC_SCF_Integral),Optional -!> The MO basis one-particle integrals. -!> \endverbatim -!> -!> \param[in] MO_ERIs -!> \verbatim -!> MO_ERIs is Type(MQC_TwoERIs),Optional -!> The MO basis two-particle integrals. -!> \endverbatim -!> !> \param[in] SubsAIn !> \verbatim !> SubsA is Type(MQC_Vector),Dimension(:),Optional @@ -21782,59 +22411,307 @@ End Subroutine MQC_Build_CI_Hamiltonian !> substitutions are equal to row basis substitutions. !> \endverbatim !> -!> \param[in] doS2 +!> \param[out] SymmFlagOut !> \verbatim -!> DoS2 is Logical,Optional -!> Logical flag to build the S^2 matrix instead of using -!> Slater Condon rules. The spatial molecular orbital basis -!> overlap (i.e. do not multiply by spin overlap which -!> always leads to zero alpha-beta blocks) must be provided -!> in MO_Core_Ham argument if this option is used) must be provided -!> in MO_Core_Ham argument if this option is used.. +!> SymmFlagOut is Logical,Optional +!> Logical flag to flag whether the left and right strings +!> correspond to a symmetric determinant basis. !> \endverbatim !> ! Authors: ! ======== !> \author L. M. Thompson -!> \date 2017,2021,2024 +!> \date 2025 ! - Subroutine MQC_Build_CI_Operator(IOut,IPrint,Operator_Matrix,Determinants,Dets2,& - MO_Core_Ham,MO_ERIs,SubsAIn,SubsBIn,Subs2AIn,Subs2BIn,doS2) + Subroutine MQC_Build_Index_Determinant_Strings(L_Alpha_Strings,R_Alpha_Strings,& + L_Beta_Strings,R_Beta_Strings,IOut,IPrint,Determinants,Dets2,SubsAIn,SubsBIn,Subs2AIn,& + Subs2BIn,SymmFlagOut) ! Implicit None + Type(MQC_Bits),Dimension(:),Allocatable,Intent(Out)::L_Alpha_Strings,R_Alpha_Strings,& + L_Beta_Strings,R_Beta_Strings Integer(kind=int64),Intent(In)::IOut,IPrint - Type(MQC_TwoERIs),Optional,Intent(In)::MO_ERIs - Type(MQC_SCF_Integral),Optional,Intent(In)::MO_Core_Ham Type(MQC_Determinant),Dimension(:),Intent(In)::Determinants Type(MQC_Determinant),Dimension(:),Optional,Intent(In)::Dets2 Type(MQC_Vector),Dimension(:),Optional,Intent(In)::SubsAIn,SubsBIn,Subs2AIn,Subs2BIn - Logical,Optional,Intent(In)::doS2 - Type(MQC_Matrix),Intent(Out)::Operator_Matrix - Integer(kind=int64)::I,J,L_Index,R_Index,NDets1T,NDets2T,R_Index_End,LDetSize,RDetSize + Logical,Optional,Intent(Out)::SymmFlagOut + Integer(kind=int64)::I,J,L_Index,R_Index,NDets1T,NDets2T Integer(kind=int64),Dimension(:),Allocatable::SubsA,SubsB,Subs2A,Subs2B,NDets1,NDets2,& NDets1A,NDets1B,NDets2A,NDets2B - Type(MQC_Bits)::Alpha_String_1,Alpha_String_2,Beta_String_1,Beta_String_2,alphaStrTemp,& - betaStrTemp Type(MQC_Vector),Dimension(:),Allocatable::lSubsT,lSubsA,lSubsB,rSubsT,rSubsA,rSubsB Type(MQC_Determinant),Dimension(:),Allocatable::Determinants2 - Logical::SymmFlag,S2Flag - Character(len=12)::SymmStr + Logical::SymmFlag ! - If(present(doS2)) then - S2Flag = doS2 - else - S2Flag = .false. - EndIf -! - If(.not.present(MO_Core_Ham).and..not.present(MO_ERIs)) & - call mqc_error('No integrals given to MQC_Build_CI_Operator') ! If(present(Dets2)) then Determinants2 = Dets2 - SymmFlag = .False. + SymmFlag = .True. Else Determinants2 = Determinants + SymmFlag = .False. + EndIf +! +! Determine different substitutions requested from input, and number of strings +! of different types. +! + Call MQC_Get_Substring_CI_Parameters(LSubsT,RSubsT,LSubsA,RSubsA,LSubsB,RSubsB,& + NDets1,NDets2,NDets1A,NDets2A,NDets1B,NDets2B,NDets1T,NDets2T,SymmFlag,Determinants,& + Dets2,SubsAIn,SubsBIn,Subs2AIn,Subs2BIn) +! +! Allocate and fill determinant string arrays +! + Allocate(L_Alpha_Strings(NDets1T),L_Beta_Strings(NDets1T),R_Alpha_Strings(NDets2T),& + R_Beta_Strings(NDets2T)) + + Do L_Index = 1, NDets1T + call mqc_get_string_from_substring(iout,iprint,l_index,determinants,& + nDets1,lsubsA,lSubsB,lSubsT,l_alpha_strings(l_index),l_beta_strings(l_index)) + EndDo + + if(symmFlag) then + r_alpha_strings = l_alpha_strings + r_beta_strings = l_beta_strings + else + Do R_Index = 1, NDets2T + call mqc_get_string_from_substring(iout,iprint,r_index,determinants2,& + nDets2,rsubsA,rSubsB,rSubsT,r_alpha_strings(r_index),r_beta_strings(r_index)) + EndDo + endIf + + if(present(symmFlagOut)) symmFlagOut = symmFlag + + End Subroutine MQC_Build_Index_Determinant_Strings +! +! +! +!===================================================================== +! +! PROCEDURE MQC_GET_SUBSTRING_CI_PARAMETERS +! +!> \brief MQC_GET_SUBSTRING_CI_PARAMETERS is a subroutine that +!> returns information about substitution levels and number of +!> determinants based on input arguments +!> +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_GET_SUBSTRING_CI_PARAMETERS is a subroutine that returns information +!> about substitution levels and number of determinants based on input +!> arguments. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[out] LSubsT +!> \verbatim +!> LSubsT is Type(MQC_Vector),Dimension(:),Allocatable +!> Array of total substitutions allowed in each orbital subset +!> of the left determinant set. +!> \endverbatim +!> +!> \param[out] RSubsT +!> \verbatim +!> RSubsT is Type(MQC_Vector),Dimension(:),Allocatable +!> Array of total substitutions allowed in each orbital subset +!> of the right determinant set. +!> \endverbatim +!> +!> \param[out] LSubsA +!> \verbatim +!> LSubsA is Type(MQC_Vector),Dimension(:),Allocatable +!> Array of alpha substitutions allowed in each orbital subset +!> of the left determinant set. +!> \endverbatim +!> +!> \param[out] RSubsA +!> \verbatim +!> RSubsA is Type(MQC_Vector),Dimension(:),Allocatable +!> Array of alpha substitutions allowed in each orbital subset +!> of the right determinant set. +!> \endverbatim +!> +!> \param[out] LSubsB +!> \verbatim +!> LSubsB is Type(MQC_Vector),Dimension(:),Allocatable +!> Array of beta substitutions allowed in each orbital subset +!> of the left determinant set. +!> \endverbatim +!> +!> \param[out] RSubsB +!> \verbatim +!> RSubsB is Type(MQC_Vector),Dimension(:),Allocatable +!> Array of beta substitutions allowed in each orbital subset +!> of the right determinant set. +!> \endverbatim +!> +!> \param[in] NDets1 +!> \verbatim +!> NDets1 is Integer(kind=int64),Dimension(:) +!> The number of total strings in each orbital subset +!> given the permitted substitutions for the left +!> determinant set. +!> \endverbatim +!> +!> \param[in] NDets2 +!> \verbatim +!> NDets2 is Integer(kind=int64),Dimension(:) +!> The number of total strings in each orbital subset +!> given the permitted substitutions for the right +!> determinant set. +!> \endverbatim +!> +!> \param[in] NDets1A +!> \verbatim +!> NDets1A is Integer(kind=int64),Dimension(:) +!> The number of alpha strings in each orbital subset +!> given the permitted substitutions for the left +!> determinant set. +!> \endverbatim +!> +!> \param[in] NDets2A +!> \verbatim +!> NDets2A is Integer(kind=int64),Dimension(:) +!> The number of alpha strings in each orbital subset +!> given the permitted substitutions for the right +!> determinant set. +!> \endverbatim +!> +!> \param[in] NDets1B +!> \verbatim +!> NDets1B is Integer(kind=int64),Dimension(:) +!> The number of beta strings in each orbital subset +!> given the permitted substitutions for the left +!> determinant set. +!> \endverbatim +!> +!> \param[in] NDets2B +!> \verbatim +!> NDets2B is Integer(kind=int64),Dimension(:) +!> The number of beta strings in each orbital subset +!> given the permitted substitutions for the right +!> determinant set. +!> \endverbatim +!> +!> \param[out] SymmFlag +!> \verbatim +!> SymmFlagOut is Logical +!> Logical flag to flag whether the left and right strings +!> correspond to a symmetric determinant basis. +!> \endverbatim +!> +!> \param[in] NDets1T +!> \verbatim +!> NDets1A is Integer(kind=int64) +!> The total number of determinants given the permitted +!> substitutions for the left determinant set. +!> \endverbatim +!> +!> \param[in] NDets2T +!> \verbatim +!> NDets2T is Integer(kind=int64),Dimension(:) +!> The total number of determinants given the permitted +!> substitutions for the right determinant set. +!> \endverbatim +!> +!> \param[in] Determinants +!> \verbatim +!> Determinants is Type(MQC_Determinant),Dimension(:) +!> The binary string occupation number vectors. If +!> opional argument Dets2 is provided Determinants +!> variable will contain the row basis strings. One +!> determinant object should be provided for each substring +!> set. +!> \endverbatim +!> +!> \param[in] Dets2 +!> \verbatim +!> Dets2 is Type(MQC_Determinant),Dimension(:),Optional +!> If provided, gives the binary string occupation number +!> vectors for the column basis, where Determinants variable +!> will continue to provide the binary string occupation +!> number vectors for the row basis. One determinant object +!> should be provided for each substring set. +!> \endverbatim +!> +!> \param[in] SubsAIn +!> \verbatim +!> SubsA is Type(MQC_Vector),Dimension(:),Optional +!> If provided, gives the permitted substitutions to include +!> in the Hamiltonian. One mqc_vector object should be +!> provided for each of the determinant objects given. +!> If not present, all substitutions are included. If SubsB +!> not present, SubsA containes the the total (alpha+beta) +!> substitutions that are permitted. If both SubsA and SubsB +!> provided, SubsA contains the alpha substitutions that are +!> permitted. If Subs2A are provided, Subs1A gives substitutions +!> for the row basis only. +!> \endverbatim +!> +!> \param[in] SubsBIn +!> \verbatim +!> SubsB is Type(MQC_Vector),Dimension(:),Optional +!> If provided, gives the permitted beta substitutions to +!> include in the Hamiltonian. One mqc_vector object should be +!> provided for each of the determinant objects given. +!> \endverbatim +!> +!> \param[in] Subs2AIn +!> \verbatim +!> Subs2A is Type(MQC_Vector),Dimension(:),Optional +!> If provided, gives the permitted substitutions to include +!> in the Hamiltonian for the column basis. One mqc_vector +!> object should be provided for each of the determinant +!> objects given. If not present, column basis substitutions +!> are equal to row basis substitutions provided by variable +!> SubsA and/or SubsB. If Subs2A provided but not Subs2B, +!> Subs2A contains the total (alpha+beta) substitutions that +!> are permitted. If both Subs2A and Subs2B provided, Subs2A +!> containes the alpha substitutions that are permitted. +!> \endverbatim +!> +!> \param[in] Subs2BIn +!> \verbatim +!> Subs2B is Type(MQC_Vector),Dimension(:),Optional +!> If provided, gives the permitted beta substitutions to +!> include in the Hamiltonian for the column basis. One +!> mqc_vector object should be provided for each of the +!> determinant objects given. If not present, column basis +!> substitutions are equal to row basis substitutions. +!> \endverbatim +!> +! Authors: +! ======== +!> \author L. M. Thompson +!> \date 2025 +! + Subroutine MQC_Get_Substring_CI_Parameters(LSubsT,RSubsT,LSubsA,RSubsA,LSubsB,RSubsB,& + NDets1,NDets2,NDets1A,NDets2A,NDets1B,NDets2B,NDets1T,NDets2T,SymmFlag,Determinants,& + Dets2,SubsAIn,SubsBIn,Subs2AIn,Subs2BIn) +! + Implicit None + Type(MQC_Vector),Dimension(:),Allocatable,Intent(Out)::lSubsT,lSubsA,lSubsB,rSubsT,rSubsA,& + rSubsB + Integer(kind=int64),Dimension(:),Allocatable,Intent(Out)::NDets1,NDets2,NDets1A,NDets1B,& + NDets2A,NDets2B + Integer(kind=int64),Intent(Out)::NDets1T,NDets2T + Logical,Intent(Out)::SymmFlag + Type(MQC_Determinant),Dimension(:),Intent(In)::Determinants + Type(MQC_Determinant),Dimension(:),Optional,Intent(In)::Dets2 + Type(MQC_Vector),Dimension(:),Optional,Intent(In)::SubsAIn,SubsBIn,Subs2AIn,Subs2BIn + Integer(kind=int64)::I,LDetSize,RDetSize + Integer(kind=int64),Dimension(:),Allocatable::SubsA,SubsB,Subs2A,Subs2B + Type(MQC_Determinant),Dimension(:),Allocatable::Determinants2 +! +! + If(present(Dets2)) then + Determinants2 = Dets2 SymmFlag = .True. + Else + Determinants2 = Determinants + SymmFlag = .False. EndIf LDetSize = Size(Determinants) @@ -21867,8 +22744,16 @@ Subroutine MQC_Build_CI_Operator(IOut,IPrint,Operator_Matrix,Determinants,Dets2, NDets2T = 1 do i = 1, RDetSize if(.not.present(subs2AIn).and..not.present(subs2BIn)) then - subs2A = subsA - subs2B = subsB + if(present(subsAIn)) then + subs2A = subsAIn(i) + else + subs2A = [integer::] + endIf + if(present(subsBIn)) then + subs2B = subsBIn(i) + else + subs2B = [integer::] + endIf elseIf(.not.present(subs2BIn)) then subs2A = subs2AIn(i) subs2B = [integer::] @@ -21887,6 +22772,309 @@ Subroutine MQC_Build_CI_Operator(IOut,IPrint,Operator_Matrix,Determinants,Dets2, NDets2(i),NDets2A(i),NDets2B(i)) NDets2T = NDets2T*NDets2(i) endDo + + End Subroutine MQC_Get_Substring_CI_Parameters +! +!===================================================================== +! +! PROCEDURE MQC_GET_STRING_FROM_SUBSTRING +! +!> \brief MQC_GET_STRING_FROM_SUBSTRING is a subroutine that returns +!> an alpha and beta determinant string, given substrings and substitution +!> level information +!> +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_GET_STRING_FROM_SUBSTRING is a subroutine that returns an alpha and +!> beta determinant string, given substrings and substitution level +!> information. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] IOut +!> \verbatim +!> IOut is Integer(kind=int64) +!> The FORTRAN file to print to. +!> \endverbatim +!> +!> \param[in] IPrint +!> \verbatim +!> IPrint is Integer(kind=int64) +!> The print level for the subroutine. +!> \endverbatim +!> +!> \param[in] Ind +!> \verbatim +!> Ind is Integer(kind=int64) +!> The Index in the determinant basis. +!> \endverbatim +!> +!> \param[in] Determinants +!> \verbatim +!> Determinants is Type(MQC_Determinant),Dimension(:) +!> The binary string occupation number vectors. +!> \endverbatim +!> +!> \param[in] NDetsArr +!> \verbatim +!> NDetsArr is Integer(kind=int64),Dimension(:) +!> The number of total strings in each orbital subset +!> given the permitted substitutions. +!> \endverbatim +!> +!> \param[in] ASubsArr +!> \verbatim +!> ASubsArr is Type(MQC_Vector),Dimension(:) +!> The array of the alpha substitutions allowed in each +!> orbital subset. +!> \endverbatim +!> +!> \param[in] BSubsArr +!> \verbatim +!> BSubsArr is Type(MQC_Vector),Dimension(:) +!> The array of the beta substitutions allowed in each +!> orbital subset. +!> \endverbatim +!> +!> \param[in] TSubsArr +!> \verbatim +!> TSubsArr is Type(MQC_Vector),Dimension(:) +!> The array of the total substitutions allowed in each +!> orbital subset. +!> \endverbatim +!> +!> \param[out] Alpha_String +!> \verbatim +!> Alpha_String is Type(MQC_Bits) +!> The composite alpha sting at the specified index. +!> \endverbatim +!> +!> \param[out] Beta_String +!> \verbatim +!> Beta_String is Type(MQC_Bits) +!> The composite beta string at the specified index. +!> \endverbatim +!> +! Authors: +! ======== +!> \author L. M. Thompson +!> \date 2025 +! + Subroutine MQC_Get_String_From_Substring(iout,iprint,ind,determinants,& + NDetsArr,ASubsArr,BSubsArr,TSubsArr,alpha_string,beta_string) +! + Implicit None + Integer(kind=int64),Intent(In)::IOut,IPrint,Ind + Type(MQC_Determinant),Dimension(:),Intent(In)::Determinants + Integer(kind=int64),Dimension(:),Intent(In)::NDetsArr + Type(MQC_Vector),Dimension(:),Intent(In)::TSubsArr,ASubsArr,BSubsArr + Type(MQC_Bits),Intent(Out)::Alpha_String,Beta_String + Integer(kind=int64)::I,J + Type(MQC_Bits)::alphaStrTemp,betaStrTemp +! + alpha_string = mqc_bits(0) + beta_string = mqc_bits(0) + Do I = 1, size(Determinants) + J = mod((Ind-1)/product(NDetsArr(1:i-1)),NDetsArr(i))+1 + call mqc_strings_at_index(iout,iprint,j,Determinants(i),ASubsArr(i),& + BSubsArr(i),TSubsArr(i),AlphaStrTemp,BetaStrTemp) + alpha_string = MQC_concatenate_bits(alphaStrTemp,alpha_string) + beta_string = MQC_concatenate_bits(betaStrTemp,beta_string) + EndDo + + End Subroutine MQC_Get_String_From_Substring +! +!===================================================================== +! +! PROCEDURE MQC_BUILD_CI_OPERATOR +! +!> \brief MQC_BUILD_CI_OPERATOR is a subroutine that builds an +!> operator matrix given MO integrals and determinant strings +!> +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_BUILD_CI_OPERATOR is a subroutine that builds an operator matrix +!> given MO integrals and determinant strings. It can be used for any one, +!> two, or one and two particle operator. The routine can also be used to +!> build the S^2 operator matrix by passing the orbital overlap matrix and +!> setting the appropriate input flag. In contrast to the similar routine +!> MQC_BUILD_CI_HAMILTONIAN, MQC_BUILD_CI_OPERATOR can accept multiple +!> determinant strings corresponsing to substring blocks. Due to the +!> requirement to compute the determinant strings corresponding to a given +!> index, which has no analytic form and must be determined by searching +!> through all strings of a given substitution level, the routine is less +!> computationally efficient for a given orbital string length than +!> MQC_BUILD_CI_HAMILTONIAN. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] IOut +!> \verbatim +!> IOut is Integer(kind=int64) +!> The FORTRAN file to print to. +!> \endverbatim +!> +!> \param[in] IPrint +!> \verbatim +!> IPrint is Integer(kind=int64) +!> The print level for the subroutine. +!> \endverbatim +!> +!> \param[out] Operator_Matrix +!> \verbatim +!> Operator_Matrix is Type(MQC_Matrix) +!> The CI operator matrix returned. +!> \endverbatim +!> +!> \param[in] Determinants +!> \verbatim +!> Determinants is Type(MQC_Determinant),Dimension(:) +!> The binary string occupation number vectors. If +!> opional argument Dets2 is provided Determinants +!> variable will contain the row basis strings. One +!> determinant object should be provided for each substring +!> set. +!> \endverbatim +!> +!> \param[in] Dets2 +!> \verbatim +!> Dets2 is Type(MQC_Determinant),Dimension(:),Optional +!> If provided, gives the binary string occupation number +!> vectors for the column basis, where Determinants variable +!> will continue to provide the binary string occupation +!> number vectors for the row basis. One determinant object +!> should be provided for each substring set. +!> \endverbatim +!> +!> \param[in] MO_Core_Ham +!> \verbatim +!> MO_Core_Ham is Type(MQC_SCF_Integral),Optional +!> The MO basis one-particle integrals. +!> \endverbatim +!> +!> \param[in] MO_ERIs +!> \verbatim +!> MO_ERIs is Type(MQC_TwoERIs),Optional +!> The MO basis two-particle integrals. +!> \endverbatim +!> +!> \param[in] SubsAIn +!> \verbatim +!> SubsA is Type(MQC_Vector),Dimension(:),Optional +!> If provided, gives the permitted substitutions to include +!> in the Hamiltonian. One mqc_vector object should be +!> provided for each of the determinant objects given. +!> If not present, all substitutions are included. If SubsB +!> not present, SubsA containes the the total (alpha+beta) +!> substitutions that are permitted. If both SubsA and SubsB +!> provided, SubsA contains the alpha substitutions that are +!> permitted. If Subs2A are provided, Subs1A gives substitutions +!> for the row basis only. +!> \endverbatim +!> +!> \param[in] SubsBIn +!> \verbatim +!> SubsB is Type(MQC_Vector),Dimension(:),Optional +!> If provided, gives the permitted beta substitutions to +!> include in the Hamiltonian. One mqc_vector object should be +!> provided for each of the determinant objects given. +!> \endverbatim +!> +!> \param[in] Subs2AIn +!> \verbatim +!> Subs2A is Type(MQC_Vector),Dimension(:),Optional +!> If provided, gives the permitted substitutions to include +!> in the Hamiltonian for the column basis. One mqc_vector +!> object should be provided for each of the determinant +!> objects given. If not present, column basis substitutions +!> are equal to row basis substitutions provided by variable +!> SubsA and/or SubsB. If Subs2A provided but not Subs2B, +!> Subs2A contains the total (alpha+beta) substitutions that +!> are permitted. If both Subs2A and Subs2B provided, Subs2A +!> containes the alpha substitutions that are permitted. +!> \endverbatim +!> +!> \param[in] Subs2BIn +!> \verbatim +!> Subs2B is Type(MQC_Vector),Dimension(:),Optional +!> If provided, gives the permitted beta substitutions to +!> include in the Hamiltonian for the column basis. One +!> mqc_vector object should be provided for each of the +!> determinant objects given. If not present, column basis +!> substitutions are equal to row basis substitutions. +!> \endverbatim +!> +!> \param[in] doS2 +!> \verbatim +!> DoS2 is Logical,Optional +!> Logical flag to build the S^2 matrix instead of using +!> Slater Condon rules. The spatial molecular orbital basis +!> overlap (i.e. do not multiply by spin overlap which +!> always leads to zero alpha-beta blocks) must be provided +!> in MO_Core_Ham argument if this option is used) must be provided +!> in MO_Core_Ham argument if this option is used.. +!> \endverbatim +!> +! Authors: +! ======== +!> \author L. M. Thompson +!> \date 2017,2021,2024 +! + Subroutine MQC_Build_CI_Operator(IOut,IPrint,Operator_Matrix,Determinants,Dets2,& + MO_Core_Ham,MO_ERIs,SubsAIn,SubsBIn,Subs2AIn,Subs2BIn,doS2) +! + Implicit None + Integer(kind=int64),Intent(In)::IOut,IPrint + Type(MQC_TwoERIs),Optional,Intent(In)::MO_ERIs + Type(MQC_SCF_Integral),Optional,Intent(In)::MO_Core_Ham + Type(MQC_Determinant),Dimension(:),Intent(In)::Determinants + Type(MQC_Determinant),Dimension(:),Optional,Intent(In)::Dets2 + Type(MQC_Vector),Dimension(:),Optional,Intent(In)::SubsAIn,SubsBIn,Subs2AIn,Subs2BIn + Logical,Optional,Intent(In)::doS2 + Type(MQC_Matrix),Intent(Out)::Operator_Matrix + Integer(kind=int64)::I,J,L_Index,R_Index,NDets1T,NDets2T,R_Index_End + Integer(kind=int64),Dimension(:),Allocatable::SubsA,SubsB,Subs2A,Subs2B,NDets1,NDets2,& + NDets1A,NDets1B,NDets2A,NDets2B + Type(MQC_Bits)::Alpha_String_1,Alpha_String_2,Beta_String_1,Beta_String_2,alphaStrTemp,& + betaStrTemp + Type(MQC_Vector),Dimension(:),Allocatable::lSubsT,lSubsA,lSubsB,rSubsT,rSubsA,rSubsB + Type(MQC_Determinant),Dimension(:),Allocatable::Determinants2 + Logical::SymmFlag,S2Flag + Character(len=12)::SymmStr +! + If(present(doS2)) then + S2Flag = doS2 + else + S2Flag = .false. + EndIf +! + If(.not.present(MO_Core_Ham).and..not.present(MO_ERIs)) & + call mqc_error('No integrals given to MQC_Build_CI_Operator') +! + If(present(Dets2)) then + Determinants2 = Dets2 + SymmFlag = .False. + Else + Determinants2 = Determinants + SymmFlag = .True. + EndIf +! +! Determine different substitutions requested from input, and number of strings +! of different types. +! + Call MQC_Get_Substring_CI_Parameters(LSubsT,RSubsT,LSubsA,RSubsA,LSubsB,RSubsB,& + NDets1,NDets2,NDets1A,NDets2A,NDets1B,NDets2B,NDets1T,NDets2T,SymmFlag,Determinants,& + Dets2,SubsAIn,SubsBIn,Subs2AIn,Subs2BIn) ! ! Allocate output operator matrix ! @@ -21906,25 +23094,11 @@ Subroutine MQC_Build_CI_Operator(IOut,IPrint,Operator_Matrix,Determinants,Dets2, Else R_Index_End = NDets2T EndIf - alpha_string_1 = mqc_bits(0) - beta_string_1 = mqc_bits(0) - Do I = 1, LDetSize - J = mod((L_Index-1)/product(NDets1(1:i-1)),NDets1(i))+1 - call mqc_strings_at_index(iout,iprint,j,Determinants(i),lsubsA(i),lSubsB(i),& - lSubsT(i),AlphaStrTemp,BetaStrTemp) - alpha_string_1 = MQC_concatenate_bits(alphaStrTemp,alpha_string_1) - beta_string_1 = MQC_concatenate_bits(betaStrTemp,beta_string_1) - EndDo + call mqc_get_string_from_substring(iout,iprint,l_index,determinants,& + nDets1,lsubsA,lSubsB,lSubsT,alpha_string_1,beta_string_1) Do R_Index = 1, R_Index_End - alpha_string_2 = mqc_bits(0) - beta_string_2 = mqc_bits(0) - Do I = 1, RDetSize - J = mod((R_Index-1)/product(NDets2(1:i-1)),NDets2(i))+1 - call mqc_strings_at_index(iout,iprint,j,Determinants2(i),rsubsA(i),rSubsB(i),& - rSubsT(i),AlphaStrTemp,BetaStrTemp) - alpha_string_2 = MQC_concatenate_bits(alphaStrTemp,alpha_string_2) - beta_string_2 = MQC_concatenate_bits(betaStrTemp,beta_string_2) - EndDo + call mqc_get_string_from_substring(iout,iprint,r_index,determinants2,& + nDets2,rsubsA,rSubsB,rSubsT,alpha_string_2,beta_string_2) ! write(*,*) '-------------------------------------------' ! write(*,*) ' L_Index: ',L_Index,' R_Index: ',R_Index ! write(*,*) '-------------------------------------------' @@ -21957,9 +23131,9 @@ End Subroutine MQC_Build_CI_Operator ! !===================================================================== ! -! PROCEDURE MQC_GET_STRINGS_AT_INDEX +! PROCEDURE MQC_STRINGS_AT_INDEX ! -!> \brief MQC_GET_STRINGS_AT_INDEX is a subroutine that returns the strings +!> \brief MQC_STRINGS_AT_INDEX is a subroutine that returns the strings !> corresponding to the index of matrices built by MQC_Build_CI_Hamiltonian !> !> \par Purpose: @@ -21967,7 +23141,7 @@ End Subroutine MQC_Build_CI_Operator !> !> \verbatim !> -!> MQC_GET_STRINGS_AT_INDEX is a subroutine that returns the strings corresponding +!> MQCSTRINGS_AT_INDEX is a subroutine that returns the strings corresponding !> to the index of matrices built by MQC_Build_CI_Hamiltonian. !> !> \endverbatim @@ -22001,19 +23175,19 @@ End Subroutine MQC_Build_CI_Operator !> !> \param[in] SubsA !> \verbatim -!> SubsA is Integer(kind=int64),Dimension(:) +!> SubsA is type(mqc_vector) !> Permitted alpha substitutions. !> \endverbatim !> !> \param[in] SubsB !> \verbatim -!> SubsB is Integer(kind=int64),Dimension(:) +!> SubsB is type(mqc_vector) !> Permitted beta substitutions. !> \endverbatim !> !> \param[in] SubsT !> \verbatim -!> SubsT is Integer(kind=int64),Dimension(:) +!> SubsT is type(mqc_vector) !> Permitted total substitutions. !> \endverbatim !>