Work done for this code was performed in the National University of Singapore (NUS).

This README file includes instructions for the calculation of g-factors at the DFT and GW (COHSEX) levels according to Phys Rev Research 2, 033256 (2020).
(An example is provided in the folder example_MoS2_inputOnly/.)

Please note that an updated approach to incorporate dynamical effects into the calculation of the g-factors is presented in our recent paper
- npj Computational Materials, 7, 198 (2021). Here, dynamical effects are included in the GW g-factors. 
(An example is provided in the folder example_WSe2_renormalized_gfactor/)

We ask that if you use or refer to these codes for calculating g-factors (or Berry curvatures), that you please cite the following reference:

Phys Rev Research 2, 033256 (2020) (mean field, static COHSEX)
npj Computational Materials, 7, 198 (2021) (dynamical GW)


QuantumESPRESSO code is used in conjunction with this g-factor code. 
QuantumESPRESSO is licensed under GNU GENERAL PUBLIC LICENSE v2.0.

We also include code that is used in conjunction with the BerkeleyGW code. 
The BerkeleyGW code is licensed under a free, open source, and permissive 3-clause modified BSD license. 
BerkeleyGW, Copyright (c) 2011, The Regents of the University of California, through Lawrence Berkeley National Laboratory (subject to receipt of any required approvals from the U.S. Dept. of Energy). All rights reserved.

Author: Xuan Fengyuan (c2dxf@nus.edu.sg)
Contact: Su Ying Quek (phyqsy@nus.edu.sg)

======================================================================================
Instructions
======================================================================================

The calculation procedures are very simple and you may find that some parts of this algorithm, code or script are not effecient.
You can change these accordingly after you have understood the purpose of each step.

1) DFT level using QuantumESPRESSO-6.1

Our target is to calculate the orbital magnetic moment defined in Eq. (10) of Ref [Phys Rev Research 2, 033256 (2020)], where du/dk is evaluated as in Eq. (12):
<du/dk|x[Hk-Enk]|du/dk>
|du/dk> = (exp^{-i theta}|k+q> -|k>)/q
exp^{i theta} = <k|k+q>/|<k|k+q>|

So we need to calculate the DFT wave functions at K+qx, K+qy, K (K is the point at which g factor will be evaluated)

You can read "QuantumESPRESSO_codestructure.txt" to get a better understanding how QE works and why we modifiy c_bands.f90 
Do the following to get the g-factor:

   a, Set k list in the QE 'bands' calculation input: K+qx, K+qy, K. (qx = qy = 0.0001 or 0.00001)
      (be careful that in this step we need this parameter to be high    conv_thr = 1.0d-10)

      K_POINTS crystal
        3
      0.333353333  0.333323333  0.000000000   1.0
      0.333333333  0.333343333  0.000000000   1.0
      0.333333333  0.333333333  0.000000000   1.0

   b, Modify PW/src/c_bands.f90, which is the subroutine that controls the calculation of the k list in input. Set 'iiiv' to be the band index. Compile QE.
      remember to set qx in c_bands.f90 according to the input k list.
      the modification of c_bands.f90 is very simple. Just calcualte du/dkx and du/dky first and then evaluate <du/dk|x[Hk-Enk]|du/dk>
   c, Run pw.x < in > out.   In the "out", the g factor and Berry curvature are output for band "iiiv" at "K".


2) GW(COHSEX) level using QE-6.1 and BGW- (> = 1.2)

