Commit 84c19986 authored by Victor's avatar Victor
Browse files

VY: Added Tomato tight-binding matrices generation.

parent bd721618
......@@ -11,6 +11,12 @@
/CB_toy_code/src/cb_davidson.x
/CB_toy_code/src/cb_davidson_rci.x
# /TB_toy_code/src
/TB_toy_code/src/*.a
/TB_toy_code/src/*.mod
/TB_toy_code/src/*.o
/TB_toy_code/src/tb_main.x
# /clib
/clib/*.a
/clib/*.o
......
......@@ -26,6 +26,10 @@ cb : bindir libfft libdavid libdavid_rci libcg libla liblapack
if test -d CB_toy_code ; then \
( cd CB_toy_code ; $(MAKE) TLDEPS= all || exit 1) ; fi
tb : bindir
if test -d TB_toy_code ; then \
( cd TB_toy_code ; $(MAKE) TLDEPS= all || exit 1) ; fi
###########################################################
# Auxiliary targets used by main targets:
# compile modules, libraries, directory for binaries, etc
......@@ -112,7 +116,7 @@ test-suite: touch-dummy
clean :
touch make.inc
for dir in \
CB_toy_code LAXlib FFTXlib KS_Solvers/Davidson KS_Solvers/CG clib \
CB_toy_code TB_toy_code LAXlib FFTXlib KS_Solvers/Davidson KS_Solvers/CG clib \
; do \
if test -d $$dir ; then \
( cd $$dir ; \
......
Fabiano Corsetti, Imperial College London, UK
* Email: fabiano.corsetti08 {at} imperial.ac.uk
* Homepage: <http://www.cmth.ph.ic.ac.uk/people/f.corsetti/>
Copyright (c) 2016, Fabiano Corsetti
All rights reserved.
Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
# Makefile for tomato_TB
sinclude ../make.inc
default: all
all: tb
tb:
( cd src ; $(MAKE) all || exit 1 )
tblib:
( cd src ; $(MAKE) tblib.a || exit 1 )
clean : examples_clean
( cd src ; $(MAKE) clean )
examples_clean:
if test -d examples ; then \
( cd examples ; ./clean_all ) ; fi
Tomato generates tight-binding toy matrices.
It works with the MatrixSwitch library.
Author of Tomato and MatrixSwitch: Fabiano Corsetti.
In the ./src directory you may find the source code of the libraries, and a simple program to generate toy matriecs using Tomato.
Currently the code only works in serial. The generated matrices are real.
Contact Victor Yu (wy29@duke.edu) for any question.
# Makefile for TB
include ../../make.inc
# MODFLAGS=$(MOD_FLAG)../../FFTXlib $(MOD_FLAG)../../LAXlib $(MOD_FLAG)../../KS_Solvers/Davidson $(MOD_FLAG)../../KS_Solvers/Davidson_RCI
TBOBJS = \
TB_main.o
TBLIBS = \
tomato_TB.o \
MatrixSwitch.o \
ms_mm_multiply.o \
ms_m_add.o \
ms_m_set.o \
ms_m_copy.o \
ms_m_register.o \
MatrixSwitch_ops.o
# QEMODS=../../FFTXlib/libqefft.a ../../KS_Solvers/Davidson_RCI/libdavid_rci.a ../../KS_Solvers/Davidson/libdavid.a ../../LAXlib/libqela.a ../../clib/clib.a
TLDEPS=bindir liblapack libblas
all : tb_main.x
tb_main.x : $(TBOBJS) tblib.a $(QEMODS)
$(LD) $(LDFLAGS) -o $@ \
$(TBOBJS) tblib.a $(QEMODS) $(LIBS)
- ( cd ../../bin; ln -fs ../TB_toy_code/src/$@ . ; )
tblib.a : $(TBLIBS)
$(AR) $(ARFLAGS) $@ $?
$(RANLIB) $@
tldeps :
if test -n "$(TLDEPS)" ; then \
( cd ../.. ; $(MAKE) $(TLDEPS) || exit 1 ) ; fi
clean :
- /bin/rm -f *.x *.o *.a *~ *_tmp.f90 *.d *.mod *.i *.L
include make.depend
# DO NOT DELETE
This diff is collapsed.
!==============================================================================!
!> @brief Storage and auxiliary operations for MatrixSwitch.
!==============================================================================!
module MatrixSwitch_ops
implicit none
public
!**** PARAMS ************************************!
integer, parameter :: dp=selected_real_kind(15,300) !< Double precision.
complex(dp), parameter :: cmplx_1=(1.0_dp,0.0_dp) !< Complex 1.
complex(dp), parameter :: cmplx_i=(0.0_dp,1.0_dp) !< Complex i.
complex(dp), parameter :: cmplx_0=(0.0_dp,0.0_dp) !< Complex 0.
!**** VARIABLES *********************************!
#if defined(__MPI)
character(1), save :: ms_lap_order !< Ordering of the BLACS process grid.
integer, save :: ms_mpi_comm !< MPI communicator.
integer, save :: ms_mpi_size !< Number of MPI processes.
integer, save :: ms_mpi_rank !< MPI process rank.
integer, save :: ms_lap_nprow !< Number of rows in the BLACS process grid.
integer, save :: ms_lap_npcol !< Number of columns in the BLACS process grid.
integer, save :: ms_lap_bs_def !< Default block size.
integer, save :: ms_lap_bs_num !< Number of block size exceptions.
integer, save :: ms_lap_icontxt !< BLACS context handle used by MatrixSwitch.
integer, allocatable, save :: ms_lap_bs_list(:,:) !< Block size exception list.
#endif
!**** TYPES *************************************!
!============================================================================!
!> @brief MatrixSwitch matrix type.
!!
!! This is the derived type that encapsulates all matrix storage
!! possibilities and hides the details from the user. Typically, the elements
!! below will never need to be accessed directly.
!============================================================================!
type matrix
character(3) :: str_type !< Label identifying the storage format.
logical :: is_initialized=.false. !< Has the matrix been initialized?
logical :: is_serial !< Is the matrix serial or parallel distributed?
logical :: is_real !< Is the matrix real or complex (both kind \p dp)?
logical :: is_square !< Is the matrix square?
logical :: is_sparse !< Is the matrix sparse?
logical :: iaux1_is_allocated=.false. !< Is iaux1 directly allocated or it is a pointer?
logical :: iaux2_is_allocated=.false. !< Is iaux2 directly allocated or it is a pointer?
logical :: iaux3_is_allocated=.false. !< Is iaux3 directly allocated or it is a pointer?
logical :: iaux4_is_allocated=.false. !< Is iaux4 directly allocated or it is a pointer?
logical :: dval_is_allocated=.false. !< Is dval directly allocated or it is a pointer?
logical :: zval_is_allocated=.false. !< Is zval directly allocated or it is a pointer?
integer :: dim1 !< Row dimension size of the matrix.
integer :: dim2 !< Column dimension size of the matrix.
integer, pointer :: iaux1(:) => null() !< Auxiliary information for certain storage formats.
integer, pointer :: iaux2(:) => null() !< Auxiliary information for certain storage formats.
integer, pointer :: iaux3(:) => null() !< Auxiliary information for certain storage formats.
integer, pointer :: iaux4(:) => null() !< Auxiliary information for certain storage formats.
real(dp), pointer :: dval(:,:) => null() !< Matrix elements for a real matrix.
complex(dp), pointer :: zval(:,:) => null() !< Matrix elements for a complex matrix.
end type matrix
!**** INTERFACES ********************************!
!============================================================================!
!> @brief \p opM parameter converter.
!!
!! Converts the input parameters \p opA and \p opB in the subroutines
!! \a mm_multiply and \a m_add from a character to a logical (real version)
!! or integer (complex version) for internal use:
!! \arg \c n / \c N mapped to \c .false. (real version) or \c 0 (complex
!! version)
!! \arg \c t / \c T mapped to \c .true. (real version) or \c 2 (complex
!! version)
!! \arg \c c / \c C mapped to \c .true. (real version) or \c 1 (complex
!! version)
!!
!! @param[in] seM Parameter
!! @param[out] luM
!============================================================================!
interface process_opM
module procedure process_lopM
module procedure process_iopM
end interface process_opM
!************************************************!
contains
!============================================================================!
!> @brief \p opM parameter converter (real version).
!============================================================================!
subroutine process_lopM(opM,trM)
implicit none
!**** INPUT ***********************************!
character(1), intent(in) :: opM
!**** INOUT ***********************************!
logical, intent(inout) :: trM
!**********************************************!
if ((opM .eq. 'T') .or. &
(opM .eq. 't') .or. &
(opM .eq. 'C') .or. &
(opM .eq. 'c')) then
trM=.true.
else if ((opM .eq. 'N') .or. &
(opM .eq. 'n')) then
trM=.false.
else
call die('process_lopM: invalid opM')
end if
end subroutine process_lopM
!============================================================================!
!> @brief \p opM parameter converter (complex version).
!============================================================================!
subroutine process_iopM(opM,tcM)
implicit none
!**** INPUT ***********************************!
character(1), intent(in) :: opM
!**** INOUT ***********************************!
integer, intent(inout) :: tcM
!**********************************************!
if ((opM .eq. 'T') .or. &
(opM .eq. 't')) then
tcM=2
else if ((opM .eq. 'C') .or. &
(opM .eq. 'c')) then
tcM=1
else if ((opM .eq. 'N') .or. &
(opM .eq. 'n')) then
tcM=0
else
call die('process_iopM: invalid opM')
end if
end subroutine process_iopM
!============================================================================!
!> @brief \p seC parameter converter.
!!
!! Converts the input parameter \p seC in the subroutine \a m_set from a
!! character to an integer for internal use:
!! \arg \c l / \c L mapped to \c 2
!! \arg \c u / \c U mapped to \c 1
!! \arg other: mapped to \c 0
!!
!! @param[in] seM Parameter
!! @param[out] luM
!============================================================================!
subroutine process_seM(seM,luM)
implicit none
!**** INPUT ***********************************!
character(1), intent(in) :: seM
!**** INOUT ***********************************!
integer, intent(inout) :: luM
!**********************************************!
if ((seM .eq. 'L') .or. &
(seM .eq. 'l')) then
luM=2
else if ((seM .eq. 'U') .or. &
(seM .eq. 'u')) then
luM=1
else
luM=0
end if
end subroutine process_seM
!============================================================================!
!> @brief Code termination from error.
!!
!! Outputs a custom error message to file and then stops execution with a
!! non-zero exit status.
!!
!! @param[in] message Error mesage to output before stopping.
!============================================================================!
subroutine die(message)
implicit none
!**** INPUT ***********************************!
character(*), intent(in), optional :: message
!**** INTERNAL ********************************!
integer :: err_unit
!**********************************************!
open(unit=err_unit,file='MatrixSwitch.err',status='replace')
write(err_unit,'(a)') 'FATAL ERROR in matrix_switch!'
if (present(message)) write(err_unit,'(a)') message
#if defined(__MPI)
write(err_unit,'(a,1x,i5)') 'MPI rank:', ms_mpi_rank
#endif
close(err_unit)
stop 1
end subroutine die
end module MatrixSwitch_ops
program TB_main
use MatrixSwitch ! Only for test matrices generation
implicit none
#if defined(__MPI) && defined(__SCALAPACK)
include "mpif.h"
#endif
! MatrixSwitch operation
character(5) :: m_storage
character(3) :: m_operation
character(128) :: arg1
! Test system info
integer :: n_basis
integer :: n_states
integer :: supercell(3)
real*8 :: frac_occ
real*8 :: sparsity
real*8 :: orb_r_cut
real*8 :: k_point(3)
! K-point sampling
integer :: i_k_point = 1
integer :: n_k_point = 1
integer :: i_coord
real*8 :: k_coord
! Matrix size
integer :: matrix_size
integer :: l_rows
integer :: l_cols
! Toy matrices
! complex*16, allocatable :: h(:,:)
! complex*16, allocatable :: s(:,:)
real*8, allocatable :: h(:,:)
real*8, allocatable :: s(:,:)
type(matrix) :: h_tmp
type(matrix) :: s_tmp
! MPI
integer :: n_proc
integer :: myid
#if defined(__MPI) && defined(__SCALAPACK)
integer :: nprow
integer :: npcol
integer :: myprow
integer :: mypcol
integer :: mpi_comm_global,mpierr
integer :: blk
integer :: BLACS_CTXT
integer, external :: numroc
#endif
#if defined(__MPI) && defined(__SCALAPACK)
! Initialize MPI
call MPI_Init(mpierr)
mpi_comm_global = MPI_COMM_WORLD
call MPI_Comm_size(mpi_comm_global,n_proc,mpierr)
call MPI_Comm_rank(mpi_comm_global,myid,mpierr)
#else
n_proc = 1
myid = 0
#endif
if(COMMAND_ARGUMENT_COUNT() == 1) then
call GET_COMMAND_ARGUMENT(1,arg1)
else
stop
endif
#if defined(__MPI) && defined(__SCALAPACK)
! Set up square-like processor grid
do npcol = nint(sqrt(real(n_proc))),2,-1
if(mod(n_proc,npcol) == 0) exit
enddo
nprow = n_proc/npcol
! Set block size
blk = 16
! Set up BLACS
BLACS_CTXT = mpi_comm_global
call BLACS_Gridinit(BLACS_CTXT,'r',nprow,npcol)
call BLACS_Gridinfo(BLACS_CTXT,nprow,npcol,myprow,mypcol)
call ms_scalapack_setup(mpi_comm_global,nprow,'r',blk,icontxt=BLACS_CTXT)
#endif
! Set parameters
#if defined(__MPI) && defined(__SCALAPACK)
m_storage = 'pddbc'
#else
m_storage = 'sdden'
#endif
m_operation = 'lap'
n_basis = 22
supercell = (/1,1,1/)
orb_r_cut = 0.5d0
k_point(1:3) = (/0.0d0,0.0d0,0.0d0/)
! Generate test matrices at this k-point
call tomato_TB(arg1,'silicon',.false.,frac_occ,n_basis,.false.,matrix_size,&
supercell,.false.,sparsity,orb_r_cut,n_states,.true.,k_point,&
.true.,0.0d0,h_tmp,s_tmp,m_storage,.true.)
#if defined(__MPI) && defined(__SCALAPACK)
l_rows = numroc(matrix_size,blk,myprow,0,nprow)
l_cols = numroc(matrix_size,blk,mypcol,0,npcol)
#else
l_rows = matrix_size
l_cols = matrix_size
#endif
allocate(h(l_rows,l_cols))
allocate(s(l_rows,l_cols))
h = h_tmp%dval
s = s_tmp%dval
if(myid == 0) then
write(*,*) "K-point", i_k_point, "generated."
endif
call m_deallocate(h_tmp)
call m_deallocate(s_tmp)
! Call the Davidson solver
! call Davidson(...)
deallocate(h)
deallocate(s)
#if defined(__MPI) && defined(__SCALAPACK)
call MPI_Finalize(mpierr)
#endif
end program TB_main
ms_mm_multiply.o : MatrixSwitch_ops.o
ms_m_add.o : MatrixSwitch_ops.o
ms_m_set.o : MatrixSwitch_ops.o
ms_m_copy.o : MatrixSwitch_ops.o
ms_m_register.o : MatrixSwitch_ops.o
MatrixSwitch.o : MatrixSwitch_ops.o
MatrixSwitch.o : ms_mm_multiply.o
MatrixSwitch.o : ms_m_add.o
MatrixSwitch.o : ms_m_set.o
MatrixSwitch.o : ms_m_copy.o
MatrixSwitch.o : ms_m_register.o
tomato_TB.o : MatrixSwitch.o
tomato_TB.o : MatrixSwitch_ops.o
TB_main.o : MatrixSwitch.o
!==============================================================================!
!> @brief Implementations of \a m_add.
!==============================================================================!
module MS_m_add
use MatrixSwitch_ops
implicit none
contains
!============================================================================!
!> @brief Matrix addition (simple dense, serial distribution, reference
!! implementation, real version).
!============================================================================!
subroutine m_add_sddenref(A,trA,C,alpha,beta)
implicit none
!**** INPUT ***********************************!
logical, intent(in) :: trA
real(dp), intent(in) :: alpha
real(dp), intent(in) :: beta
type(matrix), intent(in) :: A
!**** INOUT ***********************************!
type(matrix), intent(inout) :: C
!**** INTERNAL ********************************!
integer :: i, j
!**********************************************!
C%dval=beta*C%dval
if (.not. trA) then
C%dval=C%dval+alpha*A%dval
else if (trA) then
do i=1,C%dim1
do j=1,C%dim2
C%dval(i,j)=C%dval(i,j)+alpha*A%dval(j,i)
end do
end do
end if
end subroutine m_add_sddenref
!============================================================================!
!> @brief Matrix addition (simple dense, serial distribution, reference
!! implementation, complex version).
!============================================================================!
subroutine m_add_szdenref(A,tcA,C,alpha,beta)
implicit none
!**** INPUT ***********************************!
integer, intent(in) :: tcA
complex(dp), intent(in) :: alpha
complex(dp), intent(in) :: beta
type(matrix), intent(in) :: A
!**** INOUT ***********************************!
type(matrix), intent(inout) :: C
!**** INTERNAL ********************************!
integer :: i, j
!**********************************************!
C%zval=beta*C%zval
if (tcA==0) then