hamiltonian_cartesian_8band.c++

Go to the documentation of this file.

/*
 *An 8-band Hamiltonian on a Cartesian grid.

    Copyright (C) 2010 Joseph Pingenot

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU Affero General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Affero General Public License for more details.

    You should have received a copy of the GNU Affero General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.

*/

#include <string.h>
#include <liblapackwrap/lapackwrap_real_double.h>
#include <liblapackwrap/lapackwrap_complex_double.h>
#include <liblapackwrap/blaswrap_complex_double.h>
#include <liblapackwrap/blaswrap_real_double.h>
#include <libkdotp/hamiltonian_cartesian_8band.h++>
#include <libmodelxx/grid.h++>
#include <libkdotp/physical_constants.h++>
#include <libmodelxx/cartesian_grid.h++>
#include <libkdotp/physical_constants.h++>
#include <libmodelxx/material.h++>

using namespace kdp::constants;

#define N_BANDS 8

namespace kdotp {
  namespace libkdotp {
    namespace hamiltonian {

      typedef libmodelxx::grid::cartesian_grid<libmodelxx::structure::material*,1,1> matgrid;
      typedef libmodelxx::grid::cartesian_grid<double,1,1> potgrid;

      /**\note This matrix contains all of the *constants* of the 8-band
       *k-dot-p matrix. Any non-constants across the nanostructure (e.g.
       *the coupling constants P, gamma_i, m', etc.) are dependent upon
       *position and, therefore, not included (they're not constant across
       *the entire structure (nor can they be determined at compile time)
       */
      const static struct llw_complex_double consts_8band[] = {
  {-hbar*hbar/2.0, 0},
  {0, 0},
  {-1/__builtin_sqrt(2), 0},
  {__builtin_sqrt(2/3), 0},
  {1/__builtin_sqrt(6), 0},
  {0, 0},
  {-1/__builtin_sqrt(3), 0},
  {-1/__builtin_sqrt(3), 0},
  
  {0, 0},
  {-hbar*hbar/2.0, 0},
  {0, 0,},
  {-1/__builtin_sqrt(2), 0},
  {__builtin_sqrt(2/3), 0},
  {1/__builtin_sqrt(2), 0},
  {-1/__builtin_sqrt(3), 0},
  {1/__builtin_sqrt(3), 0},

  {-1/__builtin_sqrt(2), 0},
  {0, 0},
  /**\notehbar^2/2 for both k_|| and k_z*/
  {-hbar*hbar/(2*m0), 0},
  {2*__builtin_sqrt(3)*hbar*hbar/(2*m0), 0},
  {__builtin_sqrt(3)*hbar*hbar/(2*m0), 2*__builtin_sqrt(3)*hbar*hbar/(2*m0)},
  {0, 0},
  {-__builtin_sqrt(6)*hbar*hbar/(2*m0), 0},
  {-__builtin_sqrt(6)*hbar*hbar/(2*m0), 2*__builtin_sqrt(6)*hbar*hbar/(2*m0)},
  
  {__builtin_sqrt(2/3), 0},
  {-__builtin_sqrt(6), 0},
  {2*__builtin_sqrt(3)*hbar*hbar/(2*m0), 0},
  {-hbar*hbar/(2*m0), 0},
  {0, 0},
  {__builtin_sqrt(3)*hbar*hbar/(2*m0), 2*__builtin_sqrt(3)*hbar*hbar/(2*m0)},
  {-__builtin_sqrt(2)*hbar*hbar/(2*m0), 0},
  {-3*__builtin_sqrt(2)*hbar*hbar/(2*m0), 0},
  
  {1/__builtin_sqrt(6), 0},
  {__builtin_sqrt(2/3), 0},
  {__builtin_sqrt(3)*hbar*hbar/(2*m0), -2*__builtin_sqrt(3)*hbar*hbar/(2*m0)},
  {0, 0},
  {-hbar*hbar/(2*m0), 0},
  {-2*__builtin_sqrt(3)*hbar*hbar/(2*m0), 0},
  {3*__builtin_sqrt(2)*hbar*hbar/(2*m0), 0},
  {__builtin_sqrt(2)*hbar*hbar/(2*m0), 0},

  {0, 0},
  {1/__builtin_sqrt(2), 0},
  {0, 0},
  {__builtin_sqrt(3)*hbar*hbar/(2*m0), 2*__builtin_sqrt(3)*hbar*hbar/(2*m0)},
  {-2*__builtin_sqrt(3)*hbar*hbar/(2*m0), 0},
  {__builtin_sqrt(6)*hbar*hbar/(2*m0), -2*__builtin_sqrt(6)*hbar*hbar/(2*m0)},
  {-__builtin_sqrt(6)*hbar*hbar/(2*m0), 0},
  
  {-1/__builtin_sqrt(3), 0},
  {-1/__builtin_sqrt(3), 0},
  {-__builtin_sqrt(6)*hbar*hbar/(2*m0), 0},
  {-__builtin_sqrt(2)*hbar*hbar/(2*m0), 0},
  {3*__builtin_sqrt(2)*hbar*hbar/(2*m0), 0},
  {__builtin_sqrt(6)*hbar*hbar/(2*m0), 2*__builtin_sqrt(6)*hbar*hbar/(2*m0)},
  {-hbar*hbar/(2*m0), 0},
  {0, 0},
  
  {-1/__builtin_sqrt(3), 0},
  {1/__builtin_sqrt(3), 0},
  {-__builtin_sqrt(6)*hbar*hbar/(2*m0), 2*__builtin_sqrt(6)*hbar*hbar/(2*m0)},
  {3*__builtin_sqrt(2)*hbar*hbar/(2*m0), 0},
  {__builtin_sqrt(2)*hbar*hbar/(2*m0), 0},
  {-__builtin_sqrt(6)*hbar*hbar/(2*m0), 0},
  {0, 0},
  {-hbar*hbar/(2*m0), 0}
      };

