Added zmat construction via graph solver and other routines

src/prg_graphsolver_mod.F90

@@ -14,14 +14,15 @@ module prg_graphsolver_mod
   use prg_extras_mod
   use prg_densitymatrix_mod
   use prg_subgraphloop_mod
+  use prg_genz_mod
 
   implicit none
 
   private  !Everything is private by default
 
   integer, parameter :: dp = kind(1.0d0)
 
-  public :: prg_build_densityGP_T0
+  public :: prg_build_densityGP_T0, prg_build_zmatGP
 
 contains
 
@@ -31,13 +32,15 @@ module prg_graphsolver_mod
   !! \param rho_bml Density matrix.
   !! \param threshold Threshold for sparse matrix algebra.
   !! \param bndfil Filing factor.
+  !! \param Ef Fermi level or chemical potential.
+  !! \param nparts Number of parts to be used in graph partitioning.
+  !! \param verbose Verbosity level.
   !!
   subroutine prg_build_densityGP_T0(ham_bml, g_bml, rho_bml, threshold, bndfil, &
        & Ef, nparts, verbose)
 
     integer, parameter ::  dp = kind(1.0d0)
     type(bml_matrix_t), intent(in) :: ham_bml, g_bml
-    type(bml_matrix_t) :: aux_bml
     type(bml_matrix_t), intent(out) :: rho_bml
     real(dp) :: bndfil, threshold
     integer ::  myRank, norbs, nnodes, tnnz, i, vsize(2), inorbs
