hamiltonian_cartesian_effective_mass.c++

Go to the documentation of this file.

/*
** hamiltonain_effective_mass.c++
**
** Made by Johnny Q. Hacker
** Login   <solarion@nathan>
**

Copyright (C) 2009 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/>.
** Started on  Fri Jul  3 06:18:08 2009 Johnny Q. Hacker
// Last update Thu Oct  1 21:44:48 2009 Johnny Q. Hacker
*/

#include <string.h>
#include <liblapackwrap/blaswrap_real_double.h>
#include <liblapackwrap/lapackwrap_real_double.h>
#include <libkdotp/hamiltonian_cartesian_effective_mass.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;

namespace kdotp {
  namespace libkdotp {
    namespace hamiltonian {

      typedef libmodelxx::grid::cartesian_grid<libmodelxx::structure::material*,1,1> matgrid_t;
      typedef libmodelxx::grid::cartesian_grid<double,1,1> potgrid_t;

      cartesian_effective_mass::~cartesian_effective_mass() {
  /*We do a bunch of these frees.  Cluster them.  If they're not NULL, free them.
   *Note that these are pointers to pointers; we keep going until the array entry is NULL (if it's a real member of this class, it won't be NULL
   */
  void** things_to_free[] = {
    (void**)&upper_diagonal,
    (void**)&diagonal,
    (void**)&matrix,
    (void**)&Q,
    (void**)&internal_potential,
    /*don't free matgrid or potgrid; they are still in use!*/
    NULL
  };
  for(int i=0; things_to_free[i] != NULL; i++) {
    if(*things_to_free[i] != NULL) {
      free(*things_to_free[i]);
      *things_to_free[i] = NULL;
    }
  }
      }

      cartesian_effective_mass::cartesian_effective_mass(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_), 
    internal_potential(NULL),
    matgrid(NULL),
    potgrid(NULL)
      {
  n_neighbors=second_neighbor_?2:1;
  matgrid = (matgrid_t*)(gs->material_grid);
  potgrid = (potgrid_t*)(gs->potential_grid);
  L = matgrid->L[0];
  /*Allocate and set the internal potential when we're constructed.*/
  internal_potential = (double*)malloc(sizeof(double)*L);
  /*The internal potential is initially set to 0*/
  for(unsigned int i=0; i<L; i++) internal_potential[i] = 0;
        if(bc == HARDWALL) {
    if(second_neighbor) {
      matrix = (double*)malloc(sizeof(double)*(n_neighbors+1)*L);
      if(matrix == NULL) {
        fprintf(stderr, "Error allocating matrix.\n");
        return;
      }
    }else{
      upper_diagonal = (double*)malloc(sizeof(double)*(L-1));
      if(upper_diagonal == NULL) {
        fprintf(stderr, "Error allocating upper diagonal.\n");
        return;
      }
      diagonal = (double*)malloc(sizeof(double)*L);
      if(diagonal == NULL) {
        fprintf(stderr, "Error allocating diagonal.\n");
        return;
      }
    }
  }else{
    matrix = (double*)malloc(sizeof(double)*L*L);
    if(matrix == NULL) {
      fprintf(stderr, "Error allocating matrix.\n");
      return;
    }
  }
  dx = matgrid->a[0];
  properties = (struct effective_mass_properties*)malloc(sizeof(struct effective_mass_properties)*L);
  if(properties == NULL) {
    fprintf(stderr, "Error allocating properties.\n");
    return;
  }
  int e2=0;
  double Eg, Ev;
  Eg = Ev = 0;
  for(unsigned int i=0; i<L; i++) {
    properties[i].ooms = 1.0/matgrid->data[i]->get_property("m_e_Gamma", &e2);
      if(e2 != 0) {
        fprintf(stderr, "hamiltonians::cartesian_effective_mass: ERROR getting m_e_Gamma (err=%d).\n", e2);
        e2=0;
      }
      properties[i].epsilon = matgrid->data[i]->get_property("epsilon_0", &e2);
      if(e2 != 0) {
        fprintf(stderr, "hamiltonians::cartesian_effective_mass: ERROR getting epsilon_0 (err=%d).\n", e2);
        e2=0;
      }
      Eg = matgrid->data[i]->get_property("E_g_Gamma", &e2);
      if(e2 != 0) {
        fprintf(stderr, "hamiltonians::cartesian_effective_mass: ERROR getting E_g_Gamma (err=%d).\n", e2);
        e2=0;
      }
      Ev = matgrid->data[i]->get_property("E_v", &e2);
      if(e2 != 0) {
        fprintf(stderr, "hamiltonians::cartesian_effective_mass: ERROR getting E_v (err=%d).\n", e2);
        e2=0;
      }
      properties[i].Ec = Ev + Eg;
      properties[i].potential = potgrid->data[i];
  }
  set_up_matrix();
      }