      /**Simple helper function to get the appropriate constant out of consts_8band[]
       *\note entries are stored in column-major order.
       *\note rows, cols are in the order of page 215 in Winkler's spin-orbit book
       */
      static inline short get_consts_index(short unsigned row, short unsigned col) {
  return row*8 + col;
      }

      /**does a quick check to make sure the constants matrix is truly Hermitean
       *\param eye pointer to a short unsigned int to store the row index of the invalid constant.
       *\return -j if a real portion was not identical across a transpose
       *        +j if a complex portion was not negated across a transpose
       *        0 if identical
       *\note this is a debugging/safety aid, but costs basically no time.
       */
      static int check_constants_hermiticity(unsigned short *eye) {
  for(unsigned short i=0; i<8; i++) {
    for(unsigned short j=0; j<8; j++) {
      /*
       *This is buggy and this file currently unused, so it's been commented out so that things will build.
      if(consts_8band[get_consts_index(i, j)].re 
         != consts_8band[get_consts_index(i, j)].re) {
        *eye = i;
        return -j;
      }
      if(consts_8band[get_consts_index(i, j)].im 
         != -consts_8band[get_consts_index(i, j)].im) {
        *eye = i;
        return j;
      }
      */
    }
  }
  return 0;
      }

      static GError* store_reduced_material_parameters(struct parameters_8band_nostrain *params, matgrid *mg, potgrid *pg) {
  int e;
  char *s;
  const char* paramstring[] = {"E_g_Gamma", "E_v", "Delta", "m_e_Gamma", "gamma1", "gamma2", "gamma3", "Ep"};
  double temp_params[8];
  for(long unsigned i=0; i<mg->L[0]; i++) {
    for(short unsigned p=0; p<8; p++) {
      temp_params[i] = mg->data[i]->get_property(paramstring[p], &e);
      if(e != 0) {
        return g_error_new(HAMILTONIAN_CARTESIAN_8BAND_GERROR_DOMAIN, 10+e, "Got error %d getting property (%s) at site %lu (material %s).\n", e, paramstring[p], i, s=mg->data[i]->get_name());
      }
    }
    params[i].Ec = temp_params[0]+temp_params[1];
    params[i].Ev = temp_params[1];
    params[i].Eso = temp_params[1]-temp_params[2];
    params[i].P = sqrt(hbar*hbar/(2*m0)*temp_params[7]);
    /*For the reduced param terms, see Winkler's Spin-orbit book page 216*/
    params[i].oo_mstar = 1.0/temp_params[4] - temp_params[7]/(3*m0)*(2/(temp_params[0]+temp_params[1]) + 1/(temp_params[0]+temp_params[1]+temp_params[2]));
    params[i].gamma1 = temp_params[4] - temp_params[7]/(3*(temp_params[0] + temp_params[1]));
    params[i].gamma2 = temp_params[5] - temp_params[7]/(6*(temp_params[0] + temp_params[1]));
    params[i].gamma3 = temp_params[6] - temp_params[7]/(6*(temp_params[0] + temp_params[1]));
  }
  return NULL;
      }


