Skip to content

Commit

Permalink
Version 1.0.0 is now available, WOPP can now be correctly calculated.
Browse files Browse the repository at this point in the history
  • Loading branch information
@Chengcheng-Xiao
Chengcheng-Xiao committed Jan 8, 2019
1 parent 5eea9b7 commit f56ee21
Show file tree
Hide file tree
Showing 28 changed files with 866,977 additions and 639,641 deletions.
1 change: 1 addition & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
**/.DS_Store
relevant_document
Calculations
old_version
example/01_PbTiO3_example/CENT/wannier_data/wannier_AO
example/01_PbTiO3_example/CENT/wannier_data/wannier_MO
example/01_PbTiO3_example/FE/wannier_data/wannier_AO
Expand Down
17 changes: 7 additions & 10 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,8 +9,6 @@
```
# Wannier Orbital Overlap Population tools (WOOPs)

[THIS PACKAGE IS STILL UNDER DEVELOPMENT]

A post-processing tool written in python to get Wannier Orbital Overlap Population (WOOP), Wannier Orbital Position Population (WOPP)* from [Wannier90](https://github.com/wannier-developers/wannier90) package.


Expand Down Expand Up @@ -64,12 +62,12 @@ File need for WOOPs:

1. `wannier90_u_AO.mat` (Unitary matrix of atomic orbitals from valence band and conduction interpolation)
2. `wannier90_u_MO.mat` (Unitary matrix of molecular orbitals from only valence band interpolation)
3. `wannier90_r.dat` (Use master branch for compatibility)
3. `wannier90_tb.dat` (Use master branch for compatibility)
4. `input.woops`

Detailed preparation for Wannier90 generated files can be found in Wanneir90's user guid or in the `example` folder.
Detailed preparation for Wannier90 generated files can be found in Wanneir90's user guide, or in the `example` folder.

The format of input.woops follows:
The format of `input.woops` follows:
```
readmode = XX #[must be True or False]
cal = XX #[Things you want to do]
Expand All @@ -80,6 +78,7 @@ cell_param = XX XX XX #[cell_param]
cell_dim = XX #[cell_dim: 0D, 1D, 2D or 3D] NOTE: currently 1D->x, 2D->xy
cprec = XX #controls the PRINTING precision of C_matrix, lower means accurate. default=1e-4
bprec = XX #controls the PRINTING precision of WOOP_matrix, lower means accurate default=1e-4
pprec = XX #controls the PRINTING precision of WOPP_matrix, lower means accurate default=1e-2
```
*IMPORTANT*: Every tag must be written in lowercase, full length without abbreviation.

Expand All @@ -96,11 +95,11 @@ bprec = XX #controls the PRINTING precision of WOOP_matrix, low

When the input file is properly prepared:
```
./WOOPs.py
./WOOPs.py > woops.log &
```
or if you put it in your path:
```
WOOPs.py
WOOPs.py > woops.log &
```
This may take a few minutes, consider submitting it to the queue system.

Expand All @@ -112,13 +111,11 @@ These commands will produce a copy of data stored in hdf5 format and a text file
Go check the description in `example` folder.

## Notes
*WOPP is still under development, the problem now I'm facing is that the phase factor in Wannier90's r matrix are not "correct" (they are arbitrary). As a result, the ferroelectric polarization decomposition cannot be correctly calculated.
WOPP is still under Alpha testing phase, the problem now I'm facing is that the 'position matrix', or, Wannier90's r_matrix ( position in `hamiltonian.F90`) do not agree with Wannier Function Center (rave in `wannierise.F90`). As a result, the ferroelectric polarization decomposition cannot be correctly calculated. To avoid this, I am now using <0n|r|Rm> written in the `wanneir90_tb.dat` file.

## TODO_list
1. Complete a simplified description of WOOP and its capability.

2. Fix phase factor error in WOPP interface.

3. Increase efficiency.

4. Add progress bar.
Expand Down
66 changes: 40 additions & 26 deletions WOOP_v0.1.3.py → WOOPs.py
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
#!python
#!/usr/bin/env python
from __future__ import print_function
import numpy as np
import linecache
import warnings
with warnings.catch_warnings():
warnings.simplefilter("ignore")
Expand All @@ -18,6 +19,7 @@ def read_input():
cal = "get_WOOP"
cprec = 1e-4
bprec = 1e-4
pprec = 1e-2
for line in data:
key, value = line.split("=")
dataset[key.strip()] = value.strip()
Expand All @@ -33,7 +35,9 @@ def read_input():
cprec = float(dataset["cprec"])
if "bprec" in dataset.keys():
bprec = float(dataset["bprec"])
return number_MO, number_AO, number_kpts, cell_dim, dim, cal, readmode, cprec, bprec
if "pprec" in dataset.keys():
pprec = float(dataset["pprec"])
return number_MO, number_AO, number_kpts, cell_dim, dim, cal, readmode, cprec, bprec, pprec

def get_u_matrix(file_name,dimension_fix,dimension,num_kpoints):

Expand Down Expand Up @@ -149,21 +153,26 @@ def get_r_matrix(file_name, num_wann):
''' a subroutine to get <ul|r|um> matrix from seedname_r.mat '''
''' num_wann -> number of MO'''

for num_f in range(len(open(file_name,'r').readlines())):
if linecache.getline(file_name,num_f) == "\n":
start_line = num_f
#print(num_f)
break
file = open(file_name, "r")
content = [x.rstrip("\n") for x in file]
nrpts = int(content[2])
data = [x.split()[:22] for x in content[3:]]
nrpts = int(content[5])
mat_j = np.zeros((nrpts,num_wann,num_wann,3),dtype=complex)
R_place = np.zeros((nrpts,3))
data_start = start_line-1+((num_wann*num_wann+2)*int(nrpts))-1
for nrpt in range(nrpts):
R_place[nrpt][0] = data[nrpt*(num_wann**2)][0]
R_place[nrpt][1] = data[nrpt*(num_wann**2)][1]
R_place[nrpt][2] = data[nrpt*(num_wann**2)][2]
for i in range(num_wann):
for j in range(num_wann):
mat_j[nrpt][i][j][0]=(float(data[nrpt*(num_wann**2)+i+num_wann*j][5])+1j*float(data[nrpt*(num_wann**2)+i+num_wann*j][6])) #x
mat_j[nrpt][i][j][1]=(float(data[nrpt*(num_wann**2)+i+num_wann*j][7])+1j*float(data[nrpt*(num_wann**2)+i+num_wann*j][8])) #y
mat_j[nrpt][i][j][2]=(float(data[nrpt*(num_wann**2)+i+num_wann*j][9])+1j*float(data[nrpt*(num_wann**2)+i+num_wann*j][10])) #z
R_place[nrpt][0] = content[data_start+nrpt*(num_wann*num_wann+2)+2].split()[:6][0]
R_place[nrpt][1] = content[data_start+nrpt*(num_wann*num_wann+2)+2].split()[:6][1]
R_place[nrpt][2] = content[data_start+nrpt*(num_wann*num_wann+2)+2].split()[:6][2]
for ao_home in range(num_wann):
for ao_other in range(num_wann):
mat_j[nrpt][ao_home][ao_other][0]=float(content[data_start+nrpt*(num_wann*num_wann+2)+3+ao_home+ao_other*num_wann].split()[:16][2])+1j*float(content[data_start+nrpt*(num_wann*num_wann+2)+3+ao_home+ao_other*num_wann].split()[:16][3])
mat_j[nrpt][ao_home][ao_other][1]=float(content[data_start+nrpt*(num_wann*num_wann+2)+3+ao_home+ao_other*num_wann].split()[:16][4])+1j*float(content[data_start+nrpt*(num_wann*num_wann+2)+3+ao_home+ao_other*num_wann].split()[:16][5])
mat_j[nrpt][ao_home][ao_other][2]=float(content[data_start+nrpt*(num_wann*num_wann+2)+3+ao_home+ao_other*num_wann].split()[:16][6])+1j*float(content[data_start+nrpt*(num_wann*num_wann+2)+3+ao_home+ao_other*num_wann].split()[:16][7])
# mat_j is in mat_j[nrpt][l][m] = <ul_origin|r|um_nrpt>
return mat_j, R_place, nrpts

Expand All @@ -173,6 +182,7 @@ def get_WOPP(rmat, kpts,C_ij, AO_wann, MO_wann, use_nrpt, R_place, cell_dim):

#initialize data array
D_imn=np.zeros((3,MO_wann,AO_wann,AO_wann,len(use_nrpt),len(use_nrpt)),dtype=complex)
rmat_tmp = np.zeros(3,dtype=complex)
#Calculating WOOP
for dir in range(3):
for MO in range(MO_wann):
Expand All @@ -181,15 +191,19 @@ def get_WOPP(rmat, kpts,C_ij, AO_wann, MO_wann, use_nrpt, R_place, cell_dim):
for nrpt_2 in range(len(use_nrpt)):
for n in range(AO_wann):
#if nrpt_2 is not home cell and nrpt_2 == nrpt_1 then both AO are in cell other than home, translation neede.
rmat_tmp = rmat[R_place.tolist().index((-R_place[int(use_nrpt[nrpt_1])]+R_place[int(use_nrpt[nrpt_2])]).tolist())][m][n][dir]# + R_place[int(use_nrpt[nrpt_1])][dir]*cell_dim[dir]
if nrpt_1 != 0 and nrpt_1 == nrpt_2 and m == n and np.real(C_ij[nrpt_1][MO][m]) > 0.1:
rmat_tmp += R_place[int(use_nrpt[nrpt_1])][dir]*cell_dim[dir]
D_imn[dir][MO][m][n][nrpt_1][nrpt_2] = np.dot(np.dot(np.conj(C_ij[nrpt_1][MO][m]),rmat_tmp),C_ij[nrpt_2][MO][n])
rmat_tmp_0 = rmat[R_place.tolist().index((-R_place[int(use_nrpt[nrpt_1])]+R_place[int(use_nrpt[nrpt_2])]).tolist())][m][n][:]#[dir] + R_place[int(use_nrpt[nrpt_1])][dir]*cell_dim[dir]
if nrpt_1 != 0 and nrpt_1 == nrpt_2 and m == n:# and np.real(C_ij[nrpt_1][MO][m]) > 0.1:
rmat_tmp[0] = rmat_tmp_0[0] + R_place[int(use_nrpt[nrpt_1])][0]*cell_dim[0]
rmat_tmp[1] = rmat_tmp_0[1] + R_place[int(use_nrpt[nrpt_1])][1]*cell_dim[1]
rmat_tmp[2] = rmat_tmp_0[2] + R_place[int(use_nrpt[nrpt_1])][2]*cell_dim[2]
else:
rmat_tmp[0] = rmat_tmp_0[0]
rmat_tmp[1] = rmat_tmp_0[1]
rmat_tmp[2] = rmat_tmp_0[2]
D_imn[dir][MO][m][n][nrpt_1][nrpt_2] = np.dot(np.dot(np.conj(C_ij[nrpt_1][MO][m]),rmat_tmp[dir]),C_ij[nrpt_2][MO][n])
return D_imn


#-----------------------------------------------------------------------------------

#===================================================================================#
# #
# Here comes the main program #
Expand All @@ -201,14 +215,14 @@ def get_WOPP(rmat, kpts,C_ij, AO_wann, MO_wann, use_nrpt, R_place, cell_dim):
"| |/\\| || | | || | | | __/ __|\n"
"\\ /\\ /\\ \\_/ /\\ \\_/ / | \\__ \\\n"
" \\/ \\/ \\___/ \\___/\\_| |___/\n"
" version 0.1.1\n")
" version 1.0.0\n")


# User input value start
number_MO, number_AO, number_kpts, cell_dim, dim, cal, readmo, cprec, bprec=read_input()
number_MO, number_AO, number_kpts, cell_dim, dim, cal, readmo, cprec, bprec, pprec=read_input()
MO_filename="wannier90_u_MO.mat" #sys.argv[4]
AO_filename="wannier90_u_AO.mat" #sys.argv[5]
AO_r_filename="wannier90_r.dat" #sys.argv[6]
AO_r_filename="wannier90_tb.dat" #sys.argv[6]
# User input value end
cal_orb, cal_comple, cal_c, cal_charg, cal_woop, cal_wopp, readmode = False, False, False, False, False, False, False
if cal == "WOPP":
Expand Down Expand Up @@ -531,7 +545,7 @@ def get_WOPP(rmat, kpts,C_ij, AO_wann, MO_wann, use_nrpt, R_place, cell_dim):
print("# #",file=f)
print("*######################################################*",file=f)

for dir in range(2):
for dir in range(dim):
total_moment = 0
for i in range(number_MO):
total_part = 0
Expand All @@ -541,9 +555,9 @@ def get_WOPP(rmat, kpts,C_ij, AO_wann, MO_wann, use_nrpt, R_place, cell_dim):
for n in range(number_AO):
total_part += np.real(WOPP[dir][i][m][n][nrpt_1][nrpt_2])
total_moment += np.real(WOPP[dir][i][m][n][nrpt_1][nrpt_2])
if np.abs(2*np.real(WOPP[dir][i][m][n][nrpt_1][nrpt_2])) > 1e-2:
print("Direction= {0:2d} MO_number_i= {1:2d} AO_number_m= {2:2d} [{3:2.0f},{4:2.0f},{5:2.0f}] AO_number_n= {6:2d} [{7:2.0f},{8:2.0f},{9:2.0f}] WOPP= {10:10f}".format(dir,i,m,R_place[int(use_nrpt[nrpt_1])][0],R_place[int(use_nrpt[nrpt_1])][1],R_place[int(use_nrpt[nrpt_1])][2],n,R_place[int(use_nrpt[nrpt_2])][0],R_place[int(use_nrpt[nrpt_2])][1],R_place[int(use_nrpt[nrpt_2])][2],2*np.real(WOPP[dir][i][m][n][nrpt_1][nrpt_2])))
print("number of MO {0:2d} total_moment= {1:10f} ".format(i,total_part*-2),file=f)
print("Direction {0:2d} total_moment= {1:10f} ".format(dir,total_moment*-2),file=f)
if np.abs(2*np.real(WOPP[dir][i][m][n][nrpt_1][nrpt_2])) > pprec:
print("Direction= {0:2d} MO_number_i= {1:2d} AO_number_m= {2:2d} [{3:2.0f},{4:2.0f},{5:2.0f}] AO_number_n= {6:2d} [{7:2.0f},{8:2.0f},{9:2.0f}] WOPP= {10:10f}".format(dir,i,m,R_place[int(use_nrpt[nrpt_1])][0],R_place[int(use_nrpt[nrpt_1])][1],R_place[int(use_nrpt[nrpt_1])][2],n,R_place[int(use_nrpt[nrpt_2])][0],R_place[int(use_nrpt[nrpt_2])][1],R_place[int(use_nrpt[nrpt_2])][2],-2.*np.real(WOPP[dir][i][m][n][nrpt_1][nrpt_2])),file=f)
print("number of MO {0:2d} total_moment= {1:10f} ".format(i,total_part*-2.),file=f)
print("Direction {0:2d} total_moment= {1:10f} ".format(dir,total_moment*-2.),file=f)

print("All calculations done, see you next time :)")
Loading

0 comments on commit f56ee21

Please sign in to comment.