      void cartesian_effective_mass::update_matrix_with_internal_potential(double* old_potential) {
  /*Remember when updating that the pot for an electron is -1 V[i].
   *Therefore, when updating, we *add* the old one to the previous diagonal terms to remove it,
   *  and the *subtract* the new one (or, subtract the difference, as below; saves an op)
   */
  for(unsigned int i=0; i<L; i++) {
    if(bc == HARDWALL) {
      if(!second_neighbor) {
        diagonal[i] -= internal_potential[i] - old_potential[i];
      }else{
        matrix[get_array_index(i, i)] -=  internal_potential[i] - old_potential[i];
      }
    }else{
      matrix[get_matrix_index(i,i)] -= internal_potential[i] - old_potential[i];
    }
  }
      }

      void cartesian_effective_mass::replace_internal_potential(double* new_potential) {
  double *old_potential = internal_potential;
        internal_potential = new_potential;
  update_matrix_with_internal_potential(old_potential);
  free(old_potential);
      }

      void cartesian_effective_mass::set_up_matrix() {
  /*Set the values of our matrix members to 0 to ensure they're properly initialized before we go setting specific numbers*/
        if(bc == HARDWALL) {
    if(second_neighbor) {
      for(long unsigned i = 0; i<(n_neighbors+1)*L; i++) matrix[i] = 0;
    }else{
      for(unsigned int i=0; i<(L-1); i++) {
        upper_diagonal[i] = diagonal[i] = 0;
      }
      diagonal[L-1]=0;
    }
  }else{
    for(long unsigned int i=0; i<(L*L); i++) matrix[i] = 0;
  }
  int e2=0;
  /*hbar^2/d^2 in the notebook on April 29*/
        double dconsts = -hbar*hbar/(m0*dx*dx)*1e18;
        //double dconsts = -1/(nm_to_bohr(matgrid->a[0])*nm_to_bohr(matgrid->a[0]));
  /*Set up the constants for the diagonal, superdiagonal, and second superdiagonal for first nearest neighbor*/
  //double diag_const=-dconsts;
  double diag_const=-dconsts/4.0;
  double sd1_const=dconsts/4.0;
  double sd2_const=0;
  if(second_neighbor) {
    /*Modify our consts for second-nearest-neighbor*/
    dconsts /= 48;
    diag_const=-60*dconsts;
    sd1_const=16*dconsts;
    sd2_const=-dconsts;
  }
  
  //fprintf(stderr, "dconsts: hbar=%g m0=%g a[0]=%g dconsts=%g\n", hbar, m0, dx, dconsts);
  double ooms=0, ooms_prev=0, ooms_next=0, ooms_prev2=0, ooms_next2=0;
  /*Initialize the current ooms.*/
  ooms = properties[0].ooms;
  ooms_next = properties[1].ooms;
  for(unsigned int i=0; i<L; i++, ooms_prev2 = ooms_prev, ooms_prev=ooms, ooms=ooms_next, ooms_next = ooms_next2) {
    if(i < L-2) {
      ooms_next2 = properties[i+2].ooms;
    }
    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, "(?)", ooms_prev2, ooms_prev, ooms, ooms_next, ooms_next2, properties[i].Ec);
    if(bc == HARDWALL) {
      if(!second_neighbor) {
        /*Note the minus sign, because the potential is acting on an electron.*/
        //diagonal[i] = Eg + Ev - potgrid->data[i] + diag_const*ooms;
        diagonal[i] = properties[i].Ec - properties[i].potential - internal_potential[i] + diag_const*(2*ooms + ooms_next + ooms_prev);
        //fprintf(stderr, "%05d: Eg+Ev=%g ooms_prev=%g ooms=%g d=%g ", i, Eg+Ev, ooms_prev, ooms, diagonal[i]);
        if(i != L-1) {
    /*note that ooms_next will not get updated on the last site (address would be out of the array)*/
    upper_diagonal[i] = sd1_const*(ooms + ooms_next);
    //fprintf(stderr, "ooms_next=%g ud=%g\n", ooms_next, upper_diagonal[i]);
              //}else{
    //fprintf(stderr, "(last; no ud)\n");
        }
      }else{
        /////////////////////////////////////////////////////////////////
        //TODO: Need to fill in the appropriate constants.
        /////////////////////////////////////////////////////////////////
        /*First, the diagonal*/
        matrix[get_array_index(i, i)] = properties[i].Ec - properties[i].potential - internal_potential[i] + diag_const*ooms;
        if(i != L-1) {
    if(i != L-2) {
      /*Second superdiagonal*/
      matrix[get_array_index(i, i+2)] = sd2_const*(ooms + ooms_next2);
    }
    /*Superdiagonal*/
    matrix[get_array_index(i, i+1)] = sd1_const*(ooms + ooms_next);
        }
        //fprintf(stderr, "%05u: d=%g u=%g u2=%g\n", i, matrix[get_array_index(i, i)], matrix[get_array_index(i, i+1)], matrix[get_array_index(i, i+2)]);
      }
    }else{
      /////////////////////////////////////////////////////////////////
      //TODO: Need to fill in the appropriate constants.
      /////////////////////////////////////////////////////////////////
      /*Mostly same as for hardwall, but now we have an explicit matrix*/
      matrix[get_matrix_index(i,i)] = properties[i].Ec - properties[i].potential - internal_potential[i] + diag_const*ooms;
      if(i == 0) {
        /*Previous ooms is the final gridsite, since we have periodic BCs*/
        ooms_prev = properties[L-1].ooms;
        matrix[get_matrix_index(i, i+1)] = sd1_const*(ooms + ooms_next);
        matrix[get_matrix_index(i, L-1)] = sd1_const*(ooms + ooms_prev);
        if(second_neighbor) {
    ooms_prev2 = properties[L-2].ooms;
    matrix[get_matrix_index(i, i+2)] = sd2_const*(ooms + ooms_next2);
    matrix[get_matrix_index(i, L-2)] = sd2_const*(ooms + ooms_prev2);
        }
      }else if(i == L-1) {
        /*Next ooms is the *first* gridsite
         *We've already gotten the previous value (in the for loop)*/
        ooms_next = properties[0].ooms;
        /*Previous value is the final gridsite, since we have periodic BCs*/
        matrix[get_matrix_index(i, 0)] = sd1_const*(ooms + ooms_next);
        matrix[get_matrix_index(i, i-1)] = sd1_const*(ooms + ooms_prev);
        if(second_neighbor) {
    ooms_next2 = properties[1].ooms;
    matrix[get_matrix_index(i, 1)] = sd2_const*(ooms + ooms_next);
    matrix[get_matrix_index(i, i-2)] = sd2_const*(ooms + ooms_prev);
        }
      }else{
        /*We already have the prev ooms, but we need to fetch the next one.*/
        matrix[get_matrix_index(i, i+1)] = sd1_const*(ooms + ooms_next);
        matrix[get_matrix_index(i, i-1)] = sd1_const*(ooms + ooms_prev);
        if(second_neighbor) {
    if(i == L-2) {
      ooms_next2 = properties[0].ooms;
      matrix[get_matrix_index(i, 0)] = sd2_const*(ooms + ooms_next2);
      matrix[get_matrix_index(i, i-2)] = sd2_const*(ooms + ooms_prev2);
    } else if(i == 1) {
      ooms_prev2 = properties[L-1].ooms;
      matrix[get_matrix_index(i, i+2)] = sd2_const*(ooms + ooms_next2);
      matrix[get_matrix_index(i, L-1)] = sd2_const*(ooms + ooms_prev2);
    } else {
      matrix[get_matrix_index(i, i+2)] = sd2_const*(ooms + ooms_next2);
      matrix[get_matrix_index(i, i-2)] = sd2_const*(ooms + ooms_prev2);
    }
        }
      }
    }
  }
      }
      