      cartesian_8band::cartesian_8band(struct libmodelxx::structure::generic_structure *gs, boundary_condition bc_, bool second_neighbor_) : upper_diagonal(NULL), diagonal(NULL), matrix(NULL), Q(NULL), bc(bc_), second_neighbor(second_neighbor_), params(NULL) {
  n_neighbors=second_neighbor_?2:1;
  /*Check to make sure our consts matrix is Hermitean*/
  {
    int col;
    short unsigned row;
    if((col=check_constants_hermiticity(&row)) != 0) {
      fprintf(stderr, "Hermiticty check of 8-band constants FAILED at row=%hu col=%d.\n", row, col);
    }else{
      fprintf(stderr, "Hermiticty check of 8-band constants PASSED.\n");
    }
  }
  matgrid *mg = (matgrid*)(gs->material_grid);
  potgrid *pg = (potgrid*)(gs->potential_grid);
  L = mg->L[0];
  /*Load the reduced material parameters into the array.*/
  params = (struct parameters_8band_nostrain*)malloc(L*sizeof(struct parameters_8band_nostrain));
  if(params == NULL) return;
  GError *e = NULL;
  if((e = store_reduced_material_parameters(params, mg, pg)) != NULL) {
    fprintf(stderr, "hamiltonian_cartesian_8band: ERROR calculating reduced band parameters: %s\n", e->message);
    g_error_free(e);
    free(params);
    params =NULL;
    return;
  }
        if(bc == HARDWALL) {
    matrix = new llw_complex_double[(n_neighbors+1)*L*N_BANDS];
      for(long unsigned i = 0; i<(n_neighbors+1)*L*N_BANDS; i++) {
        matrix[i].re = 0;
        matrix[i].im = 0;
      }
  }else{
    matrix = new llw_complex_double[L*L*N_BANDS*N_BANDS];
    for(long unsigned int i=0; i<(L*L*N_BANDS*N_BANDS); i++) {
      matrix[i].re = 0;
      matrix[i].im = 0;
    }
  }
  int e2=0;
        //double dconsts = -1/(nm_to_bohr(mg->a[0])*nm_to_bohr(mg->a[0]));
  /*Set up the constants for the diagonal, superdiagonal, and second superdiagonal for first nearest neighbor*/
  
  double ooms=0, ooms_prev=0, ooms_next=0, Eg=0, Ev=0, ooms_prev2=0, ooms_next2=0;
  /*Initialize the current ooms.*/
  ooms = 1.0/mg->data[0]->get_property("m_e_Gamma", &e2);
  if(e2 != 0) {
    fprintf(stderr, "hamiltonians::cartesian_8band: ERROR getting m_e_Gamma (err=%d).\n", e2);
    e2=0;
  }
  ooms_next = 1.0/mg->data[1]->get_property("m_e_Gamma", &e2);
  if(e2 != 0) {
    fprintf(stderr, "hamiltonians::cartesian_8band: ERROR getting m_e_Gamma (err=%d).\n", e2);
    e2=0;
  }
  for(unsigned int i=0; i<L; i++, ooms_prev2 = ooms_prev, ooms_prev=ooms, ooms=ooms_next, ooms_next = ooms_next2) {
    //Eg = ev_to_hartree(mg->data[i]->get_property("E_g_Gamma", &e2));
    Eg = mg->data[i]->get_property("E_g_Gamma", &e2);
    if(e2 != 0) {
      fprintf(stderr, "hamiltonians::cartesian_8band: ERROR getting E_g_Gamma (err=%d).\n", e2);
      e2=0;
    }
    //Ev = ev_to_hartree(mg->data[i]->get_property("E_v", &e2));
    Ev = mg->data[i]->get_property("E_v", &e2);
    if(e2 != 0) {
      fprintf(stderr, "hamiltonians::cartesian_8band: ERROR getting E_v (err=%d).\n", e2);
      e2=0;
    }
    if(i < L-2) {
      ooms_next2 = 1.0/mg->data[i+2]->get_property("m_e_Gamma", &e2);
      if(e2 != 0) {
        fprintf(stderr, "hamiltonians::cartesian_8band: ERROR getting m_e_Gamma (err=%d).\n", e2);
        e2=0;
      }
    }
    char *s;
    fprintf(stderr, "%05u(%s): ooms_prev2=%g ooms_prev=%g ooms=%g ooms_next=%g ooms_next2=%g Ev+Eg=%g\n", i, s=mg->data[i]->get_name(), ooms_prev2, ooms_prev, ooms, ooms_next, ooms_next2, Ev+Eg);
    g_free(s);
    ////////////////////////////////////////////////////////////////
    // NOTHING BELOW THIS LINE HAS BEEN FULLY IMPLEMENTED FOR 8BAND;
    // THIS IS ALL PLACEHOLDER/FORM CODE WHICH HAS BEEN MADE TO
    // COMPILE.
    ////////////////////////////////////////////////////////////////
    if(bc == HARDWALL) {
    }else{
    }
  }
      }
      
