This thread is locked.Only browsing is available.
Top Page > Browsing
Diagonal terms in overlap matrix not equalled to 1
Date: 2015/06/04 17:09
Name: KZ

Hi, OpenMX users and developers,
I am using openmx to get the hamiltonian and overlap for further uses. And I am trying to use H,S to reproduce the E-k relations of monolayer MoS2.
But as I see the overlap matrix, the diagonal terms of overlap matrix I got are not 1. Could you help explain why it is? Or I just did it wrong in some steps?


Page: [1]

Re: Diagonal terms in overlap matrix not equalled to 1 ( No.1 )
Date: 2015/06/04 17:49
Name: Artem Pulkin


Could you please explain, why do you expect the diagonal terms to be unity?

Re: Diagonal terms in overlap matrix not equalled to 1 ( No.2 )
Date: 2015/06/05 04:42
Name: KZ

Hi, Artem,
I believe only if the basis sets are normolized, the overlap matrix between themselves should be unity, right?
So I am checking the basis sets OpenMX is using.

Moreover, because I also got Hamiltonian, Overlap matrix that are not strictly hermitian only because some very tiny elements in the order of 10^-15. So the bands I reproduced are still strange.
I am wondering which step is wrong for my calculations.

Re: Diagonal terms in overlap matrix not equalled to 1 ( No.3 )
Date: 2015/06/05 10:09
Name: KZ

I checked the bands and hamiltonian and overlap matrix of graphene as well. The diagonal terms are also not unity.
So could anyone help point out: am I just thinking about this in a wrong way?
Re: Diagonal terms in overlap matrix not equalled to 1 ( No.4 )
Date: 2015/06/05 18:44
Name: Artem Pulkin

Actually I do not remember if the basis set is normalized in OpenMX. I just want to emphasize that this is not important anyway. 10^-15 is below machine epsilon. To reproduce the bands you need to solve

[ H_0 + (H_x exp(ikx) + H_y exp(iky) + H_z exp(ikz) + H.c.) ] psi = E(k) *
[ S_0 + (S_x exp(ikx) + S_y exp(iky) + S_z exp(ikz) + H.c.) ] psi

You do not really care about the basis set in the above equation. It worked for me for a particular case of MoS2.


Re: Diagonal terms in overlap matrix not equalled to 1 ( No.6 )
Date: 2015/06/06 10:36
Name: KZ

I see what you were saying here. I think I am doing that part as you did.
But the bands are not correct for monolayer graphene case yet. So I am checking the elements. It seems the coupling terms between C_pz C_pz are about 0.08. Which are quite small compared to 2.7 eV which it should be.
So am I missing some units here? Or which part might be wrong?


Re: Diagonal terms in overlap matrix not equalled to 1 ( No.7 )
Date: 2015/06/08 13:44
Name: T. Ozaki


I am wondering that you got diagonal overlap integrals such as 0.9999964, which is
a bit different from unity. This is due to difference of integration scheme between
generation and use of those basis functions. OpenMX itself never assumes that the
basis functions are exactly normalized, and just solves the generalized eigenvalue
problem with overlap matrix elements which are evaluated in OpenMX.

As for the reproduction of band structure, if you properly solve the generalized eigenvalue
problem, you should obtain the same band structure as what OpenMX calculates.
If not, this is an indication that you did not solve the problem properly.


Re: Diagonal terms in overlap matrix not equalled to 1 ( No.8 )
Date: 2015/06/10 21:46
Name: Artem Pulkin

Dear KZ,

I hope that you solved your problem. Otherwise I can send you some example code in C and python that reproduces the band structure. When I was implementing it I found Hamiltonian_Band.c and Hamiltonian_Band_NC.c particularly helpful. Also try doing simpler systems without spin-orbit coupling (it is way easier for a start).

Unfortunately, without looking into the code there is no way to guess, for example, that the imaginary part of the Hamiltonian is spread over 4 different structures Hks[3] and iHks[0,1,2]. I understood it as a sort of 'compressing' Hamiltonian (still, it could be compressed further). And I have to say that I do not like the way it is done.


Re: Diagonal terms in overlap matrix not equalled to 1 ( No.9 )
Date: 2015/06/11 03:48
Name: KZ

Thanks for your response. I am still looking into the codes. I am trying to find a good place to add some lines to output the hamiltonian. If you can send me some example codes, that will be great.
My email address is:;
Thanks for your help. I appreciate it.


Re: Diagonal terms in overlap matrix not equalled to 1 ( No.11 )
Date: 2015/06/12 06:21
Name: KZ

Hi, Artem,
I am now testing the hamiltonian on monolayer graphene system without SOC. And for some reason, it still have problems. I will check more.
And I checked with the imaginary part of Hamiltonian iHks in read_scfout.c. I believe for this simple case, it will not come into the problem at all, right?
And I think all the data I need for reproducing the bands are from using read_scfout.c file to read the binary file, right?


Re: Diagonal terms in overlap matrix not equalled to 1 ( No.12 )
Date: 2015/06/12 19:15
Name: Artem Pulkin

Yep, as far as I understand iHks are (imaginary) spin-orbit and +U terms. And yes, it looks like read_scfout.c readf .scfout file with all the data.