The strategy is straight forward. Just use the GW Hamiltonian and quasiparticle energies in <du/dk|x[Hk-Enk]|du/dk>.
<du/dk|x[Hk(GW)-Enk(GW)]|du/dk> = <du/dk|x[Hk(DFT)-Enk(GW)]|du/dk> + <du/dk|x[Sigma-Vxc]|du/dk> 
The second term is reduced to:
2Im <du/dkx |[Sigma-Vxc]| du/dky>

   a, Prepare WFN and WFNq. Perform the standard GW calculation to get Enk(GW). 
      Calculate <du/dk|x[Hk(DFT)-Enk(GW)]|du/dk> using the QE code (remember to replace Enk with GW values in c_bands.f90, or we can simply use the Berry curvature multiplied by Enk(GW))
      Next, we calculate <du/dk|x[Sigma-Vxc]|du/dk> using the part of BGW that computes the off-diagonal elements of the sigma matrix.
   b, Prepare a WFN_kkxky.h5. Be careful that the k-list order is different here : K, K+qx, K+qy

      K_POINTS crystal
        3
      0.333333333  0.333333333  0.000000000   1.0
      0.333353333  0.333323333  0.000000000   1.0
      0.333333333  0.333343333  0.000000000   1.0

      Remember: After pw2bgw, transform the BIN format to hdf5 format by: wfn2hdf.x BIN WFN.complex WFN_kkxky.h5.
   c, Prepare a WFN_after_C(V)BM.complex using python script 'writedudk.py'. 
      This script replaces the wave function coeffecient of band index 13 and 14 band at K by du/dkx and du/dky.
      The readin WFN_before.h5 is on the standard GW coarse grid (12x12 for TMD with subsample for example)
      13 14 is for TMD where VBM = 13. For other materials, change it accordingly.
      Set ib to be Valence and Conduction band in 'writedudk.py' to get WFN_after_C(V)BM.complex.
      (ib = 12 for valence and 13 for conduction, this is due to python array start from 0)
      You may run the python script as :

./writedudk.py WFN_kkxky.h5 WFN_before.h5 WFN_tempout.h5 WFN_after.h5 > CBM_write.out
rm WFN_tempout.h5
mv WFN_after.h5 WFN_after_CBM.h5
hdf2wfn.x BIN WFN_after_CBM.h5 WFN_after_CBM.complex % this has "band index 13 and 14 at K : du/dkx and du/dky"

Please note that there is a quantity called the unit of k ('unitk') in writedudk.py, which needs to be modified for each system. Here is also where you can decide on the units of length. 

% We need to run this again to get WFN_after_VBM.complex % this has "band index 12 and 13 at K : du/dkx and du/dky"

      Be careful that in the output ' CBM_write.out', there is a normalization factor at the end (last line). 
      The BGW code expects all wavefunctions to be normalized. So here, we have normalized du/dkx and du/dky and the product of the normalization factors is given in the output.
      We will need to take to multiply the results by the normalization factor in the final step of 2Im <du/dkx |[Sigma-Vxc]| du/dky>, since du/dkx and du/dky are not actually normalized to 1.

   d, Perform offdiag calculation using BGW for 2Im <du/dkx |[Sigma-Vxc]| du/dky>. 
      Link WFN_after_C(V)BM.complex to WFN_outer
      Remember to multiply by the normalization factor mentioned above.
      The GW g factor is then obtained by summing the two contributions in the equation:

<du/dk|x[Hk(GW)-Enk(GW)]|du/dk> = <du/dk|x[Hk(DFT)-Enk(GW)]|du/dk> + <du/dk|x[Sigma-Vxc]|du/dk>

NOTE: In all of the above, Rydberg units are used, except that the unit of length was in Angstroms. We need to divide by 0.53^2. Propose to change the code in future. 

3) GW g-factors

The intrinsic g-factor, using the GPP approximation to the energy depedence of the self-energy, is calculated in a similar way as the static COHSEX g-factor, except that the Hamiltonian matrix elements and energies used are obtained with dynamical GW.
Basically we need to calculate g^intrinsic = <du/dk|x[Hk(DFT)-Enk(GW)]|du/dk> + <du/dk|x[Sigma-Vxc]|du/dk>
The first part "<du/dk|x[Hk(DFT)-Enk(GW)]|du/dk>" can be calculated using our modified QE code "c_bands.f90" with the GW quasiparticle energy.
The second part "<du/dk|x[Sigma-Vxc]|du/dk>" can be calculated using the "outer corrections" flag in BerkeleyGW package where we put the du/dk in the WFN_outer.

The renormalized g-factor, which includes the dynamical effect from the QP self-energy, is defined in the Eq 8 of "npj Computational Materials, 7, 198 (2021)". 
This renormalized g-factor accounts for the effect of changes in the self-energy with energy. 
g^intrinsic / g^renormalized = 1 - d Sig(E) / d E