From 46d9f3c8474b335954f71cdff06c6358c3d9876c Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 18 Mar 2021 00:55:37 +0100 Subject: [PATCH] Optimizations, maybe wrong --- src/csf/configuration_CI_sigma_helpers.irp.f | 385 ++++++++++--------- src/csf/sigma_vector.irp.f | 111 +++--- 2 files changed, 248 insertions(+), 248 deletions(-) diff --git a/src/csf/configuration_CI_sigma_helpers.irp.f b/src/csf/configuration_CI_sigma_helpers.irp.f index 5cf8cb64..951882b2 100644 --- a/src/csf/configuration_CI_sigma_helpers.irp.f +++ b/src/csf/configuration_CI_sigma_helpers.irp.f @@ -1,6 +1,6 @@ use bitmasks - BEGIN_PROVIDER [ integer(bit_kind), alphasIcfg_list , (N_int,2,N_configuration,mo_num*mo_num)] + BEGIN_PROVIDER [ integer(bit_kind), alphasIcfg_list , (N_int,2,N_configuration,elec_num)] &BEGIN_PROVIDER [ integer, NalphaIcfg_list, (N_configuration) ] implicit none !use bitmasks @@ -10,41 +10,43 @@ use bitmasks ! the input configuration Icfg. END_DOC - integer :: idxI ! The id of the Ith CFG - integer(bit_kind) :: Icfg(N_int,2) - integer :: NalphaIcfg + integer :: idxI ! The id of the Ith CFG + integer(bit_kind) :: Icfg(N_int,2) + integer :: NalphaIcfg logical,dimension(:,:),allocatable :: tableUniqueAlphas - integer :: listholes(mo_num) - integer :: holetype(mo_num) ! 1-> SOMO 2->DOMO - integer :: nholes - integer :: nvmos - integer :: listvmos(mo_num) - integer :: vmotype(mo_num) ! 1 -> VMO 2 -> SOMO - integer*8 :: Idomo - integer*8 :: Isomo - integer*8 :: Jdomo - integer*8 :: Jsomo - integer*8 :: diffSOMO - integer*8 :: diffDOMO - integer*8 :: xordiffSOMODOMO - integer :: ndiffSOMO - integer :: ndiffDOMO - integer :: nxordiffSOMODOMO - integer :: ndiffAll - integer :: i - integer :: j - integer :: k - integer :: kstart - integer :: kend - integer :: Nsomo_I - integer :: hole - integer :: p - integer :: q - integer :: countalphas - logical :: pqAlreadyGenQ - logical :: pqExistsQ - logical :: ppExistsQ + integer :: listholes(mo_num) + integer :: holetype(mo_num) ! 1-> SOMO 2->DOMO + integer :: nholes + integer :: nvmos + integer :: listvmos(mo_num) + integer :: vmotype(mo_num) ! 1 -> VMO 2 -> SOMO + integer*8 :: Idomo + integer*8 :: Isomo + integer*8 :: Jdomo + integer*8 :: Jsomo + integer*8 :: diffSOMO + integer*8 :: diffDOMO + integer*8 :: xordiffSOMODOMO + integer :: ndiffSOMO + integer :: ndiffDOMO + integer :: nxordiffSOMODOMO + integer :: ndiffAll + integer :: i + integer :: j + integer :: k + integer :: kstart + integer :: kend + integer :: Nsomo_I + integer :: hole + integer :: p + integer :: q + integer :: countalphas + logical :: pqAlreadyGenQ + logical :: pqExistsQ + logical :: ppExistsQ + double precision :: t0, t1 + call wall_time(t0) allocate(tableUniqueAlphas(mo_num,mo_num)) NalphaIcfg_list = 0 @@ -60,19 +62,19 @@ use bitmasks nholes = 0 ! holes in SOMO do i = 1,mo_num - if(POPCNT(IAND(Isomo,IBSET(0_8,i-1))) .EQ. 1) then - nholes += 1 - listholes(nholes) = i - holetype(nholes) = 1 - endif + if(POPCNT(IAND(Isomo,IBSET(0_8,i-1))) .EQ. 1) then + nholes += 1 + listholes(nholes) = i + holetype(nholes) = 1 + endif end do ! holes in DOMO do i = 1,mo_num - if(POPCNT(IAND(Idomo,IBSET(0_8,i-1))) .EQ. 1) then - nholes += 1 - listholes(nholes) = i - holetype(nholes) = 2 - endif + if(POPCNT(IAND(Idomo,IBSET(0_8,i-1))) .EQ. 1) then + nholes += 1 + listholes(nholes) = i + holetype(nholes) = 2 + endif end do ! find vmos @@ -80,16 +82,17 @@ use bitmasks vmotype = -1 nvmos = 0 do i = 1,mo_num - !print *,i,IBSET(0,i-1),POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))), POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) - if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 0 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) .EQ. 0) then + if(IAND(Idomo,(IBSET(0_8,i-1))) .EQ. 0) then + if(IAND(Isomo,(IBSET(0_8,i-1))) .EQ. 0) then nvmos += 1 listvmos(nvmos) = i vmotype(nvmos) = 1 - else if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 1 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) .EQ. 0 ) then + else if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 1) then nvmos += 1 listvmos(nvmos) = i vmotype(nvmos) = 2 - end if + end if + end if end do tableUniqueAlphas = .FALSE. @@ -106,157 +109,161 @@ use bitmasks kend = idxI-1 do i = 1,nholes - p = listholes(i) - do j = 1,nvmos - q = listvmos(j) - if(p .EQ. q) cycle - if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then - ! SOMO -> VMO - Jsomo = IBCLR(Isomo,p-1) - Jsomo = IBSET(Jsomo,q-1) - Jdomo = Idomo - kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-2))) - kend = idxI-1 - else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then - ! SOMO -> SOMO - Jsomo = IBCLR(Isomo,p-1) - Jsomo = IBCLR(Jsomo,q-1) - Jdomo = IBSET(Idomo,q-1) - kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-4))) - kend = idxI-1 - else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then - ! DOMO -> VMO - Jsomo = IBSET(Isomo,p-1) - Jsomo = IBSET(Jsomo,q-1) - Jdomo = IBCLR(Idomo,p-1) - kstart = cfg_seniority_index(Nsomo_I) - kend = idxI-1 - else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then - ! DOMO -> SOMO - Jsomo = IBSET(Isomo,p-1) - Jsomo = IBCLR(Jsomo,q-1) - Jdomo = IBCLR(Idomo,p-1) - Jdomo = IBSET(Jdomo,q-1) - kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-2))) - kend = idxI-1 - else - print*,"Something went wrong in obtain_associated_alphaI" + p = listholes(i) + do j = 1,nvmos + q = listvmos(j) + if(p .EQ. q) cycle + if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then + ! SOMO -> VMO + Jsomo = IBCLR(Isomo,p-1) + Jsomo = IBSET(Jsomo,q-1) + Jdomo = Idomo + kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-2))) + kend = idxI-1 + else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then + ! SOMO -> SOMO + Jsomo = IBCLR(Isomo,p-1) + Jsomo = IBCLR(Jsomo,q-1) + Jdomo = IBSET(Idomo,q-1) + kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-4))) + kend = idxI-1 + else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then + ! DOMO -> VMO + Jsomo = IBSET(Isomo,p-1) + Jsomo = IBSET(Jsomo,q-1) + Jdomo = IBCLR(Idomo,p-1) + kstart = cfg_seniority_index(Nsomo_I) + kend = idxI-1 + else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then + ! DOMO -> SOMO + Jsomo = IBSET(Isomo,p-1) + Jsomo = IBCLR(Jsomo,q-1) + Jdomo = IBCLR(Idomo,p-1) + Jdomo = IBSET(Jdomo,q-1) + kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-2))) + kend = idxI-1 + else + print*,"Something went wrong in obtain_associated_alphaI" + endif + + ! Again, we don't have to search from 1 + ! we just use seniority to find the + ! first index with NSOMO - 2 to NSOMO + 2 + ! this is what is done in kstart, kend + + pqAlreadyGenQ = .FALSE. + ! First check if it can be generated before + do k = kstart, kend + diffSOMO = IEOR(Jsomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k))) + ndiffSOMO = POPCNT(diffSOMO) + if((ndiffSOMO .NE. 0) .AND. (ndiffSOMO .NE. 2)) cycle + diffDOMO = IEOR(Jdomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k))) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffDOMO = POPCNT(diffDOMO) + nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) + nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO + !if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then + if((ndiffSOMO+ndiffDOMO) .EQ. 0) then + pqAlreadyGenQ = .TRUE. + ppExistsQ = .TRUE. + EXIT endif - - ! Again, we don't have to search from 1 - ! we just use seniorty to find the - ! first index with NSOMO - 2 to NSOMO + 2 - ! this is what is done in kstart, kend - - pqAlreadyGenQ = .FALSE. - ! First check if it can be generated before - do k = kstart, kend - diffSOMO = IEOR(Jsomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k))) - ndiffSOMO = POPCNT(diffSOMO) - if((ndiffSOMO .NE. 0) .AND. (ndiffSOMO .NE. 2)) cycle - diffDOMO = IEOR(Jdomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k))) - xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) - ndiffDOMO = POPCNT(diffDOMO) - nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) - nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO - !if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then - if((ndiffSOMO+ndiffDOMO) .EQ. 0) then - pqAlreadyGenQ = .TRUE. - ppExistsQ = .TRUE. - EXIT - endif - if((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then - pqAlreadyGenQ = .TRUE. - EXIT - endif - end do - - if(pqAlreadyGenQ) cycle - - pqExistsQ = .FALSE. - - if(.NOT. pqExistsQ) then - tableUniqueAlphas(p,q) = .TRUE. + if((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then + pqAlreadyGenQ = .TRUE. + EXIT endif - end do + end do + + if(pqAlreadyGenQ) cycle + + pqExistsQ = .FALSE. + + if(.NOT. pqExistsQ) then + tableUniqueAlphas(p,q) = .TRUE. + endif + end do end do - !print *,tableUniqueAlphas(:,:) + !print *,tableUniqueAlphas(:,:) - ! prune list of alphas - Isomo = Icfg(1,1) - Idomo = Icfg(1,2) - Jsomo = Icfg(1,1) - Jdomo = Icfg(1,2) - NalphaIcfg = 0 - do i = 1, nholes - p = listholes(i) - do j = 1, nvmos + ! prune list of alphas + Isomo = Icfg(1,1) + Idomo = Icfg(1,2) + Jsomo = Icfg(1,1) + Jdomo = Icfg(1,2) + NalphaIcfg = 0 + do i = 1, nholes + p = listholes(i) + do j = 1, nvmos q = listvmos(j) if(p .EQ. q) cycle if(tableUniqueAlphas(p,q)) then - if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then - ! SOMO -> VMO - Jsomo = IBCLR(Isomo,p-1) - Jsomo = IBSET(Jsomo,q-1) - Jdomo = Idomo - else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then - ! SOMO -> SOMO - Jsomo = IBCLR(Isomo,p-1) - Jsomo = IBCLR(Jsomo,q-1) - Jdomo = IBSET(Idomo,q-1) - else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then - ! DOMO -> VMO - Jsomo = IBSET(Isomo,p-1) - Jsomo = IBSET(Jsomo,q-1) - Jdomo = IBCLR(Idomo,p-1) - else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then - ! DOMO -> SOMO - Jsomo = IBSET(Isomo,p-1) - Jsomo = IBCLR(Jsomo,q-1) - Jdomo = IBCLR(Idomo,p-1) - Jdomo = IBSET(Jdomo,q-1) - else - print*,"Something went wrong in obtain_associated_alphaI" - endif + if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then + ! SOMO -> VMO + Jsomo = IBCLR(Isomo,p-1) + Jsomo = IBSET(Jsomo,q-1) + Jdomo = Idomo + else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then + ! SOMO -> SOMO + Jsomo = IBCLR(Isomo,p-1) + Jsomo = IBCLR(Jsomo,q-1) + Jdomo = IBSET(Idomo,q-1) + else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then + ! DOMO -> VMO + Jsomo = IBSET(Isomo,p-1) + Jsomo = IBSET(Jsomo,q-1) + Jdomo = IBCLR(Idomo,p-1) + else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then + ! DOMO -> SOMO + Jsomo = IBSET(Isomo,p-1) + Jsomo = IBCLR(Jsomo,q-1) + Jdomo = IBCLR(Idomo,p-1) + Jdomo = IBSET(Jdomo,q-1) + else + print*,"Something went wrong in obtain_associated_alphaI" + endif - ! SOMO - NalphaIcfg += 1 - !print *,i,j,"|",NalphaIcfg - alphasIcfg_list(1,1,idxI,NalphaIcfg) = Jsomo - alphasIcfg_list(1,2,idxI,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1) - NalphaIcfg_list(idxI) = NalphaIcfg + ! SOMO + NalphaIcfg += 1 + !print *,i,j,"|",NalphaIcfg + alphasIcfg_list(1,1,idxI,NalphaIcfg) = Jsomo + alphasIcfg_list(1,2,idxI,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1) + NalphaIcfg_list(idxI) = NalphaIcfg endif - end do - end do + end do + end do - ! Check if this Icfg has been previously generated as a mono - ppExistsQ = .False. - Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1)) - Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2)) - do k = 1, idxI-1 - diffSOMO = IEOR(Isomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k))) - diffDOMO = IEOR(Idomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k))) - xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) - ndiffSOMO = POPCNT(diffSOMO) - ndiffDOMO = POPCNT(diffDOMO) - nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) - if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then + ! Check if this Icfg has been previously generated as a mono + ppExistsQ = .False. + Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1)) + Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2)) + kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-2))) + do k = kstart, idxI-1 + diffSOMO = IEOR(Isomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k))) + ndiffSOMO = POPCNT(diffSOMO) + if (ndiffSOMO /= 2) cycle + diffDOMO = IEOR(Idomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k))) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffDOMO = POPCNT(diffDOMO) + nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) + if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4)) then ppExistsQ = .TRUE. EXIT - endif - end do - ! Diagonal part (pp,qq) - if(nholes > 0 .AND. (.NOT. ppExistsQ))then - ! SOMO - NalphaIcfg += 1 - alphasIcfg_list(1,1,idxI,NalphaIcfg) = Icfg(1,1) - alphasIcfg_list(1,2,idxI,NalphaIcfg) = Icfg(1,2) - NalphaIcfg_list(idxI) = NalphaIcfg - endif + endif + end do + ! Diagonal part (pp,qq) + if(nholes > 0 .AND. (.NOT. ppExistsQ))then + ! SOMO + NalphaIcfg += 1 + alphasIcfg_list(1,1,idxI,NalphaIcfg) = Icfg(1,1) + alphasIcfg_list(1,2,idxI,NalphaIcfg) = Icfg(1,2) + NalphaIcfg_list(idxI) = NalphaIcfg + endif - NalphaIcfg = 0 + NalphaIcfg = 0 enddo ! end loop idxI + call wall_time(t1) + print *, 'Preparation : ', t1 - t0 END_PROVIDER @@ -416,7 +423,7 @@ END_PROVIDER endif ! Again, we don't have to search from 1 - ! we just use seniorty to find the + ! we just use seniortiy to find the ! first index with NSOMO - 2 to NSOMO + 2 ! this is what is done in kstart, kend @@ -588,10 +595,10 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod !integer(bit_kind) :: Idomo(N_int) !integer(bit_kind) :: Jsomo(N_int) !integer(bit_kind) :: Jdomo(N_int) - integer*8 :: Isomo - integer*8 :: Idomo - integer*8 :: Jsomo - integer*8 :: Jdomo + integer*8 :: Isomo + integer*8 :: Idomo + integer*8 :: Jsomo + integer*8 :: Jdomo integer*8 :: mask integer :: iint, ipos !integer(bit_kind) :: Isomotmp(N_int) diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f index 01a2544a..a2ec9cbf 100644 --- a/src/csf/sigma_vector.irp.f +++ b/src/csf/sigma_vector.irp.f @@ -144,6 +144,7 @@ end subroutine get_phase_qp_to_cfg integer :: nt + norm_det1 = 0.d0 MS = elec_alpha_num - elec_beta_num print *,"Maxbfdim=",NBFMax @@ -444,7 +445,6 @@ end subroutine get_phase_qp_to_cfg print *,"Rowsmax=",rowsmax," Colsmax=",colsmax END_PROVIDER - !BEGIN_PROVIDER [ real*8, AIJpqContainer, (NSOMOMin:NSOMOMax,4,NSOMOMax+1,NSOMOMax+1,NBFMax,NBFMax)] BEGIN_PROVIDER [ real*8, AIJpqContainer, (NBFMax,NBFmax,NSOMOMax+1,NSOMOMax+1,4,NSOMOMin:NSOMOMax)] use cfunctions implicit none @@ -476,79 +476,72 @@ end subroutine get_phase_qp_to_cfg cols = -1 integer*8 MS MS = 0 - touch AIJpqMatrixDimsList real*8,dimension(:,:),allocatable :: meMatrix integer maxdim - !maxdim = max(rowsmax,colsmax) - ! allocate matrix - !print *,"rowsmax =",rowsmax," colsmax=",colsmax - !print *,"NSOMOMax = ",NSOMOMax + ! Type ! 1. SOMO -> SOMO - !print *,"Doing SOMO -> SOMO" - !AIJpqContainer(NSOMOMin,1,1,1,1,1) = 1.0d0 + AIJpqContainer = 0.d0 AIJpqContainer(1,1,1,1,1,NSOMOMin) = 1.0d0 + integer :: rows_old, cols_old + rows_old = -1 + cols_old = -1 + allocate(meMatrix(1,1)) do i = NSOMOMin+2, NSOMOMax, 2 Isomo = ISHFT(1_8,i)-1 - do j = i-2,i-2, 2 - if(j .GT. NSOMOMax .OR. j .LT. 0) cycle - !print *,"i,j=",i,j - do k = 1,i - do l = 1,i + j=i-2 + if(j .GT. NSOMOMax .OR. j .LT. 0) cycle + nsomoi = i + do k = 1,i + orbp = k + do l = 1,i - ! Define Jsomo - if(k .NE. l) then - Jsomo = IBCLR(Isomo, k-1) - Jsomo = IBCLR(Jsomo, l-1) - nsomoi = i - nsomoj = j - else - Isomo = ISHFT(1_8,i)-1 - Jsomo = ISHFT(1_8,i)-1 - nsomoi = i - nsomoj = i - endif + ! Define Jsomo + if(k .NE. l) then + Jsomo = IBCLR(Isomo, k-1) + Jsomo = IBCLR(Jsomo, l-1) + nsomoj = j + else + Isomo = ISHFT(1_8,i)-1 + Jsomo = ISHFT(1_8,i)-1 + nsomoj = i + endif - !print *,"k,l=",k,l - !call debug_spindet(Jsomo,1) - !call debug_spindet(Isomo,1) + call getApqIJMatrixDims(Isomo, & + Jsomo, & + MS, & + rows, & + cols) - !AIJpqContainer(nsomoi,1,k,l,:,:) = 0.0d0 - AIJpqContainer(:,:,k,l,1,nsomoi) = 0.0d0 - call getApqIJMatrixDims(Isomo, & - Jsomo, & - MS, & - rows, & - cols) - - orbp = k - orbq = l - allocate(meMatrix(rows,cols)) - meMatrix = 0.0d0 - ! fill matrix - call getApqIJMatrixDriver(Isomo, & - Jsomo, & - orbp, & - orbq, & - MS, & - NMO, & - meMatrix, & - rows, & - cols) - !print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols,">",rowsmax,colsmax - !call printMatrix(meMatrix,rows,cols) - ! i -> j - do ri = 1,rows - do ci = 1,cols - !AIJpqContainer(nsomoi,1,k,l,ri,ci) = meMatrix(ri, ci) - AIJpqContainer(ri,ci,k,l,1,nsomoi) = meMatrix(ri, ci) - end do + orbq = l + if ((rows /= rows_old).or.(cols /= cols_old)) then + deallocate(meMatrix) + allocate(meMatrix(rows,cols)) + rows_old = rows + cols_old = cols + endif + meMatrix = 0.0d0 + ! fill matrix + call getApqIJMatrixDriver(Isomo, & + Jsomo, & + orbp, & + orbq, & + MS, & + NMO, & + meMatrix, & + rows, & + cols) + ! i -> j + do ri = 1,rows + do ci = 1,cols + AIJpqContainer(ri,ci,k,l,1,nsomoi) = meMatrix(ri, ci) end do - deallocate(meMatrix) end do end do end do end do + deallocate(meMatrix) + ! Type ! 2. DOMO -> VMO !print *,"Doing DOMO -> VMO"