From 808faf1adacd0ef27ad14e4eb08cadf303ab6c83 Mon Sep 17 00:00:00 2001 From: Matheus Moraes Date: Tue, 18 Feb 2025 15:14:14 -0500 Subject: [PATCH 1/9] Matheus' updade --- src/mqc_algebra.F03 | 684 +++++++++++- src/mqc_est.F03 | 2469 ++++++++++++++++++++++++++++++++++++++++--- 2 files changed, 2996 insertions(+), 157 deletions(-) diff --git a/src/mqc_algebra.F03 b/src/mqc_algebra.F03 index c6e8e8a8..7918076b 100644 --- a/src/mqc_algebra.F03 +++ b/src/mqc_algebra.F03 @@ -87,6 +87,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 @@ -234,6 +236,7 @@ Module MQC_Algebra Procedure, Public::exp => MQC_Matrix_Exp !> \brief Returns the sum of elements in the MQC_Matrix Procedure, Public::sum => MQC_Matrix_Sum + Procedure, Public::psum => MQC_Matrix_Partial_Sum !> \brief Returns the MQC_Matrix raised to a specified power Procedure, Public::power => MQC_Matrix_Power Procedure, Private::at_sca => MQC_Matrix_Scalar_At_Scalar @@ -356,6 +359,7 @@ Module MQC_Algebra module procedure mqc_scalar_havereal Module procedure mqc_vector_havereal module procedure mqc_matrix_havereal + module procedure mqc_r4tensor_havereal end interface ! !> \brief Determines in an array is integer type @@ -363,6 +367,7 @@ Module MQC_Algebra module procedure mqc_scalar_haveinteger module procedure mqc_vector_haveinteger module procedure mqc_matrix_haveinteger + module procedure mqc_r4tensor_haveinteger end interface ! !> \brief Determines in an array is complex type @@ -370,6 +375,7 @@ Module MQC_Algebra module procedure mqc_scalar_havecomplex module procedure mqc_vector_havecomplex module procedure mqc_matrix_havecomplex + module procedure mqc_r4tensor_havecomplex end interface ! !> \brief Sets an array to integer type @@ -382,6 +388,7 @@ Module MQC_Algebra interface mqc_cast_real module procedure mqc_vector_cast_real module procedure mqc_matrix_cast_real + module procedure mqc_r4tensor_cast_real end interface ! !> \brief Sets an array to complex type @@ -620,6 +627,16 @@ Module MQC_Algebra module procedure MQC_Matrix_Scalar_Put_Int end interface ! +!> \brief Do the contraction between an R4Tensor and a Matrix + interface MQC_R4Tensor_Matrix_Contraction + module procedure MQC_R4Tensor_Matrix_Contraction + end interface +! +!> \brief Do a partial contraction between two R4Tensor R4Tensor + interface MQC_R4Tensor_R4Tensor_Partial_Contraction + module procedure MQC_R4Tensor_R4Tensor_Partial_Contraction + end interface +! !> \brief Puts a value in the R4tensor at the specied index interface MQC_R4Tensor_Put module procedure MQC_R4Tensor_Put_Scalar @@ -899,6 +916,10 @@ Module MQC_Algebra Interface Operator (.outer.) Module Procedure MQC_Matrix_Matrix_Outer End Interface +! + Interface Operator (.crossouter.) + Module Procedure MQC_Matrix_Matrix_Cross_Outer + End Interface ! ! Rank 4 tensor operators... ! @@ -2712,6 +2733,61 @@ Function MQC_Scalar_Get_ABS_Value(Scalar) Result(Output) END FUNCTION MQC_Scalar_Get_ABS_Value ! ! +! PROCEDURE MQC_Scalar_Get_Phase +! +!> \brief MQC_Scalar_Get_Phase is a function that returns the +!> phase angle in radian of MQC_scalar variable +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_Scalar_Get_Phase is a function that returns the phase angle (p) +!> in radian of MQC_scalar variable. Such as: +!> Scalar=exp(-ip)*abs(Scalar)=abs(Scalar)(cos(p)-i sin(p)) +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] Scalar +!> \verbatim +!> Scalar is Class(MQC_Scalar) +!> The MQC_Scalar to be tested. +!> \endverbatim +! +! Authors: +! ======== +!> \author M.M.F. Moraes +!> \date 2025 +! + Function MQC_Scalar_Get_Phase(Scalar,thresh_inp) Result(Output) +! +! Variable Declarations + Implicit None + Class(MQC_Scalar),Intent(In)::Scalar + Type(MQC_Scalar),optional::thresh_inp + Type(MQC_Scalar)::Output + real::thresh + + if(present(thresh_inp)) then + thresh=thresh_inp%rval() + else + thresh=1.0E-8 + end if + + if (abs(Scalar).lt.thresh) then + output = 0.0 + else + output = acos(Scalar%rval()/abs(Scalar)) + write(6,*) 'test', output%rval() + end if + return +! + END FUNCTION MQC_Scalar_Get_Phase +! +! ! PROCEDURE MQC_Scalar_Log ! !> \brief MQC_Scalar_Log is a function that returns the @@ -12887,6 +12963,52 @@ Function MQC_Matrix_Cast_Real(MA) Result(MB) End Function MQC_Matrix_Cast_Real ! +! PROCEDURE MQC_R4Tensor_Cast_Real +! +!> \brief MQC_R4Tensor_Cast_Real is a function that converts an MQC r4tensor to +!> its real space +!> +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_R4Tensor_Cast_Real is a function that converts an MQC r4tensor to its real +!> space. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] MA +!> \verbatim +!> RA is Type(MQC_R4Tensor) +!> The MQC r4tensor to convert. +!> \endverbatim +! +! Authors: +! ======== +!> \author M. M. F. Moraes +!> \date 2024 +! + Function MQC_R4Tensor_Cast_Real(RA) Result(RB) +! + Implicit None + Type(MQC_R4Tensor),Intent(In)::RA + Type(MQC_R4Tensor)::RB + + Call MQC_Allocate_R4Tensor(MQC_R4Tensor_Size(RA,1),MQC_R4Tensor_Size(RA,2),& + MQC_R4Tensor_Size(RA,3),MQC_R4Tensor_Size(RA,4),& + RB,'Real',RA%Storage) + if(mqc_have_int(RA)) then + RB%RTen = RA%ITen + elseif(mqc_have_complex(RA)) then + RB%RTen = RA%CTen + endif + + End Function MQC_R4Tensor_Cast_Real +! +! ! ! PROCEDURE MQC_Matrix_Cast_Integer ! @@ -23102,7 +23224,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 @@ -26294,6 +26416,7 @@ Function MQC_Givens_Matrix(m_size,angle,p,q) !> Matrix2 is type(mqc_matrix) !> The second matrix in the outer product. !> \endverbatim +!> ! ! Authors: ! ======== @@ -26312,9 +26435,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 @@ -26326,6 +26449,74 @@ Function MQC_Matrix_Matrix_Outer(matrix1,matrix2) result(outer) endDo end function MQC_Matrix_Matrix_Outer + + +! +! +! PROCEDURE MQC_Matrix_Matrix_Cross_Outer +! +!> \brief MQC_Matrix_Matrix_Cross_Outer is similar to +!> MQC_Matrix_Matrix_Outer but returns a rank-4 with swapped indexes +!> M1(i,l)xM(k,j) instead of M1(i,j)xM2(k,l) +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_Matrix_Matrix_Cross_Outer is a function that returns the rank-4 tensor containing +!> the M1(i,l)xM(k,j) outer product of two MQC matrices . +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] Matrix1 +!> \verbatim +!> Matrix1 is type(mqc_matrix) +!> The first matrix in the outer product. +!> \endverbatim +!> +!> \param[in] Matrix2 +!> \verbatim +!> Matrix2 is type(mqc_matrix) +!> The second matrix in the outer product. +!> \endverbatim +!> +! +! Authors: +! ======== +!> \author M. M. F. Moraes +!> \date 2024 +! + Function MQC_Matrix_Matrix_Cross_Outer(matrix1,matrix2) result(outer) +! +! Variable Declarations. + Implicit None + type(mqc_r4tensor)::outer + type(mqc_matrix),intent(in)::matrix1,matrix2 + integer(kind=int64)::M,N,O,P,i,j,k,l + + M = size(matrix1,1) + N = size(matrix1,2) + O = size(matrix2,1) + P = size(matrix2,2) + + call outer%init(M,P,O,N) + + do i = 1,M + do j = 1,N + do k = 1,O + do l = 1,P + call outer%put(matrix1%at(i,j)*matrix2%at(k,l),i,l,k,j) + endDo + endDo + endDo + endDo + + end function MQC_Matrix_Matrix_Cross_Outer + + ! ! ! PROCEDURE MQC_Matrix_Complex_RealPart @@ -26824,6 +27015,75 @@ function mqc_matrix_sum(A) result(output) endDo ! end function mqc_matrix_sum + +! +! +! PROCEDURE MQC_Matrix_Partial_Sum +! +!> \brief MQC_Matrix_Partial_Sum is a function that returns the sum of all +!> rows or columns elemente of an MQC matrix and returns a MQC vector. +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_Matrix_Partial_Sum is a function that returns the sum of all +!> rows or columns elemente of an MQC matrix and returns a MQC vector. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] A +!> \verbatim +!> A is Class(MQC_Matrix) +!> The name of the MQC_Matrix variable. +!> \endverbatim +!> \param[in] indexnum +!> \verbatim +!> indexnum is integer(kind=int64),optional +!> The index that will be summed over (1 or 2). +!> default=1 +!> \endverbatim +! +! Authors: +! ======== +!> \author M. M. F. Moraes +!> \date 2024 +! + function mqc_matrix_partial_sum(A,indexnum) result(output) +! + Implicit None + class(mqc_matrix),intent(in)::A + integer(kind=int64),optional::indexnum + type(mqc_vector)::output + integer(kind=int64)::i,j,ind + + ind = 1 + if(present(indexnum)) ind = indexnum + + if (ind.eq.1) then + call output%init(size(A,2)) + do j=1,size(A,2) + do i=1,size(A,1) + call output%put(output%at(j) + A%at(i,j) , j) + endDo + endDo + else if (ind.eq.2) then + call output%init(size(A,1)) + do j=1,size(A,1) + do i=1,size(A,2) + call output%put(output%at(j) + A%at(j,i) , j) + endDo + endDo + else + Call MQC_Error_I('Index out of range MQC_Matrix_Partial_Sum', 6, & + 'index', ind ) + end if +! + end function mqc_matrix_partial_sum + ! ! ! PROCEDURE MQC_Matrix_Power @@ -31336,4 +31596,422 @@ function mqc_R4Tensor_complex_imagPart(A) result(output) end function mqc_R4Tensor_complex_imagPart ! ! + +! +! +! PROCEDURE MQC_R4Tensor_Matrix_Contraction +! +!> \brief MQC_R4Tensor_Matrix_Contraction +!> +!> +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_R4Tensor_Matrix_Contraction +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] RA +!> \verbatim +!> RA is Type(MQC_R4Tensor) +!> \endverbatim +!> +!> \param[in] MB +!> \verbatim +!> MB is Type(MQC_Matrix) +!> The MQC matrix multiply. +!> \endverbatim +!> +!> \param[in] inorder +!> \verbatim +!> inorder is Integer(kind=int64),dimension(2),optional +!> Contain the indexes to be contracted. +!> Ex: if inorder=[3,2], the thrid and second indexes of RA is contacted +!> with the first and second of MB, respectivelly. +!> \endverbatim +!> +!> \endverbatim +! +! Authors: +! ======== +!> \author M. M. F. Moraes +!> \date 2024 +! + Function MQC_R4Tensor_Matrix_Contraction(RA,MB,inorder) Result(MC) +! + Implicit None + Type(MQC_R4Tensor),Intent(In)::RA + Type(MQC_Matrix),Intent(In)::MB + Integer(kind=int64),dimension(2),Intent(In),optional::inorder + Type(MQC_R4Tensor)::RAreal + Type(MQC_Matrix)::MC,MCreal,MBreal + Type(MQC_Matrix)::tmpM,tmpR + Integer(kind=int64)::ia,ib,ic,id,i,j,k + Integer(kind=int64)::dima,dimb,dimc,dimd + Integer(kind=int64),dimension(2)::va,vb + + if(present(inorder)) then + ia=inorder(1) + ib=inorder(2) + ic=0 + id=0 + if (ia.eq.ib) Call MQC_Error_I('Duplicated index in MQC_R4Tensor_Matrix_Contraction', 6, & + 'ia', ia,'ib',ib) + if (ia.gt.4.or.ia.le.0.or.ib.gt.4.or.ib.le.0) Call MQC_Error_I(& + 'Index out of range in MQC_R4Tensor_Matrix_Contraction', 6, & + 'ia', ia,'ib',ib) + do i=1,4 + if (i.ne.ia.and.i.ne.ib) then + if(ic.eq.0) then + ic=i + else + id=i + endif + endif + end do + else + ia=1 + ib=2 + ic=3 + id=4 + endif + dima = MQC_R4Tensor_Size(RA,ia) + dimb = MQC_R4Tensor_Size(RA,ib) + dimc = MQC_R4Tensor_Size(RA,ic) + dimd = MQC_R4Tensor_Size(RA,id) + If (MB%NRow.ne.dima.or.MB%NCol.ne.dimb) then + Call MQC_Error_I('R4Tensor and Matrix are not conformable for contraction', 6, & + 'MB%NRow', MB%NRow,'MQC_R4Tensor_Size(RA%ia)',dima,& + 'MB%NCol', MB%NCol,'MQC_R4Tensor_Size(RA%ib)',dimb) + EndIf + + If (MQC_R4Tensor_HaveInteger(RA)) RAreal = MQC_Cast_Real(RA) + If (MQC_Matrix_HaveInteger(MB)) MBreal = MQC_Cast_Real(MB) + If (MQC_R4Tensor_HaveComplex(RA) .or. MQC_Matrix_HaveComplex(MB)) then + Call MQC_Allocate_Matrix(dimc,dimd,MC,'Complex','StorFull') + ElseIf (MQC_R4Tensor_HaveReal(RA) .or. MQC_Matrix_HaveReal(MB)) then + Call MQC_Allocate_Matrix(dimc,dimd,MC,'Real','StorFull') + !!TODO:Consider symm + Else + Call MQC_Allocate_Matrix(dimc,dimd,MC,'Integer','StorFull') + EndIf + + va=[1,dima] + vb=[1,dimb] + do i=1,dimc + do j=1,dimd + select case(ia) + case(1) + select case(ib) + case(2) + tmpR = RA%mat(va,vb,[i],[j]) + case(3) + tmpR = RA%mat(va,[i],vb,[j]) + case(4) + tmpR = RA%mat(va,[i],[j],vb) + end select + case(2) + select case(ib) + case(1) + tmpR = transpose(RA%mat(vb,va,[i],[j])) + case(3) + tmpR = RA%mat([i],va,vb,[j]) + case(4) + tmpR = RA%mat([i],va,[j],vb) + end select + case(3) + select case(ib) + case(1) + tmpR = transpose(RA%mat(vb,[i],va,[j])) + case(2) + tmpR = transpose(RA%mat([i],vb,va,[j])) + case(4) + tmpR = RA%mat([i],[j],va,vb) + end select + case(4) + select case(ib) + case(1) + tmpR = transpose(RA%mat(vb,[i],[j],va)) + case(2) + tmpR = transpose(RA%mat([i],vb,[j],va)) + case(3) + tmpR = transpose(RA%mat([i],[j],vb,va)) + end select + end select + call MC%put(Contraction(tmpR,MB),i,j) + enddo + enddo + + + EndFunction MQC_R4Tensor_Matrix_Contraction + +! +! +! PROCEDURE MQC_R4Tensor_R4Tensor_Partial_Contraction +! +!> \brief MQC_R4Tensor_R4Tensor_Partial_Contraction +!> +!> +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_R4Tensor_R4Tensor_Partial_Contraction +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] RA +!> \verbatim +!> RA is Type(MQC_R4Tensor) +!> The MQC R4Tensor to be partially contracted. +!> \endverbatim +!> +!> \param[in] RB +!> \verbatim +!> RB is Type(MQC_R4Tensor) +!> The MQC R4Tensor to be partially contracted. +!> \endverbatim +!> +!> \param[in] orderA +!> \verbatim +!> inorder is Integer(kind=int64),dimension(4),optional +!> Contain the indexes of RA to be contracted. +!> The two first integers in the object are the ones to be contracted, +!> while the two last are the one to be multiplied. +!> Ex: If orderA=[3,2,4,1], the thrid and second indexes of RA are contacted. +!> The fourth and first are multiplied to form the first and second indexes +!> of the result matrix, respectivelly. +!> \endverbatim +!> +!> \param[in] orderB +!> \verbatim +!> inorder is Integer(kind=int64),dimension(B),optional +!> Contain the indexes of RB to be contracted. +!> The two first integers in the object are the ones to be contracted, +!> while the two last are the one to be multiplied. +!> \endverbatim +!> +!> Ex: orderA=[3,2,4,1] and orderB=[2,1,4,3] +!> sum_{i,j} RA_{ljik} RB_{jilk} = MC_{kl} +!> +!> +!> \endverbatim +! +! Authors: +! ======== +!> \author M. M. F. Moraes +!> \date 2024 +! + Function MQC_R4Tensor_R4Tensor_Partial_Contraction(RA,RB,orderA,orderB) Result(MC) +! + Implicit None + Type(MQC_R4Tensor),Intent(In)::RA,RB + Integer(kind=int64),dimension(4),Intent(In),optional::orderA,orderB + Type(MQC_R4Tensor)::RAreal,RBreal + Type(MQC_Matrix)::MC,MCreal + Type(MQC_Matrix)::tmpRB,tmpRA + Integer(kind=int64)::ia,ja,ka,la,ib,jb,kb,lb,i,j,k,l + Integer(kind=int64)::dimAi,dimAj,dimAk,dimAl,dimBi,dimBj,dimBk,dimBl + Integer(kind=int64),dimension(2)::vi,vj + + if(present(orderA)) then + ia=orderA(1) + ja=orderA(2) + ka=orderA(3) + la=orderA(4) + if (ia.eq.ja) Call MQC_Error_I('Duplicated index in MQC_R4Tensor_R4Tensor_Contraction',& + 6,'ia', ia,'ja',ja) + if (ia.eq.ka) Call MQC_Error_I('Duplicated index in MQC_R4Tensor_R4Tensor_Contraction',& + 6,'ia', ia,'ka',ka) + if (ia.eq.la) Call MQC_Error_I('Duplicated index in MQC_R4Tensor_R4Tensor_Contraction',& + 6,'ia', ia,'la',la) + if (ja.eq.ka) Call MQC_Error_I('Duplicated index in MQC_R4Tensor_R4Tensor_Contraction',& + 6,'ja', ja,'ka',ka) + if (ja.eq.la) Call MQC_Error_I('Duplicated index in MQC_R4Tensor_R4Tensor_Contraction',& + 6,'ja', ja,'la',la) + if (ka.eq.la) Call MQC_Error_I('Duplicated index in MQC_R4Tensor_R4Tensor_Contraction',& + 6,'ka', ka,'la',la) + else + ia=1 + ja=2 + ka=3 + la=4 + endif + if(present(orderB)) then + ib=orderB(1) + jb=orderB(2) + kb=orderB(3) + lb=orderB(4) + if (ib.eq.jb) Call MQC_Error_I('Duplicated index in MQC_R4Tensor_R4Tensor_Partial_Contraction',& + 6,'ib', ib,'jb',jb) + if (ib.eq.kb) Call MQC_Error_I('Duplicated index in MQC_R4Tensor_R4Tensor_Partial_Contraction',& + 6,'ib', ib,'kb',kb) + if (ib.eq.lb) Call MQC_Error_I('Duplicated index in MQC_R4Tensor_R4Tensor_Partial_Contraction',& + 6,'ib', ib,'lb',lb) + if (jb.eq.kb) Call MQC_Error_I('Duplicated index in MQC_R4Tensor_R4Tensor_Partial_Contraction',& + 6,'jb', jb,'kb',kb) + if (jb.eq.lb) Call MQC_Error_I('Duplicated index in MQC_R4Tensor_R4Tensor_Partial_Contraction',& + 6,'jb', jb,'lb',lb) + if (kb.eq.lb) Call MQC_Error_I('Duplicated index in MQC_R4Tensor_R4Tensor_Partial_Contraction',& + 6,'kb', kb,'lb',lb) + else + ib=1 + jb=2 + kb=3 + lb=4 + endif + if (ia.gt.4.or.ia.le.0.or.ib.gt.4.or.ib.le.0.or.& + ja.gt.4.or.ja.le.0.or.jb.gt.4.or.jb.le.0.or.& + ka.gt.4.or.ka.le.0.or.kb.gt.4.or.kb.le.0.or.& + la.gt.4.or.la.le.0.or.lb.gt.4.or.lb.le.0.) Call MQC_Error(& + 'Index out of range in MQC_R4Tensor_R4Tensor_Partial_Contraction', 6) + dimAi = MQC_R4Tensor_Size(RA,ia) + dimAj = MQC_R4Tensor_Size(RA,ja) + dimAk = MQC_R4Tensor_Size(RA,ka) + dimAl = MQC_R4Tensor_Size(RA,la) + + dimBi = MQC_R4Tensor_Size(RB,ib) + dimBj = MQC_R4Tensor_Size(RB,jb) + dimBk = MQC_R4Tensor_Size(RB,kb) + dimBl = MQC_R4Tensor_Size(RB,lb) + + If (dimAi.ne.dimBi.or.dimAj.ne.dimBj.or.& + dimAk.ne.dimBk.or.dimAl.ne.dimBl) Call MQC_Error(& + 'R4Tensor are not conformable for partial contraction', 6) + + If (MQC_R4Tensor_HaveInteger(RA)) RAreal = MQC_Cast_Real(RA) + If (MQC_R4Tensor_HaveInteger(RB)) RBreal = MQC_Cast_Real(RB) + If (MQC_R4Tensor_HaveComplex(RA) .or. MQC_R4Tensor_HaveComplex(RB)) then + Call MQC_Allocate_Matrix(dimAk,dimAl,MC,'Complex','StorFull') + ElseIf (MQC_R4Tensor_HaveReal(RA) .or. MQC_R4Tensor_HaveReal(RB)) then + Call MQC_Allocate_Matrix(dimAk,dimAl,MC,'Real','StorFull') + !!TODO:Consider symm + Else + Call MQC_Allocate_Matrix(dimAk,dimAl,MC,'Integer','StorFull') + EndIf + + vi=[1,dimAi] + vj=[1,dimAj] +!$OMP PARALLEL DO DEFAULT(NONE), & +!$OMP SHARED(RA,RB,vi,vj,ia,ja,la,ka,& +!$OMP ib,jb,lb,kb,MC,dimAl,dimAk),& +!$OMP PRIVATE(tmpRA,tmpRB,k,l) + do k=1,dimAl + do l=1,dimAk + select case(ia) + case(1) + select case(ja) + case(2) + if(ka.eq.3) tmpRA = RA%mat(vi,vj,[k],[l]) + if(ka.eq.4) tmpRA = RA%mat(vi,vj,[l],[k]) + case(3) + if(ka.eq.2) tmpRA = RA%mat(vi,[k],vj,[l]) + if(ka.eq.4) tmpRA = RA%mat(vi,[l],vj,[k]) + case(4) + if(ka.eq.2) tmpRA = RA%mat(vi,[k],[l],vj) + if(ka.eq.3) tmpRA = RA%mat(vi,[l],[k],vj) + end select + case(2) + select case(ja) + case(1) + if(ka.eq.3) tmpRA = transpose(RA%mat(vj,vi,[k],[l])) + if(ka.eq.4) tmpRA = transpose(RA%mat(vj,vi,[l],[k])) + case(3) + if(ka.eq.1) tmpRA = RA%mat([k],vi,vj,[l]) + if(ka.eq.4) tmpRA = RA%mat([l],vi,vj,[k]) + case(4) + if(ka.eq.1) tmpRA = RA%mat([k],vi,[l],vj) + if(ka.eq.3) tmpRA = RA%mat([l],vi,[k],vj) + end select + case(3) + select case(ja) + case(1) + if(ka.eq.2) tmpRA = transpose(RA%mat(vj,[k],vi,[l])) + if(ka.eq.4) tmpRA = transpose(RA%mat(vj,[l],vi,[k])) + case(2) + if(ka.eq.1) tmpRA = transpose(RA%mat([k],vj,vi,[l])) + if(ka.eq.4) tmpRA = transpose(RA%mat([l],vj,vi,[k])) + case(4) + if(ka.eq.1) tmpRA = RA%mat([k],[l],vi,vj) + if(ka.eq.2) tmpRA = RA%mat([l],[k],vi,vj) + end select + case(4) + select case(ib) + case(1) + if(ka.eq.2) tmpRA = transpose(RA%mat(vj,[k],[l],vi)) + if(ka.eq.3) tmpRA = transpose(RA%mat(vj,[l],[k],vi)) + case(2) + if(ka.eq.1) tmpRA = transpose(RA%mat([k],vj,[l],vi)) + if(ka.eq.3) tmpRA = transpose(RA%mat([l],vj,[k],vi)) + case(3) + if(ka.eq.1) tmpRA = transpose(RA%mat([k],[l],vj,vi)) + if(ka.eq.2) tmpRA = transpose(RA%mat([l],[k],vj,vi)) + end select + end select + select case(ib) + case(1) + select case(jb) + case(2) + if(kb.eq.3) tmpRB = RB%mat(vi,vj,[k],[l]) + if(kb.eq.4) tmpRB = RB%mat(vi,vj,[l],[k]) + case(3) + if(kb.eq.2) tmpRB = RB%mat(vi,[k],vj,[l]) + if(kb.eq.4) tmpRB = RB%mat(vi,[l],vj,[k]) + case(4) + if(kb.eq.2) tmpRB = RB%mat(vi,[k],[l],vj) + if(kb.eq.3) tmpRB = RB%mat(vi,[l],[k],vj) + end select + case(2) + select case(ja) + case(1) + if(kb.eq.3) tmpRB = transpose(RB%mat(vj,vi,[k],[l])) + if(kb.eq.4) tmpRB = transpose(RB%mat(vj,vi,[l],[k])) + case(3) + if(kb.eq.1) tmpRB = RB%mat([k],vi,vj,[l]) + if(kb.eq.4) tmpRB = RB%mat([l],vi,vj,[k]) + case(4) + if(kb.eq.1) tmpRB = RB%mat([k],vi,[l],vj) + if(kb.eq.3) tmpRB = RB%mat([l],vi,[k],vj) + end select + case(3) + select case(ja) + case(1) + if(kb.eq.2) tmpRB = transpose(RB%mat(vj,[k],vi,[l])) + if(kb.eq.4) tmpRB = transpose(RB%mat(vj,[l],vi,[k])) + case(2) + if(kb.eq.1) tmpRB = transpose(RB%mat([k],vj,vi,[l])) + if(kb.eq.4) tmpRB = transpose(RB%mat([l],vj,vi,[k])) + case(4) + if(kb.eq.1) tmpRB = RB%mat([k],[l],vi,vj) + if(kb.eq.2) tmpRB = RB%mat([l],[k],vi,vj) + end select + case(4) + select case(ib) + case(1) + if(kb.eq.2) tmpRB = transpose(RB%mat(vj,[k],[l],vi)) + if(kb.eq.3) tmpRB = transpose(RB%mat(vj,[l],[k],vi)) + case(2) + if(kb.eq.1) tmpRB = transpose(RB%mat([k],vj,[l],vi)) + if(kb.eq.3) tmpRB = transpose(RB%mat([l],vj,[k],vi)) + case(3) + if(kb.eq.1) tmpRB = transpose(RB%mat([k],[l],vj,vi)) + if(kb.eq.2) tmpRB = transpose(RB%mat([l],[k],vj,vi)) + end select + end select + call MC%put(Contraction(tmpRA,tmpRB),k,l) + enddo + enddo +!$OMP END PARALLEL DO + + + + EndFunction MQC_R4Tensor_R4Tensor_Partial_Contraction + End Module MQC_Algebra diff --git a/src/mqc_est.F03 b/src/mqc_est.F03 index 53d99040..78858c26 100644 --- a/src/mqc_est.F03 +++ b/src/mqc_est.F03 @@ -89,7 +89,7 @@ Module MQC_EST Procedure, Private::hasAlphaBeta => mqc_integral_has_alphabeta Procedure, Private::hasBetaAlpha => mqc_integral_has_betaalpha Procedure, Public::type => mqc_integral_array_type - Procedure, Private::blockSize => mqc_integral_dimension + Procedure, Public::blockSize => mqc_integral_dimension !> \brief Return the label of the MQC_SCF_Integral object Procedure, Public::getLabel => mqc_integral_array_name !> \brief Set the label of the MQC_SCF_Integral object @@ -119,10 +119,14 @@ Module MQC_EST Procedure, Public::power => mqc_scf_integral_power !> \brief Returns an MQC_SCF_Integral object raised to a power (function) Procedure, Public::powerf => mqc_scf_integral_power_func + Procedure, Public::put => mqc_integral_put 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 +!> \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 !> \brief Returns a subset of an MQC_SCF_Integral object Procedure, Public::orbitals => mqc_integral_output_orbitals !> \brief Swaps rows/columns of an MQC_SCF_Integral object @@ -131,6 +135,12 @@ Module MQC_EST Procedure, Public::combine => mqc_integral_combine_orbitals !> \brief Swaps ab and ba blocks of an MQC_SCF_Integral objects Procedure, Public::swapODB => mqc_integral_swap_offDiagBlocks +!> \brief Outer product between 2 integrals + Procedure, Public::outer => mqc_integral_outer +!> \brief Outer Product of two variables into a two ERIs + Procedure, Public::outer2ERI => mqc_integral_integral_outer_product +!> \brief Outer Product of two variables into a two ERIs + Procedure, Public::crossouter2ERI => mqc_integral_integral_crossouter_product End Type ! !> \brief EST intermediate diagonal matrix variable @@ -157,6 +167,8 @@ Module MQC_EST !> \brief Returns the value of a specified element in an MQC_SCF_Eigenvalues !> object Procedure, Public::at => mqc_eigenvalues_at +!> \brief Outer Product of two variables into a SCF Integral + Procedure, Public::outer => mqc_eigen_eigen_outer_product End Type ! !> \brief EST intermediates of single determinant wavefunctions @@ -284,6 +296,12 @@ Module MQC_EST Module Procedure mqc_scf_integral_contraction Module Procedure mqc_eri_integral_contraction Module Procedure mqc_eri_r4tensor_contraction + Module Procedure mqc_eri_eri_contraction + End Interface +! +!> \brief Contracts two objects + Interface PartialContraction + Module Procedure mqc_eri_eri_partial_contraction End Interface ! !> \brief Converts object from alpha-beta order to energy order @@ -316,12 +334,19 @@ Module MQC_EST !> \brief Sums two variables Interface Operator (+) Module Procedure MQC_Integral_Sum + Module Procedure MQC_ERI_Sum + End Interface +! +!> \brief Sums two variables dimensions presevend + Interface Operator (.sum.) + Module Procedure MQC_Integral_Merge_Sum End Interface ! !> \brief Subtracts two variables Interface Operator (-) Module Procedure MQC_Integral_Difference Module Procedure MQC_Matrix_Integral_Difference + Module Procedure MQC_ERI_Difference End Interface ! !> \brief Multiplies two variables @@ -2606,7 +2631,7 @@ subroutine mqc_twoeris_transform(twoERIs,integralType,mo_coefficients,mo_coeff_2 call mqc_twoeris_allocate(twoERIs,'symm','raffenetti3',aaaa) return elseIf(my_integralType.eq.'molecular') then - if(.not.present(mo_coefficients)) call mqc_error_l('MO oefficients myst be present to& + if(.not.present(mo_coefficients)) call mqc_error_l('MO oefficients must be present to& & transform AO 2ERIs to MO 2ERIs',6,'present(mo_coefficients)',present(mo_coefficients)) if(present(mo_coeff_2)) then call twoERI_Trans(iOut,iPrint,mo_coefficients,twoERIs,temp,mo_coeff_2) @@ -2903,16 +2928,20 @@ end subroutine mqc_integral_initialize !> = 'full': outputs the full MQC_SCF_Integral. !> = 'alpha': outputs the alpha-alpha spin block. !> = 'alpha-alpha': outputs the alpha-alpha spin block. +!> = 'aa': outputs the alpha-alpha spin block. !> = 'beta': outputs the beta-beta spin block. !> = 'beta-beta': outputs the beta-beta spin block. -!> = 'alpha-beta' outputs the alpha-beta spin block. -!> = 'beta-alpha' outputs the beta-alpha spin block. +!> = 'bb': outputs the beta-beta spin block. +!> = 'alpha-beta': outputs the alpha-beta spin block. +!> = 'ab': outputs the alpha-beta spin block. +!> = 'beta-alpha': outputs the beta-alpha spin block. +!> = 'ba': outputs the beta-alpha spin block. !> \endverbatim ! ! Authors: ! ======== -!> \author L. M. Thompson -!> \date 2016 +!> \author L. M. Thompson, M. M. F. Moraes +!> \date 2016, 2024 ! function mqc_integral_output_block(integral,blockName) result(matrixOut) ! @@ -3009,14 +3038,14 @@ function mqc_integral_output_block(integral,blockName) result(matrixOut) ! 'integral%hasBetaAlpha()', integral%hasBetaAlpha() ) ! endIf endIf - case('alpha') + case('alpha','alpha-alpha','aa') if (integral%hasAlpha()) then matrixOut = integral%alpha else call mqc_error_L('block does not exist in mqc_integral_output_block', 6, & 'integral%hasAlpha()', integral%hasAlpha() ) endIf - case('beta') + case('beta','beta-beta','bb') if (integral%hasBeta()) then matrixOut = integral%beta elseIf (integral%type().eq.'space') then @@ -3030,14 +3059,7 @@ function mqc_integral_output_block(integral,blockName) result(matrixOut) call mqc_error_A('block does not exist in mqc_integral_output_block', 6, & 'integral%type()', integral%type() ) endIf - case('alpha-alpha') - if (integral%hasAlpha()) then - matrixOut = integral%alpha - else - call mqc_error_L('block does not exist in mqc_integral_output_block', 6, & - 'integral%hasAlpha()', integral%hasAlpha() ) - endIf - case('alpha-beta') + case('alpha-beta','alphabeta','ab') if (integral%hasAlphaBeta()) then matrixOut = integral%alphaBeta else @@ -3064,7 +3086,7 @@ function mqc_integral_output_block(integral,blockName) result(matrixOut) endIf endIf endIf - case('beta-alpha') + case('beta-alpha','betaalpha','ba') if (integral%hasBetaAlpha()) then matrixOut = integral%betaAlpha else @@ -3090,20 +3112,6 @@ function mqc_integral_output_block(integral,blockName) result(matrixOut) endIf endIf endIf - case('beta-beta') - if (integral%hasBeta()) then - matrixOut = integral%beta - elseIf (integral%type().eq.'space') then - if (integral%hasAlpha()) then - matrixOut = integral%alpha - else - call mqc_error_L('block does not exist in mqc_integral_output_block', 6, & - 'integral%hasAlpha()', integral%hasAlpha() ) - endIf - else - call mqc_error_A('block does not exist in mqc_integral_output_block', 6, & - 'integral%type()', integral%type() ) - endIf case default call mqc_error_A('unrecognised block name in mqc_integral_output_block', 6, & 'myBlockName', myBlockName ) @@ -4192,8 +4200,8 @@ end function mqc_detstring_print !> \verbatim !> BlockName is character(len=*),optional !> = 'full': outputs the full MQC_SCF_Integral. -!> = 'alpha': outputs the alpha spin block. -!> = 'beta': outputs the beta spin block. +!> = 'alpha','a': outputs the alpha spin block. +!> = 'beta','b': outputs the beta spin block. !> \endverbatim ! ! Authors: @@ -4260,14 +4268,14 @@ function mqc_eigenvalues_output_block(eigenvalues,blockName) result(vectorOut) eigenvalues%hasBeta()) endIf endIf - case('alpha') + case('alpha','a') if (eigenvalues%hasAlpha()) then vectorOut = eigenvalues%alpha else call mqc_error_L('block does not exist in mqc_eigenvalues_output_block', 6, & 'eigenvalues%hasAlpha()', eigenvalues%hasAlpha() ) endIf - case('beta') + case('beta','b') if (eigenvalues%hasBeta()) then vectorOut = eigenvalues%beta elseIf (eigenvalues%type().eq.'space') then @@ -5431,6 +5439,410 @@ function mqc_integral_sum(integralA,integralB) result(integralOut) end function mqc_integral_sum ! ! +! PROCEDURE MQC_Integral_Merge_Sum +! +!> \brief MQC_Integral_Merge_Sum is used to sum two MQC integral type variables +!> the result preserves the size of the largest integral object +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_Integral_Merge_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 M. M. F. Moraes +!> \date 2024 +! + function mqc_integral_merge_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 +! + select case(integralA%array_type) + case('space') + select case (integralB%array_type) + case('space') + if(integralA%hasAlpha().and.integralB%hasAlpha()) tmpMatrixAlpha = integralA%alpha+integralB%alpha + call mqc_integral_allocate(integralOut,'','space',tmpMatrixAlpha) + case('spin') + if(integralA%hasAlpha().and.integralB%hasAlpha()) then + tmpMatrixAlpha = integralA%alpha+integralB%alpha + elseif(integralB%hasAlpha()) then + tmpMatrixAlpha = integralB%alpha + elseif(integralA%hasAlpha()) then + tmpMatrixAlpha = integralA%alpha + end if + if(integralA%hasAlpha().and.integralB%hasBeta()) then + tmpMatrixBeta = integralA%alpha+integralB%beta + elseif(integralB%hasBeta()) then + tmpMatrixBeta = integralB%beta + elseif(integralA%hasAlpha()) then + tmpMatrixBeta = integralA%alpha + end if + call mqc_integral_allocate(integralOut,'','spin',tmpMatrixAlpha,tmpMatrixBeta) + case('general') + if(integralA%hasAlpha().and.integralB%hasAlpha()) then + tmpMatrixAlpha = integralA%alpha+integralB%alpha + elseif(integralB%hasAlpha()) then + tmpMatrixAlpha = integralB%alpha + elseif(integralA%hasAlpha()) then + tmpMatrixAlpha = integralA%alpha + end if + if(integralA%hasAlpha().and.integralB%hasBeta()) then + tmpMatrixBeta = integralA%alpha+integralB%beta + elseif(integralB%hasBeta()) then + tmpMatrixBeta = integralB%beta + elseif(integralA%hasAlpha()) then + tmpMatrixBeta = integralA%alpha + end if + if(integralB%hasAlphaBeta()) then + tmpMatrixAlphaBeta = integralB%alphaBeta + end if + if(integralB%hasBetaAlpha()) then + tmpMatrixBetaAlpha = integralB%betaAlpha + end if + call mqc_integral_allocate(integralOut,'','general',tmpMatrixAlpha,tmpMatrixBeta, & + tmpMatrixAlphaBeta,tmpMatrixBetaAlpha) + case default + call mqc_error_A('Unknown integral type in mqc_integral_merge_sum', 6, & + 'integralB%array_type', integralB%array_type ) + end select + case('spin') + select case (integralB%array_type) + case('space') + if(integralA%hasAlpha().and.integralB%hasAlpha()) then + tmpMatrixAlpha = integralA%alpha+integralB%alpha + elseif(integralB%hasAlpha()) then + tmpMatrixAlpha = integralB%alpha + elseif(integralA%hasAlpha()) then + tmpMatrixAlpha = integralA%alpha + end if + if(integralA%hasBeta().and.integralB%hasAlpha()) then + tmpMatrixBeta = integralB%alpha+integralA%beta + elseif(integralA%hasBeta()) then + tmpMatrixBeta = integralA%beta + elseif(integralB%hasAlpha()) then + tmpMatrixBeta = integralB%alpha + end if + call mqc_integral_allocate(integralOut,'','spin',tmpMatrixAlpha,tmpMatrixBeta) + case('spin') + if(integralA%hasAlpha().and.integralB%hasAlpha()) then + tmpMatrixAlpha = integralA%alpha+integralB%alpha + elseif(integralB%hasAlpha()) then + tmpMatrixAlpha = integralB%alpha + elseif(integralA%hasAlpha()) then + tmpMatrixAlpha = integralA%alpha + end if + if(integralA%hasBeta().and.integralB%hasBeta()) then + tmpMatrixBeta = integralA%beta+integralB%beta + elseif(integralB%hasBeta()) then + tmpMatrixBeta = integralB%beta + elseif(integralA%hasBeta()) then + tmpMatrixBeta = integralA%beta + end if + call mqc_integral_allocate(integralOut,'','spin',tmpMatrixAlpha,tmpMatrixBeta) + case('general') + if(integralA%hasAlpha().and.integralB%hasAlpha()) then + tmpMatrixAlpha = integralA%alpha+integralB%alpha + elseif(integralB%hasAlpha()) then + tmpMatrixAlpha = integralB%alpha + elseif(integralA%hasAlpha()) then + tmpMatrixAlpha = integralA%alpha + end if + if(integralA%hasBeta().and.integralB%hasBeta()) then + tmpMatrixBeta = integralA%beta+integralB%beta + elseif(integralB%hasBeta()) then + tmpMatrixBeta = integralB%beta + elseif(integralA%hasBeta()) then + tmpMatrixBeta = integralA%beta + end if + if(integralB%hasAlphaBeta()) then + tmpMatrixAlphaBeta = integralB%alphaBeta + end if + if(integralB%hasBetaAlpha()) then + tmpMatrixBetaAlpha = integralB%betaAlpha + end if + call mqc_integral_allocate(integralOut,'','general',tmpMatrixAlpha,tmpMatrixBeta, & + tmpMatrixAlphaBeta,tmpMatrixBetaAlpha) + case default + call mqc_error_A('Unknown integral type in mqc_integral_merge_sum', 6, & + 'integralB%array_type', integralB%array_type ) + end select + case('general') + select case (integralB%array_type) + case('space') + if(integralA%hasAlpha().and.integralB%hasAlpha()) then + tmpMatrixAlpha = integralA%alpha+integralB%alpha + elseif(integralB%hasAlpha()) then + tmpMatrixAlpha = integralB%alpha + elseif(integralA%hasAlpha()) then + tmpMatrixAlpha = integralA%alpha + end if + if(integralB%hasAlpha().and.integralA%hasBeta()) then + tmpMatrixBeta = integralA%alpha+integralB%beta + elseif(integralA%hasBeta()) then + tmpMatrixBeta = integralA%beta + elseif(integralB%hasAlpha()) then + tmpMatrixBeta = integralB%alpha + end if + if(integralA%hasAlphaBeta()) then + tmpMatrixAlphaBeta = integralA%alphaBeta + end if + if(integralA%hasBetaAlpha()) then + tmpMatrixBetaAlpha = integralA%betaAlpha + end if + call mqc_integral_allocate(integralOut,'','general',tmpMatrixAlpha,tmpMatrixBeta, & + tmpMatrixAlphaBeta,tmpMatrixBetaAlpha) + case('spin') + if(integralA%hasAlpha().and.integralB%hasAlpha()) then + tmpMatrixAlpha = integralA%alpha+integralB%alpha + elseif(integralB%hasAlpha()) then + tmpMatrixAlpha = integralB%alpha + elseif(integralA%hasAlpha()) then + tmpMatrixAlpha = integralA%alpha + end if + if(integralA%hasBeta().and.integralB%hasBeta()) then + tmpMatrixBeta = integralA%beta+integralB%beta + elseif(integralB%hasBeta()) then + tmpMatrixBeta = integralB%beta + elseif(integralA%hasBeta()) then + tmpMatrixBeta = integralA%beta + end if + if(integralA%hasAlphaBeta()) then + tmpMatrixAlphaBeta = integralA%alphaBeta + end if + if(integralA%hasBetaAlpha()) then + tmpMatrixBetaAlpha = integralA%betaAlpha + end if + call mqc_integral_allocate(integralOut,'','general',tmpMatrixAlpha,tmpMatrixBeta, & + tmpMatrixAlphaBeta,tmpMatrixBetaAlpha) + case('general') + if(integralA%hasAlpha().and.integralB%hasAlpha()) then + tmpMatrixAlpha = integralA%alpha+integralB%alpha + elseif(integralB%hasAlpha()) then + tmpMatrixAlpha = integralB%alpha + elseif(integralA%hasAlpha()) then + tmpMatrixAlpha = integralA%alpha + end if + if(integralA%hasBeta().and.integralB%hasBeta()) then + tmpMatrixBeta = integralA%beta+integralB%beta + elseif(integralB%hasBeta()) then + tmpMatrixBeta = integralB%beta + elseif(integralA%hasBeta()) then + tmpMatrixBeta = integralA%beta + end if + if(integralA%hasAlphaBeta().and.integralB%hasAlphaBeta()) then + tmpMatrixAlphaBeta = integralA%alphaBeta+integralB%alphaBeta + elseif(integralA%hasAlphaBeta()) then + tmpMatrixAlphaBeta = integralA%alphaBeta + elseif(integralB%hasAlphaBeta()) then + tmpMatrixAlphaBeta = integralB%alphaBeta + end if + if(integralA%hasBetaAlpha().and.integralB%hasBetaAlpha()) then + tmpMatrixBetaAlpha = integralA%betaAlpha+integralB%betaAlpha + elseif(integralA%hasBetaAlpha()) then + tmpMatrixBetaAlpha = integralA%betaAlpha + elseif(integralB%hasBetaAlpha()) then + tmpMatrixBetaAlpha = integralB%betaAlpha + endif + call mqc_integral_allocate(integralOut,'','general',tmpMatrixAlpha,tmpMatrixBeta, & + tmpMatrixAlphaBeta,tmpMatrixBetaAlpha) + case default + call mqc_error_A('Unknown integral type in mqc_integral_merge_sum', 6, & + 'integralB%array_type', integralB%array_type ) + end select + case default + call mqc_error_A('Unknown integral type in mqc_integral_merge_sum', 6, & + 'integralA%array_type', integralA%array_type ) + end select + + end function mqc_integral_merge_sum +! +! +! 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 ) + end select + case('spin') + select case (ERIB%integraltype) + case('space') + aaaa = ERIA%getBlock('aaaa')+ERIB%getBlock('aaaa') + bbbb = ERIA%getBlock('bbbb')+ERIB%getBlock('aaaa') + aabb = ERIA%getBlock('aabb')+ERIB%getBlock('aaaa') + bbaa = ERIA%getBlock('bbaa')+ERIB%getBlock('aaaa') + call mqc_twoeris_allocate(ERIOut,'full','spin',aaaa,bbbb,& + aabb,bbaa) + case('spin') + aaaa = ERIA%getBlock('aaaa')+ERIB%getBlock('aaaa') + bbbb = ERIA%getBlock('bbbb')+ERIB%getBlock('bbbb') + aabb = ERIA%getBlock('aabb')+ERIB%getBlock('aabb') + bbaa = ERIA%getBlock('bbaa')+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('bbbb')+ERIB%getBlock('bbbb') + aabb = ERIA%getBlock('aabb')+ERIB%getBlock('aabb') + bbaa = ERIA%getBlock('bbaa')+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 ) + end select + case('general') + select case (ERIB%integraltype) + case('space') + aaaa = ERIA%getBlock('aaaa')+ERIB%getBlock('aaaa') + bbbb = ERIA%getBlock('bbbb')+ERIB%getBlock('aaaa') + aabb = ERIA%getBlock('aabb')+ERIB%getBlock('aaaa') + bbaa = ERIA%getBlock('bbaa')+ERIB%getBlock('aaaa') + call mqc_twoeris_allocate(ERIOut,'full','general',aaaa,bbbb,& + aabb,bbaa,& + ERIA%getBlock('abab'), ERIA%getBlock('abba'),& + ERIA%getBlock('baab'), ERIA%getBlock('baba'),& + ERIA%getBlock('aaab'), ERIA%getBlock('aaba'),& + ERIA%getBlock('abaa'), ERIA%getBlock('baaa'),& + ERIA%getBlock('bbba'), ERIA%getBlock('bbab'),& + ERIA%getBlock('babb'), ERIA%getBlock('abbb')) + case('spin') + aaaa = ERIA%getBlock('aaaa')+ERIB%getBlock('aaaa') + bbbb = ERIA%getBlock('bbbb')+ERIB%getBlock('bbbb') + aabb = ERIA%getBlock('aabb')+ERIB%getBlock('aabb') + bbaa = ERIA%getBlock('bbaa')+ERIB%getBlock('bbaa') + call mqc_twoeris_allocate(ERIOut,'full','general',aaaa,bbbb,& + aabb,bbaa,& + ERIA%getBlock('abab'), ERIA%getBlock('abba'),& + ERIA%getBlock('baab'), ERIA%getBlock('baba'),& + ERIA%getBlock('aaab'), ERIA%getBlock('aaba'),& + ERIA%getBlock('abaa'), ERIA%getBlock('baaa'),& + ERIA%getBlock('bbba'), ERIA%getBlock('bbab'),& + ERIA%getBlock('babb'), ERIA%getBlock('abbb')) + case('general') + aaaa = ERIA%getBlock('aaaa')+ERIB%getBlock('aaaa') + bbbb = ERIA%getBlock('bbbb')+ERIB%getBlock('bbbb') + aabb = ERIA%getBlock('aabb')+ERIB%getBlock('aabb') + bbaa = ERIA%getBlock('bbaa')+ERIB%getBlock('bbaa') + aaab = ERIA%getBlock('aaab')+ERIB%getBlock('aaab') + aaba = ERIA%getBlock('aaba')+ERIB%getBlock('aaba') + abaa = ERIA%getBlock('abaa')+ERIB%getBlock('abaa') + baaa = ERIA%getBlock('baaa')+ERIB%getBlock('baaa') + bbba = ERIA%getBlock('bbba')+ERIB%getBlock('bbba') + bbab = ERIA%getBlock('bbab')+ERIB%getBlock('bbab') + babb = ERIA%getBlock('babb')+ERIB%getBlock('babb') + abbb = ERIA%getBlock('abbb')+ERIB%getBlock('abbb') + abba = ERIA%getBlock('abba')+ERIB%getBlock('abba') + abab = ERIA%getBlock('abab')+ERIB%getBlock('abab') + baab = ERIA%getBlock('baab')+ERIB%getBlock('baab') + baba = ERIA%getBlock('baba')+ERIB%getBlock('baba') + call mqc_twoeris_allocate(ERIOut,'full','general',aaaa,bbbb,& + aabb,bbaa,abab,abba,baab,baba,aaab,aaba,abaa,baaa,bbba,bbab,babb,abbb) + case default + call mqc_error_A('Unknown twoERI type in mqc_eri_difference', 6, & + 'ERIB%integraltype', ERIB%integraltype ) + end select + case default + call mqc_error_A('Unknown twoERI type in mqc_eri_difference', 6, & + 'ERIA%integraltype', ERIA%integraltype ) + end select + end function mqc_eri_sum +! +! +! ! PROCEDURE MQC_Integral_Difference ! !> \brief MQC_Integral_Difference is used to subtract two MQC integral type @@ -5664,6 +6076,170 @@ function mqc_matrix_integral_difference(matrixA,integralB) result(integralOut) end function mqc_matrix_integral_difference ! ! +! PROCEDURE MQC_ERI_Difference +! +!> \brief MQC_ERI_Difference is used to subtract two MQC twoERIs type +!> variables +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_ERI_Difference is used to subtract two MQC twoERIs type variables. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] ERIA +!> \verbatim +!> ERIA is type(mqc_twoERIs) +!> The MQC twoERIs variable to subtract from. +!> \endverbatim +!> +!> \param[in] ERIB +!> \verbatim +!> ERIB is type(mqc_twoERIs) +!> The MQC twoERIs variable to subtract. +!> \endverbatim +! +! Authors: +! ======== +!> \author M. M. F. Moraes +!> \date 2024 +! + function mqc_eri_difference(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,& + (-1)*ERIB%getBlock('abab'),(-1)*ERIB%getBlock('abba'),& + (-1)*ERIB%getBlock('baab'),(-1)*ERIB%getBlock('baba'),& + (-1)*ERIB%getBlock('aaab'),(-1)*ERIB%getBlock('aaba'),& + (-1)*ERIB%getBlock('abaa'),(-1)*ERIB%getBlock('baaa'),& + (-1)*ERIB%getBlock('bbba'),(-1)*ERIB%getBlock('bbab'),& + (-1)*ERIB%getBlock('babb'),(-1)*ERIB%getBlock('abbb')) + case default + call mqc_error_A('Unknown twoERI type in mqc_eri_difference', 6, & + 'ERIB%integraltype', ERIB%integraltype ) + end select + case('spin') + select case (ERIB%integraltype) + case('space') + aaaa = ERIA%getBlock('aaaa')-ERIB%getBlock('aaaa') + bbbb = ERIA%getBlock('bbbb')-ERIB%getBlock('aaaa') + aabb = ERIA%getBlock('aabb')-ERIB%getBlock('aaaa') + bbaa = ERIA%getBlock('bbaa')-ERIB%getBlock('aaaa') + call mqc_twoeris_allocate(ERIOut,'full','spin',aaaa,bbbb,& + aabb,bbaa) + case('spin') + aaaa = ERIA%getBlock('aaaa')-ERIB%getBlock('aaaa') + bbbb = ERIA%getBlock('bbbb')-ERIB%getBlock('bbbb') + aabb = ERIA%getBlock('aabb')-ERIB%getBlock('aabb') + bbaa = ERIA%getBlock('bbaa')-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('bbbb')-ERIB%getBlock('bbbb') + aabb = ERIA%getBlock('aabb')-ERIB%getBlock('aabb') + bbaa = ERIA%getBlock('bbaa')-ERIB%getBlock('bbaa') + call mqc_twoeris_allocate(ERIOut,'full','general',aaaa,bbbb,& + aabb,bbaa,& + (-1)*ERIB%getBlock('abab'),(-1)*ERIB%getBlock('abba'),& + (-1)*ERIB%getBlock('baab'),(-1)*ERIB%getBlock('baba'),& + (-1)*ERIB%getBlock('aaab'),(-1)*ERIB%getBlock('aaba'),& + (-1)*ERIB%getBlock('abaa'),(-1)*ERIB%getBlock('baaa'),& + (-1)*ERIB%getBlock('bbba'),(-1)*ERIB%getBlock('bbab'),& + (-1)*ERIB%getBlock('babb'),(-1)*ERIB%getBlock('abbb')) + case default + call mqc_error_A('Unknown twoERI type in mqc_eri_difference', 6, & + 'ERIB%integraltype', ERIB%integraltype ) + end select + case('general') + select case (ERIB%integraltype) + case('space') + aaaa = ERIA%getBlock('aaaa')-ERIB%getBlock('aaaa') + bbbb = ERIA%getBlock('bbbb')-ERIB%getBlock('aaaa') + aabb = ERIA%getBlock('aabb')-ERIB%getBlock('aaaa') + bbaa = ERIA%getBlock('bbaa')-ERIB%getBlock('aaaa') + call mqc_twoeris_allocate(ERIOut,'full','general',aaaa,bbbb,& + aabb,bbaa,& + ERIA%getBlock('abab'), ERIA%getBlock('abba'),& + ERIA%getBlock('baab'), ERIA%getBlock('baba'),& + ERIA%getBlock('aaab'), ERIA%getBlock('aaba'),& + ERIA%getBlock('abaa'), ERIA%getBlock('baaa'),& + ERIA%getBlock('bbba'), ERIA%getBlock('bbab'),& + ERIA%getBlock('babb'), ERIA%getBlock('abbb')) + case('spin') + aaaa = ERIA%getBlock('aaaa')-ERIB%getBlock('aaaa') + bbbb = ERIA%getBlock('bbbb')-ERIB%getBlock('bbbb') + aabb = ERIA%getBlock('aabb')-ERIB%getBlock('aabb') + bbaa = ERIA%getBlock('bbaa')-ERIB%getBlock('bbaa') + call mqc_twoeris_allocate(ERIOut,'full','general',aaaa,bbbb,& + aabb,bbaa,& + ERIA%getBlock('abab'), ERIA%getBlock('abba'),& + ERIA%getBlock('baab'), ERIA%getBlock('baba'),& + ERIA%getBlock('aaab'), ERIA%getBlock('aaba'),& + ERIA%getBlock('abaa'), ERIA%getBlock('baaa'),& + ERIA%getBlock('bbba'), ERIA%getBlock('bbab'),& + ERIA%getBlock('babb'), ERIA%getBlock('abbb')) + case('general') + aaaa = ERIA%getBlock('aaaa')-ERIB%getBlock('aaaa') + bbbb = ERIA%getBlock('bbbb')-ERIB%getBlock('bbbb') + aabb = ERIA%getBlock('aabb')-ERIB%getBlock('aabb') + bbaa = ERIA%getBlock('bbaa')-ERIB%getBlock('bbaa') + aaab = ERIA%getBlock('aaab')-ERIB%getBlock('aaab') + aaba = ERIA%getBlock('aaba')-ERIB%getBlock('aaba') + abaa = ERIA%getBlock('abaa')-ERIB%getBlock('abaa') + baaa = ERIA%getBlock('baaa')-ERIB%getBlock('baaa') + bbba = ERIA%getBlock('bbba')-ERIB%getBlock('bbba') + bbab = ERIA%getBlock('bbab')-ERIB%getBlock('bbab') + babb = ERIA%getBlock('babb')-ERIB%getBlock('babb') + abbb = ERIA%getBlock('abbb')-ERIB%getBlock('abbb') + abba = ERIA%getBlock('abba')-ERIB%getBlock('abba') + abab = ERIA%getBlock('abab')-ERIB%getBlock('abab') + baab = ERIA%getBlock('baab')-ERIB%getBlock('baab') + baba = ERIA%getBlock('baba')-ERIB%getBlock('baba') + call mqc_twoeris_allocate(ERIOut,'full','general',aaaa,bbbb,& + aabb,bbaa,abab,abba,baab,baba,aaab,aaba,abaa,baaa,bbba,bbab,babb,abbb) + case default + call mqc_error_A('Unknown twoERI type in mqc_eri_difference', 6, & + 'ERIB%integraltype', ERIB%integraltype ) + end select + case default + call mqc_error_A('Unknown twoERI type in mqc_eri_difference', 6, & + 'ERIA%integraltype', ERIA%integraltype ) + end select + + end function mqc_eri_difference +! +! +! ! PROCEDURE MQC_Integral_Integral_Multiply ! !> \brief MQC_Integral_Integral_Multiply is used to multiply two MQC integral @@ -5908,6 +6484,9 @@ function mqc_scalar_integral_multiply(scalarIn,integral) result(integralOut) scalar = scalarIn type is (mqc_scalar) scalar = scalarIn + type is (mqc_scf_integral) + integralOut = MQC_Element_Integral_Multiply(scalarIn,integral) + return class default call mqc_error_I('scalarIn type not defined in mqc_scalar_integral_multiply',6) end select @@ -6002,6 +6581,9 @@ function mqc_integral_scalar_multiply(integral,scalarIn) result(integralOut) scalar = scalarIn type is (mqc_scalar) scalar = scalarIn + type is (mqc_scf_integral) + integralOut = mqc_element_integral_multiply(integral,scalarIn) + return class default call mqc_error_I('scalarIn type not defined in mqc_integral_scalar_multiply',6) end select @@ -6030,64 +6612,518 @@ function mqc_integral_scalar_multiply(integral,scalarIn) result(integralOut) tmpMatrixAlphaBeta,tmpMatrixBetaAlpha) case default call mqc_error_A('Unknown integral type in mqc_integral_scalar_multiply', 6, & - 'integralA%array_type', integral%array_type ) + 'integral%array_type', integral%array_type ) end select ! end function mqc_integral_scalar_multiply + ! +! PROCEDURE MQC_Element_Integral_Multiply ! -! PROCEDURE MQC_Integral_Eigenvalues_Multiply -! -!> \brief MQC_Integral_Eigenvalues_Multiply is used to multiply an MQC -!> integral type variable with an MQC eigenvalues type variable (a diagonal -!> matrix) and returns an MQC integral +!> \brief MQC_Element_Integral_Multiply is used to multiply term-by-term two MQC integral +!> an MQC integral ! !> \par Purpose: ! ============= !> !> \verbatim !> -!> MQC_Integral_Eigenvalues_Multiply is used to multiply an MQC integral type -!> variable with an MQC eigenvalues type variable (a diagonal matrix) and returns -!> an MQC integral. +!> MQC_Element_Integral_Multiply is used to multiply term-by-term two MQC integral type variable +!> and returns an MQC integral. !> !> \endverbatim ! ! Arguments: ! ========== -!> \param[in] IntegralA +!> \param[in] Integral1 !> \verbatim -!> IntegralA is type(mqc_scf_integral) +!> Integral is type(mqc_scf_integral) !> The MQC integral variable to multiply. !> \endverbatim !> -!> \param[in] EigenvaluesB -!> \verbatim -!> EigenvaluesB is type(mqc_scf_eigenvalues) -!> The MQC eigenvalues type variable (a diagonal matrix) -!> to multiply. -!> \endverbatim -!> -!> \param[in] Label +!> \param[in] Integral2 !> \verbatim -!> Label is Character(Len=*),optional -!> Description tag of output MQC integral variable. +!> Integral is type(mqc_scf_integral) +!> The MQC integral variable to multiply. !> \endverbatim ! ! Authors: ! ======== -!> \author L. M. Thompson -!> \date 2018 +!> \author M. M. F. Moraes +!> \date 2025 ! - function mqc_integral_eigenvalues_multiply(integralA,eigenvaluesB,label) result(integralOut) + function mqc_element_integral_multiply(integral1,integral2) result(integralOut) ! implicit none - type(mqc_scf_integral),intent(in)::integralA - type(mqc_scf_eigenvalues),intent(in)::eigenvaluesB - Character(Len=*),optional,intent(in)::label + type(mqc_scf_integral),intent(in)::integral1,integral2 +! Character(Len=*),optional,intent(in)::label type(mqc_scf_integral)::integralOut - type(mqc_matrix)::tmpMatrix1,tmpMatrix2,tmpMatrix3,tmpMatrix4 - type(mqc_vector)::tmpVector1,tmpVector2 + type(mqc_matrix)::tmpMatrixAlpha,tmpMatrixBeta,tmpMatrixAlphaBeta,tmpMatrixBetaAlpha + type(mqc_scalar)::scalar + Character(Len=64)::myLabel +! +! if(present(label)) then +! call string_change_case(label,'l',myLabel) +! else + myLabel = '' +! endIf + + + select case (integral1%array_type) + case('space') + select case (integral2%array_type) + case('space') + tmpMatrixAlpha=MQC_ElementMatrixProduct(integral1%getBlock('alpha'),integral2%getBlock('alpha')) + call mqc_integral_allocate(integralOut,myLabel,'space',tmpMatrixAlpha) + case('spin','general') + tmpMatrixAlpha=MQC_ElementMatrixProduct(integral1%getBlock('alpha'),integral2%getBlock('alpha')) + tmpMatrixBeta=MQC_ElementMatrixProduct(integral1%getBlock('alpha'),integral2%getBlock('beta')) + call mqc_integral_allocate(integralOut,myLabel,'spin',tmpMatrixAlpha,tmpMatrixBeta) + case default + call mqc_error_A('Unknown integral type in mqc_element_integral_multiply', 6, & + 'integral2%array_type', integral2%array_type ) + end select + case('spin') + select case (integral2%array_type) + case('space') + tmpMatrixAlpha=MQC_ElementMatrixProduct(integral1%getBlock('alpha'),integral2%getBlock('alpha')) + tmpMatrixBeta=MQC_ElementMatrixProduct(integral1%getBlock('beta'),integral2%getBlock('alpha')) + call mqc_integral_allocate(integralOut,myLabel,'space',tmpMatrixAlpha) + case('spin','general') + tmpMatrixAlpha=MQC_ElementMatrixProduct(integral1%getBlock('alpha'),integral2%getBlock('alpha')) + tmpMatrixBeta=MQC_ElementMatrixProduct(integral1%getBlock('beta'),integral2%getBlock('beta')) + call mqc_integral_allocate(integralOut,myLabel,'spin',tmpMatrixAlpha,tmpMatrixBeta) + case default + call mqc_error_A('Unknown integral type in mqc_element_integral_multiply', 6, & + 'integral2%array_type', integral2%array_type ) + end select + case('general') + select case (integral2%array_type) + case('space') + tmpMatrixAlpha=MQC_ElementMatrixProduct(integral1%getBlock('alpha'),integral2%getBlock('alpha')) + tmpMatrixBeta=MQC_ElementMatrixProduct(integral1%getBlock('beta'),integral2%getBlock('alpha')) + call mqc_integral_allocate(integralOut,myLabel,'space',tmpMatrixAlpha) + case('spin') + tmpMatrixAlpha=MQC_ElementMatrixProduct(integral1%getBlock('alpha'),integral2%getBlock('alpha')) + tmpMatrixBeta=MQC_ElementMatrixProduct(integral1%getBlock('beta'),integral2%getBlock('beta')) + call mqc_integral_allocate(integralOut,myLabel,'spin',tmpMatrixAlpha,tmpMatrixBeta) + case('general') + tmpMatrixAlpha=MQC_ElementMatrixProduct(integral1%getBlock('alpha'),integral2%getBlock('alpha')) + tmpMatrixBeta=MQC_ElementMatrixProduct(integral1%getBlock('beta'),integral2%getBlock('beta')) + tmpMatrixAlphaBeta=MQC_ElementMatrixProduct(integral1%getBlock('alpha'),integral2%getBlock('beta')) + tmpMatrixBetaAlpha=MQC_ElementMatrixProduct(integral1%getBlock('beta'),integral2%getBlock('alpha')) + call mqc_integral_allocate(integralOut,myLabel,'general',tmpMatrixAlpha,tmpMatrixBeta, & + tmpMatrixAlphaBeta,tmpMatrixBetaAlpha) + case default + call mqc_error_A('Unknown integral type in mqc_element_integral_multiply', 6, & + 'integral2%array_type', integral2%array_type ) + end select + case default + call mqc_error_A('Unknown integral type in mqc_element_integral_multiply', 6, & + 'integral1%array_type', integral1%array_type ) + end select +! + end function mqc_element_integral_multiply +! +! +! PROCEDURE MQC_Integral_Integral_Outer_Product +! +!> \brief MQC_Integral_Integral_Outer_Product is used to calculate the +!> outer product between two MQC integral type variabel and returns an MQC ERIs +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_Integral_Integral_Outer_Product is used to calculate the outer product between +!> two MQC integral and returns an MQC ERIs. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] IntegralA +!> \verbatim +!> IntegralA is type(mqc_scf_integral) +!> \endverbatim +!> +!> \param[in] IntegralB +!> \verbatim +!> IntegralB is type(mqc_scf_integral) +!> \endverbatim +! +! +! Authors: +! ======== +!> \author M. M. F. Moraes +!> \date 2024 +! + function mqc_integral_integral_outer_product(integralA,integralB) result(ERIsOut) + implicit none + class(mqc_scf_integral),intent(in)::integralA + type(mqc_scf_integral),intent(in)::integralB + type(mqc_twoERIs)::ERIsOut + + Character(len=15)::intType + Integer(kind=int64)::NBasis + Type(mqc_matrix)::AlphaA,AlphaB,BetaA,BetaB,AlBeA,AlBeB,BeAlA,BeAlB + Type(mqc_r4tensor)::Alpha,Beta,AlphaBeta,BetaAlpha,aaab,aaba,abaa,baaa,abab, & + baab,baba,abba,bbba,bbab,babb,abbb + + if(integralA%array_type.eq.'general'.or.integralB%array_type.eq.'general') then + intType='general' + elseif(integralA%array_type.eq.'spin'.or.integralB%array_type.eq.'spin') then + intType='spin' + elseif(integralA%array_type.eq.'space'.or.integralB%array_type.eq.'space') then + intType='space' + else + call mqc_error_A('Unknown integral type in mqc_integral_integral_outer_product', 6, & + 'integralA%array_type', integralA%array_type,& + 'integralB%array_type', integralB%array_type ) + end if + if(integralA%array_type.eq.'general'.or.integralA%array_type.eq.'spin') then + NBasis = max(MQC_Matrix_Rows(integralA%Alpha),MQC_Matrix_Rows(integralA%Beta)) + else + NBasis = MQC_Matrix_Rows(integralA%Alpha) + end if + + call Alpha%init(NBasis,NBasis,NBasis,NBasis) + if(intType.ne.'space') then + call Beta%init(NBasis,NBasis,NBasis,NBasis) + call AlphaBeta%init(NBasis,NBasis,NBasis,NBasis) + call BetaAlpha%init(NBasis,NBasis,NBasis,NBasis) + endif + if(intType.ne.'general') then + call aaab%init(NBasis,NBasis,NBasis,NBasis) + call aaba%init(NBasis,NBasis,NBasis,NBasis) + call abaa%init(NBasis,NBasis,NBasis,NBasis) + call baaa%init(NBasis,NBasis,NBasis,NBasis) + call abab%init(NBasis,NBasis,NBasis,NBasis) + call baab%init(NBasis,NBasis,NBasis,NBasis) + call baba%init(NBasis,NBasis,NBasis,NBasis) + call abba%init(NBasis,NBasis,NBasis,NBasis) + call bbba%init(NBasis,NBasis,NBasis,NBasis) + call bbab%init(NBasis,NBasis,NBasis,NBasis) + call babb%init(NBasis,NBasis,NBasis,NBasis) + call abbb%init(NBasis,NBasis,NBasis,NBasis) + endif + + AlphaA=integralA%getBlock('alpha') + AlphaB=integralB%getBlock('alpha') + Alpha = AlphaA.outer.AlphaB + if(intType.ne.'space') then + if(integralA%array_type.eq.'space') then + BetaA=integralA%getBlock('alpha') + else + BetaA=integralA%getBlock('beta') + end if + if(integralB%array_type.eq.'space') then + BetaB=integralB%getBlock('alpha') + else + BetaB=integralB%getBlock('beta') + end if + Beta = BetaA.outer.BetaB + AlphaBeta = AlphaA.outer.BetaB + BetaAlpha = BetaA.outer.AlphaB + endif + if(intType.eq.'general') then + if(integralA%array_type.eq.'general') then + AlBeA=integralA%getBlock('alphabeta') + BeAlA=integralA%getBlock('betaalpha') + abaa = AlBeA.outer.AlphaB + baaa = BeAlA.outer.AlphaB + abbb = AlBeA.outer.BetaB + babb = BeAlA.outer.BetaB + end if + if(integralB%array_type.eq.'general') then + AlBeB=integralB%getBlock('alphabeta') + BeAlB=integralB%getBlock('betaalpha') + aaab = AlphaA.outer.AlBeB + aaba = AlphaA.outer.BeAlB + bbab = BetaA.outer.AlBeB + bbba = BetaA.outer.BeAlB + end if + if(integralA%array_type.eq.'general'.and.integralB%array_type.eq.'general') then + abab = AlBeA.outer.AlBeB + abba = AlBeA.outer.BeAlB + baab = BeAlA.outer.AlBeB + baba = BeAlA.outer.BeAlB + end if + end if + if (intType.eq.'space') then + call mqc_twoeris_allocate(ERIsOut,'full','space',Alpha) + elseif (intType.eq.'spin') then + call mqc_twoeris_allocate(ERIsOut,'full','spin',Alpha,Beta,AlphaBeta,BetaAlpha) + elseif (intType.eq.'general') then + call mqc_twoeris_allocate(ERIsOut,'full','general',Alpha,Beta,AlphaBeta,BetaAlpha,& + abab,abba,baab,baba,aaab,aaba,abaa,baaa,bbba,bbab,babb,abbb) + endif + end function mqc_integral_integral_outer_product + +! +! +! PROCEDURE MQC_Integral_Integral_CrossOuter_Product +! +!> \brief MQC_Integral_Integral_CrossOuter_Product is used to calculate the +!> crossouter product between two MQC integral type variabel and returns an MQC ERIs +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_Integral_Integral_CrossOuter_Product is used to calculate the crossouter product between +!> two MQC integral and returns an MQC ERIs. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] IntegralA +!> \verbatim +!> IntegralA is type(mqc_scf_integral) +!> \endverbatim +!> +!> \param[in] IntegralB +!> \verbatim +!> IntegralB is type(mqc_scf_integral) +!> \endverbatim +! +! +! Authors: +! ======== +!> \author M. M. F. Moraes +!> \date 2024 +! + function mqc_integral_integral_crossouter_product(integralA,integralB) result(ERIsOut) + implicit none + class(mqc_scf_integral),intent(in)::integralA + type(mqc_scf_integral),intent(in)::integralB + type(mqc_twoERIs)::ERIsOut + + Character(len=15)::intType + Integer(kind=int64)::NBasis + Type(mqc_matrix)::AlphaA,AlphaB,BetaA,BetaB,AlBeA,AlBeB,BeAlA,BeAlB + Type(mqc_r4tensor)::Alpha,Beta,AlphaBeta,BetaAlpha,aaab,aaba,abaa,baaa,abab, & + baab,baba,abba,bbba,bbab,babb,abbb + + if(integralA%array_type.eq.'general'.or.integralB%array_type.eq.'general') then + intType='general' + elseif(integralA%array_type.eq.'spin'.or.integralB%array_type.eq.'spin') then + intType='general' + elseif(integralA%array_type.eq.'space'.or.integralB%array_type.eq.'space') then + intType='general' + else + call mqc_error_A('Unknown integral type in mqc_integral_integral_outer_product', 6, & + 'integralA%array_type', integralA%array_type,& + 'integralB%array_type', integralB%array_type ) + end if + if(integralA%array_type.eq.'general'.or.integralA%array_type.eq.'spin') then + NBasis = max(MQC_Matrix_Rows(integralA%Alpha),MQC_Matrix_Rows(integralA%Beta)) + else + NBasis = MQC_Matrix_Rows(integralA%Alpha) + end if + + call Alpha%init(NBasis,NBasis,NBasis,NBasis) + call Beta%init(NBasis,NBasis,NBasis,NBasis) + call AlphaBeta%init(NBasis,NBasis,NBasis,NBasis) + call BetaAlpha%init(NBasis,NBasis,NBasis,NBasis) + call aaab%init(NBasis,NBasis,NBasis,NBasis) + call aaba%init(NBasis,NBasis,NBasis,NBasis) + call abaa%init(NBasis,NBasis,NBasis,NBasis) + call baaa%init(NBasis,NBasis,NBasis,NBasis) + call abab%init(NBasis,NBasis,NBasis,NBasis) + call baab%init(NBasis,NBasis,NBasis,NBasis) + call baba%init(NBasis,NBasis,NBasis,NBasis) + call abba%init(NBasis,NBasis,NBasis,NBasis) + call bbba%init(NBasis,NBasis,NBasis,NBasis) + call bbab%init(NBasis,NBasis,NBasis,NBasis) + call babb%init(NBasis,NBasis,NBasis,NBasis) + call abbb%init(NBasis,NBasis,NBasis,NBasis) + + AlphaA=integralA%getBlock('alpha') + AlphaB=integralB%getBlock('alpha') + Alpha = AlphaA.crossouter.AlphaB + if(integralA%array_type.eq.'space') then + BetaA=integralA%getBlock('alpha') + else + BetaA=integralA%getBlock('beta') + end if + if(integralB%array_type.eq.'space') then + BetaB=integralB%getBlock('alpha') + else + BetaB=integralB%getBlock('beta') + end if + Beta = BetaA.crossouter.BetaB + abba = AlphaA.crossouter.BetaB + baab = BetaA.crossouter.AlphaB + if(integralA%array_type.eq.'general') then + AlBeA=integralA%getBlock('alphabeta') + BeAlA=integralA%getBlock('betaalpha') + aaab = AlBeA.crossouter.AlphaB + baaa = BeAlA.crossouter.AlphaB + abbb = AlBeA.crossouter.BetaB + bbba = BeAlA.crossouter.BetaB + end if + if(integralB%array_type.eq.'general') then + AlBeB=integralB%getBlock('alphabeta') + BeAlB=integralB%getBlock('betaalpha') + abaa = AlphaA.crossouter.AlBeB + aaba = AlphaA.crossouter.BeAlB + bbab = BetaA.crossouter.AlBeB + babb = BetaA.outer.BeAlB + end if + if(integralA%array_type.eq.'general'.and.integralB%array_type.eq.'general') then + abab = AlBeA.crossouter.AlBeB + AlphaBeta = AlBeA.crossouter.BeAlB + BetaAlpha = BeAlA.crossouter.AlBeB + baba = BeAlA.crossouter.BeAlB + end if + call mqc_twoeris_allocate(ERIsOut,'full','general',Alpha,Beta,AlphaBeta,BetaAlpha,& + abab,abba,baab,baba,aaab,aaba,abaa,baaa,bbba,bbab,babb,abbb) + end function mqc_integral_integral_crossouter_product + +! +! +! PROCEDURE MQC_Eigen_Eigen_Outer_Product +! +!> \brief MQC_Eigen_Eigen_Outer_Product is used to calculate the +!> outer product between two MQC SCF eigenvalue type variabel and returns +!> an MQC SCF Integral +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_Eigen_Eigen_Outer_Product is used to calculate the +!> outer product between two MQC SCF eigenvalue type variabel and returns +!> an MQC SCF Integral. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] eigA +!> \verbatim +!> eigA is type(mqc_scf_eigenvalues) +!> \endverbatim +!> +!> \param[in] eigB +!> \verbatim +!> eigB is type(mqc_scf_eigenvalues) +!> \endverbatim +! +! +! Authors: +! ======== +!> \author M. M. F. Moraes +!> \date 2024 +! + function mqc_eigen_eigen_outer_product(eigA,eigB) result(IntOut) + implicit none + class(mqc_scf_eigenvalues),intent(in)::eigA + type(mqc_scf_eigenvalues),intent(in)::eigB + type(mqc_scf_integral)::IntOut + + Integer(kind=int64)::NBasis + Type(mqc_vector)::AlphaA,AlphaB,BetaA,BetaB + Type(mqc_matrix)::aa,bb,ab,ba + + AlphaA=eigA%getBlock('a') + AlphaB=transpose(eigB%getBlock('a')) + aa = AlphaA.outer.AlphaB + select case(eigA%type()) + case('space') + select case(eigB%type()) + case('space') + call mqc_integral_allocate(IntOut,'','general',aa,aa,aa,aa) + case('spin','general') + BetaB=transpose(eigB%getBlock('b')) + ab = AlphaA.outer.BetaB + call mqc_integral_allocate(IntOut,'','general',aa,ab,ab,aa) + case default + call mqc_error_A('Unknown Eigenvalue type in mqc_eigen_eigen_outer_product', 6, & + 'eigB%type', eigB%type() ) + end select + case('spin','general') + select case(eigB%type()) + case('space') + BetaA=eigA%getBlock('b') + ba = BetaA.outer.AlphaB + call mqc_integral_allocate(IntOut,'','general',aa,ba,aa,ba) + case('spin','general') + BetaA=eigA%getBlock('b') + BetaB=transpose(eigB%getBlock('b')) + ab = AlphaA.outer.BetaB + ba = BetaA.outer.AlphaB + bb = BetaA.outer.BetaB + call mqc_integral_allocate(IntOut,'','general',aa,bb,ab,ba) + case default + call mqc_error_A('Unknown Eigenvalue type in mqc_eigen_eigen_outer_product', 6, & + 'eigB%type', eigB%type() ) + end select + case default + call mqc_error_A('Unknown Eigenvalue type in mqc_eigen_eigen_outer_product', 6, & + 'eigA%type', eigA%type() ) + end select + end function mqc_eigen_eigen_outer_product + +! +! +! PROCEDURE MQC_Integral_Eigenvalues_Multiply +! +!> \brief MQC_Integral_Eigenvalues_Multiply is used to multiply an MQC +!> integral type variable with an MQC eigenvalues type variable (a diagonal +!> matrix) and returns an MQC integral +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_Integral_Eigenvalues_Multiply is used to multiply an MQC integral type +!> variable with an MQC eigenvalues type variable (a diagonal matrix) and returns +!> an MQC integral. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] IntegralA +!> \verbatim +!> IntegralA is type(mqc_scf_integral) +!> The MQC integral variable to multiply. +!> \endverbatim +!> +!> \param[in] EigenvaluesB +!> \verbatim +!> EigenvaluesB is type(mqc_scf_eigenvalues) +!> The MQC eigenvalues type variable (a diagonal matrix) +!> to multiply. +!> \endverbatim +!> +!> \param[in] Label +!> \verbatim +!> Label is Character(Len=*),optional +!> Description tag of output MQC integral variable. +!> \endverbatim +! +! Authors: +! ======== +!> \author L. M. Thompson +!> \date 2018 +! + function mqc_integral_eigenvalues_multiply(integralA,eigenvaluesB,label) result(integralOut) +! + implicit none + type(mqc_scf_integral),intent(in)::integralA + type(mqc_scf_eigenvalues),intent(in)::eigenvaluesB + Character(Len=*),optional,intent(in)::label + type(mqc_scf_integral)::integralOut + type(mqc_matrix)::tmpMatrix1,tmpMatrix2,tmpMatrix3,tmpMatrix4 + type(mqc_vector)::tmpVector1,tmpVector2 Character(Len=64)::myLabel ! if(present(label)) then @@ -8230,22 +9266,22 @@ function mqc_eri_r4tensor_contraction(eris,r4tensor,label) result(output) end function mqc_eri_r4tensor_contraction ! ! -! PROCEDURE MQC_ERI_Interaction +! PROCEDURE MQC_ERI_ERI_Contraction ! -!> \brief MQC_ERI_Interaction is used to return the two-electron integrals of -!> the requested two-body interaction type +!> \brief MQC_ERI_ERI_Contraction is used to return the contraction of +!> two-ERIs with an MQC rank-4 tensor ! !> \par Purpose: ! ============= !> !> \verbatim !> -!> MQC_ERI_Interaction is used to return the two-electron integrals of the -!> requested two-body interaction type. The following options are available: +!> MQC_ERI_ERI_Contraction is used to return the contraction of two-ERIs with +!> an MQC two-ERIs. !> -!> 1. Label = 'coulomb' returns the 2ERIs with single bar coulomb integrals (default). -!> 2. Label = 'exchange' returns the 2ERIs with single bar exchange integrals. -!> 3. Label = 'doublebar' returns the 2ERIs with double-bar integrals. +!> 1. Label = 'coulomb' returns the contraction with single bar coulomb integrals. +!> 2. Label = 'exchange' returns the contraction with single bar exchange integrals. +!> 3. Label = 'doublebar' returns the contraction with double-bar integrals (default). !> !> \endverbatim ! @@ -8254,21 +9290,625 @@ end function mqc_eri_r4tensor_contraction !> \param[in] ERIs !> \verbatim !> ERIs is type(mqc_twoERIs) -!> The 2ERIs from which to construct the interaction type. +!> The 2ERIs to be contracted. +!> \endverbatim +!> +!> \param[in] R4Tensor +!> \verbatim +!> R4Tensor is type(mqc_R4Tensor) +!> The MQC rank-4 tensor to be contracted. !> \endverbatim !> !> \param[in] Label !> \verbatim !> Label is character(len=*),optional -!> = 'coulomb': returns coulomb integrals. -!> = 'exchange': returns exchange integrals. -!> = 'doublebar': returns double bar integrals (default). +!> = 'coulomb': contracts with coulomb integrals. +!> = 'exchange': contracts with exchange integrals. +!> = 'doublebar': contracts with double bar integrals +!> (default). !> \endverbatim ! ! Authors: ! ======== -!> \author L. M. Thompson -!> \date 2019 +!> \author M. M. F. Moraes +!> \date 2024 +! + function mqc_eri_eri_contraction(eris1,eris2,label) result(output) +! + implicit none + type(mqc_twoERIs),intent(in)::eris1,eris2 + character(len=*),optional,intent(in)::label + type(mqc_scalar)::output + character(len=64)::mylabel + integer(kind=int64)::r4t1,r4t2,r4t3,r4t4,nB1,nB2,nB3,nB4 + type(mqc_twoERIs)::temp +! integer(kind=int64)::k,l,m,n,nBasis +! character(len=64)::mylabel +! real(kind=real64)::zero=1.0d-10 +! + if(present(label)) then + call string_change_case(label,'L',myLabel) + else + myLabel = 'doublebar' + endIf +! + if (eris1%blockSize('alpha',1).ne.eris2%blockSize('alpha',1).or.& + eris1%blockSize('alpha',2).ne.eris2%blockSize('alpha',2).or.& + eris1%blockSize('alpha',3).ne.eris2%blockSize('alpha',3).or.& + eris1%blockSize('alpha',4).ne.eris2%blockSize('alpha',4))& + call mqc_error_i('2ERIs not correctly dimensioned in mqc_eri_eri_contraction',6,& + 'dim1',eris1%blockSize('alpha',1)-eris2%blockSize('alpha',1),& + 'dim2',eris1%blockSize('alpha',2)-eris2%blockSize('alpha',2),& + 'dim3',eris1%blockSize('alpha',3)-eris2%blockSize('alpha',3),& + 'dim4',eris1%blockSize('alpha',4)-eris2%blockSize('alpha',4)) + + select case (myLabel) + case ('coulomb') + temp = eris1 + case ('exchange') + temp = mqc_eri_interaction(eris1,'exchange') + case ('doublebar') + temp = mqc_eri_interaction(eris1,'doublebar') + case default + call mqc_error_a('Unrecognized 2ERI integral interaction specified',& + 6,'myLabel',myLabel) + end select + + output = 0.0 + select case (eris1%type()) + case ('space','regular') + output = output + contraction(temp%getBlock('aaaa'),eris2%getBlock('aaaa')) + select case (eris2%type()) + case ('space','regular') + if(myLabel.eq.'coulomb') output = 4*output + if(myLabel.eq.'exchange') output = 2*output + if(myLabel.eq.'doublebar') & + output = 2*output + 2*contraction(eris1%getBlock('aaaa'),eris2%getBlock('aaaa')) + case ('spin') + output = output + contraction(temp%getBlock('aaaa'),eris2%getBlock('bbbb')) + if(myLabel.ne.'exchange') then + output = output + contraction(eris1%getBlock('aaaa'),eris2%getBlock('aabb')) + output = output + contraction(eris1%getBlock('aaaa'),eris2%getBlock('bbaa')) + end if + case ('general') + output = output + contraction(temp%getBlock('aaaa'),eris2%getBlock('bbbb')) + if(myLabel.ne.'exchange') then + output = output + contraction(temp%getBlock('aaaa'),eris2%getBlock('aabb')) + output = output + contraction(temp%getBlock('aaaa'),eris2%getBlock('bbaa')) + end if + if(myLabel.ne.'coulomb') then + output = output + contraction(eris1%getBlock('abba'),eris2%getBlock('abba')) + output = output + contraction(eris1%getBlock('abba'),eris2%getBlock('baab')) + end if + end select + case ('spin') + output = output + contraction(temp%getBlock('aaaa'),eris2%getBlock('aaaa')) + select case (eris2%type()) + case ('space','regular') + output = output + contraction(temp%getBlock('bbbb'),eris2%getBlock('aaaa')) + if(myLabel.ne.'exchange') then + output = output + contraction(temp%getBlock('aabb'),eris2%getBlock('aaaa')) + output = output + contraction(temp%getBlock('bbaa'),eris2%getBlock('aaaa')) + end if + case ('spin') + output = output + contraction(temp%getBlock('bbbb'),eris2%getBlock('bbbb')) + if(myLabel.ne.'exchange') then + output = output + contraction(temp%getBlock('aabb'),eris2%getBlock('aabb')) + output = output + contraction(temp%getBlock('bbaa'),eris2%getBlock('bbaa')) + end if + case ('general') + output = output + contraction(temp%getBlock('bbbb'),eris2%getBlock('bbbb')) + if(myLabel.ne.'exchange') then + output = output + contraction(temp%getBlock('aabb'),eris2%getBlock('aabb')) + output = output + contraction(temp%getBlock('bbaa'),eris2%getBlock('bbaa')) + end if + if(myLabel.ne.'coulomb') then + output = output + contraction(temp%getBlock('abba'),eris2%getBlock('abba')) + output = output + contraction(temp%getBlock('baab'),eris2%getBlock('baab')) + end if + end select + case ('general') + output = output + contraction(temp%getBlock('aaaa'),eris2%getBlock('aaaa')) + select case (eris2%type()) + case ('space','regular') + output = output + contraction(temp%getBlock('bbbb'),eris2%getBlock('aaaa')) + output = output + contraction(temp%getBlock('aabb'),eris2%getBlock('aaaa')) + output = output + contraction(temp%getBlock('bbaa'),eris2%getBlock('aaaa')) + case ('spin') + output = output + contraction(temp%getBlock('bbbb'),eris2%getBlock('bbbb')) + output = output + contraction(temp%getBlock('aabb'),eris2%getBlock('aabb')) + output = output + contraction(temp%getBlock('bbaa'),eris2%getBlock('bbaa')) + case ('general') + output = output + contraction(temp%getBlock('bbbb'),eris2%getBlock('bbbb')) + output = output + contraction(temp%getBlock('aabb'),eris2%getBlock('aabb')) + output = output + contraction(temp%getBlock('bbaa'),eris2%getBlock('bbaa')) + output = output + contraction(temp%getBlock('baab'),eris2%getBlock('baab')) + output = output + contraction(temp%getBlock('abba'),eris2%getBlock('abba')) + output = output + contraction(temp%getBlock('baba'),eris2%getBlock('baba')) + output = output + contraction(temp%getBlock('abab'),eris2%getBlock('abab')) + output = output + contraction(temp%getBlock('aaab'),eris2%getBlock('aaab')) + output = output + contraction(temp%getBlock('aaba'),eris2%getBlock('aaba')) + output = output + contraction(temp%getBlock('abaa'),eris2%getBlock('abaa')) + output = output + contraction(temp%getBlock('baaa'),eris2%getBlock('baaa')) + output = output + contraction(temp%getBlock('bbba'),eris2%getBlock('bbba')) + output = output + contraction(temp%getBlock('bbab'),eris2%getBlock('bbab')) + output = output + contraction(temp%getBlock('babb'),eris2%getBlock('babb')) + output = output + contraction(temp%getBlock('abbb'),eris2%getBlock('abbb')) + end select + end select +! + end function mqc_eri_eri_contraction +! +! +! PROCEDURE MQC_ERI_ERI_Partial_Contraction +! +!> \brief MQC_ERI_ERI_Contraction is used to return the contraction of +!> two-ERIs with an MQC rank-4 tensor +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_ERI_ERI_Contraction is used to return the contraction of two-ERIs with +!> an MQC two-ERIs. +!> +!> 1. Label = 'coulomb' returns the contraction with single bar coulomb integrals. +!> 2. Label = 'exchange' returns the contraction with single bar exchange integrals. +!> 3. Label = 'doublebar' returns the contraction with double-bar integrals (default). +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] ERIs +!> \verbatim +!> ERIs is type(mqc_twoERIs) +!> The 2ERIs to be contracted. +!> \endverbatim +!> +!> \param[in] R4Tensor +!> \verbatim +!> R4Tensor is type(mqc_R4Tensor) +!> The MQC rank-4 tensor to be contracted. +!> \endverbatim +!> +!> TODO: finish doc +!> \endverbatim +! +! Authors: +! ======== +!> \author M. M. F. Moraes +!> \date 2024 +! + function mqc_eri_eri_partial_contraction(eris1,eris2,index1,index2) result(IntOut) +! + implicit none + type(mqc_twoERIs),intent(in)::eris1,eris2 + integer(kind=int64),dimension(4),intent(in),optional::index1,index2 + type(mqc_scf_integral)::IntOut + + type(mqc_matrix)::aa,ab,ba,bb + character(len=64)::mylabel,intType + integer(kind=int64)::r4t1,r4t2,r4t3,r4t4,nB1,nB2,nB3,nB4 + type(mqc_twoERIs)::temp + Integer(kind=int64)::ia,ja,ka,la,ib,jb,kb,lb,i,j,k,l + Integer(kind=int64)::dimAi,dimAj,dimAk,dimAl,dimBi,dimBj,dimBk,dimBl + character(len=4)::str1,str2,active_conf + character(len=4),dimension(4)::allowed_spin,a_bbb_array,aaa_b_array + character(len=4),dimension(6)::aa_bb_array + + + aa_bb_array=['aabb','abab','abba','baab','baba','bbaa'] + a_bbb_array=['abbb','babb','bbab','bbba'] + aaa_b_array=['baaa','abaa','aaba','aaab'] + +! + if(present(index1)) then + ia=index1(1) + ja=index1(2) + ka=index1(3) + la=index1(4) + if (ia.eq.ja) Call MQC_Error_I('Duplicated index in MQC_ERI_ERI_Contraction',& + 6,'ia', ia,'ja',ja) + if (ia.eq.ka) Call MQC_Error_I('Duplicated index in MQC_ERI_ERI_Contraction',& + 6,'ia', ia,'ka',ka) + if (ia.eq.la) Call MQC_Error_I('Duplicated index in MQC_ERI_ERI_Contraction',& + 6,'ia', ia,'la',la) + if (ja.eq.ka) Call MQC_Error_I('Duplicated index in MQC_ERI_ERI_Contraction',& + 6,'ja', ja,'ka',ka) + if (ja.eq.la) Call MQC_Error_I('Duplicated index in MQC_ERI_ERI_Contraction',& + 6,'ja', ja,'la',la) + if (ka.eq.la) Call MQC_Error_I('Duplicated index in MQC_ERI_ERI_Contraction',& + 6,'ka', ka,'la',la) + else + ia=1 + ja=2 + ka=3 + la=4 + endif + if(present(index2)) then + ib=index2(1) + jb=index2(2) + kb=index2(3) + lb=index2(4) + if (ib.eq.jb) Call MQC_Error_I('Duplicated index in MQC_ERI_ERI_Partial_Contraction',& + 6,'ib', ib,'jb',jb) + if (ib.eq.kb) Call MQC_Error_I('Duplicated index in MQC_ERI_ERI_Partial_Contraction',& + 6,'ib', ib,'kb',kb) + if (ib.eq.lb) Call MQC_Error_I('Duplicated index in MQC_ERI_ERI_Partial_Contraction',& + 6,'ib', ib,'lb',lb) + if (jb.eq.kb) Call MQC_Error_I('Duplicated index in MQC_ERI_ERI_Partial_Contraction',& + 6,'jb', jb,'kb',kb) + if (jb.eq.lb) Call MQC_Error_I('Duplicated index in MQC_ERI_ERI_Partial_Contraction',& + 6,'jb', jb,'lb',lb) + if (kb.eq.lb) Call MQC_Error_I('Duplicated index in MQC_ERI_ERI_Partial_Contraction',& + 6,'kb', kb,'lb',lb) + else + ib=1 + jb=2 + kb=3 + lb=4 + endif + if (ia.gt.4.or.ia.le.0.or.ib.gt.4.or.ib.le.0.or.& + ja.gt.4.or.ja.le.0.or.jb.gt.4.or.jb.le.0.or.& + ka.gt.4.or.ka.le.0.or.kb.gt.4.or.kb.le.0.or.& + la.gt.4.or.la.le.0.or.lb.gt.4.or.lb.le.0.) Call MQC_Error(& + 'Index out of range in MQC_ERI_ERI_Partial_Contraction', 6) + + if (eris1%blockSize('alpha',ia).ne.eris2%blockSize('alpha',ib).or.& + eris1%blockSize('alpha',ja).ne.eris2%blockSize('alpha',jb).or.& + eris1%blockSize('alpha',ka).ne.eris2%blockSize('alpha',kb).or.& + eris1%blockSize('alpha',la).ne.eris2%blockSize('alpha',lb))& + call mqc_error_i('2ERIs not correctly dimensioned in mqc_eri_eri_contraction',6,& + 'dim1',eris1%blockSize('alpha',ia)-eris2%blockSize('alpha',ib),& + 'dim2',eris1%blockSize('alpha',ja)-eris2%blockSize('alpha',jb),& + 'dim3',eris1%blockSize('alpha',ka)-eris2%blockSize('alpha',kb),& + 'dim4',eris1%blockSize('alpha',la)-eris2%blockSize('alpha',lb)) + + call aa%init(eris1%blockSize('alpha',ka),eris1%blockSize('alpha',la)) + call ab%init(eris1%blockSize('alpha',ka),eris1%blockSize('beta',la)) + call ba%init(eris1%blockSize('beta',ka),eris1%blockSize('alpha',la)) + call bb%init(eris1%blockSize('beta',ka),eris1%blockSize('beta',la)) + + aa = aa + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock('aaaa'),eris2%getBlock('aaaa'),index1,index2) + select case (eris1%type()) + case ('space','regular') + select case (eris2%type()) + case ('space','regular') + intType='space' + case ('spin') + intType='spin' + bb = bb + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock('aaaa'),eris2%getBlock('bbbb'),index1,index2) + do i=1,6 + active_conf=aa_bb_array(i) + str1=mqc_index_to_string([ia,ja,ka,la],active_conf(:)) + str2=mqc_index_to_string([ib,jb,kb,lb],active_conf(:)) + if(str1.eq.'aabb'.or.str1.eq.'bbaa') then + if(str2.eq.'aabb'.or.str1.eq.'bbaa') then + if (active_conf(3:4).eq.'aa') then + aa = aa + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + else if (active_conf(3:4).eq.'bb') then + bb = bb + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + else if (active_conf(3:4).eq.'ab') then + intType='general' + ab = ab + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + else if (active_conf(3:4).eq.'ba') then + intType='general' + ba = ba + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + end if + end if + end if + end do + case ('general') + intType='spin' + bb = bb + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock('aaaa'),eris2%getBlock('bbbb'),index1,index2) + do i=1,6 + active_conf=aa_bb_array(i) + str1=mqc_index_to_string([ia,ja,ka,la],active_conf(:)) + str2=mqc_index_to_string([ib,jb,kb,lb],active_conf(:)) + if(str1.eq.'aabb'.or.str1.eq.'bbaa') then + if (active_conf(3:4).eq.'aa') then + aa = aa + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + else if (active_conf(3:4).eq.'bb') then + bb = bb + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + else if (active_conf(3:4).eq.'ab') then + intType='general' + ab = ab + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + else if (active_conf(3:4).eq.'ba') then + intType='general' + ba = ba + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + end if + end if + end do + case default + call mqc_error_a('Unrecognized 2ERI integral type',6,& + 'eris2%type',eris2%type()) + end select + case ('spin') + select case (eris2%type()) + case ('space','regular','spin') + intType='spin' + bb = bb + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock('bbbb'),eris2%getBlock('bbbb'),index1,index2) + do i=1,6 + active_conf=aa_bb_array(i) + str1=mqc_index_to_string([ia,ja,ka,la],active_conf(:)) + str2=mqc_index_to_string([ib,jb,kb,lb],active_conf(:)) + if(str1.eq.'aabb'.or.str1.eq.'bbaa') then + if(str2.eq.'aabb'.or.str1.eq.'bbaa') then + if (active_conf(3:4).eq.'aa') then + aa = aa + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + else if (active_conf(3:4).eq.'bb') then + bb = bb + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + else if (active_conf(3:4).eq.'ab') then + intType='general' + ab = ab + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + else if (active_conf(3:4).eq.'ba') then + intType='general' + ba = ba + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + end if + end if + end if + end do + case ('general') + intType='spin' + bb = bb + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock('bbbb'),eris2%getBlock('bbbb'),index1,index2) + do i=1,6 + active_conf=aa_bb_array(i) + str1=mqc_index_to_string([ia,ja,ka,la],active_conf(:)) + str2=mqc_index_to_string([ib,jb,kb,lb],active_conf(:)) + if(str1.eq.'aabb'.or.str1.eq.'bbaa') then + if (active_conf(3:4).eq.'aa') then + aa = aa + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + else if (active_conf(3:4).eq.'bb') then + bb = bb + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + else if (active_conf(3:4).eq.'ab') then + intType='general' + ab = ab + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + else if (active_conf(3:4).eq.'ba') then + intType='general' + ba = ba + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + end if + end if + end do + case default + call mqc_error_a('Unrecognized 2ERI integral type',6,& + 'eris2%type',eris2%type()) + end select + case ('general') + select case (eris2%type()) + case ('space','regular','spin') + intType='spin' + bb = bb + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock('bbbb'),eris2%getBlock('bbbb'),index1,index2) + do i=1,6 + active_conf=aa_bb_array(i) + str1=mqc_index_to_string([ia,ja,ka,la],active_conf(:)) + str2=mqc_index_to_string([ib,jb,kb,lb],active_conf(:)) + if(str2.eq.'aabb'.or.str2.eq.'bbaa') then + if (active_conf(3:4).eq.'aa') then + aa = aa + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + else if (active_conf(3:4).eq.'bb') then + bb = bb + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + else if (active_conf(3:4).eq.'ab') then + intType='general' + ab = ab + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + else if (active_conf(3:4).eq.'ba') then + intType='general' + ba = ba + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + end if + end if + end do + case ('general') + intType='spin' + bb = bb + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock('bbbb'),eris2%getBlock('bbbb'),index1,index2) + do i=1,6 + active_conf=aa_bb_array(i) + str1=mqc_index_to_string([ia,ja,ka,la],active_conf(:)) + str2=mqc_index_to_string([ib,jb,kb,lb],active_conf(:)) + if (active_conf(3:4).eq.'aa') then + aa = aa + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + else if (active_conf(3:4).eq.'bb') then + bb = bb + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + else if (active_conf(3:4).eq.'ab') then + intType='general' + ab = ab + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + else if (active_conf(3:4).eq.'ba') then + intType='general' + ba = ba + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + end if + end do + do i=1,4 + active_conf=a_bbb_array(i) + str1=mqc_index_to_string([ia,ja,ka,la],active_conf(:)) + str2=mqc_index_to_string([ib,jb,kb,lb],active_conf(:)) + if (active_conf(3:4) == 'aa') then + aa = aa + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + else if (active_conf(3:4) == 'bb') then + bb = bb + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + else if (active_conf(3:4) == 'ab') then + intType='general' + ab = ab + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + else if (active_conf(3:4) == 'ba') then + intType='general' + ba = ba + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + end if + end do + do i=1,4 + active_conf=aaa_b_array(i) + str1=mqc_index_to_string([ia,ja,ka,la],active_conf) + str2=mqc_index_to_string([ib,jb,kb,lb],active_conf) + if (active_conf(3:4) == 'aa') then + aa = aa + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + else if (active_conf(3:4) == 'bb') then + bb = bb + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + else if (active_conf(3:4) == 'ab') then + intType='general' + ab = ab + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + else if (active_conf(3:4) == 'ba') then + intType='general' + ba = ba + MQC_R4Tensor_R4Tensor_Partial_Contraction(& + eris1%getBlock(str1),eris2%getBlock(str2),index1,index2) + end if + end do + case default + call mqc_error_a('Unrecognized 2ERI integral type',6,& + 'eris2%type',eris2%type()) + end select + case default + call mqc_error_a('Unrecognized 2ERI integral type',6,& + 'eris1%type',eris1%type()) + end select + if (intType.eq.'space') then + call mqc_integral_allocate(IntOut,'full','space',aa) + elseif (intType.eq.'spin') then + call mqc_integral_allocate(IntOut,'full','spin',aa,bb) + elseif (intType.eq.'general') then + call mqc_integral_allocate(IntOut,'full','general',aa,bb,ab,ba) + endif + +! + end function mqc_eri_eri_partial_contraction +! +! +! PROCEDURE MQC_Index_to_String +! +!> \brief MQC_Index_to_String is a support function to +!> MQC_ERI_ERI_Partial_Contraction that returns an alpha/beta string +!> representing the ERI block to be contracted. +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_Index_to_String is used to return an alpha/beta string based on a target +!> and a reference set of indexes. +!> +!> +!> Example: +!> ref: (4,3,1,2) +!> target: (a,b,a,a) +!> ^ ^ +!> to be contracted +!> output: (a,a,b,a) +!> +!> A partial contraction (via MQC_ERI_ERI_Partial_Contraction) with +!> (4,3,1,2) indexes (ref_index). +!> The desired result is the contactions of one alpha and one +!> beta type orbtials (a,b) and return a alpha,alpha matrix (a,a) (target) +!> The ERI block symmetri to be used is 'aaba' (output) +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] ref_index +!> \verbatim +!> ref_index is integer(kind=int64),dimension(4) +!> An arbitrary permutation of (1,2,3,4) +!> \endverbatim +!> +!> \param[in] targ_string +!> \verbatim +!> target_string character(len=4) +!> An ERI alpha/beta string, ex: 'aaaa' or 'abaa' +!> \endverbatim +!> + + function mqc_index_to_string(ref_index,target_string) result(string) + implicit none + integer(kind=int64),dimension(4),intent(in)::ref_index + character(len=4),intent(in)::target_string + character(len=4)::string + integer(kind=int64)::i,refi + + string=' ' + do i=1,4 + refi=ref_index(i) + string=trim(string(1:i))//target_string(refi:refi) + end do + + end function mqc_index_to_string +! +! +! +! PROCEDURE MQC_ERI_Interaction +! +!> \brief MQC_ERI_Interaction is used to return the two-electron integrals of +!> the requested two-body interaction type +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_ERI_Interaction is used to return the two-electron integrals of the +!> requested two-body interaction type. The following options are available: +!> +!> 1. Label = 'coulomb' returns the 2ERIs with single bar coulomb integrals (default). +!> 2. Label = 'exchange' returns the 2ERIs with single bar exchange integrals. +!> 3. Label = 'doublebar' returns the 2ERIs with double-bar integrals. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] ERIs +!> \verbatim +!> ERIs is type(mqc_twoERIs) +!> The 2ERIs from which to construct the interaction type. +!> \endverbatim +!> +!> \param[in] Label +!> \verbatim +!> Label is character(len=*),optional +!> = 'coulomb': returns coulomb integrals. +!> = 'exchange': returns exchange integrals. +!> = 'doublebar': returns double bar integrals (default). +!> \endverbatim +! +! Authors: +! ======== +!> \author L. M. Thompson +!> \date 2019 ! function mqc_eri_interaction(eris,label) result(output) ! @@ -9471,94 +11111,405 @@ subroutine mqc_integral_set_energy_list(integral,elist) if(allocated(integral%Energy_List)) deallocate(integral%Energy_List) - integral%Energy_List = elist -! - end subroutine mqc_integral_set_energy_list + integral%Energy_List = elist +! + end subroutine mqc_integral_set_energy_list +! +! +! PROCEDURE MQC_Integral_Get_Energy_List +! +!> \brief MQC_Integral_Get_Energy_List is a subroutine that returns the energy +!> list object of an MQC integral +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_Integral_Get_Energy_List is a subroutine that returns the energy list object +!> of an MQC integral +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] Integral +!> \verbatim +!> Integral is Class(MQC_SCF_Integral) +!> The MQC integral from which EList will be returned. +!> \endverbatim +! +! Authors: +! ======== +!> \author A. D. Mahler +!> \date 2018 +! + function mqc_integral_get_energy_list(integral) result(elist) +! + implicit none + class(mqc_scf_integral)::integral + integer(kind=int64),dimension(:),allocatable::elist + + if(.not.allocated(integral%Energy_List)) then + call mqc_error('Energy_List is not set in MQC_Integral object') + else + elist = integral%Energy_List + endif +! + end function mqc_integral_get_energy_list +! +! +! PROCEDURE MQC_Integral_Delete_Energy_List +! +!> \brief MQC_Integral_Delete_Energy_List is a subroutine that deletes the +!> energy list object of an MQC integral +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_Integral_Delete_Energy_List is a subroutine that deletes the energy list +!> object of an MQC integral +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] Integral +!> \verbatim +!> Integral is Class(MQC_SCF_Integral) +!> The MQC integral from which EList will be deleted. +!> \endverbatim +! +! Authors: +! ======== +!> \author A. D. Mahler +!> \date 2018 +! + subroutine mqc_integral_delete_energy_list(integral) +! + implicit none + class(mqc_scf_integral)::integral + + if(.not.allocated(integral%Energy_List)) then + call mqc_error('Energy_List was not set in MQC_Integral object') + else + deallocate(integral%Energy_List) + end if +! + end subroutine mqc_integral_delete_energy_list +! +! +! PROCEDURE MQC_Integral_Sum_Elements +! +!> \brief MQC_Integral_Sum_Elements is a function that returns the sum of all elements +!> of an MQC SCF Integral +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_Integral_Sum is a function that returns the sum of all elements of an MQC +!> SCF Integral. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] A +!> \verbatim +!> A is Class(MQC_scf_integral) +!> The name of the MQC_SCF_Integral variable. +!> \endverbatim +! +! Authors: +! ======== +!> \author M. M. F. Moraes +!> \date 2024 +! + function mqc_integral_sum_elements(A) result(output) +! + Implicit None + class(mqc_scf_integral),intent(in)::A + type(mqc_scalar)::output + + output = 0.0 + select case(A%array_type) + case('space') + output = output + 2*sum(A%getBlock('aa')) + case('spin') + output = output + sum(A%getBlock('aa')) + output = output + sum(A%getBlock('bb')) + case('general') + output = output + sum(A%getBlock('aa')) + output = output + sum(A%getBlock('bb')) + output = output + sum(A%getBlock('ab')) + output = output + sum(A%getBlock('ba')) + case default + call mqc_error_a('Unrecognized data type in MQC_Integral_Sum_Elements',6,& + 'integral%array_type',A%array_type) + end select + end function mqc_integral_sum_elements +! +! +! PROCEDURE MQC_Integral_Partial_Sum_Elements +! +!> \brief MQC_Integral_Partial_Sum_Elements is a function that returns +!> the sum of all elements in one index of an MQC SCF Integral +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_Integral_Partial_Sum_Elements is a function that returns +!> the sum of all elements in one index of an MQC SCF Integral. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] A +!> \verbatim +!> A is Class(MQC_scf_integral) +!> The name of the MQC_SCF_Integral variable. +!> \endverbatim +! +! Authors: +! ======== +!> \author M. M. F. Moraes +!> \date 2024 +! + function mqc_integral_partial_sum_elements(A,indexnum) result(output) +! + Implicit None + class(mqc_scf_integral),intent(in)::A + integer(kind=int64),optional::indexnum + type(mqc_matrix)::activeM,activeM2,activeM3,activeM4 + type(mqc_vector)::aa,bb + type(mqc_scf_eigenvalues)::output + integer(kind=int64)::i,j,ind + + ind = 1 + if(present(indexnum)) ind = indexnum + + select case(A%array_type) + case('space') + activeM = A%getBlock('aa') + if (ind.eq.1) then + call aa%init(size(activeM,2)) + do j=1,size(activeM,2) + do i=1,size(activeM,1) + call aa%put(aa%at(j) + activeM%at(i,j) , j) + endDo + endDo + else if (ind.eq.2) then + call aa%init(size(activeM,1)) + do j=1,size(activeM,1) + do i=1,size(activeM,2) + call aa%put(aa%at(j) + activeM%at(j,i) , j) + endDo + endDo + else + Call MQC_Error_I('Index out of range MQC_Integral_Partial_Sum_Element', & + 6,'index', ind ) + end if + call mqc_eigenvalues_allocate(output,'','space',aa) + case('spin') + activeM = A%getBlock('aa') + activeM2 = A%getBlock('bb') + if (ind.eq.1) then + call aa%init(size(activeM,2)) + do j=1,size(activeM,2) + do i=1,size(activeM,1) + call aa%put(aa%at(j) + activeM%at(i,j) , j) + endDo + endDo + call bb%init(size(activeM2,2)) + do j=1,size(activeM2,2) + do i=1,size(activeM2,1) + call bb%put(bb%at(j) + activeM2%at(i,j) , j) + endDo + endDo + else if (ind.eq.2) then + call aa%init(size(activeM,1)) + do j=1,size(activeM,1) + do i=1,size(activeM,2) + call aa%put(aa%at(j) + activeM%at(j,i) , j) + endDo + endDo + call bb%init(size(activeM2,1)) + do j=1,size(activeM2,1) + do i=1,size(activeM2,2) + call bb%put(bb%at(j) + activeM2%at(j,i) , j) + endDo + endDo + else + Call MQC_Error_I('Index out of range MQC_Integral_Partial_Sum_Element', & + 6,'index', ind ) + end if + call mqc_eigenvalues_allocate(output,'','spin',aa,bb) + case('general') + activeM = A%getBlock('aa') + activeM2 = A%getBlock('bb') + activeM3 = A%getBlock('ab') + activeM4 = A%getBlock('ba') + if (ind.eq.1) then + call aa%init(size(activeM,2)) + call bb%init(size(activeM2,2)) + do j=1,size(activeM,2) + do i=1,size(activeM,1) + call aa%put(aa%at(j) + activeM%at(i,j) , j) + endDo + do i=1,size(activeM4,1) + call aa%put(aa%at(j) + activeM4%at(i,j) , j) + endDo + endDo + do j=1,size(activeM2,2) + do i=1,size(activeM2,1) + call bb%put(bb%at(j) + activeM2%at(i,j) , j) + endDo + do i=1,size(activeM3,1) + call bb%put(bb%at(j) + activeM3%at(i,j) , j) + endDo + endDo + else if (ind.eq.2) then + call aa%init(size(activeM,1)) + call bb%init(size(activeM2,1)) + do j=1,size(activeM,1) + do i=1,size(activeM,2) + call aa%put(aa%at(j) + activeM%at(j,i) , j) + endDo + do i=1,size(activeM3,2) + call aa%put(aa%at(j) + activeM3%at(j,i) , j) + endDo + endDo + do j=1,size(activeM2,1) + do i=1,size(activeM2,2) + call bb%put(bb%at(j) + activeM2%at(j,i) , j) + endDo + do i=1,size(activeM4,2) + call bb%put(bb%at(j) + activeM4%at(j,i) , j) + endDo + endDo + endif + call mqc_eigenvalues_allocate(output,'','spin',aa,bb) + case default + call mqc_error_a('Unrecognized data type in MQC_Integral_Partial_Sum_Elements',6,& + 'integral%array_type',A%array_type) + end select + end function mqc_integral_partial_sum_elements + ! ! -! PROCEDURE MQC_Integral_Get_Energy_List +! PROCEDURE MQC_Integral_Multiply ! -!> \brief MQC_Integral_Get_Energy_List is a subroutine that returns the energy -!> list object of an MQC integral +!> \brief MQC_Integral_Multiply is a function that returns +!> the product of each elements between two MQC SCF Integral ! !> \par Purpose: ! ============= !> !> \verbatim !> -!> MQC_Integral_Get_Energy_List is a subroutine that returns the energy list object -!> of an MQC integral +!> MQC_Integral_Multiply is a function that returns the product of each elements +!> between two MQC SCF Integral. !> !> \endverbatim ! ! Arguments: ! ========== -!> \param[in] Integral -!> \verbatim -!> Integral is Class(MQC_SCF_Integral) -!> The MQC integral from which EList will be returned. -!> \endverbatim -! -! Authors: -! ======== -!> \author A. D. Mahler -!> \date 2018 -! - function mqc_integral_get_energy_list(integral) result(elist) -! - implicit none - class(mqc_scf_integral)::integral - integer(kind=int64),dimension(:),allocatable::elist - - if(.not.allocated(integral%Energy_List)) then - call mqc_error('Energy_List is not set in MQC_Integral object') - else - elist = integral%Energy_List - endif -! - end function mqc_integral_get_energy_list -! -! -! PROCEDURE MQC_Integral_Delete_Energy_List -! -!> \brief MQC_Integral_Delete_Energy_List is a subroutine that deletes the -!> energy list object of an MQC integral -! -!> \par Purpose: -! ============= -!> +!> \param[in] A !> \verbatim -!> -!> MQC_Integral_Delete_Energy_List is a subroutine that deletes the energy list -!> object of an MQC integral -!> +!> A is Class(MQC_scf_integral) +!> The name of the MQC_SCF_Integral variable. !> \endverbatim ! -! Arguments: -! ========== -!> \param[in] Integral +!> \param[in] B !> \verbatim -!> Integral is Class(MQC_SCF_Integral) -!> The MQC integral from which EList will be deleted. +!> B is Class(MQC_scf_integral) +!> The name of the MQC_SCF_Integral variable. !> \endverbatim ! ! Authors: ! ======== -!> \author A. D. Mahler -!> \date 2018 +!> \author M. M. F. Moraes +!> \date 2024 ! - subroutine mqc_integral_delete_energy_list(integral) + function mqc_integral_multiply(A,B) result(integralOut) ! - implicit none - class(mqc_scf_integral)::integral + Implicit None + type(mqc_scf_integral),intent(in)::A,B + type(mqc_scf_integral)::integralOut - if(.not.allocated(integral%Energy_List)) then - call mqc_error('Energy_List was not set in MQC_Integral object') - else - deallocate(integral%Energy_List) - end if -! - end subroutine mqc_integral_delete_energy_list + type(mqc_matrix)::Aaa,Abb,Oaa,Obb,Oab,Oba + Character(Len=64)::OutType + + select case(A%array_type) + case('space') + Aaa = A%getBlock('aa') + select case(B%array_type) + case('space') + Oaa = B%getBlock('aa')*Aaa + OutType = 'space' + case('spin','general') + Oaa = B%getBlock('aa')*Aaa + Obb = B%getBlock('bb')*Aaa + OutType = 'spin' + case default + call mqc_error_a('Unrecognized data type in MQC_Integral_Multiply',6,& + 'integralB%array_type',B%array_type) + end select + case('spin') + Aaa = A%getBlock('aa') + Abb = A%getBlock('bb') + select case(B%array_type) + case('space','spin','general') + Oaa = B%getBlock('aa')*Aaa + Obb = B%getBlock('bb')*Abb + OutType = 'spin' + case default + call mqc_error_a('Unrecognized data type in MQC_Integral_Multiply',6,& + 'integralB%array_type',B%array_type) + end select + case('general') + Aaa = A%getBlock('aa') + Abb = A%getBlock('bb') + select case(B%array_type) + case('space','spin') + Oaa = B%getBlock('aa')*Aaa + Obb = B%getBlock('bb')*Abb + OutType = 'spin' + case('general') + Oaa = B%getBlock('aa')*Aaa + Obb = B%getBlock('bb')*Abb + Oab = B%getBlock('ab')*A%getBlock('ab') + Oba = B%getBlock('ba')*A%getBlock('ba') + OutType = 'general' + case default + call mqc_error_a('Unrecognized data type in MQC_Integral_Multiply',6,& + 'integralB%array_type',B%array_type) + end select + case default + call mqc_error_a('Unrecognized data type in MQC_Integral_Multiply',6,& + 'integralA%array_type',A%array_type) + end select + + select case(OutType) + case('space') + call mqc_integral_allocate(integralOut,'','space',Oaa) + case('spin') + call mqc_integral_allocate(integralOut,'','spin',Oaa,Obb) + case('general') + call mqc_integral_allocate(integralOut,'','general',Oaa,Obb,Oab,Oba) + case default + call mqc_error_a('Unrecognized data type in MQC_Integral_Multiply',6,& + 'OuType',OutType) + end select + + end function mqc_integral_multiply ! ! ! PROCEDURE MQC_SCF_Eigenvalues_Power @@ -10854,6 +12805,175 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen endIf end function mqc_twoERIs_At + +! +! +! PROCEDURE MQC_Integral_Put +! +!> \brief MQC_Integral_Put is a function that updates the value of an element +!> of a MQC integral variable +!> +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_Integral_Put is a function that updates the value of (I,J)th element of a +!> MQC integral variable as an MQC scalar. If I or J is negative, the (N-I+1)th +!> index is selected. Optional argument SpinBlock specifies the spin block of MQC +!> integrals. If I or J 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. +!> 3. SpinBlock = 'alphaBeta' uses the alpha-beta spin block. +!> 4. SpinBlock = 'betaAlpha' uses the beta-alpha spin block. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] Integral +!> \verbatim +!> Integral is Class(MQC_SCF_Integral) +!> The MQC integral to return the value of the (I,J)th +!> element. +!> \endverbatim +!> +!> \param[in] I +!> \verbatim +!> I is Integer(kind=int64) +!> The first dimension of the element in MQC integral +!> If I>0 row count is from first index +!> If I<0 row count is from last index. +!> \endverbatim +!> +!> \param[in] J +!> \verbatim +!> J is Integer(kind=int64) +!> The second dimension of the element in MQC integral +!> If J>0 row count is from first index +!> If J<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 +!> = 'alphaBeta': alpha-beta spin block +!> = 'betaAlpha': beta-alpha spin block. +!> \endverbatim +! +! Authors: +! ======== +!> \author M. M. F. Moraes +!> \date 2024 +! + subroutine mqc_integral_put(integral,element,i,j,spinBlockIn) +! + implicit none + class(mqc_scf_integral)::integral + class(*),intent(in)::i,j + integer(kind=int64)::myI,myJ + 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 + select type (j) + type is (integer) + myJ = j + type is (mqc_scalar) + myJ = j + class default + call mqc_error('index type not defined in MQC_twoERIs_put') + end select + + if(myI.gt.integral%blockSize('alpha',1)+integral%blockSize('beta',1)) & + call mqc_error('Requested I dimension out of bounds in mqc_integral_put') + if(myJ.gt.integral%blockSize('alpha',2)+integral%blockSize('beta',2)) & + call mqc_error('Requested J dimension out of bounds in mqc_integral_put') + + 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(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) + 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) + endIf + elseIf(integral%array_type.eq.'general') 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.eq.'alphabeta') then + if(integral%hasAlphaBeta()) call integral%alphaBeta%put(element,myI,myJ) + elseif(spinblock.eq.'betaalpha') then + if(integral%hasBetaAlpha()) call integral%betaAlpha%put(element,myI,myJ) + 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_Integral_Put',6,'Array_Type',integral%array_type) + endIf + else + if(integral%array_type.eq.'space') then + if(myI.le.integral%blockSize('alpha',1).and.myJ.le.integral%blockSize('alpha',2)) then + if(integral%hasAlpha()) call integral%alpha%put(element,myI,myJ) + elseIf(myI.gt.integral%blockSize('alpha',1).and.myJ.gt.integral%blockSize('alpha',2)) then + myI = myI - integral%blockSize('alpha',1) + myJ = myJ - integral%blockSize('alpha',2) + if(integral%hasAlpha()) call integral%alpha%put(element,myI,myJ) + endIf + elseIf(integral%array_type.eq.'spin') then + if(myI.le.integral%blockSize('alpha',1).and.myJ.le.integral%blockSize('alpha',2)) then + if(integral%hasAlpha()) call integral%alpha%put(element,myI,myJ) + elseIf(myI.gt.integral%blockSize('alpha',1).and.myJ.gt.integral%blockSize('alpha',2)) then + myI = myI - integral%blockSize('alpha',1) + myJ = myJ - integral%blockSize('alpha',2) + if(integral%hasBeta()) call integral%beta%put(element,myI,myJ) + endIf + elseIf(integral%array_type.eq.'general') then + if(myI.le.integral%blockSize('alpha',1).and.myJ.le.integral%blockSize('alpha',2)) then + if(integral%hasAlpha()) call integral%alpha%put(element,myI,myJ) + elseIf(myI.le.integral%blockSize('alpha',1).and.myJ.gt.integral%blockSize('alpha',2)) then + myJ = myJ - integral%blockSize('alpha',2) + if(integral%hasAlphaBeta()) call integral%alphaBeta%put(element,myI,myJ) + elseIf(myI.gt.integral%blockSize('alpha',1).and.myJ.le.integral%blockSize('alpha',2)) then + myI = myI - integral%blockSize('alpha',1) + if(integral%hasBetaAlpha()) call integral%betaAlpha%put(element,myI,myJ) + elseIf(myI.gt.integral%blockSize('alpha',1).and.myJ.gt.integral%blockSize('alpha',2)) then + myI = myI - integral%blockSize('alpha',1) + myJ = myJ - integral%blockSize('alpha',2) + if(integral%hasBeta()) call integral%beta%put(element,myI,myJ) + endIf + else + call mqc_error_A('Array type not valid in MQC_Integral_Put',6,'Array_Type',integral%array_type) + endIf + endIf + + end subroutine mqc_integral_put + + ! ! ! PROCEDURE MQC_Integral_At @@ -12588,7 +14708,7 @@ Function Slater_Condon(IOut,IPrint,NBasisIn,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' ! ! This section computes the Coulomb contribution to the matrix ! element @@ -12603,8 +14723,7 @@ Function Slater_Condon(IOut,IPrint,NBasisIn,Alpha_String_1,Beta_String_1, & ! ! This section computes the value of the matrix element ! - If(TwoPartInts%type().eq.'space'.or.(Spin(1).eq.0.and.Spin(2).eq.0.and.& - Spin(3).eq.0.and.Spin(4).eq.0)) then + If(Spin(1).eq.0.and.Spin(2).eq.0.and.Spin(3).eq.0.and.Spin(4).eq.0) then If(Det(1).eq.1) ERI1 = TwoPartInts%at(Orbs(1),Orbs(3),Orbs(2),Orbs(4),'alpha') If(Det(1).eq.2) ERI1 = TwoPartInts%at(Orbs(3),Orbs(1),Orbs(4),Orbs(2),'alpha') ElseIf(Spin(1).eq.0.and.Spin(2).eq.1.and.Spin(3).eq.0.and.Spin(4).eq.1) then @@ -12669,8 +14788,7 @@ Function Slater_Condon(IOut,IPrint,NBasisIn,Alpha_String_1,Beta_String_1, & ERI1 = Zero EndIf - If(TwoPartInts%type().eq.'space'.or.(Spin(1).eq.0.and.Spin(2).eq.0.and. & - Spin(3).eq.0.and.Spin(4).eq.0)) then + If(Spin(1).eq.0.and.Spin(2).eq.0.and.Spin(3).eq.0.and.Spin(4).eq.0) then If(Det(1).eq.1) ERI2 = TwoPartInts%at(Orbs(1),Orbs(4),Orbs(2),Orbs(3),'alpha') If(Det(1).eq.2) ERI2 = TwoPartInts%at(Orbs(3),Orbs(2),Orbs(4),Orbs(1),'alpha') ElseIf(Spin(1).eq.0.and.Spin(2).eq.1.and.Spin(3).eq.1.and.Spin(4).eq.0) then @@ -12741,8 +14859,7 @@ Function Slater_Condon(IOut,IPrint,NBasisIn,Alpha_String_1,Beta_String_1, & If(Spin(1).gt.Spin(3)) ISgn = ISgn+1 If(Spin(2).gt.Spin(4)) ISgn = ISgn+1 - If(TwoPartInts%type().eq.'space'.or.(Spin(1).eq.0.and.Spin(2).eq.0.and.& - Spin(3).eq.0.and.Spin(4).eq.0)) then + If((Spin(1).eq.0.and.Spin(2).eq.0.and.Spin(3).eq.0.and.Spin(4).eq.0)) then If(Det(1).eq.1) ERI1 = TwoPartInts%at(Orbs(1),Orbs(2),Orbs(3),Orbs(4),'alpha') If(Det(1).eq.2) ERI1 = TwoPartInts%at(Orbs(2),Orbs(1),Orbs(4),Orbs(3),'alpha') ElseIf(Spin(1).eq.0.and.Spin(2).eq.0.and.Spin(3).eq.1.and.Spin(4).eq.1) then @@ -12807,8 +14924,7 @@ Function Slater_Condon(IOut,IPrint,NBasisIn,Alpha_String_1,Beta_String_1, & ERI1 = Zero EndIf - If(TwoPartInts%type().eq.'space'.or.(Spin(1).eq.0.and.Spin(2).eq.0.and.& - Spin(3).eq.0.and.Spin(4).eq.0)) then + If(Spin(1).eq.0.and.Spin(2).eq.0.and.Spin(3).eq.0.and.Spin(4).eq.0) then If(Det(1).eq.1) ERI2 = TwoPartInts%at(Orbs(1),Orbs(4),Orbs(3),Orbs(2),'alpha') If(Det(1).eq.2) ERI2 = TwoPartInts%at(Orbs(2),Orbs(3),Orbs(4),Orbs(1),'alpha') ElseIf(Spin(1).eq.0.and.Spin(2).eq.1.and.Spin(3).eq.1.and.Spin(4).eq.0) then @@ -12879,8 +14995,7 @@ Function Slater_Condon(IOut,IPrint,NBasisIn,Alpha_String_1,Beta_String_1, & If(Spin(1).gt.Spin(4)) ISgn = ISgn+1 If(Spin(2).gt.Spin(3)) ISgn = ISgn+1 - If(TwoPartInts%type().eq.'space'.or.(Spin(1).eq.0.and.Spin(2).eq.0.and. & - Spin(3).eq.0.and.Spin(4).eq.0)) then + If(Spin(1).eq.0.and.Spin(2).eq.0.and.Spin(3).eq.0.and.Spin(4).eq.0) then If(Det(1).eq.1) ERI1 = TwoPartInts%at(Orbs(1),Orbs(2),Orbs(4),Orbs(3),'alpha') If(Det(1).eq.2) ERI1 = TwoPartInts%at(Orbs(2),Orbs(1),Orbs(3),Orbs(4),'alpha') ElseIf(Spin(1).eq.0.and.Spin(2).eq.0.and.Spin(3).eq.1.and.Spin(4).eq.1) then @@ -12945,8 +15060,7 @@ Function Slater_Condon(IOut,IPrint,NBasisIn,Alpha_String_1,Beta_String_1, & ERI1 = Zero EndIf - If(TwoPartInts%type().eq.'space'.or.(Spin(1).eq.0.and.Spin(2).eq.0.and.& - Spin(3).eq.0.and.Spin(4).eq.0)) then + If(Spin(1).eq.0.and.Spin(2).eq.0.and.Spin(3).eq.0.and.Spin(4).eq.0) then If(Det(1).eq.1) ERI2 = TwoPartInts%at(Orbs(1),Orbs(3),Orbs(4),Orbs(2),'alpha') If(Det(1).eq.2) ERI2 = TwoPartInts%at(Orbs(2),Orbs(4),Orbs(3),Orbs(1),'alpha') ElseIf(Spin(1).eq.0.and.Spin(2).eq.1.and.Spin(3).eq.0.and.Spin(4).eq.1) then @@ -16282,6 +18396,53 @@ Subroutine TwoERI_Trans(IOut,IPrint,MO_Coeff,ERIs,MO_ERIs,Right_MOs) ! End Subroutine TwoERI_Trans +! +! PROCEDURE MQC_Integral_Outer +! +!> \brief MQC_Print_Integral is a subroutine used to print an MQC +!> EST SCF integral object +!> +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> lines before or after output. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] Integral +!> \endverbatim +! +! Authors: +! ======== +!> \author M. M. F. Moraes +!> \date 2024 +! + function mqc_integral_outer(integral,integral2,cross) result(R4out) +! + implicit none + class(mqc_scf_integral),intent(in)::integral,integral2 + logical::cross + type(mqc_r4tensor)::R4out + + type(mqc_matrix)::tmp1,tmp2 + +! +!! TODO:check if the sizes are consistent + tmp1=integral%getBlock('full') + tmp2=integral2%getBlock('full') + if(cross) then + R4out=tmp1.crossouter.tmp2 + else + R4out=tmp1.outer.tmp2 + end if +! + end function mqc_integral_outer +! +! ! !===================================================================== From e1f2b55ad62d634fa25890cc8e7ec3a00872b31c Mon Sep 17 00:00:00 2001 From: Lee Thompson Date: Sun, 18 May 2025 17:46:36 -0400 Subject: [PATCH 2/9] Minor bug fixes --- src/mqc_algebra.F03 | 107 +++++++++++++++++++++----------------------- src/mqc_est.F03 | 4 +- 2 files changed, 53 insertions(+), 58 deletions(-) diff --git a/src/mqc_algebra.F03 b/src/mqc_algebra.F03 index 31ebc69f..7b5c4cac 100644 --- a/src/mqc_algebra.F03 +++ b/src/mqc_algebra.F03 @@ -26838,71 +26838,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 diff --git a/src/mqc_est.F03 b/src/mqc_est.F03 index b01b25f5..efa87497 100644 --- a/src/mqc_est.F03 +++ b/src/mqc_est.F03 @@ -11702,8 +11702,8 @@ subroutine mqc_scf_integral_generalized_eigensystem(integralA,integralB,eVals,rE call tmpMat1%eigensys(eigenvals=tmpVec1,reigenvecs=tmpMat3,leigenvecs=tmpMat4) call tmpMat2%eigensys(eigenvals=tmpVec2,reigenvecs=tmpMat5,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) From 133f3e50d8bb7ab7ddfdc797450e90e67f705456 Mon Sep 17 00:00:00 2001 From: Lee Thompson Date: Mon, 19 May 2025 16:56:52 -0400 Subject: [PATCH 3/9] Added eigenvalues put routine (via Matheus) --- src/mqc_est.F03 | 126 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 126 insertions(+) diff --git a/src/mqc_est.F03 b/src/mqc_est.F03 index efa87497..bc7bd66d 100644 --- a/src/mqc_est.F03 +++ b/src/mqc_est.F03 @@ -125,7 +125,9 @@ 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 @@ -179,6 +181,9 @@ Module MQC_EST !> \brief Returns the value of a specified element in an MQC_SCF_Eigenvalues !> object Procedure, Public::at => mqc_eigenvalues_at +!> \brief Updates the value of a specified element in an MQC_SCF_Eigenvalues +!> object + Procedure, Public::put => mqc_eigenvalues_put !> \brief Outer Product of two variables into a SCF Integral Procedure, Public::outer => mqc_eigen_eigen_outer_product !> \brief Returns the minimum value in MQC_SCF_Eigenvalues @@ -14612,6 +14617,127 @@ function mqc_eigenvalues_at(eigenvalues,i,spinBlockIn) result(element) end function mqc_eigenvalues_at ! ! +! PROCEDURE MQC_Eigenvalues_Put +! +!> \brief MQC_Eigenvalues_Put is a function that updates the value of an element +!> of a MQC integral variable +!> +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> 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] 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 Type(mqc_scalar) +!> The value to put in the integral. +!> \endverbatim +!> +!> \param[in] I +!> \verbatim +!> 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 M. M. F. Moraes +!> \date 2025 +! + subroutine mqc_eigenvalues_put(eigenvalues,element,i,spinBlockIn) +! + implicit none + 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 + end subroutine mqc_eigenvalues_put +! +! ! PROCEDURE MQC_SCF_Transformation_Matrix ! !> \brief MQC_SCF_Transformation_Matrix is a subroutine that returns the From 5383e05491809b60e3ef39df3b0a5d44424500a6 Mon Sep 17 00:00:00 2001 From: Matheus Moraes Date: Wed, 21 May 2025 10:30:01 -0400 Subject: [PATCH 4/9] eigenvector%put and minor fixes in eigensys --- src/mqc_algebra.F03 | 88 ++++++++++++++++++++++- src/mqc_est.F03 | 172 +++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 257 insertions(+), 3 deletions(-) diff --git a/src/mqc_algebra.F03 b/src/mqc_algebra.F03 index 97c270bb..1c258646 100644 --- a/src/mqc_algebra.F03 +++ b/src/mqc_algebra.F03 @@ -264,7 +264,9 @@ Module MQC_Algebra !> \brief Updates the specified submatrix of the MQC Matrix to the specified !> matrix Procedure, Public::mput => MQC_Matrix_Matrix_Put -!> \brief Updates the specified element of the MQC Matrix to the specified +!> \brief Increases a MQC Matrix dimension and fill the additional elements with +!> with zeros. + Procedure, Public::mpad => MQC_Matrix_Pad !> value End Type MQC_Matrix ! @@ -15714,6 +15716,90 @@ Subroutine MQC_Matrix_SymmMatrix_Put_Complex(mat,symmMatrixIn,label) end subroutine MQC_Matrix_SymmMatrix_Put_Complex ! ! +! PROCEDURE MQC_Matrix_Pad +! +!> \brief MQC_Matrix_Pad is a function that increases the dimension of a +!> MQC matrix and set the new elements to zero. +!> +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_Matrix_Pad is a function that increases the dimension of a +!> MQC matrix and set the new elements to zero. Each dimension of the +!> final MQC matrix cannot be smaller that the input one. +!> To reduce the dimension of a MQC matrix use MQC_Matrix_Matrix_At. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] Matrix +!> \verbatim +!> Mat is Class(MQC_Matrix) +!> The MQC matrix from which the dimension will be increased. +!> \endverbatim +!> +!> \param[in] Rows +!> \verbatim +!> Rows is Integer(kind=int64) +!> The specification of the number of rows after padding. +!> The final number of rows must be larger or equal to the initial. +!> Otherwise the rows of the final matrix is the same as the initial +!> one. To reduce the dimension of a MQC matrix use MQC_Matrix_Matrix_At. +!> If Rows > Matrix%NRow +!> The MatrixOut%Row = Rows +!> else +!> The MatrixOut%Row = Matrix%Row +!> \endverbatim +!> +!> \param[in] Cols +!> \verbatim +!> Cols is Integer(kind=int64) +!> The specification of the number of columns after padding. +!> The final number of columns must be larger or equal to the initial. +!> Otherwise the columns of the final matrix is the same as the initial +!> one. To reduce the dimension of a MQC matrix use MQC_Matrix_Matrix_At. +!> If Cols > Matrix%NCol +!> The MatrixOut%Col = Cols +!> else +!> The MatrixOut%Col = Matrix%Col +!> \endverbatim +! +! Authors: +! ======== +!> \author M. M. F. Moraes +!> \date 2025 +! + function MQC_Matrix_Pad(Matrix,Rows,Cols) result(MatrixOut) +! + Implicit None + Class(MQC_Matrix),Intent(In)::Matrix + Integer(kind=int64),Intent(In)::Rows,Cols + Type(MQC_Matrix)::MatrixOut + + if (Rows.le.Matrix%NRow) then + if (Cols.le.Matrix%NCol) then + MatrixOut = Matrix + else + call MatrixOut%init(Matrix%NRow,Cols) + call MatrixOut%mput(Matrix,[1,Matrix%NRow],[1,Matrix%NCol]) + end if + else + if (Cols.le.Matrix%NCol) then + call MatrixOut%init(Rows,Matrix%NCol) + call MatrixOut%mput(Matrix,[1,Matrix%NRow],[1,Matrix%NCol]) + else + call MatrixOut%init(Rows,Cols) + call MatrixOut%mput(Matrix,[1,Matrix%NRow],[1,Matrix%NCol]) + end if + end if +! + end function MQC_Matrix_Pad +! +! +! ! PROCEDURE MQC_Matrix_Matrix_Put ! !> \brief MQC_Matrix_Matrix_Put is a subroutine that writes a submatrix to the diff --git a/src/mqc_est.F03 b/src/mqc_est.F03 index d74dbd16..aee669ff 100644 --- a/src/mqc_est.F03 +++ b/src/mqc_est.F03 @@ -130,6 +130,7 @@ Module MQC_EST 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 @@ -6176,6 +6177,173 @@ function mqc_matrix_integral_multiply(matrixA,integralB,label) result(integralOu end function mqc_matrix_integral_multiply ! ! +! PROCEDURE MQC_Integral_Pad +! +!> \brief MQC_Integral_Pad is used to increase the dimension +!> of a MQC integral variable +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> 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] 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 +!> \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 M. M. F. Moraes +!> \date 2025 +! + function mqc_integral_pad(integral,axis1,axis2,blockname) result(integralOut) +! + 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 +! + 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 @@ -11446,8 +11614,8 @@ subroutine mqc_scf_integral_generalized_eigensystem(integralA,integralB,eVals,rE call tmpMat1%eigensys(eigenvals=tmpVec1,reigenvecs=tmpMat3,leigenvecs=tmpMat4) call tmpMat2%eigensys(eigenvals=tmpVec2,reigenvecs=tmpMat5,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) From a418f4d203cb1687b05967a33b949458e7dd0318 Mon Sep 17 00:00:00 2001 From: Matheus Moraes Date: Tue, 10 Jun 2025 16:32:09 -0400 Subject: [PATCH 5/9] Eigenvalues and Integral related updates Eigenvalues: Init and push functions implemented Integrals: Allows multiplication of nx0 type integrals --- src/mqc_algebra.F03 | 145 ++++++++++++++++++ src/mqc_est.F03 | 357 +++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 494 insertions(+), 8 deletions(-) diff --git a/src/mqc_algebra.F03 b/src/mqc_algebra.F03 index e4dccd93..022ad2f4 100644 --- a/src/mqc_algebra.F03 +++ b/src/mqc_algebra.F03 @@ -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 @@ -8192,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 @@ -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, & diff --git a/src/mqc_est.F03 b/src/mqc_est.F03 index b4a1f213..795f90c6 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 @@ -169,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 @@ -180,6 +184,9 @@ Module MQC_EST !> \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 @@ -1771,6 +1778,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 @@ -1994,6 +2101,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 @@ -7461,9 +7656,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') @@ -11709,8 +11916,8 @@ 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','spin',tmpMat3,tmpMat4) If(present(lEVecs)) call mqc_integral_allocate(lEVecs,'eigenvectors','spin',tmpMat5,tmpMat6) @@ -14237,18 +14444,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 @@ -14742,6 +14961,128 @@ subroutine mqc_eigenvalues_put(eigenvalues,element,i,spinBlockIn) endIf end subroutine mqc_eigenvalues_put +! +! PROCEDURE MQC_Eigenvalues_Push +! +!> \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 From 1102aab9ed7badd7b97fc67cd5ee8cf309eade00 Mon Sep 17 00:00:00 2001 From: Lee Thompson Date: Sat, 28 Jun 2025 13:49:14 -0400 Subject: [PATCH 6/9] Moved sections of MQC_Build_CI_Operator to subroutines. --- src/mqc_algebra.F03 | 4 +- src/mqc_est.F03 | 933 +++++++++++++++++++++++++++++++++++++------- 2 files changed, 786 insertions(+), 151 deletions(-) diff --git a/src/mqc_algebra.F03 b/src/mqc_algebra.F03 index 7b5c4cac..103c48e6 100644 --- a/src/mqc_algebra.F03 +++ b/src/mqc_algebra.F03 @@ -653,12 +653,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 diff --git a/src/mqc_est.F03 b/src/mqc_est.F03 index bc7bd66d..04e2093b 100644 --- a/src/mqc_est.F03 +++ b/src/mqc_est.F03 @@ -389,6 +389,11 @@ Module MQC_EST Module Procedure MQC_Eigenvalues_MinLoc End Interface ! +!> \brief Returns the identity MQC_SCF_Integral of the specified size + interface intiden + module procedure MQC_Integral_Identity_Func + end interface +! ! !---------------------------------------------------------------- ! | @@ -3162,6 +3167,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 @@ -11704,8 +11796,8 @@ 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','spin',tmpMat3,tmpMat4) If(present(lEVecs)) call mqc_integral_allocate(lEVecs,'eigenvectors','spin',tmpMat5,tmpMat6) @@ -21463,8 +21555,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: @@ -21781,37 +21872,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) @@ -21824,12 +21932,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(:) @@ -21850,18 +21952,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 @@ -21908,149 +21998,694 @@ 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 -! - If(present(doS2)) then - S2Flag = doS2 - else - S2Flag = .false. - EndIf + Logical::SymmFlag ! - 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 = .True. + SymmFlag = .False. EndIf - - LDetSize = Size(Determinants) - RDetSize = Size(Determinants2) - allocate(NDets1(LDetSize),NDets1A(LDetSize),NDets1B(LDetSize),NDets2(RDetSize),& - NDets2A(RDetSize),NDets2B(RDetSize),lSubsT(LDetSize),lSubsA(LDetSize),& - lSubsB(LDetSize),rSubsT(RDetSize),rSubsA(RDetSize),rSubsB(RDetSize)) ! ! Determine different substitutions requested from input, and number of strings ! of different types. ! - NDets1T = 1 - do i = 1, LDetSize - if(present(subsAIn)) then - subsA = subsAIn(i) - else - subsA = [integer::] - endIf - if(present(subsBIn)) then - subsB = subsBIn(i) - else - subsB = [integer::] - endIf - call get_subLevels(Determinants(i),subsA,subsB,lsubsA(i),lsubsB(i),lsubsT(i)) - call get_expansion_size(Determinants(i),lSubsT(i),lSubsA(i),lSubsB(i),& - NDets1(i),NDets1A(i),NDets1B(i)) - NDets1T = NDets1T*NDets1(i) - endDo - - NDets2T = 1 - do i = 1, RDetSize - if(.not.present(subs2AIn).and..not.present(subs2BIn)) then - subs2A = subsA - subs2B = subsB - elseIf(.not.present(subs2BIn)) then - subs2A = subs2AIn(i) - subs2B = [integer::] - symmFlag = .False. - elseIf(.not.present(subs2AIn)) then - subs2A = [integer::] - subs2B = subs2BIn(i) - symmFlag = .False. - else - subs2A = subs2AIn(i) - subs2B = subs2BIn(i) - symmFlag = .False. - endIf - call get_subLevels(Determinants2(i),subs2A,subs2B,rsubsA(i),rsubsB(i),rsubsT(i)) - call get_expansion_size(Determinants2(i),rSubsT(i),rSubsA(i),rSubsB(i),& - NDets2(i),NDets2A(i),NDets2B(i)) - NDets2T = NDets2T*NDets2(i) - endDo + 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 +! Allocate and fill determinant string arrays ! - If(SymmFlag) then - if(NDets1T.ne.NDets2T) call mqc_error_i('SymmFlag is true but Hamiltonian dimensions are not equal',& - 6,'NDets1T',NDets1T,'NDets2T',NDets2T) - Call Operator_Matrix%init(NDets1T,NDets2T,Storage='StorHerm') - SymmStr = 'hermitian' - Else - Call Operator_Matrix%init(NDets1T,NDets2T) - SymmStr = 'element' - EndIf + Allocate(L_Alpha_Strings(NDets1T),L_Beta_Strings(NDets1T),R_Alpha_Strings(NDets2T),& + R_Beta_Strings(NDets2T)) Do L_Index = 1, NDets1T - If(SymmFlag) then - R_Index_End = L_Index - 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) + 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 - 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 + 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) + RDetSize = Size(Determinants2) + allocate(NDets1(LDetSize),NDets1A(LDetSize),NDets1B(LDetSize),NDets2(RDetSize),& + NDets2A(RDetSize),NDets2B(RDetSize),lSubsT(LDetSize),lSubsA(LDetSize),& + lSubsB(LDetSize),rSubsT(RDetSize),rSubsA(RDetSize),rSubsB(RDetSize)) +! +! Determine different substitutions requested from input, and number of strings +! of different types. +! + NDets1T = 1 + do i = 1, LDetSize + if(present(subsAIn)) then + subsA = subsAIn(i) + else + subsA = [integer::] + endIf + if(present(subsBIn)) then + subsB = subsBIn(i) + else + subsB = [integer::] + endIf + call get_subLevels(Determinants(i),subsA,subsB,lsubsA(i),lsubsB(i),lsubsT(i)) + call get_expansion_size(Determinants(i),lSubsT(i),lSubsA(i),lSubsB(i),& + NDets1(i),NDets1A(i),NDets1B(i)) + NDets1T = NDets1T*NDets1(i) + endDo + + NDets2T = 1 + do i = 1, RDetSize + if(.not.present(subs2AIn).and..not.present(subs2BIn)) then + 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::] + symmFlag = .False. + elseIf(.not.present(subs2AIn)) then + subs2A = [integer::] + subs2B = subs2BIn(i) + symmFlag = .False. + else + subs2A = subs2AIn(i) + subs2B = subs2BIn(i) + symmFlag = .False. + endIf + call get_subLevels(Determinants2(i),subs2A,subs2B,rsubsA(i),rsubsB(i),rsubsT(i)) + call get_expansion_size(Determinants2(i),rSubsT(i),rSubsA(i),rSubsB(i),& + 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 +! + If(SymmFlag) then + if(NDets1T.ne.NDets2T) call mqc_error_i('SymmFlag is true but Hamiltonian dimensions are not equal',& + 6,'NDets1T',NDets1T,'NDets2T',NDets2T) + Call Operator_Matrix%init(NDets1T,NDets2T,Storage='StorHerm') + SymmStr = 'hermitian' + Else + Call Operator_Matrix%init(NDets1T,NDets2T) + SymmStr = 'element' + EndIf + + Do L_Index = 1, NDets1T + If(SymmFlag) then + R_Index_End = L_Index + Else + R_Index_End = NDets2T + EndIf + 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 + 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(*,*) '-------------------------------------------' From d6e84be9a7e2d594a8406d3df0fb6b15f1885d6d Mon Sep 17 00:00:00 2001 From: Lee Thompson Date: Fri, 25 Jul 2025 11:54:38 -0400 Subject: [PATCH 7/9] Fixed bug in integral combine routine --- src/mqc_est.F03 | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/mqc_est.F03 b/src/mqc_est.F03 index c6945724..b884936d 100644 --- a/src/mqc_est.F03 +++ b/src/mqc_est.F03 @@ -4468,7 +4468,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') @@ -4489,7 +4489,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') @@ -4510,7 +4510,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') @@ -4532,7 +4532,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') @@ -4564,7 +4564,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') @@ -4595,7 +4595,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') @@ -4627,7 +4627,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') @@ -4659,7 +4659,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') @@ -4697,7 +4697,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 From afbb257aeffae7a3481160ecaf181bed97e8d00a Mon Sep 17 00:00:00 2001 From: Lee Thompson Date: Fri, 25 Jul 2025 18:44:21 -0400 Subject: [PATCH 8/9] Added mqc_eigenvalues absolute function --- src/mqc_est.F03 | 48 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/src/mqc_est.F03 b/src/mqc_est.F03 index 7d42ef04..587dd1c7 100644 --- a/src/mqc_est.F03 +++ b/src/mqc_est.F03 @@ -200,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 @@ -395,6 +397,11 @@ 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 @@ -15619,6 +15626,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 From a72e6eb6332236009f9d677922a8f52750c31c9f Mon Sep 17 00:00:00 2001 From: Lee Thompson Date: Tue, 5 Aug 2025 18:05:31 -0400 Subject: [PATCH 9/9] Fixed bug in allocations and intializations of integral objects when alpha or beta dimension is zero (fix thanks to Matheus). --- src/mqc_est.F03 | 38 ++++++++++++++++++++++++++++---------- 1 file changed, 28 insertions(+), 10 deletions(-) diff --git a/src/mqc_est.F03 b/src/mqc_est.F03 index 587dd1c7..4daefb52 100644 --- a/src/mqc_est.F03 +++ b/src/mqc_est.F03 @@ -2770,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 & @@ -2780,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 & @@ -3538,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 ! @@ -3563,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) @@ -3571,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 ! @@ -16677,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 @@ -23113,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: @@ -23123,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 @@ -23157,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 !>