I will try to comment some parts of the code taken from Hamiltonian_Band_NC.c and analysis_example.c, probably it will help you. Please note that this is not a working code.

// The piece of code constructs an H(k) spin-orbit Hamiltonian having a size of
// 2*NUM = 2*sum(Total_NumOrbs) where 2 goes for spin index.

// Initialize arrays. The full Hamiltonian is presented as 3 blocks: H11, H22, H12
// where '1' is a spin-up index, '2' is a spin-down. As you know well, H21 is just
// a Hermitian conjugate of H12 so it is skipped in this part. 'r' index goes for
// real part, 'i' for imaginary.

for (i=1; i<=NUM; i++){
for (j=1; j<=NUM; j++){
H11r[i][j] = 0.0;
H11i[i][j] = 0.0;
H22r[i][j] = 0.0;
H22i[i][j] = 0.0;
H12r[i][j] = 0.0;
H12i[i][j] = 0.0;

// Now iterate over the DFT matrix elements as it is explained in analysis_example.c
// First we iterate over non-spin-orbit part 'Hks', then over 'iHks'.

// The load cycle (note the iteration over 'spin' removed'), taken from analysis_example.c:

for (ct_AN=1; ct_AN<=atomnum; ct_AN++){
TNO1 = Total_NumOrbs[ct_AN];
for (h_AN=0; h_AN<=FNAN[ct_AN]; h_AN++){
Gh_AN = natn[ct_AN][h_AN];
Rn = ncn[ct_AN][h_AN];
TNO2 = Total_NumOrbs[Gh_AN];

// This part is taken from Hamiltonian_Band_NC.c. It constructs a phase
// factor exp(ik*r) for the H(k) Hamiltonian. l1,l2,l3 are integers defining
// 'r' and 'k' goes from -0.5 to 0.5. In other words, all these quantities
// are dimensionless. The factor exp(ik*r) = co + i*si, i - imaginary 1.

l1 = atv_ijk[Rn][1];
l2 = atv_ijk[Rn][2];
l3 = atv_ijk[Rn][3];
kRn = k1*(double)l1 + k2*(double)l2 + k3*(double)l3;

si = sin(2.0*PI*kRn);
co = cos(2.0*PI*kRn);

for (i=0; i<TNO1; i++){
for (j=0; j<TNO2; j++){

// Now the actual load. This part is also taken (corrected) from
// Hamiltonian_Band_NC.c. It writes corresponding terms of our H(k)
// Hamiltonian. 'Anum' and 'Bnum' are corresponding index offsets
// in H(k) for atoms 'ct_AN' and 'Gh_AN'.

// The first index (spin) of Hks is the following:
// 0 corresponds to H11r
// 1 corresponds to H22r
// 2 corresponds to H12r
// 3 corresponds to H12i
// I.e. we have 3 real parts here and one imaginary.
// While the first index of iHks is the following:
// 0 corresponds to H11i
// 1 corresponds to H22i
// 2 corresponds to H12i
// I.e. all 3 parts are imaginary.
// In complex notation, this piece of code just calculates:
// (H11r + i*H11i) += exp(ikr) * (Hks[0] + i*iHks[0])
// (H22r + i*H22i) += exp(ikr) * (Hks[1] + i*iHks[1])
// (H12r + i*H12i) += exp(ikr) * (Hks[2] + i*Hks[3] + i*iHks[2])

H11r[Anum+i][Bnum+j] += co*Hks[0][MA_AN][LB_AN][i][j] - si*iHks[0][MA_AN][LB_AN][i][j];
H11i[Anum+i][Bnum+j] += si*Hks[0][MA_AN][LB_AN][i][j] + co*iHks[0][MA_AN][LB_AN][i][j];
H22r[Anum+i][Bnum+j] += co*Hks[1][MA_AN][LB_AN][i][j] - si*iHks[1][MA_AN][LB_AN][i][j];
H22i[Anum+i][Bnum+j] += si*Hks[1][MA_AN][LB_AN][i][j] + co*iHks[1][MA_AN][LB_AN][i][j];
H12r[Anum+i][Bnum+j] += co*Hks[2][MA_AN][LB_AN][i][j] - si*(Hks[3][MA_AN][LB_AN][i][j]
+ iHks[2][MA_AN][LB_AN][i][j]);
H12i[Anum+i][Bnum+j] += si*Hks[2][MA_AN][LB_AN][i][j] + co*(Hks[3][MA_AN][LB_AN][i][j]
+ iHks[2][MA_AN][LB_AN][i][j]);

// Finally you put everything into a 2D array:
// NUM = sum(Total_NumOrbs) is just a half of the size of the full Hk.

for (i=1; i<=NUM; i++){
for (j=1; j<=NUM; j++){
H[i ][j ].r = H11r[i][j];
H[i ][j ].i = H11i[i][j];
H[i+NUM][j+NUM].r = H22r[i][j];
H[i+NUM][j+NUM].i = H22i[i][j];
H[i ][j+NUM].r = H12r[i][j];
H[i ][j+NUM].i = H12i[i][j];
// H21 = hermitian_conjugate(H12).
H[j+NUM][i ].r = H[i][j+NUM].r;
H[j+NUM][i ].i = -H[i][j+NUM].i;

// H(k) is ready. Sorry for the typos.

// Regads,

// Artem

Page: [1]