Add dmconstruction3d benchmark

benchmarks/CMakeLists.txt

@@ -27,11 +27,16 @@ endfunction(progress_benchmark)
 
 progress_benchmark(dmconstruction dmconstruction/dmconstruction.F90)
 
+progress_benchmark(dmconstruction3d dmconstruction3d/dmconstruction3d.F90)
+
+progress_benchmark(dmconstruction_bio dmconstruction3d/dmconstruction_bio.F90
+                                      dmconstruction_graphBased/aux_mod.F90)
+
 progress_benchmark(dmconstruction_gp dmconstruction_graphBased/dmconstruction_graphBased.F90)
 
 progress_benchmark(dmconstruction_gpbio dmconstruction_graphBased/dmconstruction_graphBased_bio.F90 dmconstruction_graphBased/aux_mod.F90)
 
 SET(EXECUTABLE_OUTPUT_PATH ${PROJECT_BINARY_DIR})
 
-install(TARGETS dmconstruction
+install(TARGETS dmconstruction dmconstruction_bio
         DESTINATION ${CMAKE_INSTALL_PREFIX}/bin)

benchmarks/dmconstruction3d/dmconstruction3d.F90

@@ -0,0 +1,146 @@
+!> High-level program to construct a model Hamiltonian
+!!
+program hmodel
+
+  !BML lib.
+  use bml
+
+  !PROGRESS lib modes
+  use prg_modelham_mod
+  use prg_system_mod
+  use prg_densitymatrix_mod
+  use prg_dos_mod
+  use prg_sp2_mod
+  use prg_timer_mod
+  use prg_extras_mod
+  use prg_parallel_mod
+
+  implicit none
+  integer, parameter ::  dp = kind(1.0d0)
+  integer ::  norbs,seed,i,nreps
+  integer ::  verbose
+  character(20) :: filename,arg
+  type(bml_matrix_t) ::  ham_bml,rho_bml,rhos_bml,evects_bml,aux_bml
+  type(mham_type) ::  mham
+  type(system_type) ::  sys
+  real(dp) ::  threshold, bndfil, tol, n1d
+  real(dp), allocatable :: eigenvalues(:)
+  real(dp) :: ef,sparsity,dec,mlsi,mlsf,bnorm
+  character(20) :: bml_dmode
+
+  call prg_initParallel()
+
+  if (getNRanks().gt.1)then
+    if (printRank() .eq. 1)print*,'BML_DMODE_DISTRIBUTED'
+    bml_dmode = BML_DMODE_DISTRIBUTED
+  else
+    print*,'BML_DMODE_SEQUENTIAL'
+    bml_dmode = BML_DMODE_SEQUENTIAL
+  endif
+
+  !Parsing input file.
+  call getarg(1,filename)
+  call prg_parse_mham(mham,trim(adjustl(filename))) !Reads the input for modelham
+  nreps=1
+  if (iargc().gt.1)then
+    call getarg(2,arg)
+    read(arg,*)nreps
+  endif
+  !Number of orbitals/matrix size
+  norbs=mham%norbs
+
+  if (printRank() .eq. 1)print*,'norbs = ',norbs
+  n1d=norbs/2.
+  n1d=n1d**(1./3.)
+  if (printRank() .eq. 1)print*,'n1d = ',n1d
+  norbs=n1d**3
+  norbs=norbs*2
+  if (printRank() .eq. 1)print*,'Uses ',norbs,' orbitals'
+
+  !Allocating bml matrices
+  call bml_zero_matrix(mham%bml_type,bml_element_real,dp,norbs,norbs,ham_bml, &
+       &               bml_dmode)
+  call bml_zero_matrix(mham%bml_type,bml_element_real,dp,norbs,norbs,rho_bml, &
+       &               bml_dmode)
+  call bml_zero_matrix(mham%bml_type,bml_element_real,dp,norbs,norbs, &
+       &               evects_bml, bml_dmode)
+  call bml_zero_matrix(mham%bml_type,bml_element_real,dp,norbs,norbs,rhos_bml, &
+       &               bml_dmode)
+  call bml_zero_matrix(mham%bml_type,bml_element_real,dp,norbs,norbs,aux_bml, &
+       &               bml_dmode)
+
+  seed = 1000 !Seed to reproduce the Hamiltonian build
+  verbose = 1 !Verbosity level
+  threshold = 1.0d-5 !Threshold value for the matrices through the whole code
+  bndfil = 0.5d0 !Fraction of orbitals that will be filled
+
+  allocate(eigenvalues(norbs))
+
+  !Constructing the Hamiltonian
+  call prg_twolevel_model3d(mham%ea, mham%eb, mham%dab, mham%daiaj, mham%dbibj, &
+       &mham%dec, mham%rcoeff, mham%reshuffle, mham%seed, ham_bml, verbose)
+  call bml_threshold(ham_bml, threshold)
+  call bml_print_matrix("ham_bml",ham_bml,0,10,0,10)
+
+  sparsity = bml_get_sparsity(ham_bml, 1.0D-5)
+  if (printRank() .eq. 1)write(*,*)"Sparsity Ham=",sparsity
+
+  !Computing the density matrix with diagonalization
+  if (printRank() .eq. 1)print*,'prg_build_density_T0'
+  !run loop twice to have an accurate timing in 2nd call
+  do i =1, nreps
+    mlsi = mls()
+    call prg_build_density_T0(ham_bml, rho_bml, threshold, bndfil, eigenvalues)
+    mlsf = mls()
+    if (printRank() .eq. 1)write(*,*)"Time_for_prg_build_density_T0",mlsf-mlsi
+  enddo
+  print*,'Eigenvalues:'
+  do i = norbs/2-2, norbs/2+3
+    write(*,*)eigenvalues(i)
+  enddo
+
+  sparsity = bml_get_sparsity(rho_bml, 1.0D-5)
+  if (printRank() .eq. 1)write(*,*)"Sparsity Rho=",sparsity
+
+  !Getting the fermi level
+  ef = (eigenvalues(int(norbs/2)+1) + eigenvalues(int(norbs/2)))/2
+  eigenvalues = eigenvalues - ef
+
+  !Writting the total DOS
+  call prg_write_tdos(eigenvalues, 0.05d0, 10000, -20.0d0, 20.0d0, "tdos.dat")
+
+  tol = 2.0D-5*norbs*bndfil
+
+  !run loop twice to have an accurate timing in 2nd call
+  do i =1, nreps
+    !Solving for Rho using SP2
+    mlsi = mls()
+    call prg_sp2_alg1(ham_bml,rhos_bml,threshold,bndfil,15,100 &
+         ,"Rel",tol,20)
+    mlsf = mls()
+    if (printRank() .eq. 1)write(*,*)"Time_for_prg_sp2_alg1",mlsf-mlsi
+  enddo
+  call bml_print_matrix("rho_bml",rho_bml,0,10,0,10)
+  call bml_print_matrix("rhos_bml",rhos_bml,0,10,0,10)
+
+  call bml_copy(rhos_bml,aux_bml)
+  call bml_add(aux_bml,rho_bml,1.0d0,-1.0d0,threshold)
+  bnorm=bml_fnorm(aux_bml)
+  if (printRank() .eq. 1)write(*,*)"||DM_sp2-DM_diag||_F",bnorm
+
+  call bml_multiply(rhos_bml, rhos_bml, aux_bml, 0.5_dp, 0.0_dp, threshold)
+  call bml_print_matrix("rhos_bml^2",aux_bml,0,10,0,10)
+  call bml_add(aux_bml,rhos_bml,1.0d0,-1.0d0,threshold)
+  bnorm=bml_fnorm(aux_bml)
+  if (printRank() .eq. 1)write(*,*)"||DM_sp2-DM_sp2^2||_F",bnorm
+
+  call bml_multiply(ham_bml,rhos_bml,aux_bml,1.0_dp,0.0_dp,threshold)
+  call bml_multiply(rhos_bml,ham_bml,aux_bml,1.0_dp,-1.0_dp,threshold)
+  bnorm=bml_fnorm(aux_bml)
+  if (printRank() .eq. 1)write(*,*)"||DM_sp2*H-H*DM_sp2||_F",bnorm
+
+  !call bml_write_matrix(ham_bml, "hamiltonian.mtx")
+
+  call prg_shutdownParallel()
+
+end program hmodel

benchmarks/dmconstruction3d/dmconstruction_bio.F90

@@ -0,0 +1,151 @@
+!> Graph-based aproach driver.
+!!
+program biosolve
+
+  use bml
+  use aux_mod
+  use prg_sp2_mod
+
+  !PROGRESS lib modes
+  use prg_modelham_mod
+  use prg_system_mod
+  use prg_timer_mod
+  use prg_extras_mod
+  use prg_parallel_mod
+  use prg_progress_mod
+  use prg_densitymatrix_mod
+  use ham_latte_mod
+  use tbparams_latte_mod
+  use ppot_latte_mod
+  use prg_ptable_mod
+  use prg_genz_mod
+  use prg_nonortho_mod
+
+  integer, parameter ::  dp = kind(1.0d0)
+  character(2), allocatable         ::  TypeA(:,:), TypeB(:,:)
+  character(3), allocatable         ::  intKind(:)
+  integer                           ::  myRank, nel, norbs, nnz
+  integer, allocatable              ::  hindex(:,:)
+  real(dp)                          ::  bndfil, ef, mlssp2, mlsi
+  real(dp)                          ::  mlsdiag, sparsity, threshold
+  real(dp), allocatable             ::  onsitesH(:,:), onsitesS(:,:), origin(:), trace(:)
+  type(bioham_type)                 ::  bioham
+  type(bml_matrix_t)                ::  aux_bml, ham_bml, oham_bml
+  type(bml_matrix_t)                ::  over_bml, rho_bml
+  type(intpairs_type), allocatable  ::  intPairsH(:,:), intPairsS(:,:)
+  type(ppot_type), allocatable      ::  ppot(:,:)
+  type(system_type)                 ::  sy, syf
+  type(tbparams_type)               ::  tb
+  real(dp) ::  tol
+  real(dp), allocatable :: eigenvalues(:)
+
+  call prg_progress_init()
+  myRank = getMyRank() + 1
+
+  ! Parsing input file
+  call prg_parse_bioham(bioham,"input.in") !Reads the input for modelham
+
+  ! Reading the system
+  call prg_parse_system(sy,"prot.pdb")
+
+  call prg_replicate_system(sy,syf,bioham%replicatex,bioham%replicatey,bioham%replicatez)
+  call prg_destroy_system(sy)
+
+  ! Center sytem inside the box and fold it by the lattice_vectors.
+  allocate(origin(3)); origin = 0.0_dp
+  call prg_translateandfoldtobox(syf%coordinate,syf%lattice_vector,origin)
+
+  call prg_write_system(syf,"protf.pdb")
+
+  ! Constructing the Hamiltonian
+  call load_latteTBparams(tb,syf%splist,bioham%parampath)
+
+  ! Get the mapping of the Hamiltonian index with the atom index
+  allocate(hindex(2,syf%nats))
+
+
+  ! Bond integrals parameters for LATTE Hamiltonian.
+  call load_bintTBparamsH(syf%splist,tb%onsite_energ,&
+       typeA,typeB,intKind,onsitesH,onsitesS,intPairsH,intPairsS,bioham%parampath)
+
+  call write_bintTBparamsH(typeA,typeB,&
+       intKind,intPairsH,intPairsS,adjustl(trim(bioham%jobname))//"_mybondints.nonorth")
+
+  ! Load Pair potentials for LATTE TB.
+  call load_PairPotTBparams(bioham%parampath,syf%splist,ppot)
+
+  call get_hindex(syf%spindex,tb%norbi,hindex,norbs)
+
+  call bml_zero_matrix(bioham%bml_type,bml_element_real,dp,norbs,norbs,ham_bml)
+  call bml_zero_matrix(bioham%bml_type,bml_element_real,dp,norbs,norbs,over_bml)
+
+  call get_hsmat(ham_bml,over_bml,syf%coordinate,&
+       syf%lattice_vector,syf%spindex,&
+       tb%norbi,hindex,onsitesH,onsitesS,intPairsH,intPairsS,bioham%threshold)
+
+  if(bioham%mdim == 0) bioham%mdim = norbs
+  ! Get occupation based on last shell population.
+  nel = sum(element_numel(syf%atomic_number(:)),&
+       & size(syf%atomic_number,dim=1))
+  bndfil = nel/(2.0_dp*real(norbs,dp))
+
+  call bml_threshold(ham_bml,bioham%threshold)
+
+  if(myRank == 1)call bml_print_matrix("ham_bml",ham_bml,0,10,0,10)
+  sparsity = bml_get_sparsity(ham_bml,bioham%threshold)
+  if(myRank == 1)write(*,*)"Sparsity Ham=",sparsity
+
+  call bml_zero_matrix(bioham%bml_type,bml_element_real,dp,norbs,norbs,aux_bml)
+
+  call prg_buildzdiag(over_bml,aux_bml,bioham%threshold,bioham%mdim,bioham%bml_type)
+  if(myRank == 1)call bml_print_matrix("zmat",aux_bml,0,10,0,10)
+  call bml_deallocate(over_bml)
+
+  call bml_zero_matrix(bioham%bml_type,bml_element_real,dp,norbs,norbs,oham_bml)
+  call prg_orthogonalize(ham_bml,aux_bml,oham_bml,&
+       bioham%threshold,bioham%bml_type,bioham%verbose)
+  call bml_deallocate(ham_bml)
+
+  call bml_zero_matrix(bioham%bml_type,bml_element_real,dp,norbs,norbs,rho_bml)
+  ! Call SP2
+  mlsi = mls()
+  tol = 2.0D-5*norbs*bndfil
+  call prg_sp2_alg1(oham_bml, rho_bml, 1.0D-5, bndfil, 15,100, "Rel", tol, 20)
+  mlssp2 = mls()-mlsi
+
+  call bml_deallocate(aux_bml)
+  call bml_zero_matrix(bioham%bml_type,bml_element_real,dp,norbs,norbs,aux_bml)
+  allocate(eigenvalues(norbs))
+
+  ! Construct the density matrix from diagonalization of full matrix to compare with
+  mlsi = mls()
+  call prg_build_density_T0(oham_bml,aux_bml,bioham%threshold, bndfil, eigenvalues)
+  mlsdiag = mls()-mlsi
+
+  ! count number of non-zeros in DM
+  nnz = 0
+  do i=1,norbs
+    nnz = nnz + bml_get_row_bandwidth(rho_bml,i)
+  enddo
+
+  if(myRank == 1)call bml_print_matrix("rhoSP2",rho_bml,0,10,0,10)
+  if(myRank == 1)call bml_print_matrix("rhoDIAG",aux_bml,0,10,0,10)
+  call bml_add(aux_bml,rho_bml,1.0d0,-1.0d0,threshold)
+
+  if(myRank == 1)then
+    write(*,*)"Nnz in DM                = ",nnz
+    write(*,*)"DM sparsity              = ",1.0D0-real(nnz)/real(norbs*norbs)
+    write(*,*)"Threshold                = ",bioham%threshold
+    write(*,*)"Number of replicas       = ",bioham%replicatex*bioham%replicatey*bioham%replicatez
+    write(*,*)"Number of atoms          = ",syf%nats
+    write(*,*)"Number of orbitals       = ",norbs
+    write(*,*)"Time for SP2             = ",mlssp2
+    write(*,*)"Time for Diagonalization = ",mlsdiag
+    write(*,*)"Speedup                  = ",mlsdiag/mlssp2
+    write(*,*)"Error DM                 = ",bml_fnorm(aux_bml)/real(norbs*norbs,dp)
+  endif
+
+  call bml_deallocate(aux_bml)
+  call bml_deallocate(rho_bml)
+
+end program biosolve

benchmarks/dmconstruction3d/input.in.bio

@@ -0,0 +1,11 @@
+
+BIOHAM{ 
+  BMLType= Ellpack
+  ReplicateX= 1
+  ReplicateY= 1
+  ReplicateZ= 1
+  Threshold= 1.0E-6
+  SystemFileName= 'prot.pdb'
+  ParamPath= '../../examples/latteTBparams/'
+}
+

benchmarks/dmconstruction3d/input.in.semicond

@@ -0,0 +1,15 @@
+MHAM{
+  NOrbs= size
+  BMLType= format
+  EpsilonA= -10.0
+  EpsilonB=  0.0
+  DeltaAB= -2.0
+  DeltaAiAj= -0.0
+  DeltaBiBj= -1.0
+  Decay= -1000.0
+  RCoeff= 0.0
+  Seed= 100
+  Reshuffle= F
+}
+
+