.. _program_listing_file_src_qcm_Chern.cpp:

Program Listing for File Chern.cpp
==================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_qcm_Chern.cpp>` (``src_qcm/Chern.cpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   
   #include <iostream>
   #include <fstream>
   
   #define MAX_K_SIDE 5000
   
   #ifdef _OPENMP
   #include <omp.h>
   #endif
   
   #include "lattice_model_instance.hpp"
   #include "integrate.hpp"
   #include "parser.hpp"
   
   bool recursive=false;
   const double threshold = 0.1*M_PI;
   vector3D<double> increment(double d, int no, int dir);
   
   template <typename T>
   void gauge_field(matrix<T> &A, matrix<T> &B, vector<T> &u){
       for(int i=0; i<A.c; ++i){
           u[i] = 0.0;
           for(int j=0; j<A.r; ++j){
               u[i] += conjugate(A(j,i))*B(j,i);
           }
       }
   }
   
   void lattice_model_instance::Green_eigensystem(Green_function &G, const vector3D<double> &k0, vector<double> &e, matrix<Complex> &U, int opt)
   {
     vector3D<double> k = k0;
     if(opt&1) k = model->superdual.to(model->dual.from(k));
     else if(opt&2){
       // do nothing : leave k as is
     }
     else k = model->superdual.to(model->physdual.from(k));
     
     Green_function_k K(G,k);
     if(opt&2){
       set_Gcpt(K);
       K.Gcpt.eigensystem(e, U);
     }
     else{
       periodized_Green_function(K);
       K.g.eigensystem(e, U);
     }
   }
   
   vector3D<double> increment(double d, int no, int dir)
   {
     switch(dir){
       case  3:
         switch(no){
           case 0: return { d, 0.0, 0.0};
           case 1: return { 0.0, d, 0.0};
           case 2: return {-d, 0.0, 0.0};
         }
       case  -3:
         switch(no){
           case 0: return { d, 0.0, 0.0};
           case 1: return { 0.0,-d, 0.0};
           case 2: return {-d, 0.0, 0.0};
         }
       case  2:
         switch(no){
           case 0: return {-d, 0.0, 0.0};
           case 1: return { 0.0, 0.0, d};
           case 2: return { d, 0.0, 0.0};
         }
       case  -2:
         switch(no){
           case 0: return {-d, 0.0, 0.0};
           case 1: return { 0.0, 0.0, -d};
           case 2: return { d, 0.0, 0.0};
         }
       case  1:
         switch(no){
           case 0: return {0.0, d, 0.0};
           case 1: return {0.0, 0.0, d};
           case 2: return {0.0, -d, 0.0};
         }
       case  -1:
         switch(no){
           case 0: return {0.0, d, 0.0};
           case 1: return {0.0, 0.0,-d};
           case 2: return {0.0,-d, 0.0};
         }
     }
     return {0.0, 0.0, 0.0};
   }
   
   
   
   
   double lattice_model_instance::Berry_flux(const vector<vector3D<double>> &k, int orb, bool spin_down)
   {
     check_signals();
     double eta = global_double("eta");
     Green_function G = cluster_Green_function(Complex(0, eta), false, spin_down);
   
     size_t ng = model->dim_reduced_GF;
     if(orb > ng) qcm_throw("the orbital number specified in Berry curvature computations is beyond the number of orbitals in the lattice model");
       
     matrix<Complex> U(ng), U0(ng), U1(ng);
       vector<Complex> u(ng);
       vector<double> e(ng), e0(ng), e1(ng);
   
     // loop over the wavevectors of the path
     Green_eigensystem(G, k[0], e1, U1, 0); for(auto& x : e0) x = -1.0/x;
     U0 = U1;
     e0 = e1;
     Complex z(1.0);
     for(int i=1; i<k.size(); i++){
       Green_eigensystem(G, k[i], e, U, 0); for(auto& x : e) x = -1.0/x;
       gauge_field(U0, U, u);
       if (orb==0){
         for(int b=0; b<ng; b++) if(e[b] < 0.0) z *= u[b];
       }
       else z *= u[orb-1];
       U0 = U;
       e0 = e;
     }
     gauge_field(U0, U1, u);
     if (orb==0){
         for(int b=0; b<ng; b++) if(e0[b]+e[b] < 0.0) z *= u[b];
     }
     else z *= u[orb-1];
     double flux = -arg(z)/(2*M_PI);
   
     if(model->mixing == HS_mixing::normal) flux *= 2.0;
     if(model->mixing == HS_mixing::up_down and spin_down == false) flux += Berry_flux(k, true, 0);
     return flux;
   }
   
   
   
   
   double lattice_model_instance::Berry_plaquette(Green_function &G, const vector3D<double> &k1, const double deltax, const double deltay, const int opt, int dir, int orb)
   {
     check_signals();
   
     size_t ng;
     if(opt&2) ng = model->dim_GF;
     else ng = model->dim_reduced_GF;
     if(orb > ng) qcm_throw("the orbital label specified in Berry curvature computations is beyond the number of orbitals in the lattice model");
   
       matrix<Complex> U1(ng), U2(ng), U3(ng), U4(ng);
       vector<Complex> u1(ng), u2(ng), u3(ng), u4(ng);
       vector<double> e1(ng), e2(ng), e3(ng), e4(ng);
     vector3D<double> k(k1);
   
     Green_eigensystem(G, k, e1, U1, opt);
     for(auto& x : e1) x = -1.0/x;
   
     k += increment(deltax, 0, dir);
     Green_eigensystem(G, k, e2, U2, opt);
     for(auto& x : e2) x = -1.0/x;
     
     k += increment(deltay, 1, dir);
     Green_eigensystem(G, k, e3, U3, opt);
     for(auto& x : e3) x = -1.0/x;
   
     k += increment(deltax, 2, dir);
     Green_eigensystem(G, k, e4, U4, opt);
     for(auto& x : e4) x = -1.0/x;
   
       gauge_field(U1, U2, u1);
       gauge_field(U2, U3, u2);
       gauge_field(U3, U4, u3);
       gauge_field(U4, U1, u4);
       Complex z(1.0);
     if (orb==0){
       for(int b=0; b<ng; b++){
         if(e1[b]+e2[b]+e3[b]+e4[b] < 0.0) z *= u1[b]*u2[b]*u3[b]*u4[b];
       }
     }
     else{
       int b = orb-1;
       z *= u1[b]*u2[b]*u3[b]*u4[b];
     }
     double phase = arg(z);
     if(abs(phase) > threshold && recursive){
       phase = 0.0;
       double dx = deltax/2;
       double dy = deltay/2;
       k = k1;
       phase += Berry_plaquette(G, k, dx, dy, opt, dir, orb);
       k += increment(dx, 0, dir);
       phase += Berry_plaquette(G, k, dx, dy, opt, dir, orb);
       k += increment(dy, 1, dir);
       phase += Berry_plaquette(G, k, dx, dy, opt, dir, orb);
       k += increment(dx, 2, dir);
       phase += Berry_plaquette(G, k, dx, dy, opt, dir, orb);
     }
   
   
       return phase;
   }
   
   
   
   
   vector<double> lattice_model_instance::Berry_curvature(vector3D<double>& k1, vector3D<double>& k2, int nk, int orb, bool rec, int dir)
   {
     recursive = rec;
   
     if(model->spatial_dimension < 2) qcm_throw("'Berry_curvature' can only be applied to a 2D or 3D Brillouin zone!");
     if(nk > MAX_K_SIDE) qcm_throw("The wavevector grid has too many points. nk should be <= "+to_string<int>(MAX_K_SIDE));
     if(nk < 10) qcm_throw("The wavevector grid has too few points. nk should be >= 10");
     
     int opt=0;
     if(global_bool("dual_basis")) opt += 1;
     if(global_char("periodization") == 'N') opt += 2;
   
       
     vector<double> B(nk*nk, 0.0);
     vector3D<double> delta1;
     vector3D<double> delta2;
     double d1, d2;
     if(dir==1){
       if(fabs(k1.x-k2.x)>1e-10) qcm_throw("in Berry curvature, k1 and k2 must have the same 3rd component");
       d1 = (k2.y-k1.y)/nk;
       d2 = (k2.z-k1.z)/nk;
       delta1.y = d1;
       delta2.z = d2;
     }
     else if(dir==2){
       if(fabs(k1.y-k2.y)>1e-10) qcm_throw("in Berry curvature, k1 and k2 must have the same 3rd component");
       d1 = (k2.x-k1.x)/nk;
       d2 = (k2.z-k1.z)/nk;
       delta1.x = d1;
       delta2.z = d2;
     }
     else{
       if(fabs(k1.z-k2.z)>1e-10) qcm_throw("in Berry curvature, k1 and k2 must have the same 3rd component");
       d1 = (k2.x-k1.x)/nk;
       d2 = (k2.y-k1.y)/nk;
       delta1.x = d1;
       delta2.y = d2;
     }
   
     double eta = global_double("eta");
     Green_function G = cluster_Green_function(Complex(0, eta), false, false);
           
       // loop over the plaquettes
     // #pragma omp parallel for
       for(int j=0; j<nk; j++){
           for(int i=0; i<nk; i++){
               B[i+nk*j] = Berry_plaquette(G, k1 + i*delta1 + j*delta2, d1, d2, opt, dir, orb);
           }
       }
       
       if(model->mixing == HS_mixing::up_down){
           G = cluster_Green_function(Complex(0, 1e-8), false, true);
       for(int j=0; j<nk; j++){
         for(int i=0; i<nk; i++){
           B[i+nk*j] = Berry_plaquette(G, k1 + i*delta1 + j*delta2, d1, d2, opt, dir, orb);
         }
       }
     }
     else if(model->mixing == HS_mixing::normal) B *= 2.0;
     // else if(model->mixing == HS_mixing::full) B *= 0.5;
     B *= -1.0/(2*M_PI*d1*d2);
    
     return B;
   }
   
   
   
   
   
   double lattice_model_instance::monopole_part(vector3D<double>& k, double a, int nk, int orb, bool rec, int dir, bool spin_down)
   {
     recursive = rec;
   
     int opt=0;
   
     double charge = 0.0;
     double d = 2*a/nk;
     
     double eta = global_double("eta");
     // double eta = 1e-6;
     Green_function G = cluster_Green_function(Complex(0, eta), false, spin_down);
           
       // loop over the plaquettes
     // #pragma omp parallel for
       for(int j=0; j<nk; j++){
           for(int i=0; i<nk; i++){
         vector3D<double>kb = k;
         switch(dir){
           case 1 :
             kb.x += a;
             kb.y += i*d - a;
             kb.z += j*d - a;
             break;
           case -1 :
             kb.x -= a;
             kb.y -= i*d - a;
             kb.z -= j*d - a;
             break;
           case 2 :
             kb.y += a;
             kb.z += i*d - a;
             kb.x += j*d - a;
             break;
           case -2 :
             kb.y -= a;
             kb.z -= i*d - a;
             kb.x -= j*d - a;
             break;
           case 3 :
             kb.z += a;
             kb.x += i*d - a;
             kb.y += j*d - a;
             break;
           case -3 :
             kb.z -= a;
             kb.x -= i*d - a;
             kb.y -= j*d - a;
             break;
         }
               charge += Berry_plaquette(G, kb, d, d, opt, dir, orb);
           }
       } 
     return charge;
   }
   
   
   
   double lattice_model_instance::monopole(vector3D<double>& k, double a, int nk, int orb, bool rec)
   {
     recursive = rec;
   
     if(model->spatial_dimension != 3) qcm_throw("'monopole' can only be applied to a 3D model!");
     if(nk > 1000) qcm_throw("The wavevector grid has too many points. nk should be <= 1000");
     if(nk < 2) qcm_throw("The wavevector grid has too few points. nk should be >= 2");
     
       double charge = 0.0;
   
     charge += monopole_part(k, a, nk, orb, rec, 1, false);
     charge += monopole_part(k, a, nk, orb, rec,-1, false);
     charge += monopole_part(k, a, nk, orb, rec, 2, false);
     charge += monopole_part(k, a, nk, orb, rec,-2, false);
     charge += monopole_part(k, a, nk, orb, rec, 3, false);
     charge += monopole_part(k, a, nk, orb, rec,-3, false);
    
     if(model->mixing == HS_mixing::up_down){
       charge += monopole_part(k, a, nk, orb, rec, 1, true);
       charge += monopole_part(k, a, nk, orb, rec,-1, true);
       charge += monopole_part(k, a, nk, orb, rec, 2, true);
       charge += monopole_part(k, a, nk, orb, rec,-2, true);
       charge += monopole_part(k, a, nk, orb, rec, 3, true);
       charge += monopole_part(k, a, nk, orb, rec,-3, true);
     }
     else if(model->mixing == HS_mixing::normal) charge *= 2.0;
     charge *= -1.0/(2*M_PI);
     return charge;
   }