@@ -63,18 +66,16 @@ subroutine prg_build_densityGP_T0(ham_bml, g_bml, rho_bml, threshold, bndfil, &
     if(.not. bml_allocated(rho_bml))then
       call bml_zero_matrix(bml_type,bml_element_real,dp,norbs,norbs,rho_bml)
     endif
-    call bml_zero_matrix(bml_type,bml_element_real,dp,norbs,norbs,aux_bml)
 
     call bml_print_matrix("ham",ham_bml,0,10,0,10)
 
-    ! Doing graph partitioning using Metis
+    ! Using Metis for graph partitioning
     allocate(xadj(norbs+1)) ! Adjacency in csr format
 
     tnnz = 0
     do i=1,norbs
       tnnz = tnnz + bml_get_row_bandwidth(g_bml,i)
     enddo
-    write(6,*) 'tnnz=', tnnz
 
     allocate(adjncy(tnnz+1))
 
@@ -103,7 +104,6 @@ subroutine prg_build_densityGP_T0(ham_bml, g_bml, rho_bml, threshold, bndfil, &
     call prg_wait
 
     if(present(verbose) .and. (verbose > 1 .and. myRank == 1)) mlsi = mls()
-    mlsi = mls()
 
     ! Extract core halo indices from partitions
     do i=1,gpat%TotalParts
@@ -113,11 +113,11 @@ subroutine prg_build_densityGP_T0(ham_bml, g_bml, rho_bml, threshold, bndfil, &
            vsize,.true.)
       gpat%sgraph(i)%lsize = vsize(1)
       gpat%sgraph(i)%llsize = vsize(2)
-      write(6,*) 'nodeInpart size', size(gpat%sgraph(i)%nodeInPart)
-      write(6,*) 'nodeInpart(i)', gpat%nnodesInPart(i)
-
-      if(present(verbose) .and. (verbose > 1 .and. myRank == 1))&
-           write(*,*)"Core and CH sizes:",vsize(2),vsize(1)
+      if(present(verbose) .and. (verbose > 1 .and. myRank == 1))then
+        write(*,*) 'nodeInpart size', size(gpat%sgraph(i)%nodeInPart)
+        write(*,*) 'nodeInpart(i)', gpat%nnodesInPart(i)
+        write(*,*)"Core and CH sizes:",vsize(2),vsize(1)
+      endif
     enddo
 
     ! Construct DM for every part
@@ -129,7 +129,6 @@ subroutine prg_build_densityGP_T0(ham_bml, g_bml, rho_bml, threshold, bndfil, &
            vsize,.true.)
       gpat%sgraph(i)%lsize = vsize(1)
       gpat%sgraph(i)%llsize = vsize(2)
-      if(myRank == 1)write(*,*)"Core and CH sizes:",gpat%sgraph(i)%llsize,gpat%sgraph(i)%lsize
       inorbs = vsize(1)
       call bml_zero_matrix("dense",bml_element_real,dp,inorbs,inorbs,syprt(i)%estr%ham)
       if(allocated(vector))deallocate(vector)
@@ -138,8 +137,6 @@ subroutine prg_build_densityGP_T0(ham_bml, g_bml, rho_bml, threshold, bndfil, &
       call bml_matrix2submatrix(ham_bml, syprt(i)%estr%ham, vector, &
            & inorbs)
 
-      call bml_print_matrix("ham_part_bml",syprt(i)%estr%ham,0,10,0,10)
-
       !Computing the density matrix with diagonalization
       mlsii = mls()
       call bml_zero_matrix("dense",bml_element_real,dp,inorbs,inorbs,syprt(i)%estr%rho)
@@ -156,11 +153,9 @@ subroutine prg_build_densityGP_T0(ham_bml, g_bml, rho_bml, threshold, bndfil, &
     call prg_collectMatrixFromParts(gpat, rho_bml)
     call prg_wait
 
-    call bml_print_matrix("rho_bml",rho_bml,0,10,0,10)
     if(present(verbose) .and. (verbose > 1 .and. myRank == 1)) &
          write(*,*)"Total time for graph solver =",mls()-mlsi
 
-    CALL BML_DEALLOCATE(AUX_BML)
     do i = gpat%localPartMin(myRank), gpat%localPartMax(myRank)
       call bml_deallocate(syprt(i)%estr%ham)
       call bml_deallocate(syprt(i)%estr%rho)
@@ -171,4 +166,141 @@ subroutine prg_build_densityGP_T0(ham_bml, g_bml, rho_bml, threshold, bndfil, &
     if(getNRanks() == 0) call prg_progress_shutdown()
   end subroutine prg_build_densityGP_T0
 
+
+  !> Builds the inverse overlap factor matrix from \f$ S \f$ using a graph-based approach.
+  !! \param over_bml Input Overlap matrix.
+  !! \param g_bml Matrix to extract the graph from.
+  !! \param zmat_bml Overlap matrix.
+  !! \param threshold Threshold for sparse matrix algebra.
+  !! \param nparts Number of parts to be used in graph partitioning.
+  !! \param verbose Verbosity level.
+  !!
+  subroutine prg_build_zmatGP(over_bml, g_bml, zmat_bml, threshold, &
+       & nparts, verbose)
+
+    integer, parameter ::  dp = kind(1.0d0)
+    type(bml_matrix_t), intent(in) :: over_bml, g_bml
+    type(bml_matrix_t), intent(out) :: zmat_bml
+    real(dp) :: threshold
+    integer ::  myRank, norbs, nnodes, tnnz, i, vsize(2), inorbs
+    integer, intent(inout) :: nparts
+    character(20) :: bml_type
+    integer, allocatable :: xadj(:), adjncy(:), vector(:)
+    integer, allocatable  ::  part(:), core_count(:), Halo_count(:,:), CH_count(:)
+    type(graph_partitioning_t) ::  gpat
+    real(dp) :: maxCH, pnorm=6, smooth_maxCH, sumCubes, mlsi, mlsii
+    type(system_type), allocatable    ::  syprt(:)
+    integer, optional, intent(in) :: verbose
+
+    ! Initialize progress MPI and get ranks
+    if(getNRanks() == 0) call prg_progress_init()
+    !call prg_progress_init()
+    myRank = getMyRank() + 1
+
+    bml_type = bml_get_type(over_bml)
+
+    ! Allocate bml matrices
+    norbs = bml_get_N(over_bml)
+    if(.not. bml_allocated(zmat_bml))then
+      call bml_zero_matrix(bml_type,bml_element_real,dp,norbs,norbs,zmat_bml)
+    endif
+
+    ! Using Metis to do graph partitioning
+    allocate(xadj(norbs+1)) ! Adjacency in csr format
+
+    tnnz = 0
+    do i=1,norbs
+      tnnz = tnnz + bml_get_row_bandwidth(g_bml,i)
+    enddo
+    write(6,*) 'tnnz=', tnnz
+
+    allocate(adjncy(tnnz+1))
+
+    call bml_adjacency(g_bml, xadj, adjncy, 1)
+
+    nnodes = norbs
+
+    allocate(part(nnodes))
+    allocate(core_count(nparts))
+    allocate(CH_count(nparts))
+    allocate(Halo_count(nparts, nnodes))
+
+    call prg_metisPartition(gpat, nnodes, nnodes, xadj, adjncy, nparts, part, core_count,&
+         CH_count, Halo_count, sumCubes, maxCH, smooth_maxCH, pnorm)
+
+    if(present(verbose) .and. (verbose > 1 .and. myRank == 1)) then
+      write(*,*)"gpat = ",gpat%pname
+      write(*,*)"totalParts = ",gpat%totalParts
+      write(*,*)"totalNodes = ",gpat%totalNodes
+      write(*,*)"partMin = ",gpat%localPartMin(1)
+      write(*,*)"partMax = ",gpat%localPartMax(1)
+    endif
+
+    allocate(syprt(gpat%TotalParts))
+
+    call prg_wait
+
+    if(present(verbose) .and. (verbose > 1 .and. myRank == 1)) mlsi = mls()
+    mlsi = mls()
+
+    ! Extract core halo indices from partitions
+    do i=1,gpat%TotalParts
+      call bml_matrix2submatrix_index(g_bml,&
+           gpat%sgraph(i)%nodeInPart,gpat%nnodesInPart(i),&
+           gpat%sgraph(i)%core_halo_index, &
+           vsize,.true.)
+      gpat%sgraph(i)%lsize = vsize(1)
+      gpat%sgraph(i)%llsize = vsize(2)
+      if(present(verbose) .and. (verbose > 1 .and. myRank == 1))then
+        write(*,*) 'nodeInpart size', size(gpat%sgraph(i)%nodeInPart)
+        write(*,*) 'nodeInpart(i)', gpat%nnodesInPart(i)
+        write(*,*)"Core and CH sizes:",vsize(2),vsize(1)
+      endif
+    enddo
+
+    ! Construct Zmat for every part
+    do i = gpat%localPartMin(myRank), gpat%localPartMax(myRank)
+      !  do i=1,gpat%TotalParts ! If no MPI
+      call bml_matrix2submatrix_index(g_bml,&
+           gpat%sgraph(i)%nodeInPart,gpat%nnodesInPart(i),&
+           gpat%sgraph(i)%core_halo_index, &
+           vsize,.true.)
+      gpat%sgraph(i)%lsize = vsize(1)
+      gpat%sgraph(i)%llsize = vsize(2)
+      inorbs = vsize(1)
+      call bml_zero_matrix("dense",bml_element_real,dp,inorbs,inorbs,syprt(i)%estr%over)
+      if(allocated(vector))deallocate(vector)
+      allocate(vector(gpat%sgraph(i)%lsize))
+      vector(:) = gpat%sgraph(i)%core_halo_index(1:gpat%sgraph(i)%lsize)
+      call bml_matrix2submatrix(over_bml, syprt(i)%estr%over, vector, &
+           & inorbs)
+
+      !Computing the zmatrix with diagonalization
+      mlsii = mls()
+      call bml_zero_matrix("dense",bml_element_real,dp,inorbs,inorbs,syprt(i)%estr%zmat)
+      call prg_buildzdiag(syprt(i)%estr%over,syprt(i)%estr%zmat,threshold,inorbs,"dense")
+
+      if(present(verbose) .and. verbose > 1 .and. myRank == 1)&
+           write(*,*)"Time for prg_build_density_T0",mls()-mlsii
+      call bml_submatrix2matrix(syprt(i)%estr%zmat,zmat_bml,vector,gpat%sgraph(i)%lsize,gpat%sgraph(i)%llsize,threshold)
+    enddo
+
+    ! Collect the small DM matrices into its full system corresponding
+    call prg_partOrdering(gpat)
+    call prg_collectMatrixFromParts(gpat, zmat_bml)
+    call prg_wait
+
+    if(present(verbose) .and. (verbose > 1 .and. myRank == 1)) &
+         write(*,*)"Total time for graph solver =",mls()-mlsi
+
+    do i = gpat%localPartMin(myRank), gpat%localPartMax(myRank)
+      call bml_deallocate(syprt(i)%estr%over)
+      call bml_deallocate(syprt(i)%estr%zmat)
+    enddo
+    deallocate(syprt)
+
+    !call prg_progress_shutdown
+    if(getNRanks() == 0) call prg_progress_shutdown()
+  end subroutine prg_build_zmatGP
+
 end module prg_graphsolver_mod

src/prg_system_mod.F90

@@ -198,6 +198,7 @@ module prg_system_mod
   public :: prg_graph2vector, prg_vector2graph, prg_sortadj, prg_get_recip_vects, prg_translatetogeomcandfoldtobox
   public :: prg_write_trajectoryandproperty, prg_get_distancematrix, prg_parameters_to_vectors, prg_vectors_to_parameters
   public :: prg_get_dihedral, prg_wraparound, prg_centeratbox, prg_replicate, prg_cleanuprepeatedatoms
+  public :: prg_replicate_system, prg_destroy_system
 
 contains
 
@@ -635,6 +636,35 @@ subroutine prg_parse_system(system,filename,extin)
 
   end subroutine prg_parse_system
 
+
+  !> Deallocates all the arrays within a system
+  !! \param sy System type
+  !!
+  subroutine prg_destroy_system(sy)
+    implicit none
+    type(system_type), intent(inout)  ::  sy
+
+    if(allocated(sy%symbol)) deallocate(sy%symbol)
+    if(allocated(sy%atomic_number)) deallocate(sy%atomic_number)
+    if(allocated(sy%coordinate)) deallocate(sy%coordinate)
+    if(allocated(sy%velocity)) deallocate(sy%velocity)
+    if(allocated(sy%force)) deallocate(sy%force)
+    if(allocated(sy%net_charge)) deallocate(sy%net_charge)
+    if(allocated(sy%mass)) deallocate(sy%mass)
+    if(allocated(sy%lattice_vector)) deallocate(sy%lattice_vector)
+    if(allocated(sy%recip_vector)) deallocate(sy%recip_vector)
+    if(allocated(sy%spindex)) deallocate(sy%spindex)
+    if(allocated(sy%splist)) deallocate(sy%splist)
+    if(allocated(sy%spatnum)) deallocate(sy%spatnum)
+    if(allocated(sy%spmass)) deallocate(sy%spmass)
+    if(allocated(sy%userdef)) deallocate(sy%userdef)
+    if(allocated(sy%resindex)) deallocate(sy%resindex)
+    if(allocated(sy%resname)) deallocate(sy%resname)
+    if(allocated(sy%atomname)) deallocate(sy%atomname)
+
+  end subroutine prg_destroy_system
+
+
   !>  Write system in .xyz, .dat or pdb file.
   !! \param system System to be constructed.
   !! \param filename File name.
@@ -1495,6 +1525,105 @@ subroutine prg_replicate(coords,symbols,lattice_vectors,nx,ny,nz)
 
   end subroutine prg_replicate
 
+
+  !> Extend/replicate a system type along the lattice vectors.
+  !! \param sy System type.
+  !! \param syf System type output.
+  !! \param nx Number of lattice points in the v1 direction.
+  !! \param ny Number of lattice points in the v2 direction.
+  !! \param nz Number of lattice points in the v2 direction.
+  !!
+  subroutine prg_replicate_system(sy,syf,nx,ny,nz)
+    implicit none
+    integer                                  ::  i, j, k, l
+    integer                                  ::  ini, fin, cont, ii
+    integer, intent(in)                      ::  nx, ny, nz
+    type(system_type), intent(inout)         ::  sy
+    type(system_type)                        ::  syf
+
+    if(.not. allocated(sy%coordinate))then
+      write(*,*)"ERROR: System does not have coordinate"
+      stop
+    endif
+    if(.not. allocated(sy%symbol))then
+      write(*,*)"ERROR: System does not have atomic symbols"
+      stop
+    endif
+    syf%nats = sy%nats*nx*ny*nz
+    allocate(syf%coordinate(3,syf%nats))
+    allocate(syf%symbol(syf%nats))
+    allocate(syf%resindex(syf%nats))
+    allocate(syf%atomname(syf%nats))
+    allocate(syf%atomic_number(syf%nats))
+    allocate(syf%spindex(syf%nats))
+
+    cont = 0
+    do i = 1,nx
+      do j = 1,ny
+        do k = 1,nz
+          do ii = 1,sy%nats
+            cont = cont + 1
+            syf%coordinate(1,cont) = sy%coordinate(1,ii) + &
+                 & sy%lattice_vector(1,1)*i + sy%lattice_vector(2,1)*i + sy%lattice_vector(3,1)*i
+            syf%coordinate(2,cont) = sy%coordinate(2,ii) + &
+                 & sy%lattice_vector(1,2)*j + sy%lattice_vector(2,2)*j + sy%lattice_vector(3,2)*j
+            syf%coordinate(3,cont) = sy%coordinate(3,ii) + &
+                 & sy%lattice_vector(1,3)*k + sy%lattice_vector(2,3)*k + sy%lattice_vector(3,3)*k
+            syf%symbol(cont) = sy%symbol(ii)
+          enddo
+        enddo
+      enddo
+    enddo
+
+    if(allocated(sy%resindex))then
+      do i = 1,nx*ny*nz
+        ini = (i - 1)*sy%nats + 1
+        fin = i*sy%nats
+        syf%resindex(ini:fin) = sy%resindex(:)
+      enddo
+    else
+      syf%resindex = 1
+    endif
+
+    if(allocated(sy%atomname))then
+      do i = 1,nx*ny*nz
+        ini = (i - 1)*sy%nats + 1
+        fin = i*sy%nats
+        syf%atomname(ini:fin) = sy%atomname(:)
+      enddo
+    else
+      syf%atomname = "At"
+    endif
+
+    if(allocated(sy%atomic_number))then
+      do i = 1,nx*ny*nz
+        ini = (i - 1)*sy%nats + 1
+        fin = i*sy%nats
+        syf%atomic_number(ini:fin) = sy%atomic_number(:)
+      enddo
+    endif
+    if(allocated(sy%spindex))then
+      do i = 1,nx*ny*nz
+        ini = (i - 1)*sy%nats + 1
+        fin = i*sy%nats
+        syf%spindex(ini:fin) = sy%spindex(:)
+      enddo
+    endif
+
+    allocate(syf%lattice_vector(3,3))
+    syf%lattice_vector(1,:) = nx*sy%lattice_vector(1,:)
+    syf%lattice_vector(2,:) = ny*sy%lattice_vector(2,:)
+    syf%lattice_vector(3,:) = nz*sy%lattice_vector(3,:)
+
+    if(allocated(sy%splist))then
+      syf%nsp = sy%nsp
+      allocate(syf%splist(syf%nsp))
+      syf%splist = sy%splist
+    endif
+
+  end subroutine prg_replicate_system
+
+
   !> Cleanup repeated atoms we might have in the system.
   !! \param nats Number of atoms in the system.
   !! \param coords Coordinates of the system (see system_type).