From f75b418e5f46cf441237364eb7f3f14018938a51 Mon Sep 17 00:00:00 2001 From: Lee Thompson Date: Tue, 1 Mar 2022 17:22:49 -0500 Subject: [PATCH 01/10] Added code to read multiple asymmetric matrices from matrix file along with option in readESTObj to recover velocity gauge dipole matrices --- src/mqc_gaussian.F03 | 244 +++++++++++++++++++++++++++++++++++++++-- src/mqc_matwrapper.F03 | 2 +- 2 files changed, 237 insertions(+), 9 deletions(-) diff --git a/src/mqc_gaussian.F03 b/src/mqc_gaussian.F03 index 12ee59c0..d5d09bf0 100644 --- a/src/mqc_gaussian.F03 +++ b/src/mqc_gaussian.F03 @@ -1360,6 +1360,18 @@ subroutine MQC_Gaussian_Unformatted_Matrix_Read_Array(fileinfo, & 6,'Present(mqcVarOut)',Present(mqcVarOut),'Present(matrixOut)',Present(matrixOut)) endIf deallocate(arrayTmp) + case('REAL-ASYMMATRIXN') + allocate(arrayTmp(LR)) + call Rd_RBuf(fileinfo%unitNumber,NTot,LenBuf,arrayTmp) + if(Present(matrixOut)) then + call MQC_Matrix_SymmMatrix_Put(matrixOut,arrayTmp((myArrayNum-1)*(LR/N3)+1:myArrayNum*(LR/N3)),'antisymmetric') + elseIf(Present(mqcVarOut)) then + mqcVarOut = mqc_matrixSymm2Full(arrayTmp((myArrayNum-1)*(LR/N3)+1:myArrayNum*(LR/N3)),'U') + else + call mqc_error_l('Reading matrix from Gaussian matrix file, but NO MATRIX SENT to procedure.', & + 6,'Present(mqcVarOut)',Present(mqcVarOut),'Present(matrixOut)',Present(matrixOut)) + endIf + deallocate(arrayTmp) case('COMPLEX-VECTOR') allocate(complexTmp(LR)) call Rd_CBuf(fileinfo%unitNumber,NTot,LenBuf,complexTmp) @@ -1408,10 +1420,6 @@ subroutine MQC_Gaussian_Unformatted_Matrix_Read_Array(fileinfo, & allocate(complexTmp(LR)) call Rd_CBuf(fileinfo%unitNumber,NTot,LenBuf,complexTmp) call MQC_Matrix_SymmMatrix_Put(matrixOut,complexTmp,'antihermitian') -! Matrix files have either symmetric/hermitian storage or antisymmetric/ -! anthermitian storage. MQC currently has a symmetric only storage for both real -! and complex parts so make nonsymmetric matrices square. -! call mqc_matrix_symm2full(matrixOut,'antihermitian') ! Triangular matrices are stored in the order (A(J,I),J=1,I),I=1,N) on the matrix ! file, where first index is the row. Therefore, we need to transpose matrix file ! storage to the MQC lower trangular matrix after reading for correct storage. @@ -1424,6 +1432,12 @@ subroutine MQC_Gaussian_Unformatted_Matrix_Read_Array(fileinfo, & call MQC_Matrix_SymmMatrix_Put(matrixOut,complexTmp((myArrayNum-1)*(LR/N3)+1:myArrayNum*(LR/N3)),'hermitian') matrixOut = transpose(matrixOut) deallocate(complexTmp) + case('COMPLEX-ASYMMATRIXN') + allocate(complexTmp(LR)) + call Rd_CBuf(fileinfo%unitNumber,NTot,LenBuf,complexTmp) + call MQC_Matrix_SymmMatrix_Put(matrixOut,complexTmp((myArrayNum-1)*(LR/N3)+1:myArrayNum*(LR/N3)),'antihermitian') + matrixOut = transpose(matrixOut) + deallocate(complexTmp) case('MIXED') write(*,1020) write(*,1040)' LR = ',LR @@ -3348,11 +3362,14 @@ subroutine mqc_gaussian_unformatted_matrix_get_EST_object(fileinfo,label, & ! 'density' return the density matrix. ! 'scf density' return the SCF density matrix. ! 'overlap' return the overlap matrix. -! 'dipole x' return the x dipole AO integrals. -! 'dipole y' return the y dipole AO integrals. -! 'dipole z' return the z dipole AO integrals. +! 'dipole x' return the length gauge x dipole AO integrals. +! 'dipole y' return the length gauge y dipole AO integrals. +! 'dipole z' return the length gauge z dipole AO integrals. +! 'vel dipole x' return the velocity gauge x dipole AO integrals. +! 'vel dipole y' return the velocity gauge y dipole AO integrals. +! 'vel dipole z' return the velocity gauge z dipole AO integrals. ! 'wavefunction' load the wavefunction object (does not contain -! dipoles. +! dipoles). ! ! * not yet implemented ! @@ -4185,6 +4202,213 @@ subroutine mqc_gaussian_unformatted_matrix_get_EST_object(fileinfo,label, & 'fileinfo%isUnrestricted()', fileinfo%isUnrestricted(), & 'fileinfo%isGeneral()', fileinfo%isGeneral() ) endIf + case('vel dipole x') + if(fileinfo%isRestricted()) then + call fileInfo%getArray('DIP VEL INTEGRALS',tmpMatrixAlpha,foundOut=found) + if(present(foundObj)) foundObj = found + if(.not.found) then + errorMsg = 'DIP VEL INTEGRALS (x) not present on file' + if(present(foundObj)) then + write(6,'(A)') errorMsg + call fileinfo%load() + else + call mqc_error_l(trim(errorMsg),6,'found',found) + endIf + else +! if(MQC_Matrix_HaveComplex(tmpMatrixAlpha)) then +! call mqc_matrix_symm2full(tmpMatrixAlpha,'hermitian') +! tmpMatrixAlpha = transpose(tmpMatrixAlpha) +! endIf + call mqc_integral_allocate(est_integral,'vel dipole x','space',tmpMatrixAlpha) + endIf + elseIf(fileinfo%isUnrestricted()) then + call fileInfo%getArray('DIP VEL INTEGRALS',tmpMatrixAlpha,foundOut=found) + if(present(foundObj)) foundObj = found + if(.not.found) then + errorMsg = 'DIP VEL INTEGRALS (x) not present on file' + if(present(foundObj)) then + write(6,'(A)') errorMsg + call fileinfo%load() + else + call mqc_error_l(trim(errorMsg),6,'found',found) + endIf + else +! if(MQC_Matrix_HaveComplex(tmpMatrixAlpha)) then +! call mqc_matrix_symm2full(tmpMatrixAlpha,'hermitian') +! tmpMatrixAlpha = transpose(tmpMatrixAlpha) +! endIf + call mqc_integral_allocate(est_integral,'vel dipole x','spin',tmpMatrixAlpha, & + tmpMatrixAlpha) + endIf + elseIf(fileinfo%isGeneral()) then + call fileInfo%getArray('DIP VEL INTEGRALS',tmpMatrixAlpha,foundOut=found) + if(present(foundObj)) foundObj = found + if(.not.found) then + errorMsg = 'DIP VEL INTEGRALS (x) not present on file' + if(present(foundObj)) then + write(6,'(A)') errorMsg + call fileinfo%load() + else + call mqc_error_l(trim(errorMsg),6,'found',found) + endIf + else +! if(MQC_Matrix_HaveComplex(tmpMatrixAlpha)) then +! call mqc_matrix_symm2full(tmpMatrixAlpha,'hermitian') +! tmpMatrixAlpha = transpose(tmpMatrixAlpha) +! endIf + nBasis = fileInfo%getVal('nBasis') + call mqc_matrix_spinBlockGHF(tmpMatrixAlpha) + tmpMatrixBeta = tmpMatrixAlpha%mat([nBasis+1,-1],[nBasis+1,-1]) + tmpMatrixBetaAlpha = tmpMatrixAlpha%mat([1,nBasis],[nBasis+1,-1]) + tmpMatrixAlphaBeta = tmpMatrixAlpha%mat([nBasis+1,-1],[1,nBasis]) + tmpMatrixAlpha = tmpMatrixAlpha%mat([1,nBasis],[1,nBasis]) + call mqc_integral_allocate(est_integral,'vel dipole x','general',tmpMatrixAlpha, & + tmpMatrixBeta,tmpMatrixAlphaBeta,tmpMatrixBetaAlpha) + endIf + else + call mqc_error_L('Unknown wavefunction type in getESTObj', 6, & + 'fileinfo%isRestricted()', fileinfo%isRestricted(), & + 'fileinfo%isUnrestricted()', fileinfo%isUnrestricted(), & + 'fileinfo%isGeneral()', fileinfo%isGeneral() ) + endIf + case('vel dipole y') + if(fileinfo%isRestricted()) then + call fileInfo%getArray('DIP VEL INTEGRALS',tmpMatrixAlpha,arrayNum=2,foundOut=found) + if(present(foundObj)) foundObj = found + if(.not.found) then + errorMsg = 'DIP VEL INTEGRALS (y) not present on file' + if(present(foundObj)) then + write(6,'(A)') errorMsg + call fileinfo%load() + else + call mqc_error_l(trim(errorMsg),6,'found',found) + endIf + else +! if(MQC_Matrix_HaveComplex(tmpMatrixAlpha)) then +! call mqc_matrix_symm2full(tmpMatrixAlpha,'hermitian') +! tmpMatrixAlpha = transpose(tmpMatrixAlpha) +! endIf + call mqc_integral_allocate(est_integral,'vel dipole y','space',tmpMatrixAlpha) + endIf + elseIf(fileinfo%isUnrestricted()) then + call fileInfo%getArray('DIP VEL INTEGRALS',tmpMatrixAlpha,arrayNum=2,foundOut=found) + if(present(foundObj)) foundObj = found + if(.not.found) then + errorMsg = 'DIP VEL INTEGRALS (y) not present on file' + if(present(foundObj)) then + write(6,'(A)') errorMsg + call fileinfo%load() + else + call mqc_error_l(trim(errorMsg),6,'found',found) + endIf + else +! if(MQC_Matrix_HaveComplex(tmpMatrixAlpha)) then +! call mqc_matrix_symm2full(tmpMatrixAlpha,'hermitian') +! tmpMatrixAlpha = transpose(tmpMatrixAlpha) +! endIf + call mqc_integral_allocate(est_integral,'vel dipole y','spin',tmpMatrixAlpha, & + tmpMatrixAlpha) + endIf + elseIf(fileinfo%isGeneral()) then + call fileInfo%getArray('DIP VEL INTEGRALS',tmpMatrixAlpha,arrayNum=2,foundOut=found) + if(present(foundObj)) foundObj = found + if(.not.found) then + errorMsg = 'DIP VEL INTEGRALS (y) not present on file' + if(present(foundObj)) then + write(6,'(A)') errorMsg + call fileinfo%load() + else + call mqc_error_l(trim(errorMsg),6,'found',found) + endIf + else +! if(MQC_Matrix_HaveComplex(tmpMatrixAlpha)) then +! call mqc_matrix_symm2full(tmpMatrixAlpha,'hermitian') +! tmpMatrixAlpha = transpose(tmpMatrixAlpha) +! endIf + nBasis = fileInfo%getVal('nBasis') + call mqc_matrix_spinBlockGHF(tmpMatrixAlpha) + tmpMatrixBeta = tmpMatrixAlpha%mat([nBasis+1,-1],[nBasis+1,-1]) + tmpMatrixBetaAlpha = tmpMatrixAlpha%mat([1,nBasis],[nBasis+1,-1]) + tmpMatrixAlphaBeta = tmpMatrixAlpha%mat([nBasis+1,-1],[1,nBasis]) + tmpMatrixAlpha = tmpMatrixAlpha%mat([1,nBasis],[1,nBasis]) + call mqc_integral_allocate(est_integral,'vel dipole y','general',tmpMatrixAlpha, & + tmpMatrixBeta,tmpMatrixAlphaBeta,tmpMatrixBetaAlpha) + endIf + else + call mqc_error_L('Unknown wavefunction type in getESTObj', 6, & + 'fileinfo%isRestricted()', fileinfo%isRestricted(), & + 'fileinfo%isUnrestricted()', fileinfo%isUnrestricted(), & + 'fileinfo%isGeneral()', fileinfo%isGeneral() ) + endIf + case('vel dipole z') + if(fileinfo%isRestricted()) then + call fileInfo%getArray('DIP VEL INTEGRALS',tmpMatrixAlpha,arrayNum=3,foundOut=found) + if(present(foundObj)) foundObj = found + if(.not.found) then + errorMsg = 'DIP VEL INTEGRALS (z) not present on file' + if(present(foundObj)) then + write(6,'(A)') errorMsg + call fileinfo%load() + else + call mqc_error_l(trim(errorMsg),6,'found',found) + endIf + else +! if(MQC_Matrix_HaveComplex(tmpMatrixAlpha)) then +! call mqc_matrix_symm2full(tmpMatrixAlpha,'hermitian') +! tmpMatrixAlpha = transpose(tmpMatrixAlpha) +! endIf + call mqc_integral_allocate(est_integral,'vel dipole z','space',tmpMatrixAlpha) + endIf + elseIf(fileinfo%isUnrestricted()) then + call fileInfo%getArray('DIP VEL INTEGRALS',tmpMatrixAlpha,arrayNum=3,foundOut=found) + if(present(foundObj)) foundObj = found + if(.not.found) then + errorMsg = 'DIP VEL INTEGRALS (z) not present on file' + if(present(foundObj)) then + write(6,'(A)') errorMsg + call fileinfo%load() + else + call mqc_error_l(trim(errorMsg),6,'found',found) + endIf + else +! if(MQC_Matrix_HaveComplex(tmpMatrixAlpha)) then +! call mqc_matrix_symm2full(tmpMatrixAlpha,'hermitian') +! tmpMatrixAlpha = transpose(tmpMatrixAlpha) +! endIf + call mqc_integral_allocate(est_integral,'vel dipole z','spin',tmpMatrixAlpha, & + tmpMatrixAlpha) + endIf + elseIf(fileinfo%isGeneral()) then + call fileInfo%getArray('DIP VEL INTEGRALS',tmpMatrixAlpha,arrayNum=3,foundOut=found) + if(present(foundObj)) foundObj = found + if(.not.found) then + errorMsg = 'DIP VEL INTEGRALS (z) not present on file' + if(present(foundObj)) then + write(6,'(A)') errorMsg + call fileinfo%load() + else + call mqc_error_l(trim(errorMsg),6,'found',found) + endIf + else +! if(MQC_Matrix_HaveComplex(tmpMatrixAlpha)) then +! call mqc_matrix_symm2full(tmpMatrixAlpha,'hermitian') +! tmpMatrixAlpha = transpose(tmpMatrixAlpha) +! endIf + nBasis = fileInfo%getVal('nBasis') + call mqc_matrix_spinBlockGHF(tmpMatrixAlpha) + tmpMatrixBeta = tmpMatrixAlpha%mat([nBasis+1,-1],[nBasis+1,-1]) + tmpMatrixBetaAlpha = tmpMatrixAlpha%mat([1,nBasis],[nBasis+1,-1]) + tmpMatrixAlphaBeta = tmpMatrixAlpha%mat([nBasis+1,-1],[1,nBasis]) + tmpMatrixAlpha = tmpMatrixAlpha%mat([1,nBasis],[1,nBasis]) + call mqc_integral_allocate(est_integral,'vel dipole y','general',tmpMatrixAlpha, & + tmpMatrixBeta,tmpMatrixAlphaBeta,tmpMatrixBetaAlpha) + endIf + else + call mqc_error_L('Unknown wavefunction type in getESTObj', 6, & + 'fileinfo%isRestricted()', fileinfo%isRestricted(), & + 'fileinfo%isUnrestricted()', fileinfo%isUnrestricted(), & + 'fileinfo%isGeneral()', fileinfo%isGeneral() ) + endIf case('wavefunction') if(present(foundObj)) foundObj = .true. if(fileinfo%isRestricted()) then @@ -5099,6 +5323,7 @@ Function MQC_Gaussian_Unformatted_Matrix_Array_Type(NI,NR,N1,N2,N3,N4,N5,NRI,ASy ! "SYMMATRIX" A symmetric/hermitian matrix. ! "ASYMMATRIX" An antisymmetric/antihermitian matrix. ! "SYMMATRIXN" Multiple (N) symmetric/hermitian matrix. +! "ASYMMATRIXN" Multiple (N) asymmetric/antihermitian matrix. ! ! If the input flags do not uniquely identify a known array type, then this ! function returns "UNKNOWN". @@ -5156,6 +5381,9 @@ Function MQC_Gaussian_Unformatted_Matrix_Array_Type(NI,NR,N1,N2,N3,N4,N5,NRI,ASy elseIf(N1.le.-1.and.N2.gt.1.and.N3.ge.1.and.N4.eq.1.and.N5.eq.1.and..not.ASym) then MQC_Gaussian_Unformatted_Matrix_Array_Type = & TRIM(MQC_Gaussian_Unformatted_Matrix_Array_Type)//"-SYMMATRIXN" + elseIf(N1.le.-1.and.N2.gt.1.and.N3.ge.1.and.N4.eq.1.and.N5.eq.1.and.ASym) then + MQC_Gaussian_Unformatted_Matrix_Array_Type = & + TRIM(MQC_Gaussian_Unformatted_Matrix_Array_Type)//"-ASYMMATRIXN" endIf ! return diff --git a/src/mqc_matwrapper.F03 b/src/mqc_matwrapper.F03 index 46200c5d..f18a700c 100644 --- a/src/mqc_matwrapper.F03 +++ b/src/mqc_matwrapper.F03 @@ -1,4 +1,4 @@ -#define UCMGAUOPEN +!#define UCMGAUOPEN Module MQC_MatWrapper ! ! ********************************************************************** From ca941b003e070f86522d03d4bbbd046548bd1547 Mon Sep 17 00:00:00 2001 From: Lee Thompson Date: Sat, 2 Apr 2022 15:58:11 -0400 Subject: [PATCH 02/10] Changed behavior so a warning is given instead of failure if energy list is not set when writing a ghf solution, typically caused by trying to output general est objects when non general est objects were stored. There should only be a problem if you expect that mqc has set an energy list, such as when reading in a ghf solution. Also changed est objects to be able to write and read custom records off the matrix file. This extends the behavior of getarrray and writearray to est objects. --- src/mqc_est.F03 | 5 +- src/mqc_gaussian.F03 | 186 +++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 184 insertions(+), 7 deletions(-) diff --git a/src/mqc_est.F03 b/src/mqc_est.F03 index ca53f2ff..e14f03c7 100644 --- a/src/mqc_est.F03 +++ b/src/mqc_est.F03 @@ -6575,9 +6575,10 @@ subroutine mqc_matrix_undoSpinBlockGHF_integral(integralIn,matrixOut) endif END DO END DO - else if(integralIn%getLabel().eq.'mo coefficients') then - call mqc_error('Energy_List is not set in MQC_Integral object') else + if(integralIn%getLabel().eq.'mo coefficients') & + Write(6,'(A)') 'WARNING: Outputting general MO coefficients but Energy_List is not set. & + &Assuming alternating alpha and Beta orbitals.' j = 1 do i = 1,nAlpha2 call matrixOut%vput(tmpMatrix%vat([0],[i]),[0],[j]) diff --git a/src/mqc_gaussian.F03 b/src/mqc_gaussian.F03 index d5d09bf0..b979f04d 100644 --- a/src/mqc_gaussian.F03 +++ b/src/mqc_gaussian.F03 @@ -1364,7 +1364,8 @@ subroutine MQC_Gaussian_Unformatted_Matrix_Read_Array(fileinfo, & allocate(arrayTmp(LR)) call Rd_RBuf(fileinfo%unitNumber,NTot,LenBuf,arrayTmp) if(Present(matrixOut)) then - call MQC_Matrix_SymmMatrix_Put(matrixOut,arrayTmp((myArrayNum-1)*(LR/N3)+1:myArrayNum*(LR/N3)),'antisymmetric') +! call MQC_Matrix_SymmMatrix_Put(matrixOut,arrayTmp((myArrayNum-1)*(LR/N3)+1:myArrayNum*(LR/N3)),'antisymmetric') + call MQC_Matrix_SymmMatrix_Put(matrixOut,arrayTmp((myArrayNum-1)*(LR/N3)+1:myArrayNum*(LR/N3))) elseIf(Present(mqcVarOut)) then mqcVarOut = mqc_matrixSymm2Full(arrayTmp((myArrayNum-1)*(LR/N3)+1:myArrayNum*(LR/N3)),'U') else @@ -3304,8 +3305,43 @@ subroutine mqc_gaussian_unformatted_matrix_write_EST_object(fileinfo,label, & mqc_integral_array_type(est_wavefunction%fock_matrix) ) endIf case default - call mqc_error_A('Invalid label sent to %writeESTObj.', 6, & - 'mylabel', mylabel ) + call String_Change_Case(label,'u',myLabel) + If(present(est_eigenvalues)) then + if(my_integral_type.eq.'space') then + call fileInfo%writeArray('ALPHA '//trim(myLabel), & + vectorIn=est_eigenvalues%getBlock('alpha')) + elseIf(my_integral_type.eq.'spin') then + call fileInfo%writeArray('ALPHA '//trim(myLabel), & + vectorIn=est_eigenvalues%getBlock('alpha')) + call fileInfo%writeArray('BETA '//trim(myLabel), & + vectorIn=est_eigenvalues%getBlock('beta')) + elseIf(my_integral_type.eq.'general') then + call mqc_matrix_undoSpinBlockGHF(est_eigenvalues,tmpVector) + call fileInfo%writeArray('ALPHA '//trim(myLabel),vectorIn=tmpVector) + else + call mqc_error_a('Unknown wavefunction type in getESTObj', 6, & + 'my_integral_type', my_integral_type ) + endIf + elseIf(present(est_integral)) then + if(my_integral_type.eq.'space') then + call fileInfo%writeArray('ALPHA '//trim(myLabel), & + matrixIn=est_integral%getBlock('alpha'),storage='full') + elseIf(my_integral_type.eq.'spin') then + call fileInfo%writeArray('ALPHA '//trim(myLabel), & + matrixIn=est_integral%getBlock('alpha'),storage='full') + call fileInfo%writeArray('BETA '//trim(myLabel), & + matrixIn=est_integral%getBlock('beta'),storage='full') + elseIf(my_integral_type.eq.'general') then + call mqc_matrix_undoSpinBlockGHF(est_integral,tmpMatrix) + if(.not.mqc_matrix_haveComplex(tmpMatrix)) call MQC_Matrix_Copy_Real2Complex(tmpMatrix) + call fileInfo%writeArray('ALPHA '//trim(myLabel),matrixIn=tmpMatrix,storage='full') + else + call mqc_error_a('Unknown wavefunction type in writeESTObj', 6, & + 'my_integral_type', my_integral_type ) + endIf + else + call mqc_error('Unsuported open label EST object in writeESTObj') + endIf end select ! return @@ -4855,8 +4891,148 @@ subroutine mqc_gaussian_unformatted_matrix_get_EST_object(fileinfo,label, & 'fileinfo%isGeneral()', fileinfo%isGeneral() ) endIf case default - call mqc_error_A('Invalid label sent to %getESTObj.', 6, & - 'mylabel', mylabel ) + call String_Change_Case(label,'u',myLabel) + If(present(est_eigenvalues)) then + if(fileinfo%isRestricted()) then + call fileInfo%getArray('ALPHA '//trim(myLabel),vectorOut=tmpVectorAlpha,foundOut=found) + if(present(foundObj)) foundObj = found + if(.not.found) then + errorMsg = 'ALPHA '//trim(myLabel)//' not present on file' + if(present(foundObj)) then + write(6,'(A)') errorMsg + call fileinfo%load() + else + call mqc_error_l(trim(errorMsg),6,'found',found) + endIf + else + call mqc_eigenvalues_allocate(est_eigenvalues,'','space',tmpVectorAlpha) + endIf + elseIf(fileinfo%isUnrestricted()) then + call fileInfo%getArray('ALPHA '//trim(myLabel),vectorOut=tmpVectorAlpha,foundOut=found) + if(present(foundObj)) foundObj = found + if(.not.found) then + errorMsg = 'ALPHA '//trim(myLabel)//' not present on file' + if(present(foundObj)) then + write(6,'(A)') errorMsg + call fileinfo%load() + else + call mqc_error_l(trim(errorMsg),6,'found',found) + endIf + else + call fileInfo%getArray('BETA '//trim(myLabel),vectorOut=tmpVectorBeta,foundOut=found) + if(present(foundObj)) foundObj = found + if(.not.found) then + errorMsg = 'BETA '//trim(myLabel)//' not present on file' + if(present(foundObj)) then + write(6,'(A)') errorMsg + call fileinfo%load() + else + call mqc_error_l(trim(errorMsg),6,'found',found) + endIf + else + call mqc_eigenvalues_allocate(est_eigenvalues,'','spin',tmpVectorAlpha, & + tmpVectorBeta) + endIf + endIf + elseIf(fileinfo%isGeneral()) then + call fileInfo%getArray('ALPHA '//trim(myLabel),vectorOut=tmpVectorAlpha,foundOut=found) + if(present(foundObj)) foundObj = found + if(.not.found) then + errorMsg = 'ALPHA '//trim(myLabel)//' not present on file' + if(present(foundObj)) then + write(6,'(A)') errorMsg + call fileinfo%load() + else + call mqc_error_l(trim(errorMsg),6,'found',found) + endIf + else + nBasis = fileInfo%getVal('nBasis') + call mqc_matrix_spinBlockGHF(tmpVectorAlpha) + tmpVectorBeta = tmpVectorAlpha%vat(nBasis+1,-1) + tmpVectorAlpha = tmpVectorAlpha%vat(1,nBasis) + call mqc_eigenvalues_allocate(est_eigenvalues,'','general',tmpVectorAlpha, & + tmpVectorBeta) + endIf + else + call mqc_error_L('Unknown wavefunction type in getESTObj', 6, & + 'fileinfo%isRestricted()', fileinfo%isRestricted(), & + 'fileinfo%isUnrestricted()', fileinfo%isUnrestricted(), & + 'fileinfo%isGeneral()', fileinfo%isGeneral() ) + endIf + elseIf(present(est_integral)) then + if(fileinfo%isRestricted()) then + call fileInfo%getArray('ALPHA '//trim(myLabel),tmpMatrixAlpha,foundOut=found) + if(present(foundObj)) foundObj = found + if(.not.found) then + errorMsg = 'ALPHA '//trim(myLabel)//' not present on file' + if(present(foundObj)) then + write(6,'(A)') errorMsg + call fileinfo%load() + else + call mqc_error_l(trim(errorMsg),6,'found',found) + endIf + else + call mqc_integral_allocate(est_integral,'','space',tmpMatrixAlpha) + endIf + elseIf(fileinfo%isUnrestricted()) then + call fileInfo%getArray('ALPHA '//trim(myLabel),tmpMatrixAlpha,foundOut=found) + if(present(foundObj)) foundObj = found + if(.not.found) then + errorMsg = 'ALPHA '//trim(myLabel)//' not present on file' + if(present(foundObj)) then + write(6,'(A)') errorMsg + call fileinfo%load() + else + call mqc_error_l(trim(errorMsg),6,'found',found) + endIf + else + call fileInfo%getArray('BETA '//trim(myLabel),tmpMatrixBeta,foundOut=found) + if(present(foundObj)) foundObj = found + if(.not.found) then + errorMsg = 'BETA '//trim(myLabel)//' not present on file' + if(present(foundObj)) then + write(6,'(A)') errorMsg + call fileinfo%load() + else + call mqc_error_l(trim(errorMsg),6,'found',found) + endIf + else + call mqc_integral_allocate(est_integral,'','spin',tmpMatrixAlpha, & + tmpMatrixBeta) + endIf + endIf + elseIf(fileinfo%isGeneral()) then + call fileInfo%getArray('ALPHA '//trim(myLabel),tmpMatrixAlpha,foundOut=found) + if(present(foundObj)) foundObj = found + if(.not.found) then + errorMsg = 'ALPHA '//trim(myLabel)//' not present on file' + if(present(foundObj)) then + write(6,'(A)') errorMsg + call fileinfo%load() + else + call mqc_error_l(trim(errorMsg),6,'found',found) + endIf + else + nBasis = fileInfo%getVal('nBasis') + call mqc_matrix_spinBlockGHF(tmpMatrixAlpha,fileInfo%getVal('nElectrons'), & + fileInfo%getVal('multiplicity'),elist) + tmpMatrixBeta = tmpMatrixAlpha%mat([nBasis+1,-1],[nBasis+1,-1]) + tmpMatrixBetaAlpha = tmpMatrixAlpha%mat([1,nBasis],[nBasis+1,-1]) + tmpMatrixAlphaBeta = tmpMatrixAlpha%mat([nBasis+1,-1],[1,nBasis]) + tmpMatrixAlpha = tmpMatrixAlpha%mat([1,nBasis],[1,nBasis]) + call mqc_integral_allocate(est_integral,'','general',tmpMatrixAlpha, & + tmpMatrixBeta,tmpMatrixAlphaBeta,tmpMatrixBetaAlpha) + call est_integral%setEList(elist) + endIf + else + call mqc_error_L('Unknown wavefunction type in getESTObj', 6, & + 'fileinfo%isRestricted()', fileinfo%isRestricted(), & + 'fileinfo%isUnrestricted()', fileinfo%isUnrestricted(), & + 'fileinfo%isGeneral()', fileinfo%isGeneral() ) + endIf + else + call mqc_error('Unsuported open label EST object in getESTObj') + endIf end select if(allocated(elist)) deallocate(elist) From 25656c1d56339358a0acd191e5e6c9470bbadb2d Mon Sep 17 00:00:00 2001 From: Lee Thompson Date: Mon, 9 May 2022 15:46:04 -0400 Subject: [PATCH 03/10] Fixed bug in one PDM calculaiton where number of orbitals of interest requires more bit storage than the number of bits to store core plus active orbitals. --- src/mqc_est.F03 | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/src/mqc_est.F03 b/src/mqc_est.F03 index e14f03c7..7b0d5257 100644 --- a/src/mqc_est.F03 +++ b/src/mqc_est.F03 @@ -16330,7 +16330,8 @@ function get_one_gamma_matrix(iOut,iPrint,nBasisIn,determinants,ci_amplitudes,UH Type(MQC_Scalar)::Sgn Integer(kind=int64)::NBasis,nCore,nOrbs,IPos,Det_Diff,ISgn,NBit_Ints,m,n,o,mm,nn,rdet,ldet, & Alpha_Diff_Cnt,Beta_Diff_Cnt,L_A_Start,L_B_Start,L_A_End,L_B_End,R_A_Start,R_B_Start,R_A_End,& - R_B_End,L_A_Sub,L_B_Sub,R_A_Sub,R_B_Sub,L_A_Sub_End,L_B_Sub_End,R_A_Sub_End,R_B_Sub_End + R_B_End,L_A_Sub,L_B_Sub,R_A_Sub,R_B_Sub,L_A_Sub_End,L_B_Sub_End,R_A_Sub_End,R_B_Sub_End, & + string_bits integer(kind=int64),dimension(2)::orb Integer(kind=int64),Dimension(:),Allocatable::Alpha_String_1,Alpha_String_2,Beta_String_1, & Beta_String_2,Alpha_Diff,Beta_Diff @@ -16347,8 +16348,9 @@ function get_one_gamma_matrix(iOut,iPrint,nBasisIn,determinants,ci_amplitudes,UH nOrbs = nBasis endIf nBit_Ints = (nBasis/(Bit_Size(0)-1))+1 - Allocate(Alpha_String_1(NBit_Ints),Alpha_String_2(NBit_Ints),Beta_String_1(NBit_Ints),Beta_String_2(NBit_Ints)) - Allocate(Alpha_Diff(NBit_Ints),Beta_Diff(NBit_Ints)) + string_bits = ((nCore+nOrbs)/(Bit_Size(0)-1))+1 + Allocate(Alpha_String_1(string_bits),Alpha_String_2(string_bits),Beta_String_1(string_bits),Beta_String_2(string_bits)) + Allocate(Alpha_Diff(string_bits),Beta_Diff(string_bits)) call onePDMalpha%init(nOrbs,nOrbs) call onePDMbeta%init(nOrbs,nOrbs) @@ -16412,14 +16414,17 @@ function get_one_gamma_matrix(iOut,iPrint,nBasisIn,determinants,ci_amplitudes,UH else rdet = 1+(R_B_String-1)*determinants%NAlpStr+(R_A_String-1) endIf - Alpha_String_1 = Determinants%Strings%Alpha%vat([L_A_String],[1,NBit_Ints]) - Alpha_String_2 = Determinants%Strings%Alpha%vat([R_A_String],[1,NBit_Ints]) - Beta_String_1 = Determinants%Strings%Beta%vat([L_B_String],[1,NBit_Ints]) - Beta_String_2 = Determinants%Strings%Beta%vat([R_B_String],[1,NBit_Ints]) - Det_Diff = 0 + Alpha_String_1 = 0 + Alpha_String_2 = 0 + Beta_String_1 = 0 + Beta_String_2 = 0 Alpha_Diff_Cnt = 0 Beta_Diff_Cnt = 0 Do m = 1,NBit_Ints + Alpha_String_1(m) = Determinants%Strings%Alpha%at(L_A_String,m) + Alpha_String_2(m) = Determinants%Strings%Alpha%at(R_A_String,m) + Beta_String_1(m) = Determinants%Strings%Beta%at(L_B_String,m) + Beta_String_2(m) = Determinants%Strings%Beta%at(R_B_String,m) Alpha_Diff(m) = IEOR(Alpha_String_1(m),Alpha_String_2(m)) Alpha_Diff_Cnt = Alpha_Diff_Cnt + PopCnt(Alpha_Diff(m)) Beta_Diff(m) = IEOR(Beta_String_1(m),Beta_String_2(m)) From 2b459d58f171123b51cd82669fef35249bc91480 Mon Sep 17 00:00:00 2001 From: Lee Thompson Date: Wed, 25 May 2022 19:21:37 -0400 Subject: [PATCH 04/10] Bug fixes in complex routines in MQC algebra, bug fixes in twoeris_at, completion of 2ERI transformation for general spin orbitals, including parallelization. --- src/mqc_algebra.F03 | 194 ++- src/mqc_est.F03 | 3331 +++++++++++++++++++++++-------------------- 2 files changed, 1918 insertions(+), 1607 deletions(-) diff --git a/src/mqc_algebra.F03 b/src/mqc_algebra.F03 index ccb430b4..afea7f8c 100644 --- a/src/mqc_algebra.F03 +++ b/src/mqc_algebra.F03 @@ -465,6 +465,7 @@ Module MQC_Algebra Module Procedure MQC_Scalar_Complex_RealPart Module Procedure MQC_Vector_Complex_RealPart Module Procedure MQC_Matrix_Complex_RealPart + Module Procedure MQC_R4Tensor_Complex_RealPart End Interface ! !> \brief Returns the imaginary part @@ -472,6 +473,7 @@ Module MQC_Algebra Module Procedure MQC_Scalar_Complex_ImagPart Module Procedure MQC_Vector_Complex_ImagPart Module Procedure MQC_Matrix_Complex_ImagPart + Module Procedure MQC_R4Tensor_Complex_ImagPart End Interface ! !> \brief Returns the sine @@ -14059,12 +14061,11 @@ end function MQC_Matrix_DiagMatrix_Put_Vector_Func !> \author L. M. Thompson !> \date 2018 ! - Subroutine MQC_Matrix_DiagMatrix_Put_Vector(diagVectorIn,mat,iprint) + Subroutine MQC_Matrix_DiagMatrix_Put_Vector(diagVectorIn,mat) ! implicit none class(MQC_Matrix),intent(inOut)::mat class(MQC_Vector),intent(in)::diagVectorIn - integer,optional::iprint ! integer(kind=int64)::n character(len=64)::data_type @@ -14723,8 +14724,8 @@ Recursive Subroutine MQC_Matrix_Matrix_Put(Mat,MatrixIn,Rows,Cols) L = 1 Do M = 1, LenCol Do N = 1, M + Mat%MatI(IndI+M-1,IndK+N-1) = MatrixIn%MatI(L,1) Mat%MatI(IndI+N-1,IndK+M-1) = MatrixIn%MatI(L,1) - Mat%MatI(IndK+M-1,IndI+N-1) = MatrixIn%MatI(L,1) L = L + 1 endDo endDo @@ -14732,8 +14733,8 @@ Recursive Subroutine MQC_Matrix_Matrix_Put(Mat,MatrixIn,Rows,Cols) L = 1 Do M = 1, LenCol Do N = 1, M - Mat%MatI(IndI+N-1,IndK+M-1) = MatrixIn%MatI(L,1) - Mat%MatI(IndK+M-1,IndI+N-1) = (-1)*MatrixIn%MatI(L,1) + Mat%MatI(IndI+M-1,IndK+N-1) = MatrixIn%MatI(L,1) + Mat%MatI(IndI+N-1,IndK+M-1) = (-1)*MatrixIn%MatI(L,1) L = L + 1 endDo endDo @@ -14767,8 +14768,8 @@ Recursive Subroutine MQC_Matrix_Matrix_Put(Mat,MatrixIn,Rows,Cols) L = 1 Do M = 1, LenCol Do N = 1, M + Mat%MatI(IndI+M-1,IndK+N-1) = MatrixIn%MatR(L,1) Mat%MatI(IndI+N-1,IndK+M-1) = MatrixIn%MatR(L,1) - Mat%MatI(IndK+M-1,IndI+N-1) = MatrixIn%MatR(L,1) L = L + 1 endDo endDo @@ -14776,8 +14777,8 @@ Recursive Subroutine MQC_Matrix_Matrix_Put(Mat,MatrixIn,Rows,Cols) L = 1 Do M = 1, LenCol Do N = 1, M - Mat%MatI(IndI+N-1,IndK+M-1) = MatrixIn%MatR(L,1) - Mat%MatI(IndK+M-1,IndI+N-1) = (-1)*MatrixIn%MatR(L,1) + Mat%MatI(IndI+M-1,IndK+N-1) = MatrixIn%MatR(L,1) + Mat%MatI(IndI+N-1,IndK+M-1) = (-1)*MatrixIn%MatR(L,1) L = L + 1 endDo endDo @@ -14811,8 +14812,8 @@ Recursive Subroutine MQC_Matrix_Matrix_Put(Mat,MatrixIn,Rows,Cols) L = 1 Do M = 1, LenCol Do N = 1, M + Mat%MatI(IndI+M-1,IndK+N-1) = MatrixIn%MatC(L,1) Mat%MatI(IndI+N-1,IndK+M-1) = MatrixIn%MatC(L,1) - Mat%MatI(IndK+M-1,IndI+N-1) = MatrixIn%MatC(L,1) L = L + 1 endDo endDo @@ -14820,8 +14821,8 @@ Recursive Subroutine MQC_Matrix_Matrix_Put(Mat,MatrixIn,Rows,Cols) L = 1 Do M = 1, LenCol Do N = 1, M - Mat%MatI(IndI+N-1,IndK+M-1) = MatrixIn%MatC(L,1) - Mat%MatI(IndK+M-1,IndI+N-1) = (-1)*MatrixIn%MatC(L,1) + Mat%MatI(IndI+M-1,IndK+N-1) = MatrixIn%MatC(L,1) + Mat%MatI(IndI+N-1,IndK+M-1) = (-1)*MatrixIn%MatC(L,1) L = L + 1 endDo endDo @@ -14829,8 +14830,8 @@ Recursive Subroutine MQC_Matrix_Matrix_Put(Mat,MatrixIn,Rows,Cols) L = 1 Do M = 1, LenCol Do N = 1, M - Mat%MatI(IndI+N-1,IndK+M-1) = MatrixIn%MatC(L,1) - Mat%MatI(IndK+M-1,IndI+N-1) = conjg(MatrixIn%MatC(L,1)) + Mat%MatI(IndI+M-1,IndK+N-1) = MatrixIn%MatC(L,1) + Mat%MatI(IndI+N-1,IndK+M-1) = conjg(MatrixIn%MatC(L,1)) L = L + 1 endDo endDo @@ -14838,8 +14839,8 @@ Recursive Subroutine MQC_Matrix_Matrix_Put(Mat,MatrixIn,Rows,Cols) L = 1 Do M = 1, LenCol Do N = 1, M - Mat%MatI(IndI+N-1,IndK+M-1) = MatrixIn%MatC(L,1) - Mat%MatI(IndK+M-1,IndI+N-1) = (-1)*conjg(MatrixIn%MatC(L,1)) + Mat%MatI(IndI+M-1,IndK+N-1) = MatrixIn%MatC(L,1) + Mat%MatI(IndI+N-1,IndK+M-1) = (-1)*conjg(MatrixIn%MatC(L,1)) L = L + 1 endDo endDo @@ -14877,8 +14878,8 @@ Recursive Subroutine MQC_Matrix_Matrix_Put(Mat,MatrixIn,Rows,Cols) L = 1 Do M = 1, LenCol Do N = 1, M + Mat%MatR(IndI+M-1,IndK+N-1) = MatrixIn%MatI(L,1) Mat%MatR(IndI+N-1,IndK+M-1) = MatrixIn%MatI(L,1) - Mat%MatR(IndK+M-1,IndI+N-1) = MatrixIn%MatI(L,1) L = L + 1 endDo endDo @@ -14886,8 +14887,8 @@ Recursive Subroutine MQC_Matrix_Matrix_Put(Mat,MatrixIn,Rows,Cols) L = 1 Do M = 1, LenCol Do N = 1, M - Mat%MatR(IndI+N-1,IndK+M-1) = MatrixIn%MatI(L,1) - Mat%MatR(IndK+M-1,IndI+N-1) = (-1)*MatrixIn%MatI(L,1) + Mat%MatR(IndI+M-1,IndK+N-1) = MatrixIn%MatI(L,1) + Mat%MatR(IndI+N-1,IndK+M-1) = (-1)*MatrixIn%MatI(L,1) L = L + 1 endDo endDo @@ -14920,8 +14921,8 @@ Recursive Subroutine MQC_Matrix_Matrix_Put(Mat,MatrixIn,Rows,Cols) L = 1 Do M = 1, LenCol Do N = 1, M + Mat%MatR(IndI+M-1,IndK+N-1) = MatrixIn%MatR(L,1) Mat%MatR(IndI+N-1,IndK+M-1) = MatrixIn%MatR(L,1) - Mat%MatR(IndK+M-1,IndI+N-1) = MatrixIn%MatR(L,1) L = L + 1 endDo endDo @@ -14929,8 +14930,8 @@ Recursive Subroutine MQC_Matrix_Matrix_Put(Mat,MatrixIn,Rows,Cols) L = 1 Do M = 1, LenCol Do N = 1, M - Mat%MatR(IndI+N-1,IndK+M-1) = MatrixIn%MatR(L,1) - Mat%MatR(IndK+M-1,IndI+N-1) = (-1)*MatrixIn%MatR(L,1) + Mat%MatR(IndI+M-1,IndK+N-1) = MatrixIn%MatR(L,1) + Mat%MatR(IndI+N-1,IndK+M-1) = (-1)*MatrixIn%MatR(L,1) L = L + 1 endDo endDo @@ -14964,8 +14965,8 @@ Recursive Subroutine MQC_Matrix_Matrix_Put(Mat,MatrixIn,Rows,Cols) L = 1 Do M = 1, LenCol Do N = 1, M + Mat%MatR(IndI+M-1,IndK+N-1) = MatrixIn%MatC(L,1) Mat%MatR(IndI+N-1,IndK+M-1) = MatrixIn%MatC(L,1) - Mat%MatR(IndK+M-1,IndI+N-1) = MatrixIn%MatC(L,1) L = L + 1 endDo endDo @@ -14973,8 +14974,8 @@ Recursive Subroutine MQC_Matrix_Matrix_Put(Mat,MatrixIn,Rows,Cols) L = 1 Do M = 1, LenCol Do N = 1, M - Mat%MatR(IndI+N-1,IndK+M-1) = MatrixIn%MatC(L,1) - Mat%MatR(IndK+M-1,IndI+N-1) = (-1)*MatrixIn%MatC(L,1) + Mat%MatR(IndI+M-1,IndK+N-1) = MatrixIn%MatC(L,1) + Mat%MatR(IndI+N-1,IndK+M-1) = (-1)*MatrixIn%MatC(L,1) L = L + 1 endDo endDo @@ -14982,8 +14983,8 @@ Recursive Subroutine MQC_Matrix_Matrix_Put(Mat,MatrixIn,Rows,Cols) L = 1 Do M = 1, LenCol Do N = 1, M - Mat%MatR(IndI+N-1,IndK+M-1) = MatrixIn%MatC(L,1) - Mat%MatR(IndK+M-1,IndI+N-1) = conjg(MatrixIn%MatC(L,1)) + Mat%MatR(IndI+M-1,IndK+N-1) = MatrixIn%MatC(L,1) + Mat%MatR(IndI+N-1,IndK+M-1) = conjg(MatrixIn%MatC(L,1)) L = L + 1 endDo endDo @@ -14991,8 +14992,8 @@ Recursive Subroutine MQC_Matrix_Matrix_Put(Mat,MatrixIn,Rows,Cols) L = 1 Do M = 1, LenCol Do N = 1, M - Mat%MatR(IndI+N-1,IndK+M-1) = MatrixIn%MatC(L,1) - Mat%MatR(IndK+M-1,IndI+N-1) = (-1)*conjg(MatrixIn%MatC(L,1)) + Mat%MatR(IndI+M-1,IndK+N-1) = MatrixIn%MatC(L,1) + Mat%MatR(IndI+N-1,IndK+M-1) = (-1)*conjg(MatrixIn%MatC(L,1)) L = L + 1 endDo endDo @@ -15030,8 +15031,8 @@ Recursive Subroutine MQC_Matrix_Matrix_Put(Mat,MatrixIn,Rows,Cols) L = 1 Do M = 1, LenCol Do N = 1, M + Mat%MatC(IndI+M-1,IndK+N-1) = MatrixIn%MatI(L,1) Mat%MatC(IndI+N-1,IndK+M-1) = MatrixIn%MatI(L,1) - Mat%MatC(IndK+M-1,IndI+N-1) = MatrixIn%MatI(L,1) L = L + 1 endDo endDo @@ -15039,8 +15040,8 @@ Recursive Subroutine MQC_Matrix_Matrix_Put(Mat,MatrixIn,Rows,Cols) L = 1 Do M = 1, LenCol Do N = 1, M - Mat%MatC(IndI+N-1,IndK+M-1) = MatrixIn%MatI(L,1) - Mat%MatC(IndK+M-1,IndI+N-1) = (-1)*MatrixIn%MatI(L,1) + Mat%MatC(IndI+M-1,IndK+N-1) = MatrixIn%MatI(L,1) + Mat%MatC(IndI+N-1,IndK+M-1) = (-1)*MatrixIn%MatI(L,1) L = L + 1 endDo endDo @@ -15073,8 +15074,8 @@ Recursive Subroutine MQC_Matrix_Matrix_Put(Mat,MatrixIn,Rows,Cols) L = 1 Do M = 1, LenCol Do N = 1, M + Mat%MatC(IndI+M-1,IndK+N-1) = MatrixIn%MatR(L,1) Mat%MatC(IndI+N-1,IndK+M-1) = MatrixIn%MatR(L,1) - Mat%MatC(IndK+M-1,IndI+N-1) = MatrixIn%MatR(L,1) L = L + 1 endDo endDo @@ -15082,8 +15083,8 @@ Recursive Subroutine MQC_Matrix_Matrix_Put(Mat,MatrixIn,Rows,Cols) L = 1 Do M = 1, LenCol Do N = 1, M - Mat%MatC(IndI+N-1,IndK+M-1) = MatrixIn%MatR(L,1) - Mat%MatC(IndK+M-1,IndI+N-1) = (-1)*MatrixIn%MatR(L,1) + Mat%MatC(IndI+M-1,IndK+N-1) = MatrixIn%MatR(L,1) + Mat%MatC(IndI+N-1,IndK+M-1) = (-1)*MatrixIn%MatR(L,1) L = L + 1 endDo endDo @@ -15116,8 +15117,8 @@ Recursive Subroutine MQC_Matrix_Matrix_Put(Mat,MatrixIn,Rows,Cols) L = 1 Do M = 1, LenCol Do N = 1, M + Mat%MatC(IndI+M-1,IndK+N-1) = MatrixIn%MatC(L,1) Mat%MatC(IndI+N-1,IndK+M-1) = MatrixIn%MatC(L,1) - Mat%MatC(IndK+M-1,IndI+N-1) = MatrixIn%MatC(L,1) L = L + 1 endDo endDo @@ -15125,8 +15126,8 @@ Recursive Subroutine MQC_Matrix_Matrix_Put(Mat,MatrixIn,Rows,Cols) L = 1 Do M = 1, LenCol Do N = 1, M - Mat%MatC(IndI+N-1,IndK+M-1) = MatrixIn%MatC(L,1) - Mat%MatC(IndK+M-1,IndI+N-1) = (-1)*MatrixIn%MatC(L,1) + Mat%MatC(IndI+M-1,IndK+N-1) = MatrixIn%MatC(L,1) + Mat%MatC(IndI+N-1,IndK+M-1) = (-1)*MatrixIn%MatC(L,1) L = L + 1 endDo endDo @@ -15134,8 +15135,8 @@ Recursive Subroutine MQC_Matrix_Matrix_Put(Mat,MatrixIn,Rows,Cols) L = 1 Do M = 1, LenCol Do N = 1, M - Mat%MatC(IndI+N-1,IndK+M-1) = MatrixIn%MatC(L,1) - Mat%MatC(IndK+M-1,IndI+N-1) = conjg(MatrixIn%MatC(L,1)) + Mat%MatC(IndI+M-1,IndK+N-1) = MatrixIn%MatC(L,1) + Mat%MatC(IndI+N-1,IndK+M-1) = conjg(MatrixIn%MatC(L,1)) L = L + 1 endDo endDo @@ -15143,8 +15144,8 @@ Recursive Subroutine MQC_Matrix_Matrix_Put(Mat,MatrixIn,Rows,Cols) L = 1 Do M = 1, LenCol Do N = 1, M - Mat%MatC(IndI+N-1,IndK+M-1) = MatrixIn%MatC(L,1) - Mat%MatC(IndK+M-1,IndI+N-1) = (-1)*conjg(MatrixIn%MatC(L,1)) + Mat%MatC(IndI+M-1,IndK+N-1) = MatrixIn%MatC(L,1) + Mat%MatC(IndI+N-1,IndK+M-1) = (-1)*conjg(MatrixIn%MatC(L,1)) L = L + 1 endDo endDo @@ -27523,12 +27524,12 @@ Subroutine MQC_Set_Tensor2ComplexArray(ArrayOut,TensorIn) Size(TensorIn%ITen,3),Size(TensorIn%ITen,4))) ArrayOut = TensorIn%ITen ElseIf(MQC_R4Tensor_HaveReal(TensorIn)) then - allocate(ArrayOut(Size(TensorIn%ITen,1),Size(TensorIn%ITen,2),& - Size(TensorIn%ITen,3),Size(TensorIn%ITen,4))) + allocate(ArrayOut(Size(TensorIn%RTen,1),Size(TensorIn%RTen,2),& + Size(TensorIn%RTen,3),Size(TensorIn%RTen,4))) ArrayOut = TensorIn%RTen ElseIf(MQC_R4Tensor_HaveComplex(TensorIn)) then - allocate(ArrayOut(Size(TensorIn%ITen,1),Size(TensorIn%ITen,2),& - Size(TensorIn%ITen,3),Size(TensorIn%ITen,4))) + allocate(ArrayOut(Size(TensorIn%CTen,1),Size(TensorIn%CTen,2),& + Size(TensorIn%CTen,3),Size(TensorIn%CTen,4))) ArrayOut = TensorIn%CTen Else Call MQC_Error_A('TensorIn type unknown in MQC_Tensor2ComplexArray', 6, & @@ -29321,6 +29322,7 @@ Function MQC_R4Tensor_Matrix_At(TensorIn,D1,D2,D3,D4) Result(MatrixOut) if((size(D1).eq.2.or.D1(1).eq.0).and.(size(D2).eq.2.or.D2(1).eq.0)) then CompR4Tensor = MQC_R4Tensor_R4Tensor_At(TensorIn,D1,D2, & [D3(1),D3(1)],[D4(1),D4(1)]) + flush(6) matrixOut = CompR4Tensor(:,:,1,1) elseIf((size(D1).eq.2.or.D1(1).eq.0).and.(size(D3).eq.2.or.D3(1).eq.0)) then CompR4Tensor = MQC_R4Tensor_R4Tensor_At(TensorIn,D1,[D2(1),D2(1)], & @@ -29937,7 +29939,7 @@ Recursive Subroutine MQC_R4Tensor_R4Tensor_Put(Ten,TensorIn,D1,D2,D3,D4) Len1 = IndJ-IndI+1 If (Len1.le.0.or.Len1.ne.Size(TensorIn,1).or.Len1.gt.Size(Ten,1)) & Call MQC_Error_I('Dimension 1 length badly defined in MQC_R4Tensor_R4Tensor_Put', & - 6,'Len1',Len1,'Size(TensorIn,1)',Size(TensorIn,1) ) + 6,'Len1',Len1,'Size(TensorIn,1)',Size(TensorIn,1),'Size(Ten,1)',Size(Ten,1)) If (IndI.le.0.or.IndI.gt.(Size(Ten,1)-Len1+1)) Call MQC_Error_I('Index I out of & & bounds in MQC_R4Tensor_R4Tensor_Put',6,'IndI',IndI,'Size(Ten,1)', Size(Ten,1), & 'Len1', Len1 ) @@ -30621,5 +30623,105 @@ Function MQC_R4Tensor_Complex_Conjugate(R4TensorIn) Result(R4TensorOut) ! Return End Function MQC_R4Tensor_Complex_Conjugate +! +! +! PROCEDURE MQC_R4Tensor_Complex_RealPart +! +!> \brief MQC_R4Tensor_Complex_RealPart is a function that returns a MQC R4Tensor +!> with elements containing the real part of elements of another MQC R4Tensor +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_R4Tensor_Complex_RealPart is a function that returns a MQC R4Tensor with +!> elements containing the real part of elements of another MQC R4Tensor. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] A +!> \verbatim +!> A is Class(MQC_R4Tensor) +!> The name of the MQC_R4Tensor variable. +!> \endverbatim +! +! Authors: +! ======== +!> \author L. M. Thompson +!> \date 2022 +! + function mqc_R4Tensor_complex_realPart(A) result(output) +! + Implicit None + class(mqc_R4Tensor),intent(in)::A + type(mqc_R4Tensor)::output + integer(kind=int64)::i +! +! Do the work. +! + output = A + if(A%Data_Type.eq.'Complex') then + call MQC_R4Tensor_Copy_Complex2Real(output) + endIf +! + end function mqc_R4Tensor_complex_realPart +! +! +! PROCEDURE MQC_R4Tensor_Complex_ImagPart +! +!> \brief MQC_R4Tensor_Complex_ImagPart is a function that returns a MQC R4Tensor +!> with elements containing the imaginary part of elements of another MQC R4Tensor +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_R4Tensor_Complex_ImagPart is a function that returns a MQC R4Tensor with +!> elements containing the imaginary part of elements of another MQC R4Tensor. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] A +!> \verbatim +!> A is Class(MQC_R4Tensor) +!> The name of the MQC_R4Tensor variable. +!> \endverbatim +! +! Authors: +! ======== +!> \author L. M. Thompson +!> \date 2022 +! + function mqc_R4Tensor_complex_imagPart(A) result(output) +! + Implicit None + class(mqc_R4Tensor),intent(in)::A + type(mqc_R4Tensor)::output + integer(kind=int64)::i +! +! Do the work. +! + output = A + select case (A%Data_Type) + case('Integer') + output%ITen = 0 + case('Real') + output%RTen = 0.0 + case('Complex') + call MQC_R4Tensor_Copy_Complex2Real(output) + output%RTen = aimag(A%CTen) + case default + call mqc_error_a('Unrecognized MQC_R4Tensor data type in mqc_matrix_complex_realPart',& + 6,'A%Data_Type',A%Data_Type) + end select +! + end function mqc_R4Tensor_complex_imagPart +! ! End Module MQC_Algebra diff --git a/src/mqc_est.F03 b/src/mqc_est.F03 index 7b0d5257..bf187e7c 100644 --- a/src/mqc_est.F03 +++ b/src/mqc_est.F03 @@ -3275,8 +3275,8 @@ function mqc_integral_output_orbitals(integral,orbString,alphaOrbsIn,betaOrbsIn, elseIf(my_axis.eq.1) then call outMatrixAlpha%init(size(alphaOrbs),nDimAlpha2) call outMatrixBeta%init(size(betaOrbs),nDimBeta2) - call outMatrixBetaAlpha%init(size(betaOrbs),nDimAlpha1) - call outMatrixAlphaBeta%init(size(alphaOrbs),nDimBeta1) + call outMatrixBetaAlpha%init(size(betaOrbs),nDimAlpha2) + call outMatrixAlphaBeta%init(size(alphaOrbs),nDimBeta2) do i = 1,size(alphaOrbs) if(alphaOrbs(i).ge.0) then call outMatrixAlpha%vput(tmpMatrixAlpha%vat([alphaOrbs(i)],[0]),[i],[0]) @@ -3288,8 +3288,8 @@ function mqc_integral_output_orbitals(integral,orbString,alphaOrbsIn,betaOrbsIn, endDo do i = 1,size(betaOrbs) if(betaOrbs(i).ge.0) then - call outMatrixBeta%vput(tmpMatrixBeta%vat([0],[betaOrbs(i)]),[0],[i]) - call outMatrixBetaAlpha%vput(tmpMatrixBetaAlpha%vat([0],[betaOrbs(i)]),[0],[i]) + call outMatrixBeta%vput(tmpMatrixBeta%vat([betaOrbs(i)],[0]),[0],[i]) + call outMatrixBetaAlpha%vput(tmpMatrixBetaAlpha%vat([betaOrbs(i)],[0]),[0],[i]) else call outMatrixBetaAlpha%vput(tmpMatrixAlpha%vat([abs(betaOrbs(i))],[0]),[i],[0]) call outMatrixBeta%vput(tmpMatrixAlphaBeta%vat([abs(betaOrbs(i))],[0]),[i],[0]) @@ -9661,6 +9661,7 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen endIf endIf elseIf(twoERIs%integralType.eq.'spin') then + !AAAA if(myI.le.twoERIs%blockSize('alpha',1).and.myJ.le.twoERIs%blockSize('alpha',2).and. & myK.le.twoERIs%blockSize('alpha',3).and.myL.le.twoERIs%blockSize('alpha',4)) then if(twoERIs%hasSpinBlock('aaaa')) then @@ -9672,6 +9673,7 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen element = twoERIs%alpha%at(myI,myJ,myK,myL) - twoERIs%alpha%at(myI,myL,myK,myJ) endIf endIf + !BBBB elseIf(myI.gt.twoERIs%blockSize('alpha',1).and.myJ.gt.twoERIs%blockSize('alpha',2).and. & myK.gt.twoERIs%blockSize('alpha',3).and.myL.gt.twoERIs%blockSize('alpha',4)) then myI = myI - twoERIs%blockSize('alpha',1) @@ -9687,6 +9689,7 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen element = twoERIs%beta%at(myI,myJ,myK,myL) - twoERIs%beta%at(myI,myL,myK,myJ) endIf endIf + !AABB elseIf(myI.le.twoERIs%blockSize('alpha',1).and.myJ.le.twoERIs%blockSize('alpha',2).and. & myK.gt.twoERIs%blockSize('alpha',3).and.myL.gt.twoERIs%blockSize('alpha',4)) then myK = myK - twoERIs%blockSize('alpha',3) @@ -9694,13 +9697,13 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen if(interaction.eq.'coulomb'.or.interaction.eq.'doublebar') then if(twoERIs%hasSpinBlock('aabb')) element = twoERIs%alphaBeta%at(myI,myJ,myK,myL) endIf + !BBAA elseIf(myI.gt.twoERIs%blockSize('alpha',1).and.myJ.gt.twoERIs%blockSize('alpha',2).and. & myK.le.twoERIs%blockSize('alpha',3).and.myL.le.twoERIs%blockSize('alpha',4)) then myI = myI - twoERIs%blockSize('alpha',1) myJ = myJ - twoERIs%blockSize('alpha',2) if(twoERIs%storageType.eq.'symm') then if(interaction.eq.'coulomb'.or.interaction.eq.'doublebar') then -! if(twoERIs%hasSpinBlock('aabb')) element = twoERIs%alphaBeta%at(myI,myJ,myK,myL) if(twoERIs%hasSpinBlock('aabb')) element = twoERIs%alphaBeta%at(myK,myL,myI,myJ) endIf else @@ -9708,6 +9711,7 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen if(twoERIs%hasSpinBlock('bbaa')) element = twoERIs%betaAlpha%at(myI,myJ,myK,myL) endIf endIf + !ABBA elseIf(myI.le.twoERIs%blockSize('alpha',1).and.myJ.gt.twoERIs%blockSize('alpha',2).and. & myK.gt.twoERIs%blockSize('alpha',3).and.myL.le.twoERIs%blockSize('alpha',4)) then myJ = myJ - twoERIs%blockSize('alpha',2) @@ -9717,16 +9721,15 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen elseif(interaction.eq.'doublebar') then if(twoERIs%hasSpinBlock('aabb')) element = (-1)*twoERIs%alphaBeta%at(myI,myL,myK,myJ) endIf + !BAAB elseIf(myI.gt.twoERIs%blockSize('alpha',1).and.myJ.le.twoERIs%blockSize('alpha',2).and. & myK.le.twoERIs%blockSize('alpha',3).and.myL.gt.twoERIs%blockSize('alpha',4)) then myI = myI - twoERIs%blockSize('alpha',1) myL = myL - twoERIs%blockSize('alpha',4) if(twoERIs%storageType.eq.'symm') then if(interaction.eq.'exchange') then -! if(twoERIs%hasSpinBlock('aabb')) element = twoERIs%alphaBeta%at(myI,myL,myK,myJ) if(twoERIs%hasSpinBlock('aabb')) element = twoERIs%alphaBeta%at(myK,myJ,myI,myL) elseIf(interaction.eq.'doublebar') then -! if(twoERIs%hasSpinBlock('aabb')) element = (-1)*twoERIs%alphaBeta%at(myI,myL,myK,myJ) if(twoERIs%hasSpinBlock('aabb')) element = (-1)*twoERIs%alphaBeta%at(myK,myJ,myI,myL) endIf else @@ -9738,6 +9741,7 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen endIf endIf elseIf(twoERIs%integralType.eq.'general') then + !AAAA if(myI.le.twoERIs%blockSize('alpha',1).and.myJ.le.twoERIs%blockSize('alpha',2).and. & myK.le.twoERIs%blockSize('alpha',3).and.myL.le.twoERIs%blockSize('alpha',4)) then if(twoERIs%hasSpinBlock('aaaa')) then @@ -9749,6 +9753,7 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen element = twoERIs%alpha%at(myI,myJ,myK,myL) - twoERIs%alpha%at(myI,myL,myK,myJ) endIf endIf + !BBBB elseIf(myI.gt.twoERIs%blockSize('alpha',1).and.myJ.gt.twoERIs%blockSize('alpha',2).and. & myK.gt.twoERIs%blockSize('alpha',3).and.myL.gt.twoERIs%blockSize('alpha',4)) then myI = myI - twoERIs%blockSize('alpha',1) @@ -9764,6 +9769,7 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen element = twoERIs%beta%at(myI,myJ,myK,myL) - twoERIs%beta%at(myI,myL,myK,myJ) endIf endIf + !AABB elseIf(myI.le.twoERIs%blockSize('alpha',1).and.myJ.le.twoERIs%blockSize('alpha',2).and. & myK.gt.twoERIs%blockSize('alpha',3).and.myL.gt.twoERIs%blockSize('alpha',4)) then myK = myK - twoERIs%blockSize('alpha',3) @@ -9776,20 +9782,17 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen if(twoERIs%hasSpinBlock('aabb')) element = twoeris%alphaBeta%at(myI,myJ,myK,myL) if(twoERIs%hasSpinBlock('abba')) element = element - twoeris%abba%at(myI,myL,myK,myJ) endIf + !BBAA elseIf(myI.gt.twoERIs%blockSize('alpha',1).and.myJ.gt.twoERIs%blockSize('alpha',2).and. & myK.le.twoERIs%blockSize('alpha',3).and.myL.le.twoERIs%blockSize('alpha',4)) then myI = myI - twoERIs%blockSize('alpha',1) myJ = myJ - twoERIs%blockSize('alpha',2) if(twoERIs%storageType.eq.'symm') then if(interaction.eq.'coulomb') then -! if(twoERIs%hasSpinBlock('aabb')) element = twoERIs%alphaBeta%at(myI,myJ,myK,myL) if(twoERIs%hasSpinBlock('aabb')) element = twoERIs%alphaBeta%at(myK,myL,myI,myJ) elseif(interaction.eq.'exchange') then -! if(twoERIs%hasSpinBlock('abba')) element = twoeris%abba%at(myI,myL,myK,myJ) if(twoERIs%hasSpinBlock('abba')) element = twoeris%abba%at(myK,myJ,myI,myL) elseif(interaction.eq.'doublebar') then -! if(twoERIs%hasSpinBlock('aabb')) element = twoeris%alphaBeta%at(myI,myJ,myK,myL) -! if(twoERIs%hasSpinBlock('abba')) element = element - twoeris%abba%at(myI,myL,myK,myJ) if(twoERIs%hasSpinBlock('aabb')) element = twoeris%alphaBeta%at(myK,myL,myI,myJ) if(twoERIs%hasSpinBlock('abba')) element = element - twoeris%abba%at(myK,myJ,myI,myL) endIf @@ -9803,6 +9806,7 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen if(twoERIs%hasSpinBlock('baab')) element = element - twoeris%baab%at(myI,myL,myK,myJ) endIf endIf + !ABBA elseIf(myI.le.twoERIs%blockSize('alpha',1).and.myJ.gt.twoERIs%blockSize('alpha',2).and. & myK.gt.twoERIs%blockSize('alpha',3).and.myL.le.twoERIs%blockSize('alpha',4)) then myJ = myJ - twoERIs%blockSize('alpha',2) @@ -9815,22 +9819,19 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen if(twoERIs%hasSpinBlock('abba')) element = twoeris%abba%at(myI,myJ,myK,myL) if(twoERIs%hasSpinBlock('aabb')) element = element - twoeris%alphaBeta%at(myI,myL,myK,myJ) endIf + !BAAB elseIf(myI.gt.twoERIs%blockSize('alpha',1).and.myJ.le.twoERIs%blockSize('alpha',2).and. & myK.le.twoERIs%blockSize('alpha',3).and.myL.gt.twoERIs%blockSize('alpha',4)) then myI = myI - twoERIs%blockSize('alpha',1) myL = myL - twoERIs%blockSize('alpha',4) if(twoERIs%storageType.eq.'symm') then if(interaction.eq.'coulomb') then -! if(twoERIs%hasSpinBlock('abba')) element = twoERIs%abba%at(myI,myJ,myK,myL) if(twoERIs%hasSpinBlock('abba')) element = twoERIs%abba%at(myK,myL,myI,myJ) elseif(interaction.eq.'exchange') then -! if(twoERIs%hasSpinBlock('aabb')) element = twoeris%alphaBeta%at(myI,myL,myK,myJ) if(twoERIs%hasSpinBlock('aabb')) element = twoeris%alphaBeta%at(myK,myJ,myI,myL) elseif(interaction.eq.'doublebar') then -! if(twoERIs%hasSpinBlock('abba')) element = twoeris%abba%at(myK,myL,myI,myJ) -! if(twoERIs%hasSpinBlock('aabb')) element = element - twoeris%alphaBeta%at(myK,myJ,myI,myL) - if(twoERIs%hasSpinBlock('abba')) element = twoeris%abba%at(myI,myJ,myK,myL) - if(twoERIs%hasSpinBlock('aabb')) element = element - twoeris%alphaBeta%at(myI,myL,myK,myJ) + if(twoERIs%hasSpinBlock('abba')) element = twoeris%abba%at(myK,myL,myI,myJ) + if(twoERIs%hasSpinBlock('aabb')) element = element - twoeris%alphaBeta%at(myK,myJ,myI,myL) endIf else if(interaction.eq.'coulomb') then @@ -9842,6 +9843,7 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen if(twoERIs%hasSpinBlock('bbaa')) element = element - twoeris%betaAlpha%at(myI,myL,myK,myJ) endIf endIf + !ABAB elseIf(myI.le.twoERIs%blockSize('alpha',1).and.myJ.gt.twoERIs%blockSize('alpha',2).and. & myK.le.twoERIs%blockSize('alpha',3).and.myL.gt.twoERIs%blockSize('alpha',4)) then myJ = myJ - twoERIs%blockSize('alpha',2) @@ -9855,6 +9857,7 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen element = twoERIs%abab%at(myI,myJ,myK,myL) - twoERIs%abab%at(myI,myL,myK,myJ) endIf endIf + !BABA elseIf(myI.gt.twoERIs%blockSize('alpha',1).and.myJ.le.twoERIs%blockSize('alpha',2).and. & myK.gt.twoERIs%blockSize('alpha',3).and.myL.le.twoERIs%blockSize('alpha',4)) then myI = myI - twoERIs%blockSize('alpha',1) @@ -9862,13 +9865,10 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen if(twoERIs%storageType.eq.'symm') then if(twoERIs%hasSpinBlock('abab')) then if(interaction.eq.'coulomb') then -! element = conjg(twoERIs%abab%at(myI,myJ,myK,myL)) element = conjg(twoERIs%abab%at(myJ,myI,myL,myK)) elseIf(interaction.eq.'exchange') then -! element = conjg(twoeris%abab%at(myI,myL,myK,myJ)) element = conjg(twoeris%abab%at(myL,myI,myJ,myK)) elseIf(interaction.eq.'doublebar') then -! element = conjg(twoERIs%abab%at(myI,myJ,myK,myL)) - conjg(twoERIs%abab%at(myI,myL,myK,myJ)) element = conjg(twoERIs%abab%at(myJ,myI,myL,myK)) - conjg(twoERIs%abab%at(myL,myI,myJ,myK)) endIf endIf @@ -9877,12 +9877,13 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen if(interaction.eq.'coulomb') then element = twoERIs%baba%at(myI,myJ,myK,myL) elseIf(interaction.eq.'exchange') then - element = twoeris%baba%at(myI,myL,myK,myJ) + element = twoeris%baba%at(myI,myL,myK,myJ) elseIf(interaction.eq.'doublebar') then element = twoERIs%baba%at(myI,myJ,myK,myL) - twoERIs%baba%at(myI,myL,myK,myJ) endIf endIf endIf + !AAAB elseIf(myI.le.twoERIs%blockSize('alpha',1).and.myJ.le.twoERIs%blockSize('alpha',2).and. & myK.le.twoERIs%blockSize('alpha',3).and.myL.gt.twoERIs%blockSize('alpha',4)) then myL = myL - twoERIs%blockSize('alpha',4) @@ -9890,7 +9891,6 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen if(twoERIs%hasSpinBlock('aaab')) element = twoERIs%aaab%at(myI,myJ,myK,myL) elseif(interaction.eq.'exchange') then if(twoERIs%storageType.eq.'symm') then -! if(twoERIs%hasSpinBlock('aaab')) element = twoeris%aaab%at(myI,myL,myK,myJ) if(twoERIs%hasSpinBlock('aaab')) element = twoeris%aaab%at(myK,myJ,myI,myL) else if(twoERIs%hasSpinBlock('abaa')) element = twoeris%abaa%at(myI,myL,myK,myJ) @@ -9898,27 +9898,23 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen elseif(interaction.eq.'doublebar') then if(twoERIs%hasSpinBlock('aaab')) element = twoeris%aaab%at(myI,myJ,myK,myL) if(twoERIs%storageType.eq.'symm') then -! if(twoERIs%hasSpinBlock('aaab')) element = element - twoeris%aaab%at(myI,myL,myK,myJ) if(twoERIs%hasSpinBlock('aaab')) element = element - twoeris%aaab%at(myK,myJ,myI,myL) else if(twoERIs%hasSpinBlock('abaa')) element = element - twoeris%abaa%at(myI,myL,myK,myJ) endIf endIf + !AABA elseIf(myI.le.twoERIs%blockSize('alpha',1).and.myJ.le.twoERIs%blockSize('alpha',2).and. & myK.gt.twoERIs%blockSize('alpha',3).and.myL.le.twoERIs%blockSize('alpha',4)) then myK = myK - twoERIs%blockSize('alpha',3) if(twoERIs%storageType.eq.'symm') then if(interaction.eq.'coulomb') then -! if(twoERIs%hasSpinBlock('aaab')) element = conjg(twoERIs%aaab%at(myI,myJ,myK,myL)) - if(twoERIs%hasSpinBlock('aaab')) element = conjg(twoERIs%aaab%at(myI,myJ,myL,myK)) + if(twoERIs%hasSpinBlock('aaab')) element = conjg(twoERIs%aaab%at(myJ,myI,myL,myK)) elseif(interaction.eq.'exchange') then -! if(twoERIs%hasSpinBlock('aaab')) element = conjg(twoeris%aaab%at(myI,myL,myK,myJ)) - if(twoERIs%hasSpinBlock('aaab')) element = conjg(twoeris%aaab%at(myI,myL,myJ,myK)) + if(twoERIs%hasSpinBlock('aaab')) element = conjg(twoeris%aaab%at(myL,myI,myJ,myK)) elseif(interaction.eq.'doublebar') then -! if(twoERIs%hasSpinBlock('aaab')) element = conjg(twoeris%aaab%at(myI,myJ,myK,myL)) - & -! conjg(twoeris%aaab%at(myI,myL,myK,myJ)) - if(twoERIs%hasSpinBlock('aaab')) element = conjg(twoeris%aaab%at(myI,myJ,myL,myK)) - & - conjg(twoeris%aaab%at(myI,myL,myJ,myK)) + if(twoERIs%hasSpinBlock('aaab')) element = conjg(twoeris%aaab%at(myJ,myI,myL,myK)) - & + conjg(twoeris%aaab%at(myL,myI,myJ,myK)) endIf else if(interaction.eq.'coulomb') then @@ -9930,18 +9926,16 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen twoeris%aaab%at(myI,myL,myK,myJ) endIf endIf + !ABAA elseIf(myI.le.twoERIs%blockSize('alpha',1).and.myJ.gt.twoERIs%blockSize('alpha',2).and. & myK.le.twoERIs%blockSize('alpha',3).and.myL.le.twoERIs%blockSize('alpha',4)) then myJ = myJ - twoERIs%blockSize('alpha',2) if(twoERIs%storageType.eq.'symm') then if(interaction.eq.'coulomb') then -! if(twoERIs%hasSpinBlock('aaab')) element = twoERIs%aaab%at(myI,myJ,myK,myL) if(twoERIs%hasSpinBlock('aaab')) element = twoERIs%aaab%at(myK,myL,myI,myJ) elseif(interaction.eq.'exchange') then if(twoERIs%hasSpinBlock('aaab')) element = twoeris%aaab%at(myI,myL,myK,myJ) elseif(interaction.eq.'doublebar') then -! if(twoERIs%hasSpinBlock('aaab')) element = twoeris%aaab%at(myI,myJ,myK,myL) - & -! twoeris%aaab%at(myI,myL,myK,myJ) if(twoERIs%hasSpinBlock('aaab')) element = twoeris%aaab%at(myK,myL,myI,myJ) - & twoeris%aaab%at(myI,myL,myK,myJ) endIf @@ -9955,21 +9949,18 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen if(twoERIs%hasSpinBlock('aaab')) element = element - twoeris%aaab%at(myI,myL,myK,myJ) endIf endIf + !BAAA elseIf(myI.gt.twoERIs%blockSize('alpha',1).and.myJ.le.twoERIs%blockSize('alpha',2).and. & myK.le.twoERIs%blockSize('alpha',3).and.myL.le.twoERIs%blockSize('alpha',4)) then myI = myI - twoERIs%blockSize('alpha',1) if(twoERIs%storageType.eq.'symm') then if(interaction.eq.'coulomb') then -! if(twoERIs%hasSpinBlock('aaab')) element = conjg(twoERIs%aaab%at(myI,myJ,myK,myL)) - if(twoERIs%hasSpinBlock('aaab')) element = conjg(twoERIs%aaab%at(myL,myK,myI,myJ)) + if(twoERIs%hasSpinBlock('aaab')) element = conjg(twoERIs%aaab%at(myL,myK,myJ,myI)) elseif(interaction.eq.'exchange') then -! if(twoERIs%hasSpinBlock('aaab')) element = conjg(twoeris%aaab%at(myI,myL,myK,myJ)) - if(twoERIs%hasSpinBlock('aaab')) element = conjg(twoeris%aaab%at(myJ,myK,myI,myL)) + if(twoERIs%hasSpinBlock('aaab')) element = conjg(twoeris%aaab%at(myJ,myK,myL,myI)) elseif(interaction.eq.'doublebar') then -! if(twoERIs%hasSpinBlock('aaab')) element = conjg(twoeris%aaab%at(myI,myJ,myK,myL)) - & -! conjg(twoeris%aaab%at(myI,myL,myK,myJ)) - if(twoERIs%hasSpinBlock('aaab')) element = conjg(twoeris%aaab%at(myL,myK,myI,myJ)) - & - conjg(twoeris%aaab%at(myJ,myK,myI,myL)) + if(twoERIs%hasSpinBlock('aaab')) element = conjg(twoeris%aaab%at(myL,myK,myJ,myI)) - & + conjg(twoeris%aaab%at(myJ,myK,myL,myI)) endIf else if(interaction.eq.'coulomb') then @@ -9981,6 +9972,7 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen twoeris%baaa%at(myI,myL,myK,myJ) endIf endIf + !BBBA elseIf(myI.gt.twoERIs%blockSize('alpha',1).and.myJ.gt.twoERIs%blockSize('alpha',2).and. & myK.gt.twoERIs%blockSize('alpha',3).and.myL.le.twoERIs%blockSize('alpha',4)) then myI = myI - twoERIs%blockSize('alpha',1) @@ -9988,16 +9980,12 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen myK = myK - twoERIs%blockSize('alpha',3) if(twoERIs%storageType.eq.'symm') then if(interaction.eq.'coulomb') then -! if(twoERIs%hasSpinBlock('abbb')) element = conjg(twoERIs%abbb%at(myI,myJ,myK,myL)) - if(twoERIs%hasSpinBlock('abbb')) element = conjg(twoERIs%abbb%at(myK,myL,myJ,myI)) + if(twoERIs%hasSpinBlock('abbb')) element = conjg(twoERIs%abbb%at(myL,myK,myJ,myI)) elseif(interaction.eq.'exchange') then -! if(twoERIs%hasSpinBlock('abbb')) element = conjg(twoeris%abbb%at(myI,myL,myK,myJ)) - if(twoERIs%hasSpinBlock('abbb')) element = conjg(twoeris%abbb%at(myL,myI,myK,myJ)) + if(twoERIs%hasSpinBlock('abbb')) element = conjg(twoeris%abbb%at(myL,myI,myJ,myK)) elseif(interaction.eq.'doublebar') then -! if(twoERIs%hasSpinBlock('abbb')) element = conjg(twoeris%abbb%at(myI,myJ,myK,myL)) - & -! conjg(twoeris%abbb%at(myI,myL,myK,myJ)) - if(twoERIs%hasSpinBlock('abbb')) element = conjg(twoeris%abbb%at(myI,myJ,myK,myL)) - & - conjg(twoeris%abbb%at(myL,myI,myK,myJ)) + if(twoERIs%hasSpinBlock('abbb')) element = conjg(twoeris%abbb%at(myL,myK,myJ,myI)) - & + conjg(twoeris%abbb%at(myL,myI,myJ,myK)) endIf else if(interaction.eq.'coulomb') then @@ -10009,6 +9997,7 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen if(twoERIs%hasSpinBlock('babb')) element = element - twoeris%babb%at(myI,myL,myK,myJ) endIf endIf + !BBAB elseIf(myI.gt.twoERIs%blockSize('alpha',1).and.myJ.gt.twoERIs%blockSize('alpha',2).and. & myK.le.twoERIs%blockSize('alpha',3).and.myL.gt.twoERIs%blockSize('alpha',4)) then myI = myI - twoERIs%blockSize('alpha',1) @@ -10016,14 +10005,10 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen myL = myL - twoERIs%blockSize('alpha',4) if(twoERIs%storageType.eq.'symm') then if(interaction.eq.'coulomb') then -! if(twoERIs%hasSpinBlock('abbb')) element = twoERIs%abbb%at(myI,myJ,myK,myL) if(twoERIs%hasSpinBlock('abbb')) element = twoERIs%abbb%at(myK,myL,myI,myJ) elseif(interaction.eq.'exchange') then -! if(twoERIs%hasSpinBlock('abbb')) element = twoeris%abbb%at(myI,myL,myK,myJ) if(twoERIs%hasSpinBlock('abbb')) element = twoeris%abbb%at(myK,myJ,myI,myL) elseif(interaction.eq.'doublebar') then -! if(twoERIs%hasSpinBlock('abbb')) element = twoeris%abbb%at(myI,myJ,myK,myL) - & -! twoeris%abbb%at(myI,myL,myK,myJ) if(twoERIs%hasSpinBlock('abbb')) element = twoeris%abbb%at(myK,myL,myI,myJ) - & twoeris%abbb%at(myK,myJ,myI,myL) endIf @@ -10037,6 +10022,7 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen twoeris%bbab%at(myI,myL,myK,myJ) endIf endIf + !BABB elseIf(myI.gt.twoERIs%blockSize('alpha',1).and.myJ.le.twoERIs%blockSize('alpha',2).and. & myK.gt.twoERIs%blockSize('alpha',3).and.myL.gt.twoERIs%blockSize('alpha',4)) then myI = myI - twoERIs%blockSize('alpha',1) @@ -10044,16 +10030,12 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen myL = myL - twoERIs%blockSize('alpha',4) if(twoERIs%storageType.eq.'symm') then if(interaction.eq.'coulomb') then -! if(twoERIs%hasSpinBlock('abbb')) element = conjg(twoERIs%abbb%at(myI,myJ,myK,myL)) - if(twoERIs%hasSpinBlock('abbb')) element = conjg(twoERIs%abbb%at(myJ,myI,myK,myL)) + if(twoERIs%hasSpinBlock('abbb')) element = conjg(twoERIs%abbb%at(myJ,myI,myL,myK)) elseif(interaction.eq.'exchange') then -! if(twoERIs%hasSpinBlock('abbb')) element = conjg(twoeris%abbb%at(myI,myL,myK,myJ)) if(twoERIs%hasSpinBlock('abbb')) element = conjg(twoeris%abbb%at(myK,myJ,myL,myI)) elseif(interaction.eq.'doublebar') then -! if(twoERIs%hasSpinBlock('abbb')) element = conjg(twoeris%abbb%at(myI,myJ,myK,myL)) - & -! conjg(twoeris%abbb%at(myI,myL,myK,myJ)) - if(twoERIs%hasSpinBlock('abbb')) element = conjg(twoeris%abbb%at(myJ,myI,myK,myL)) - & - conjg(twoeris%abbb%at(myK,myJ,myL,myI)) + if(twoERIs%hasSpinBlock('abbb')) element = conjg(twoeris%abbb%at(myJ,myI,myL,myK)) - & + conjg(twoeris%abbb%at(myJ,myK,myL,myI)) endIf else if(interaction.eq.'coulomb') then @@ -10065,6 +10047,7 @@ function MQC_twoERIs_At(twoERIs,i,j,k,l,spinBlockIn,interactionIn) result(elemen if(twoERIs%hasSpinBlock('bbba')) element = twoeris%bbba%at(myI,myL,myK,myJ) endIf endIf + !ABBB elseIf(myI.le.twoERIs%blockSize('alpha',1).and.myJ.gt.twoERIs%blockSize('alpha',2).and. & myK.gt.twoERIs%blockSize('alpha',3).and.myL.gt.twoERIs%blockSize('alpha',4)) then myJ = myJ - twoERIs%blockSize('alpha',2) @@ -10796,8 +10779,8 @@ Subroutine Gen_Det_Str(IOut,IPrint,NBasisIn,NAlphaIn,NBetaIn, & subNum = 0 call mqc_deallocate_vector(determinants%nSubsAlpha) do while (.true.) - call determinants%nSubsAlpha%push(int(bin_coeff(nAlpha,subNum)*& - bin_coeff(nBasis-nAlpha,subNum))) + call determinants%nSubsAlpha%push(bin_coeff(nAlpha,subNum)*& + bin_coeff(nBasis-nAlpha,subNum)) if (determinants%nSubsAlpha%at(subNum+1).ne.0) then subNum = subNum + 1 else @@ -10808,8 +10791,8 @@ Subroutine Gen_Det_Str(IOut,IPrint,NBasisIn,NAlphaIn,NBetaIn, & subNum = 0 call mqc_deallocate_vector(determinants%nSubsBeta) do while (.true.) - call determinants%nSubsBeta%push(int(bin_coeff(nBeta,subNum)*& - bin_coeff(nBasis-nBeta,subNum))) + call determinants%nSubsBeta%push(bin_coeff(nBeta,subNum)*& + bin_coeff(nBasis-nBeta,subNum)) if (determinants%nSubsBeta%at(subNum+1).ne.0) then subNum = subNum + 1 else @@ -12938,8 +12921,8 @@ Subroutine TwoERI_Trans(IOut,IPrint,MO_Coeff,ERIs,MO_ERIs,Right_MOs) tmpR4TensorABAA,tmpR4TensorBAAA,tmpR4TensorBBBA,tmpR4TensorBBAB,tmpR4TensorBABB,tmpR4TensorABBB Type(MQC_R4Tensor)::tmpR4TensorLoc,tmpR4TensorLoc2 Type(MQC_R4Tensor)::tmpR4Tensor - Type(MQC_TwoERIs),Intent(In)::ERIs - Type(MQC_TwoERIs),Intent(Out)::MO_ERIs + Type(MQC_TwoERIs),Intent(InOut)::ERIs + Type(MQC_TwoERIs),Intent(InOut)::MO_ERIs Type(MQC_SCF_Integral),Intent(In)::MO_Coeff Type(MQC_SCF_Integral),Optional,Intent(In)::Right_MOs Type(MQC_SCF_Integral)::MO_Coeff_2 @@ -13036,17 +13019,17 @@ Subroutine TwoERI_Trans(IOut,IPrint,MO_Coeff,ERIs,MO_ERIs,Right_MOs) !AAAA Call tmpR4TensorAlpha%put(tmpR4TensorAlpha%at(P,Q,R,S) + & conjg(MO_Coeff%Alpha%at(Mu,P)) * MO_Coeff%Alpha%at(Nu,Q) * & - ERIs%alpha%at(Mu,Nu,Lambda,Sigma) * conjg(MO_Coeff_2%AlphaBeta%at(Lambda,R)) * & - MO_Coeff_2%AlphaBeta%at(Sigma,S),P,Q,R,S) + ERIs%alpha%at(Mu,Nu,Lambda,Sigma) * conjg(MO_Coeff_2%alphaBeta%at(Lambda,R)) * & + MO_Coeff_2%alphaBeta%at(Sigma,S),P,Q,R,S) Call tmpR4TensorAlpha%put(tmpR4TensorAlpha%at(P,Q,R,S) + & - conjg(MO_Coeff%AlphaBeta%at(Mu,P)) * MO_Coeff%AlphaBeta%at(Nu,Q) * & + conjg(MO_Coeff%alphaBeta%at(Mu,P)) * MO_Coeff%alphaBeta%at(Nu,Q) * & ERIs%alpha%at(Mu,Nu,Lambda,Sigma) * conjg(MO_Coeff_2%Alpha%at(Lambda,R)) * & MO_Coeff_2%Alpha%at(Sigma,S),P,Q,R,S) Call tmpR4TensorAlpha%put(tmpR4TensorAlpha%at(P,Q,R,S) + & - conjg(MO_Coeff%AlphaBeta%at(Mu,P)) * MO_Coeff%AlphaBeta%at(Nu,Q) * & - ERIs%alpha%at(Mu,Nu,Lambda,Sigma) * conjg(MO_Coeff_2%AlphaBeta%at(Lambda,R)) * & - MO_Coeff_2%AlphaBeta%at(Sigma,S),P,Q,R,S) - !BBBB + conjg(MO_Coeff%alphaBeta%at(Mu,P)) * MO_Coeff%alphaBeta%at(Nu,Q) * & + ERIs%alpha%at(Mu,Nu,Lambda,Sigma) * conjg(MO_Coeff_2%alphaBeta%at(Lambda,R)) * & + MO_Coeff_2%alphaBeta%at(Sigma,S),P,Q,R,S) + !BBBB Call tmpR4TensorBeta%put(tmpR4TensorBeta%at(P,Q,R,S) + & conjg(MO_Coeff%Beta%at(Mu,P)) * MO_Coeff%Beta%at(Nu,Q) * & ERIs%alpha%at(Mu,Nu,Lambda,Sigma) * conjg(MO_Coeff_2%BetaAlpha%at(Lambda,R)) * & @@ -13324,10 +13307,10 @@ Subroutine TwoERI_Trans(IOut,IPrint,MO_Coeff,ERIs,MO_ERIs,Right_MOs) Call X%init(NBasis,NBasis) Call Y%init(NBasis,NBasis) + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) !$OMP PARALLEL DEFAULT(NONE), & - !$OMP SHARED(ERIs,tmpR4TensorAlpha,MO_Coeff,MO_Coeff_2,nBasis), & - !$OMP FIRSTPRIVATE(X), & - !$OMP PRIVATE(P,Q,tmpR4TensorLoc) + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) !$OMP DO COLLAPSE(2) Do P = 1, NBasis @@ -13339,21 +13322,18 @@ Subroutine TwoERI_Trans(IOut,IPrint,MO_Coeff,ERIs,MO_ERIs,Right_MOs) EndDo !$OMP END DO !$OMP CRITICAL - tmpR4TensorAlpha = tmpR4TensorAlpha + tmpR4TensorLoc + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc !$OMP END CRITICAL !$OMP END PARALLEL - tmpR4TensorLoc2 = tmpR4TensorAlpha - call tmpR4TensorAlpha%init(nBasis,nBasis,nBasis,nBasis) !$OMP PARALLEL DEFAULT(NONE), & - !$OMP SHARED(tmpR4TensorAlpha,MO_Coeff,MO_Coeff_2,nBasis), & - !$OMP FIRSTPRIVATE(X,tmpR4TensorLoc2), & - !$OMP PRIVATE(R,S,tmpR4TensorLoc) + !$OMP SHARED(tmpR4Tensor,tmpR4TensorAlpha,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - X = tmpR4TensorLoc2%mat([0],[0],[R],[S]) + X = tmpR4Tensor%mat([0],[0],[R],[S]) X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%Alpha Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo @@ -13366,10 +13346,10 @@ Subroutine TwoERI_Trans(IOut,IPrint,MO_Coeff,ERIs,MO_ERIs,Right_MOs) If(intType.eq.'spin'.or.intType.eq.'general') then + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) !$OMP PARALLEL DEFAULT(NONE), & - !$OMP SHARED(ERIs,tmpR4TensorBeta,MO_Coeff,MO_Coeff_2,nBasis), & - !$OMP FIRSTPRIVATE(X), & - !$OMP PRIVATE(P,Q,tmpR4TensorLoc) + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) !$OMP DO COLLAPSE(2) Do P = 1, NBasis @@ -13381,21 +13361,18 @@ Subroutine TwoERI_Trans(IOut,IPrint,MO_Coeff,ERIs,MO_ERIs,Right_MOs) EndDo !$OMP END DO !$OMP CRITICAL - tmpR4TensorBeta = tmpR4TensorBeta + tmpR4TensorLoc + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc !$OMP END CRITICAL !$OMP END PARALLEL - tmpR4TensorLoc2 = tmpR4TensorBeta - call tmpR4TensorBeta%init(nBasis,nBasis,nBasis,nBasis) !$OMP PARALLEL DEFAULT(NONE), & - !$OMP SHARED(tmpR4TensorBeta,MO_Coeff,MO_Coeff_2,nBasis), & - !$OMP FIRSTPRIVATE(X,tmpR4TensorLoc2), & - !$OMP PRIVATE(R,S,tmpR4TensorLoc) + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBeta,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - X = tmpR4TensorLoc2%mat([0],[0],[R],[S]) + X = tmpR4Tensor%mat([0],[0],[R],[S]) X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%Beta Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo @@ -13406,10 +13383,10 @@ Subroutine TwoERI_Trans(IOut,IPrint,MO_Coeff,ERIs,MO_ERIs,Right_MOs) !$OMP END CRITICAL !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) !$OMP PARALLEL DEFAULT(NONE), & - !$OMP SHARED(ERIs,tmpR4TensorAlphaBeta,MO_Coeff,MO_Coeff_2,nBasis), & - !$OMP FIRSTPRIVATE(X), & - !$OMP PRIVATE(P,Q,tmpR4TensorLoc) + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) !$OMP DO COLLAPSE(2) Do P = 1, NBasis @@ -13421,21 +13398,18 @@ Subroutine TwoERI_Trans(IOut,IPrint,MO_Coeff,ERIs,MO_ERIs,Right_MOs) EndDo !$OMP END DO !$OMP CRITICAL - tmpR4TensorAlphaBeta = tmpR4TensorAlphaBeta + tmpR4TensorLoc + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc !$OMP END CRITICAL !$OMP END PARALLEL - tmpR4TensorLoc2 = tmpR4TensorAlphaBeta - call tmpR4TensorAlphaBeta%init(nBasis,nBasis,nBasis,nBasis) !$OMP PARALLEL DEFAULT(NONE), & - !$OMP SHARED(tmpR4TensorAlphaBeta,MO_Coeff,MO_Coeff_2,nBasis), & - !$OMP FIRSTPRIVATE(X,tmpR4TensorLoc2), & - !$OMP PRIVATE(R,S,tmpR4TensorLoc) + !$OMP SHARED(tmpR4Tensor,tmpR4TensorAlphaBeta,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - X = tmpR4TensorLoc2%mat([0],[0],[R],[S]) + X = tmpR4Tensor%mat([0],[0],[R],[S]) X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%Alpha Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo @@ -13448,10 +13422,10 @@ Subroutine TwoERI_Trans(IOut,IPrint,MO_Coeff,ERIs,MO_ERIs,Right_MOs) If(storage.eq.'full') then + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) !$OMP PARALLEL DEFAULT(NONE), & - !$OMP SHARED(ERIs,tmpR4TensorBetaAlpha,MO_Coeff,MO_Coeff_2,nBasis), & - !$OMP FIRSTPRIVATE(X), & - !$OMP PRIVATE(P,Q,tmpR4TensorLoc) + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) !$OMP DO COLLAPSE(2) Do P = 1, NBasis @@ -13463,21 +13437,18 @@ Subroutine TwoERI_Trans(IOut,IPrint,MO_Coeff,ERIs,MO_ERIs,Right_MOs) EndDo !$OMP END DO !$OMP CRITICAL - tmpR4TensorBetaAlpha = tmpR4TensorBetaAlpha + tmpR4TensorLoc + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc !$OMP END CRITICAL !$OMP END PARALLEL - tmpR4TensorLoc2 = tmpR4TensorBetaAlpha - call tmpR4TensorBetaAlpha%init(nBasis,nBasis,nBasis,nBasis) !$OMP PARALLEL DEFAULT(NONE), & - !$OMP SHARED(tmpR4TensorBetaAlpha,MO_Coeff,MO_Coeff_2,nBasis), & - !$OMP FIRSTPRIVATE(X,tmpR4TensorLoc2), & - !$OMP PRIVATE(R,S,tmpR4TensorLoc) + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBetaAlpha,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - X = tmpR4TensorLoc2%mat([0],[0],[R],[S]) + X = tmpR4Tensor%mat([0],[0],[R],[S]) X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%Beta Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo @@ -13494,2021 +13465,2259 @@ Subroutine TwoERI_Trans(IOut,IPrint,MO_Coeff,ERIs,MO_ERIs,Right_MOs) If(intType.eq.'general') then - call tmpR4Tensor%init(NBasis,NBasis,NBasis,NBasis) - ! AAAA + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%alphaBeta).dot.X - X = Y.dot.MO_Coeff%alpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%alphaBeta).dot.X.dot.MO_Coeff%alphaBeta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorAlpha,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%alphaBeta).dot.X - X = Y.dot.MO_Coeff%alpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorAlpha%put(tmpR4TensorAlpha%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%alpha).dot.X.dot.MO_Coeff%alpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorAlpha = tmpR4TensorAlpha + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%alpha).dot.X - X = Y.dot.MO_Coeff%alphaBeta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%alpha).dot.X.dot.MO_Coeff%alpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorAlpha,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%alpha).dot.X - X = Y.dot.MO_Coeff%alphaBeta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorAlpha%put(tmpR4TensorAlpha%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%alphaBeta).dot.X.dot.MO_Coeff%alphaBeta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorAlpha = tmpR4TensorAlpha + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%alphaBeta).dot.X - X = Y.dot.MO_Coeff%alphaBeta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%alphaBeta).dot.X.dot.MO_Coeff%alphaBeta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorAlpha,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%alphaBeta).dot.X - X = Y.dot.MO_Coeff%alphaBeta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorAlpha%put(tmpR4TensorAlpha%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%alphaBeta).dot.X.dot.MO_Coeff%alphaBeta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo - + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorAlpha = tmpR4TensorAlpha + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL +! ! BBBB + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%Beta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBeta,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%Beta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBeta%put(tmpR4TensorBeta%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBeta = tmpR4TensorBeta + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBeta,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%beta).dot.X - X = Y.dot.MO_Coeff%betaAlpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBeta%put(tmpR4TensorBeta%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%beta).dot.X.dot.MO_Coeff%beta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBeta = tmpR4TensorBeta + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%betaAlpha).dot.X - X = Y.dot.MO_Coeff%betaAlpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%betaAlpha).dot.X.dot.MO_Coeff%betaAlpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBeta,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%betaAlpha).dot.X - X = Y.dot.MO_Coeff%betaAlpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBeta%put(tmpR4TensorBeta%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%betaAlpha).dot.X.dot.MO_Coeff%betaAlpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBeta = tmpR4TensorBeta + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL ! AABB + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%betaAlpha).dot.X - X = Y.dot.MO_Coeff%betaAlpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%betaAlpha).dot.X.dot.MO_Coeff%betaAlpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorAlphaBeta,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%alpha).dot.X - X = Y.dot.MO_Coeff%alpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensoralphaBeta%put(tmpR4TensoralphaBeta%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%alpha).dot.X.dot.MO_Coeff%alpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorAlphaBeta = tmpR4TensorAlphaBeta + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%Beta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo - Do R = 1, NBasis - Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%alphaBeta).dot.X - X = Y.dot.MO_Coeff%alphaBeta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensoralphaBeta%put(tmpR4TensoralphaBeta%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo - EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorAlphaBeta,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) + Do R = 1, NBasis + Do S = 1, NBasis + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%alphaBeta).dot.X.dot.MO_Coeff%alphaBeta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) + EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorAlphaBeta = tmpR4TensorAlphaBeta + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%betaAlpha).dot.X - X = Y.dot.MO_Coeff%betaAlpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%betaAlpha).dot.X.dot.MO_Coeff%betaAlpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorAlphaBeta,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%alphaBeta).dot.X - X = Y.dot.MO_Coeff%alphaBeta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensoralphaBeta%put(tmpR4TensoralphaBeta%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%alphaBeta).dot.X.dot.MO_Coeff%alphaBeta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorAlphaBeta = tmpR4TensorAlphaBeta + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL ! BBAA If(storage.eq.'full') then + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%alphaBeta).dot.X - X = Y.dot.MO_Coeff%alphaBeta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%alphaBeta).dot.X.dot.MO_Coeff%alphaBeta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBetaAlpha,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%beta).dot.X - X = Y.dot.MO_Coeff%beta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBetaAlpha%put(tmpR4TensorBetaAlpha%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%beta).dot.X.dot.MO_Coeff%beta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBetaAlpha = tmpR4TensorBetaAlpha + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%alpha).dot.X - X = Y.dot.MO_Coeff%alpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%alpha).dot.X.dot.MO_Coeff%alpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBetaAlpha,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%betaAlpha).dot.X - X = Y.dot.MO_Coeff%betaAlpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBetaAlpha%put(tmpR4TensorBetaAlpha%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%betaAlpha).dot.X.dot.MO_Coeff%betaAlpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBetaAlpha = tmpR4TensorBetaAlpha + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%alphaBeta).dot.X - X = Y.dot.MO_Coeff%alphaBeta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%alphaBeta).dot.X.dot.MO_Coeff%alphaBeta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBetaAlpha,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%betaAlpha).dot.X - X = Y.dot.MO_Coeff%betaAlpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBetaAlpha%put(tmpR4TensorBetaAlpha%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%betaAlpha).dot.X.dot.MO_Coeff%betaAlpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBetaAlpha = tmpR4TensorBetaAlpha + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL EndIf ! ABAB + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Alpha).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorABAB,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Alpha).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorABAB%put(tmpR4TensorABAB%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorABAB = tmpR4TensorABAB + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%AlphaBeta).dot.X - X = Y.dot.MO_Coeff%Beta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%AlphaBeta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorABAB,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Alpha).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorABAB%put(tmpR4TensorABAB%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorABAB = tmpR4TensorABAB + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Alpha).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorABAB,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%AlphaBeta).dot.X - X = Y.dot.MO_Coeff%Beta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorABAB%put(tmpR4TensorABAB%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%AlphaBeta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorABAB = tmpR4TensorABAB + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%AlphaBeta).dot.X - X = Y.dot.MO_Coeff%Beta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%AlphaBeta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorABAB,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%alphaBeta).dot.X - X = Y.dot.MO_Coeff%Beta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorABAB%put(tmpR4TensorABAB%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%alphaBeta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorABAB = tmpR4TensorABAB + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL ! BABA If(storage.eq.'full') then + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%Alpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%Alpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBABA,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%Alpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBABA%put(tmpR4TensorBABA%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%Alpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBABA = tmpR4TensorBABA + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBABA,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%Alpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBABA%put(tmpR4TensorBABA%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%Alpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBABA = tmpR4TensorBABA + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%Alpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%Alpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBABA,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%AlphaBeta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBABA%put(tmpR4TensorBABA%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%AlphaBeta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBABA = tmpR4TensorBABA + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%AlphaBeta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%AlphaBeta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBABA,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBABA%put(tmpR4TensorBABA%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBABA = tmpR4TensorBABA + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL EndIf ! ABBA + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%Alpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%Alpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorABBA,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Alpha).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorABBA%put(tmpR4TensorABBA%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorABBA = tmpR4TensorABBA + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis - Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%AlphaBeta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + Do Q = 1, NBasis + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%AlphaBeta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorABBA,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Alpha).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorABBA%put(tmpR4TensorABBA%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorABBA = tmpR4TensorABBA + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%Alpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%Alpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorABBA,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%AlphaBeta).dot.X - X = Y.dot.MO_Coeff%Beta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorABBA%put(tmpR4TensorABBA%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%AlphaBeta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorABBA = tmpR4TensorABBA + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%AlphaBeta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%AlphaBeta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorABBA,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%alphaBeta).dot.X - X = Y.dot.MO_Coeff%Beta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorABBA%put(tmpR4TensorABBA%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%alphaBeta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorABBA = tmpR4TensorABBA + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL ! BAAB If(storage.eq.'full') then + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Alpha).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBAAB,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%Alpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBAAB%put(tmpR4TensorBAAB%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%Alpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBAAB = tmpR4TensorBAAB + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%AlphaBeta).dot.X - X = Y.dot.MO_Coeff%Beta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%AlphaBeta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBAAB,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%Alpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBAAB%put(tmpR4TensorBAAB%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%Alpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBAAB = tmpR4TensorBAAB + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Alpha).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBAAB,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%AlphaBeta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBAAB%put(tmpR4TensorBAAB%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%AlphaBeta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBAAB = tmpR4TensorBAAB + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%AlphaBeta).dot.X - X = Y.dot.MO_Coeff%Beta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%AlphaBeta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBAAB,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%AlphaBeta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBAAB%put(tmpR4TensorBAAB%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%AlphaBeta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBAAB = tmpR4TensorBAAB + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL EndIf ! AAAB + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Alpha).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorAAAB,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Alpha).dot.X - X = Y.dot.MO_Coeff%Alpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorAAAB%put(tmpR4TensorAAAB%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%Alpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorAAAB = tmpR4TensorAAAB + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%AlphaBeta).dot.X - X = Y.dot.MO_Coeff%Beta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%AlphaBeta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorAAAB,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Alpha).dot.X - X = Y.dot.MO_Coeff%Alpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorAAAB%put(tmpR4TensorAAAB%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%Alpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorAAAB = tmpR4TensorAAAB + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Alpha).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorAAAB,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%AlphaBeta).dot.X - X = Y.dot.MO_Coeff%AlphaBeta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorAAAB%put(tmpR4TensorAAAB%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%AlphaBeta).dot.X.dot.MO_Coeff%AlphaBeta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorAAAB = tmpR4TensorAAAB + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%AlphaBeta).dot.X - X = Y.dot.MO_Coeff%Beta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%AlphaBeta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorAAAB,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%alphaBeta).dot.X - X = Y.dot.MO_Coeff%alphaBeta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorAAAB%put(tmpR4TensorAAAB%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%alphaBeta).dot.X.dot.MO_Coeff%alphaBeta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorAAAB = tmpR4TensorAAAB + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL ! AABA If(storage.eq.'full') then + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%Alpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%Alpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorAABA,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Alpha).dot.X - X = Y.dot.MO_Coeff%Alpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorAABA%put(tmpR4TensorAABA%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%Alpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorAABA = tmpR4TensorAABA + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%AlphaBeta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%AlphaBeta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorAABA,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Alpha).dot.X - X = Y.dot.MO_Coeff%Alpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorAABA%put(tmpR4TensorAABA%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%Alpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorAABA = tmpR4TensorAABA + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%Alpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%Alpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorAABA,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%AlphaBeta).dot.X - X = Y.dot.MO_Coeff%AlphaBeta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorAABA%put(tmpR4TensorAABA%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%AlphaBeta).dot.X.dot.MO_Coeff%AlphaBeta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorAABA = tmpR4TensorAABA + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%AlphaBeta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%AlphaBeta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorAABA,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%alphaBeta).dot.X - X = Y.dot.MO_Coeff%alphaBeta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorAABA%put(tmpR4TensorAABA%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%alphaBeta).dot.X.dot.MO_Coeff%alphaBeta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorAABA = tmpR4TensorAABA + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL EndIf ! ABAA If(storage.eq.'full') then + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Alpha).dot.X - X = Y.dot.MO_Coeff%Alpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%Alpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorABAA,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Alpha).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorABAA%put(tmpR4TensorABAA%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorABAA = tmpR4TensorABAA + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%AlphaBeta).dot.X - X = Y.dot.MO_Coeff%AlphaBeta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%AlphaBeta).dot.X.dot.MO_Coeff%AlphaBeta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorABAA,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Alpha).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorABAA%put(tmpR4TensorABAA%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorABAA = tmpR4TensorABAA + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Alpha).dot.X - X = Y.dot.MO_Coeff%Alpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%Alpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorABAA,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%AlphaBeta).dot.X - X = Y.dot.MO_Coeff%Beta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorABAA%put(tmpR4TensorABAA%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%AlphaBeta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorABAA = tmpR4TensorABAA + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%AlphaBeta).dot.X - X = Y.dot.MO_Coeff%AlphaBeta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%AlphaBeta).dot.X.dot.MO_Coeff%AlphaBeta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorABAA,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%alphaBeta).dot.X - X = Y.dot.MO_Coeff%Beta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorABAA%put(tmpR4TensorABAA%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%alphaBeta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorABAA = tmpR4TensorABAA + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL EndIf ! BAAA If(storage.eq.'full') then + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Alpha).dot.X - X = Y.dot.MO_Coeff%Alpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%Alpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBAAA,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%Alpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBAAA%put(tmpR4TensorBAAA%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%Alpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBAAA = tmpR4TensorBAAA + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%AlphaBeta).dot.X - X = Y.dot.MO_Coeff%AlphaBeta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%AlphaBeta).dot.X.dot.MO_Coeff%AlphaBeta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBAAA,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%Alpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBAAA%put(tmpR4TensorBAAA%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%Alpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBAAA = tmpR4TensorBAAA + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Alpha).dot.X - X = Y.dot.MO_Coeff%Alpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%Alpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBAAA,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%AlphaBeta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBAAA%put(tmpR4TensorBAAA%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%AlphaBeta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBAAA = tmpR4TensorBAAA + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%AlphaBeta).dot.X - X = Y.dot.MO_Coeff%AlphaBeta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%AlphaBeta).dot.X.dot.MO_Coeff%AlphaBeta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBAAA,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%AlphaBeta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBAAA%put(tmpR4TensorBAAA%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%AlphaBeta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBAAA = tmpR4TensorBAAA + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL EndIf ! BBBA If(storage.eq.'full') then + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%Alpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%Alpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBBBA,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBBBA%put(tmpR4TensorBBBA%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBBBA = tmpR4TensorBBBA + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%AlphaBeta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%AlphaBeta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBBBA,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBBBA%put(tmpR4TensorBBBA%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBBBA = tmpR4TensorBBBA + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%Alpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%Alpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBBBA,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%Beta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBBBA%put(tmpR4TensorBBBA%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBBBA = tmpR4TensorBBBA + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%AlphaBeta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%AlphaBeta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBBBA,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%Beta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBBBA%put(tmpR4TensorBBBA%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBBBA = tmpR4TensorBBBA + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL EndIf ! BBAB If(storage.eq.'full') then + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Alpha).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBBAB,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBBAB%put(tmpR4TensorBBAB%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBBAB = tmpR4TensorBBAB + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%AlphaBeta).dot.X - X = Y.dot.MO_Coeff%Beta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%AlphaBeta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBBAB,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBBAB%put(tmpR4TensorBBAB%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBBAB = tmpR4TensorBBAB + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Alpha).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBBAB,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%Beta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBBAB%put(tmpR4TensorBBAB%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBBAB = tmpR4TensorBBAB + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%AlphaBeta).dot.X - X = Y.dot.MO_Coeff%Beta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%AlphaBeta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBBAB,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%Beta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBBAB%put(tmpR4TensorBBAB%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBBAB = tmpR4TensorBBAB + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL EndIf ! BABB If(storage.eq.'full') then + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBABB,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%Alpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBABB%put(tmpR4TensorBABB%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%Alpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBABB = tmpR4TensorBABB + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%Beta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBABB,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%Alpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBABB%put(tmpR4TensorBABB%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%Alpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBABB = tmpR4TensorBABB + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBABB,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%AlphaBeta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBABB%put(tmpR4TensorBABB%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%AlphaBeta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBABB = tmpR4TensorBABB + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%Beta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorBABB,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%AlphaBeta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorBABB%put(tmpR4TensorBABB%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%AlphaBeta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorBABB = tmpR4TensorBABB + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL EndIf ! ABBB + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorABBB,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Alpha).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorABBB%put(tmpR4TensorABBB%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorABBB = tmpR4TensorABBB + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%Beta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorABBB,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Alpha).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorABBB%put(tmpR4TensorABBB%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%Alpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorABBB = tmpR4TensorABBB + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%BetaAlpha).dot.X - X = Y.dot.MO_Coeff%BetaAlpha - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%BetaAlpha).dot.X.dot.MO_Coeff%BetaAlpha + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorABBB,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%AlphaBeta).dot.X - X = Y.dot.MO_Coeff%Beta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorABBB%put(tmpR4TensorABBB%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%AlphaBeta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorABBB = tmpR4TensorABBB + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + call tmpR4Tensor%init(nBasis,nBasis,nBasis,nBasis) + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(ERIs,tmpR4Tensor,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,P,Q,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do P = 1, NBasis Do Q = 1, NBasis - Do R = 1, NBasis - Do S = 1, NBasis - Call X%put(ERIs%alpha%at(P,Q,R,S),R,S) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%Beta).dot.X - X = Y.dot.MO_Coeff%Beta - Do R = 1, NBasis - Do S = 1, NBasis - Call tmpR4Tensor%put(X%at(R,S),P,Q,R,S) - EndDo - EndDo + X = ERIs%alpha%mat([P],[Q],[0],[0]) + X = Dagger(MO_Coeff_2%Beta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[P],[Q],[0],[0]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4Tensor = tmpR4Tensor + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE), & + !$OMP SHARED(tmpR4Tensor,tmpR4TensorABBB,MO_Coeff,MO_Coeff_2,nBasis), & + !$OMP PRIVATE(X,R,S,tmpR4TensorLoc) + call tmpR4TensorLoc%init(nBasis,nBasis,nBasis,nBasis) + !$OMP DO COLLAPSE(2) Do R = 1, NBasis Do S = 1, NBasis - Do P = 1, NBasis - Do Q = 1, NBasis - Call X%put(tmpR4Tensor%at(P,Q,R,S),P,Q) - EndDo - EndDo - Y = Dagger(MO_Coeff_2%alphaBeta).dot.X - X = Y.dot.MO_Coeff%Beta - Do P = 1, NBasis - Do Q = 1, NBasis - Call tmpR4TensorABBB%put(tmpR4TensorABBB%at(P,Q,R,S)+X%at(P,Q),P,Q,R,S) - EndDo - EndDo + X = tmpR4Tensor%mat([0],[0],[R],[S]) + X = Dagger(MO_Coeff_2%alphaBeta).dot.X.dot.MO_Coeff%Beta + Call tmpR4TensorLoc%mput(X,[0],[0],[R],[S]) EndDo EndDo + !$OMP END DO + !$OMP CRITICAL + tmpR4TensorABBB = tmpR4TensorABBB + tmpR4TensorLoc + !$OMP END CRITICAL + !$OMP END PARALLEL EndIf From ada292e53b7a4ec980e65936cacec8cb3c0b0222 Mon Sep 17 00:00:00 2001 From: Lee Thompson Date: Thu, 2 Jun 2022 17:45:08 -0400 Subject: [PATCH 05/10] Added routine to return character array of element symbols --- src/mqc_molecule.F03 | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/src/mqc_molecule.F03 b/src/mqc_molecule.F03 index 4e56aff7..bd0228bd 100644 --- a/src/mqc_molecule.F03 +++ b/src/mqc_molecule.F03 @@ -271,6 +271,28 @@ Function MQC_Get_Nuclear_Repulsion(Molecule_Info) Result(Vnn) End Function MQC_Get_Nuclear_Repulsion ! ! +! PROCEDURE MQC_GET_NUCLEAR_SYMBOLS + Function MQC_Get_Nuclear_Symbols(Molecule_Info) Result(Symbols) +! +! This function returns the nuclear repulsion energy given the +! molecular data object. +! +! Lee M. Thompson, 2017. +! +! + Implicit None + Class(MQC_Molecule_Data),Intent(In)::Molecule_Info + Character(len=2),Dimension(:),allocatable::Symbols + Integer::I +! + Allocate(Symbols(int(Molecule_Info%NAtoms))) + Do I = 1,int(Molecule_Info%NAtoms) + Symbols(I) = mqc_element_symbol(molecule_info%atomic_numbers%at(i)) + EndDo +! + End Function MQC_Get_Nuclear_Symbols +! +! ! PROCEDURE MQC_PRINT_NUCLEAR_GEOMETRY subroutine mqc_print_nuclear_geometry(molecule_info,iOut,bohr) ! From 9374270cca816c5ff001a30db6b0b1c098906703 Mon Sep 17 00:00:00 2001 From: Lee Thompson Date: Tue, 7 Jun 2022 14:40:05 -0400 Subject: [PATCH 06/10] Added conjugation to mqc_scf_integral type --- src/mqc_est.F03 | 66 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 66 insertions(+) diff --git a/src/mqc_est.F03 b/src/mqc_est.F03 index bf187e7c..69477825 100644 --- a/src/mqc_est.F03 +++ b/src/mqc_est.F03 @@ -279,6 +279,11 @@ Module MQC_EST Module Procedure MQC_Matrix_UndoSpinBlockGHF_Integral End Interface ! +!> \brief Returns the complex conjugate<\b> + Interface Conjg + Module Procedure mqc_scf_integral_conjg + End Interface +! !---------------------------------------------------------------- ! | ! OPERATOR INTERFACES | @@ -8792,6 +8797,67 @@ function mqc_scf_integral_determinant(integral) result(determinant) end function mqc_scf_integral_determinant ! ! +! PROCEDURE MQC_SCF_Integral_Conjg +! +!> \brief MQC_SCF_Integral_Conjg is a function that returns the complex conjugate +!> of an MQC SCF integral +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_SCF_Integral_Conjg is a function that returns the complex conjugate of an MQC SCF +!> integral. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] Integral +!> \verbatim +!> Integral is Class(MQC_SCF_Integral) +!> The MQC integral which will be evaluated. +!> \endverbatim +! +! Authors: +! ======== +!> \author L. M. Thompson +!> \date 2022 +! + function mqc_scf_integral_conjg(integral,label) result(integralOut) +! + implicit none + class(mqc_scf_integral),intent(in)::integral + Character(Len=*),optional,intent(in)::label + type(mqc_scf_integral)::integralOut + type(mqc_matrix)::tmpMat1,tmpMat2,tmpMat3,tmpMat4 + Character(Len=64)::myLabel +! + if(present(label)) then + call string_change_case(label,'l',myLabel) + else + myLabel = '' + endIf + + if(integral%type().eq.'space') then + tmpMat1 = conjg(integral%getblock('alpha')) + call mqc_integral_allocate(integralOut,myLabel,'space',tmpMat1) + elseIf(integral%type().eq.'spin') then + tmpMat1 = conjg(integral%getblock('alpha')) + tmpMat2 = conjg(integral%getblock('beta')) + call mqc_integral_allocate(integralOut,myLabel,'spin',tmpMat1,tmpMat2) + elseIf(integral%type().eq.'general') then + tmpMat1 = conjg(integral%getblock('alpha')) + tmpMat2 = conjg(integral%getblock('beta')) + tmpMat3 = conjg(integral%getblock('alpha-beta')) + tmpMat4 = conjg(integral%getblock('beta-alpha')) + call mqc_integral_allocate(integralOut,myLabel,'general',tmpMat1,tmpMat2,tmpMat3,tmpMat4) + endIf +! + end function mqc_scf_integral_conjg +! +! ! PROCEDURE MQC_Integral_Set_Energy_List ! !> \brief MQC_Integral_Set_Energy_List is a subroutine that sets the energy From 38f9fe9ef78d9c6f5cec97bc136126f52ff60439 Mon Sep 17 00:00:00 2001 From: Lee Thompson Date: Thu, 9 Jun 2022 16:03:50 -0400 Subject: [PATCH 07/10] SCF integral power subroutine added --- src/mqc_est.F03 | 104 +++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 103 insertions(+), 1 deletion(-) diff --git a/src/mqc_est.F03 b/src/mqc_est.F03 index 69477825..33183318 100644 --- a/src/mqc_est.F03 +++ b/src/mqc_est.F03 @@ -120,6 +120,8 @@ Module MQC_EST Procedure, Public::det => mqc_scf_integral_determinant !> \brief Returns the norm of an MQC_SCF_Integral object Procedure, Public::norm => mqc_integral_norm +!> \brief Returns an MQC_SCF_Integral object raised to a power + Procedure, Public::power => mqc_scf_integral_power Procedure, Public::at => mqc_integral_at Procedure, Public::setEList => mqc_integral_set_energy_list Procedure, Public::getEList => mqc_integral_get_energy_list @@ -279,7 +281,7 @@ Module MQC_EST Module Procedure MQC_Matrix_UndoSpinBlockGHF_Integral End Interface ! -!> \brief Returns the complex conjugate<\b> +!> \brief Returns the complex conjugate Interface Conjg Module Procedure mqc_scf_integral_conjg End Interface @@ -8819,6 +8821,12 @@ end function mqc_scf_integral_determinant !> Integral is Class(MQC_SCF_Integral) !> The MQC integral which will be evaluated. !> \endverbatim +!> +!> \param[in] Label +!> \verbatim +!> Label is Character(Len=*),optional +!> Label for the output SCF integral. +!> \endverbatim ! ! Authors: ! ======== @@ -9060,6 +9068,100 @@ subroutine mqc_scf_eigenvalues_power(eigenvalues,power) end subroutine mqc_scf_eigenvalues_power ! ! +! PROCEDURE MQC_SCF_Integral_Power +! +!> \brief MQC_SCF_Integral_Power is a function that returns the value of +!> an MQC integral variable raised to a power +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_SCF_Integral_Power is a function that returns an MQC +!> integral variable raised to a power. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in,out] Integral +!> \verbatim +!> Integral is Class(MQC_SCF_Integral) +!> The name of the MQC_SCF_Integral variable. +!> \endverbatim +!> +!> \param[in] Power +!> \verbatim +!> Power is Class(*) +!> The power to raise the MQC integral. +!> \endverbatim +!> +!> \param[in] Label +!> \verbatim +!> Label is Character(Len=*),optional +!> Label for the output SCF integral. +!> \endverbatim +! +! Authors: +! ======== +!> \author L. M. Thompson +!> \date 2022 +! + subroutine mqc_scf_integral_power(integral,power,label) +! + implicit none + class(mqc_scf_integral),intent(inOut)::integral + class(*)::power + type(mqc_scalar)::scalar + type(mqc_matrix)::tmpMat1,tmpMat2,tmpMat3,tmpMat4 + integer(kind=int64)::nDimAlpha=0 + Character(Len=*),optional,intent(in)::label + Character(Len=64)::myLabel +! + select type(power) + type is (integer) + scalar = power + type is (real) + scalar = power + type is (complex) + scalar = power + type is (mqc_scalar) + scalar = power + class default + call mqc_error_I('power type not defined in MQC_SCF_Integral_Power',6) + end select + + if(present(label)) then + call string_change_case(label,'l',myLabel) + else + myLabel = '' + endIf + + if(integral%type().eq.'space') then + tmpMat1 = integral%getblock('alpha') + call tmpMat1%power(scalar) + call mqc_integral_allocate(integral,myLabel,'space',tmpMat1) + elseIf(integral%type().eq.'spin') then + tmpMat1 = integral%getblock('alpha') + tmpMat2 = integral%getblock('beta') + call tmpMat1%power(scalar) + call tmpMat2%power(scalar) + call mqc_integral_allocate(integral,myLabel,'spin',tmpMat1,tmpMat2) + elseIf(integral%type().eq.'general') then + nDimAlpha = integral%blockSize('alpha') + tmpMat1 = integral%getblock('full') + call tmpMat1%power(scalar) + tmpMat2 = tmpMat1%mat([NDimAlpha+1,-1],[NDimAlpha+1,-1]) + tmpMat3 = tmpMat1%mat([NDimAlpha+1,-1],[1,NDimAlpha]) + tmpMat4 = tmpMat1%mat([1,NDimAlpha],[NDimAlpha+1,-1]) + tmpMat1 = tmpMat1%mat([1,NDimAlpha],[1,NDimAlpha]) + call mqc_integral_allocate(integral,myLabel,'general',tmpMat1,tmpMat2,tmpMat3,tmpMat4) + endIf +! + end subroutine mqc_scf_integral_power +! +! ! PROCEDURE MQC_twoERIs_At ! !> \brief MQC_twoERIs_At is a function that returns the value of an element From 0d44c184279449cdbe12ec5a9301f1eec5675e33 Mon Sep 17 00:00:00 2001 From: Lee Thompson <32308301+leethomo86@users.noreply.github.com> Date: Tue, 14 Jun 2022 17:18:08 -0500 Subject: [PATCH 08/10] Update README.md --- README.md | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/README.md b/README.md index a22cc6a5..2c7031ef 100644 --- a/README.md +++ b/README.md @@ -15,17 +15,14 @@ ** ** ** The Merced Quantum Chemistry Package ** ** (MQCPack) ** - ** Development Version ** - ** Based On: ** - ** Development Version 0.1 ** + ** Version 22.06.1 ** + ** Released: June 14, 2022 ** ** ** ** ** ** Written By: ** ** Lee M. Thompson, Xianghai Sheng, Andrew D. Mahler, Dave Mullally ** ** and Hrant P. Hratchian ** ** ** - ** Version 20.05.1 Completed ** - ** May 22, 2020 ** ** ** ************************************************************************** ************************************************************************** From 0285ac1562f2072e7ad70e2f9930c74ec694e8d0 Mon Sep 17 00:00:00 2001 From: Lee Thompson Date: Tue, 14 Jun 2022 18:25:18 -0400 Subject: [PATCH 09/10] Added scalar vector sum (order no longer matters) and small fix to matrix power routine. --- src/mqc_algebra.F03 | 108 +++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 107 insertions(+), 1 deletion(-) diff --git a/src/mqc_algebra.F03 b/src/mqc_algebra.F03 index afea7f8c..4e24e308 100644 --- a/src/mqc_algebra.F03 +++ b/src/mqc_algebra.F03 @@ -801,6 +801,7 @@ Module MQC_Algebra Interface Operator (+) Module Procedure MQC_VectorVectorSum Module Procedure MQC_ScalarVectorSum + Module Procedure MQC_VectorScalarSum End Interface ! ! Documentation in scalar operator section @@ -8702,6 +8703,110 @@ Function MQC_ScalarVectorSum(ScalarIn,VectorIn) End Function MQC_ScalarVectorSum ! ! +! PROCEDURE MQC_VectorScalarSum +! +!> \brief MQC_VectorScalarSum is a function that adds an MQC scalar to all +!> elements of an MQC vector +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> MQC_VectorScalarSum is a function that adds an MQC scalar to all elements of +!> an MQC vector. +!> +!> \endverbatim +! +! Arguments: +! ========== +!> \param[in] VectorIn +!> \verbatim +!> VectorIn is Type(MQC_Vector) +!> The MQC vector with elements to sum with ScalarIn. +!> \endverbatim +!> +!> \param[in] ScalarIn +!> \verbatim +!> ScalarIn is Type(MQC_Scalar) +!> The MQC scalar to add to the MQC vector. +!> \endverbatim +! +! Authors: +! ======== +!> \author L. M. Thompson +!> \date 2022 +! + Function MQC_VectorScalarSum(VectorIn,ScalarIn) +! + Implicit None + Type(MQC_Vector)::MQC_VectorScalarSum + Type(MQC_Vector),Intent(In)::VectorIn + Type(MQC_Scalar),Intent(In)::ScalarIn + Integer(kind=int64) :: VectorLen + + VectorLen = MQC_Length_Vector(VectorIn) + + If(VectorIn%Data_type.eq.'Real') then + If(ScalarIn%Data_type.eq.'Real') then + Call MQC_Allocate_Vector(VectorLen,MQC_VectorScalarSum,'Real') + MQC_VectorScalarSum%vecr = VectorIn%vecr + ScalarIn%scar + MQC_VectorScalarSum%column = VectorIn%column + ElseIf(ScalarIn%Data_type.eq.'Integer') then + Call MQC_Allocate_Vector(VectorLen,MQC_VectorScalarSum,'Real') + MQC_VectorScalarSum%vecr = VectorIn%vecr + ScalarIn%scai + MQC_VectorScalarSum%column = VectorIn%column + ElseIf(ScalarIn%Data_type.eq.'Complex') then + Call MQC_Allocate_Vector(VectorLen,MQC_VectorScalarSum,'Complex') + MQC_VectorScalarSum%vecc = VectorIn%vecr + ScalarIn%scac + MQC_VectorScalarSum%column = VectorIn%column + Else + Call MQC_Error_A('ScalarIn type unspecified in MQC_VectorScalarSum', 6,& + 'ScalarIn%Data_type', ScalarIn%Data_type ) + EndIf + ElseIf(VectorIn%Data_type.eq.'Integer') then + If(ScalarIn%Data_type.eq.'Real') then + Call MQC_Allocate_Vector(VectorLen,MQC_VectorScalarSum,'Real') + MQC_VectorScalarSum%vecr = VectorIn%veci + ScalarIn%scar + MQC_VectorScalarSum%column = VectorIn%column + ElseIf(ScalarIn%Data_type.eq.'Integer') then + Call MQC_Allocate_Vector(VectorLen,MQC_VectorScalarSum,'Integer') + MQC_VectorScalarSum%veci = VectorIn%veci + ScalarIn%scai + MQC_VectorScalarSum%column = VectorIn%column + ElseIf(ScalarIn%Data_type.eq.'Complex') then + Call MQC_Allocate_Vector(VectorLen,MQC_VectorScalarSum,'Complex') + MQC_VectorScalarSum%vecc = VectorIn%veci + ScalarIn%scac + MQC_VectorScalarSum%column = VectorIn%column + Else + Call MQC_Error_A('ScalarIn type unspecified in MQC_VectorScalarSum', 6, & + 'ScalarIn%Data_type', ScalarIn%Data_type ) + EndIf + ElseIf(VectorIn%Data_type.eq.'Complex') then + If(ScalarIn%Data_type.eq.'Real') then + Call MQC_Allocate_Vector(VectorLen,MQC_VectorScalarSum,'Complex') + MQC_VectorScalarSum%vecc = VectorIn%vecc + ScalarIn%scar + MQC_VectorScalarSum%column = VectorIn%column + ElseIf(ScalarIn%Data_type.eq.'Integer') then + Call MQC_Allocate_Vector(VectorLen,MQC_VectorScalarSum,'Complex') + MQC_VectorScalarSum%vecc = VectorIn%vecc + ScalarIn%scai + MQC_VectorScalarSum%column = VectorIn%column + ElseIf(ScalarIn%Data_type.eq.'Complex') then + Call MQC_Allocate_Vector(VectorLen,MQC_VectorScalarSum,'Complex') + MQC_VectorScalarSum%vecc = VectorIn%vecc + ScalarIn%scac + MQC_VectorScalarSum%column = VectorIn%column + Else + Call MQC_Error_A('ScalarIn type unspecified in MQC_VectorScalarSum', 6, & + 'ScalarIn%Data_type', ScalarIn%Data_type ) + EndIf + Else + Call MQC_Error_A('VectorIn type unspecified in MQC_VectorScalarSum', 6, & + 'VectorIn%Data_type', VectorIn%Data_type ) + EndIf + + Return + End Function MQC_VectorScalarSum +! +! ! PROCEDURE MQC_ScalarVectorDifference ! !> \brief MQC_ScalarVectorDifference is a function that subtracts an MQC @@ -26272,7 +26377,8 @@ subroutine mqc_matrix_power(A,P) call mqc_error_I('P type not defined in MQC_Vector_Power',6) end select ! - call A%eigensys(eigenvals=sVals,reigenvecs=uMatrix) +! call A%eigensys(eigenvals=sVals,reigenvecs=uMatrix) + call A%diag(evals=sVals,evecs=uMatrix) call sVals%power(scalar) call sVals%diag(A) A = matmul(matmul(uMatrix,A),dagger(uMatrix)) From 4f7f846c656800e0a7522df87112893580c77494 Mon Sep 17 00:00:00 2001 From: Lee Thompson Date: Tue, 14 Jun 2022 18:50:54 -0400 Subject: [PATCH 10/10] Updated version number for release --- src/mqc_general.F03 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mqc_general.F03 b/src/mqc_general.F03 index 943aaea3..4a6296b7 100644 --- a/src/mqc_general.F03 +++ b/src/mqc_general.F03 @@ -155,9 +155,9 @@ subroutine mqc_version(major,minor,revision,versionString) character(len=*),OPTIONAL,intent(out)::versionString ! if(PRESENT(major)) major = 22 - if(PRESENT(minor)) minor = 4 + if(PRESENT(minor)) minor = 6 if(PRESENT(revision)) revision = 1 - if(PRESENT(versionString)) versionString = '22.04.1' + if(PRESENT(versionString)) versionString = '22.06.1' ! return end subroutine mqc_version