Hi Peter, Thank you so much for sharing your code with me. I did not understand what are the modifications needed for a periodic case: I should just take the transpose of the matrix?
Thanks very much, -Doron. On Sat, March 13, 2010 3:11 am, Peter Koval wrote: > Hello Doron! > > > I think you can use my program for this, although it is not perfect. The > program reads the .HSX file. At least one thing must be obviously > corrected in the program: the row and columns are interchanged, therefore > it might be necessary to interchange them in periodic systems. > > Best regards, > > > Peter > > > !! > !! > !! > subroutine import_hsx(fname, haux, saux, iv) use modul_log use > modul_precision use modul_orbital_vars, only : norbitals, nspin, orb_occ, > Temp, > total_electronic_charge !! external > character(len=*), intent(in) :: fname real(dp) :: haux(norbitals, > norbitals, nspin), saux(norbitals, norbitals); integer, intent(in) :: iv > > > !! internal > integer, allocatable :: int_buff(:) !! buffer for pointers (to nonzero > elements) within a column real(sp), allocatable :: sp_buff(:) !! > buffer for vector values (of nonzero elements) within a column integer :: > ifile, iostat, icol, i, ispin, sum_col2nnzero, maxnnzero integer(4) :: > norbitals_in, norbitals_in_supercell, nspin_in logical(4) :: gamma > ! Gamma point or not > > > ifile = get_free_handle(); > > open(ifile,file=trim(fname),form='unformatted',action='read',status='old' > ,iostat=iostat); > if(iostat/=0) then; write(ilog,*)'import_hsx: error: file ', trim(fname), " > ?"; stop; endif; > rewind(ifile) read(ifile,iostat=iostat) norbitals_in, > norbitals_in_supercell, nspin_in, nnonzero if (norbitals /= norbitals_in) > then write(ilog,*)norbitals_in, norbitals_in_supercell, nspin_in, nnonzero > write(ilog,*)"import_hsx: norbitals, norbitals_in:", norbitals, > norbitals_in stop "import_hsx: (norbitals /= norbitals_in)" endif if > (nspin_in /= nspin) then > write(ilog,*) "import_hsx: nspin, nspin_in:", nspin, nspin_in stop > "import_hsx: (nspin /= nspin_in)" > endif if (nspin==1) then; orb_occ=2; else if (nspin==2) then; orb_occ=1; > else; write(ilog,*)'import_hsx: nspin', nspin; stop '(nspin/=1 .and. > nspin/=2)'; endif > > if(iv>1) write(ilog,*) "import_hsx: norbitals_in, norbitals_in_supercell, > nspin_in, nnonzero" if(iv>1) write(ilog,*) norbitals_in, > norbitals_in_supercell, nspin_in, nnonzero > > read(ifile,iostat=iostat) gamma if (.not. gamma)stop "import_hsx: .not. > gamma"; > > !! allocate the buffers > allocate(col2nnzero(norbitals), col2displ(norbitals), > sparse_ind2row(nnonzero)) allocate(H_sparse(nnonzero,nspin), stat=iostat); > ! Hamiltonian matrix in > packed form allocate(S_sparse(nnonzero), stat=iostat); ! Overlap > matrix in packed form > > read(ifile,iostat=iostat)col2nnzero sum_col2nnzero = sum(col2nnzero) if > (sum_col2nnzero > nnonzero) then > write(ilog,*) 'import_hsx: sum_col2nnzero > nnonzero ', sum_col2nnzero, > nnonzero; write(ilog,*) col2nnzero; stop 'import_hsx:'; endif > > maxnnzero = maxval(col2nnzero) allocate(int_buff(maxnnzero), > sp_buff(maxnnzero)); > > !! Fill the displacements (according to col2nnzero) col2displ > col2displ(1)=0 do icol=2, norbitals col2displ(icol) = col2displ(icol-1) + > col2nnzero(icol-1) enddo > > !! Fill the rows for each index in *_sparse arrays > do icol=1, norbitals read(ifile,iostat=iostat)int_buff(1:col2nnzero(icol)) > ! read set of rows > where nonzero elements reside if (iostat/=0) stop "import_hsx: (iostat/=0) > int_buff" > > do i=1, col2nnzero(icol) sparse_ind2row(col2displ(icol)+i) = int_buff(i) > enddo enddo > > !! Read the data to H_sparse array > do ispin=1,nspin do icol=1,norbitals > read(ifile,iostat=iostat)sp_buff(1:col2nnzero(icol)) > if (iostat /= 0) stop "import_hsx: (iostat/=0) Hamiltonian matrix" do > i=1,col2nnzero(icol); H_sparse(col2displ(icol)+i, ispin) = sp_buff(i); > enddo; enddo enddo > > !! Read the data to S_sparse array > do icol=1,norbitals read(ifile,iostat=iostat)sp_buff(1:col2nnzero(icol)) > if (iostat /= 0) stop "import_hsx: (iostat/=0) overlap matrix" do > i=1,col2nnzero(icol); S_sparse(col2displ(icol)+i) = sp_buff(i); enddo > enddo > > do ispin=1, nspin call sparse2full(norbitals, Haux(:,:,ispin), > H_sparse(:,ispin), > col2nnzero, col2displ, sparse_ind2row); if(ispin==1) & call > sparse2full(norbitals, Saux, S_sparse(:), col2nnzero, col2displ, > sparse_ind2row); enddo > > read(ifile,iostat=iostat) total_electronic_charge, Temp ! Total > electronic charge and Temperature if(iv>0)write(ilog,*) "import_hsx: > total_electronic_charge, Temp (Ry):", real(total_electronic_charge,4), > real(Temp,4) deallocate(int_buff, sp_buff); close(ifile); > > end subroutine !import_hsx > > !! > !! > !! > subroutine sparse2full(ndim, M_full, M_sparse, col2nnzero, col2displ, > sparse_ind2row) use modul_precision !! external > integer(4), intent(in) :: ndim real(dp), intent(out) :: M_full(ndim,ndim) > real(dp), intent(in) :: M_sparse(:) integer(4), intent(in) :: > col2nnzero(ndim), col2displ(ndim), sparse_ind2row(:) > > > !! internal > integer :: icol, i, irow, sparse_ind > > do icol=1,ndim do i=1,col2nnzero(icol); sparse_ind = col2displ(icol)+i; irow > = sparse_ind2row(sparse_ind); > M_full(irow, icol) = M_sparse(sparse_ind) > enddo enddo > > end subroutine !sparse2full > > > On Fri, Mar 12, 2010 at 7:41 PM, Doron Naveh <[email protected]> wrote: > > >> Hi, >> I'm trying to obtain the overlap matrix of basis set functions, >> does anyone know how? Thanks, >> Doron. >> >> >> > > > -- > Dr. Peter Koval > email: [email protected] > inet: http://sites.google.com/site/kovalpeter/ > >