      double* cartesian_effective_mass::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_real_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_EFFECTIVE_MASS_GERROR_DOMAIN, 1, "Got error from llw_real_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;
  int n_eigvecs;
  int* isuppz;
  info = llw_real_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_EFFECTIVE_MASS_GERROR_DOMAIN, 6, "in return from llw_real_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_EFFECTIVE_MASS_GERROR_DOMAIN, 1, "Bad argument to LAPACK call DSTEGR: %d\n", -info);
      }else if(info < 10) {
        *err = g_error_new(HAMILTONIAN_CARTESIAN_EFFECTIVE_MASS_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_EFFECTIVE_MASS_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_EFFECTIVE_MASS_GERROR_DOMAIN, 5, "Error in liblapackwrap llw_real_double_stegr_a: malloc failed (%d)\n", info/100);
      }else{
        *err = g_error_new(HAMILTONIAN_CARTESIAN_EFFECTIVE_MASS_GERROR_DOMAIN, 4, "Got error (info=%d) from DSTEGR, but don't know what it means.\n", info);
      }
    }
    fprintf(stderr, "Got error; freeing evals\n");
    free(*evals);
    free(eigvecs);
    *evals=NULL;
    eigvecs=NULL;
  }else{
    free(isuppz);
    if(Q != NULL) {
      fprintf(stderr, "Transforming eigenvectors by Q.\n");
      double* final_eigvecs = NULL;
      double one = 1.0;
      info = llw_real_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_EFFECTIVE_MASS_GERROR_DOMAIN, 4, "Got error (info=%d) from llw_real_double_gemmm_a\n", info);
        if(final_eigvecs != NULL) free(final_eigvecs);
        return eigvecs;
      }
      free(eigvecs);
      return final_eigvecs;
    }
  }
  return eigvecs;
      }
    }
  }
}