      cartesian_8band::~cartesian_8band(){
  if(matrix != NULL) delete[] matrix;
  if(diagonal != NULL) delete[] diagonal;
  if(upper_diagonal != NULL) delete[] upper_diagonal;
  if(Q != NULL) delete[] Q;
  if(matrix != NULL) delete[] matrix;
      }

      double* cartesian_8band::get_eigenvectors(char VoI, double eVmin, double eVmax, int imin, int imax, double **evals, int *n_eigs, GError **err)  {
  int info;
  if ((bc == PERIODIC) && (diagonal == NULL) && (upper_diagonal == NULL) && (Q == NULL)) {
    /*If we're not hardwall, we first have to reduce to tridiagonal.*/
    /*We're a symmetric matrix, but not tridiag.
     *Call DSYTRD to get the tri-diag matrix*/
    /*Call DORGTR to get the matrix Q to convert the eigenfunctions back to the position basis.*/
  } else if (second_neighbor && (diagonal == NULL) && (upper_diagonal == NULL) && (Q == NULL)) {
    /*We're a banded matrix.*/
    /*Call DSBTRD to get the tri-diagonal matrix and Q*/
    info = llw_complex_double_sbtrd_a('V', 'U', L, n_neighbors, matrix, n_neighbors+1, &diagonal, &upper_diagonal, &Q);
    if(info != 0) {
      *err = g_error_new(HAMILTONIAN_CARTESIAN_8BAND_GERROR_DOMAIN, 1, "Got error from llw_complex_double_sbtrd_a: %d\n", info);
      if(diagonal != NULL) {
        free(diagonal);
        diagonal = NULL;
      }
      if(upper_diagonal != NULL) {
        free(upper_diagonal);
        upper_diagonal = NULL;
      }
      if(Q != NULL) {
        free(Q);
        Q = NULL;
      }
      return NULL;
    }
  }
  double *eigvecs = NULL;
  int n_eigvecs = -1;
  int* isuppz = NULL;
  //info = llw_complex_double_stegr_a('V', VoI, L, diagonal, upper_diagonal, eVmin, eVmax, imin, imax, n_eigs, evals, &eigvecs, &n_eigvecs, &isuppz);
  if((n_eigvecs <= *n_eigs) || (info != 0)) {
    if(n_eigvecs <= *n_eigs) {
      *err = g_error_new(HAMILTONIAN_CARTESIAN_8BAND_GERROR_DOMAIN, 6, "in return from llw_complex_double_stegr_a: failed to get same number or more eigenvec storage than there were eigenvals (%d and %d, respectively)", n_eigvecs, *n_eigs);
    }else{
      if(info < 0) {
        *err = g_error_new(HAMILTONIAN_CARTESIAN_8BAND_GERROR_DOMAIN, 1, "Bad argument to LAPACK call DSTEGR: %d\n", -info);
      }else if(info < 10) {
        *err = g_error_new(HAMILTONIAN_CARTESIAN_8BAND_GERROR_DOMAIN, 2, "Error in BLAS call DLARRE called by LAPACK call DSGEGR: %d\n", info);
      }else if(info < 20) {
        *err = g_error_new(HAMILTONIAN_CARTESIAN_8BAND_GERROR_DOMAIN, 3, "Error in BLAS call DLARRV called by LAPACK call DSGEGR: %d\n", info/2);
      }else if(info == 100) {
        *err = g_error_new(HAMILTONIAN_CARTESIAN_8BAND_GERROR_DOMAIN, 5, "Error in liblapackwrap llw_complex_double_stegr_a: malloc failed (%d)\n", info/100);
      }else{
        *err = g_error_new(HAMILTONIAN_CARTESIAN_8BAND_GERROR_DOMAIN, 4, "Got error (info=%d) from DSTEGR, but don't know what it means.\n", info);
      }
    }
    free(*evals);
    free(eigvecs);
    *evals=NULL;
    eigvecs=NULL;
  }else{
    free(isuppz);
    if(Q != NULL) {
      double* final_eigvecs = NULL;
      double one = 1.0;
      //info = llw_complex_double_gemm_a('N', 'N', L, n_eigvecs, L, &one, Q, L, eigvecs, L, &final_eigvecs);
      if(info != 0) {
        *err = g_error_new(HAMILTONIAN_CARTESIAN_8BAND_GERROR_DOMAIN, 4, "Got error (info=%d) from llw_complex_double_gemmm_a\n", info);
        if(final_eigvecs != NULL) free(final_eigvecs);
        return eigvecs;
      }
      free(eigvecs);
      return final_eigvecs;
    }
  }
  return eigvecs;
      }


    }
  }
}