Update density matrix build

Do not consider basis functions removed due to ill-conditioning.

CHANGELOG.md

 # ELSI changelog
 
+## Not released
+
+### ELSI interface
+* Computation of density matrix from eigenvectors was made more robust.
+
+### Known issues
+* The ELPA code can not be compiled with the NAG Fortran compiler, due to the
+  use of GNU extensions in ELPA.
+* Depending on the choice of k-points, the complex PEXSI solver may randomly
+  fail at the inertia counting stage.
+
 ## v2.6.3 (November 2020)
 
 ### ELSI interface
@@ -14,12 +25,6 @@
 ### EigenExa
 * Interface compatible with EigenExa 2.6.
 
-### Known issues
-* The ELPA code can not be compiled with the NAG Fortran compiler, due to the
-  use of GNU extensions in ELPA.
-* Depending on the choice of k-points, the complex PEXSI solver may randomly
-  fail at the inertia counting stage.
-
 ## v2.6.2 (July 2020)
 
 ### ELPA

CMakeLists.txt

@@ -7,7 +7,7 @@ SET(elsi_URL "http://elsi-interchange.org")
 SET(elsi_EMAIL "elsi-team@duke.edu")
 SET(elsi_LICENSE "BSD 3")
 SET(elsi_DESCRIPTION "Electronic Structure Infrastructure")
-SET(elsi_DATESTAMP "20201105")
+SET(elsi_DATESTAMP "20201117")
 
 ### CMake modules ###
 LIST(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)

src/elsi_solver.f90

@@ -2188,7 +2188,8 @@ subroutine elsi_compute_edm_real(eh,eval,occ,evec,edm)
 
    call elsi_allocate(eh%bh,factor,eh%ph%n_states,"factor",caller)
 
-   factor(:) = -occ*eval
+   factor(1:eh%ph%n_states_solve) = -occ(1:eh%ph%n_states_solve)&
+      *eval(1:eh%ph%n_states_solve)
 
    call elsi_build_dm_edm(eh%ph,eh%bh,factor,evec,edm,GET_EDM)
 
@@ -2217,7 +2218,8 @@ subroutine elsi_compute_edm_complex(eh,eval,occ,evec,edm)
 
    call elsi_allocate(eh%bh,factor,eh%ph%n_states,"factor",caller)
 
-   factor(:) = -occ*eval
+   factor(1:eh%ph%n_states_solve) = -occ(1:eh%ph%n_states_solve)&
+      *eval(1:eh%ph%n_states_solve)
 
    call elsi_build_dm_edm(eh%ph,eh%bh,factor,evec,edm,GET_EDM)
 

src/elsi_util.f90

@@ -818,7 +818,7 @@ subroutine elsi_build_dm_edm_real(ph,bh,factor,evec,dm,which)
    end if
 
    ! Methfessel-Paxton or energy density matrix may have negative factors
-   if(any(factor < 0.0_r8)) then
+   if(any(factor(1:ph%n_states_solve) < 0.0_r8)) then
       do i = 1,bh%n_lcol
          call elsi_get_gid(bh%my_pcol,bh%n_pcol,bh%blk,i,gid)
 
@@ -829,8 +829,8 @@ subroutine elsi_build_dm_edm_real(ph,bh,factor,evec,dm,which)
 
       call pdgemm("N","T",ph%n_basis,ph%n_basis,ph%n_states_solve,alpha,tmp,1,&
            1,bh%desc,evec,1,1,bh%desc,0.0_r8,dm,1,1,bh%desc)
-   else if(any(factor > 0.0_r8)) then
-      do i = 1,ph%n_states
+   else if(any(factor(1:ph%n_states_solve) > 0.0_r8)) then
+      do i = 1,ph%n_states_solve
          if(factor(i) > 0.0_r8) then
             max_state = i
          end if
@@ -939,8 +939,8 @@ subroutine elsi_build_dm_edm_cmplx(ph,bh,factor,evec,dm,which)
       alpha = (-1.0_r8,0.0_r8)
    end if
 
-   if(any(factor /= 0.0_r8)) then
-      do i = 1,ph%n_states
+   if(any(factor(1:ph%n_states_solve) /= 0.0_r8)) then
+      do i = 1,ph%n_states_solve
          if(factor(i) /= 0.0_r8) then
             max_state = i
          end if