Skip to content

Commit 69b2c8e

Browse files
avmaharajavmaharaj
committedApr 1, 2016
Added Fortran implementation
1 parent c09704b commit 69b2c8e

9 files changed

+591
-0
lines changed
 

‎.DS_Store

4 KB
Binary file not shown.

‎fortran_implementation/.DS_Store

6 KB
Binary file not shown.

‎fortran_implementation/README.md

+19
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,19 @@
1+
# QO_fort
2+
### A newer Fortran implementation of the recursive calculation of Quantum Oscillations
3+
4+
This is also written in an object oriented manner, but now in the much faster (compiled)
5+
Fortran language. Usage is virtually identical to the python implementation. (see main_recursive.f90)
6+
7+
### Description copied from python implementation
8+
An object oriented version of the recursive Green's function approach to calculation Quantum Oscillations. Here we include a 6 band Hamiltonian representing the three t2g orbitals ($d_{xz}$, $d_{yx}$ and $d_{xy}$), and including spin orbit coupling.
9+
10+
Details of the recursive algorithm are given here:
11+
http://diku.dk/forskning/Publikationer/tekniske_rapporter/2010/Inversion-of-Block-Tridiagonal-Matrices.pdf
12+
13+
### General format of the code:
14+
Initialize an object which grabs input parameters from the commnad line
15+
Initialize a Hamiltonian object which has appropriate (tight binding) parameters.
16+
Pass this Hamiltonian object to the Solver object, which implements the recursive algorithm.
17+
18+
### Example usage from command line:
19+
./qo length=16394 mu=-0.5 nsteps=1000 spinOrbit=0.01 minInvField=10 maxInvField=1000

‎fortran_implementation/get_command_line_input.f90

+112
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,112 @@
1+
module get_command_line_input
2+
implicit none
3+
private
4+
integer,parameter :: dp = selected_real_kind(15,307)
5+
character(32), dimension(8) :: inchars
6+
7+
8+
9+
type, public :: InputParameters
10+
integer :: length = 16394
11+
integer :: nsteps = 2000
12+
real(kind=8) :: theta = 0.d0
13+
real(kind=8) :: delta = 0.001d0
14+
real(kind=8) :: t2 = 0.5d0
15+
real(kind=8) :: mu = -4.96d0
16+
real(kind=8) :: spinOrbit = 0.005d0
17+
real(kind=8) :: minInvField = 30.d0
18+
real(kind=8) :: maxInvField = 6000.d0
19+
contains
20+
procedure :: grabinput => get_input1
21+
end type InputParameters
22+
23+
contains
24+
25+
subroutine get_input1(this)
26+
class(InputParameters), intent(in) :: this
27+
integer :: i
28+
character(len=32) :: arg
29+
30+
i = 0
31+
do
32+
call get_command_argument(i, arg)
33+
if (len_trim(arg) == 0) EXIT
34+
inchars(i+1) = trim(arg)
35+
i = i+1
36+
enddo
37+
38+
call initialize_vars(this,i)
39+
40+
end subroutine get_input1
41+
42+
43+
44+
subroutine initialize_vars(this,imax)
45+
class(InputParameters) :: this
46+
integer, intent(in) :: imax
47+
integer :: i,intvalue
48+
real(kind=8) :: realvalue
49+
character(len=32) :: argtot,name,val
50+
51+
do i=2,imax
52+
call split_string(inchars(i),name,val,'=')
53+
54+
select case (name)
55+
case ('length')
56+
read(val,'(i10)') intvalue
57+
this%length = intvalue
58+
case ('mu')
59+
read(val, '(F7.5)') realvalue
60+
this%mu = realvalue
61+
case ('theta')
62+
read(val, '(F10.5)') realvalue
63+
this%theta = realvalue
64+
case ('nsteps')
65+
read(val,'(i10)') intvalue
66+
this%nsteps = intvalue
67+
case ('delta')
68+
read(val,'(F10.5)') realvalue
69+
this%delta = realvalue
70+
case ('t2')
71+
read(val, '(F4.2)') realvalue
72+
this%t2 = realvalue
73+
case ('spinOrbit')
74+
read(val, '(F7.5)') realvalue
75+
this%spinOrbit = realvalue
76+
case ('minInvField')
77+
read(val,'(F12.5)') realvalue
78+
this%minInvField = realvalue
79+
case ('maxInvField')
80+
read(val,'(F12.5)') realvalue
81+
this%maxInvField = realvalue
82+
case default
83+
write(6,*) trim(inchars(i))," is an invalid Input argument!!"
84+
write(6,*) "Input must be of the form arg=xxx"
85+
write(6,*) "With arg= length, theta, nsteps, delta,t2, spinOrbit minInvField, maxInvField"
86+
write(6,*) "Exiting the program ....."
87+
call exit
88+
end select
89+
enddo
90+
91+
end subroutine initialize_vars
92+
93+
94+
95+
96+
subroutine split_string(instring, string1, string2, delim)
97+
! split a string into 2 either side of a delimiter token
98+
! taken from https://gist.github.com/ponderomotion/
99+
character(32) :: instring
100+
character :: delim
101+
character(32),INTENT(OUT):: string1,string2
102+
integer :: index
103+
104+
instring = TRIM(instring)
105+
106+
index = SCAN(instring,delim)
107+
string1 = instring(1:index-1)
108+
string2 = instring(index+1:)
109+
110+
end subroutine split_string
111+
112+
end module get_command_line_input

‎fortran_implementation/linalg.f90

+54
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,54 @@
1+
module linalg
2+
3+
Implicit none
4+
private
5+
integer,parameter :: dp = selected_real_kind(15,307)
6+
complex(kind=8),allocatable,dimension(:,:)::A
7+
complex(kind=8),allocatable,dimension(:,:)::res
8+
complex(kind=8),allocatable,dimension(:)::WORK
9+
integer,allocatable,dimension(:)::IPIV
10+
integer i,j,info,error, M
11+
public :: findinv
12+
13+
contains
14+
15+
function findinv(Ain, dim) result(res)
16+
integer :: dim
17+
complex(kind=8), dimension(dim,dim)::Ain
18+
complex(kind=8),dimension(dim,dim)::res
19+
!definition of the test matrix A
20+
M= dim
21+
22+
allocate(A(M,M),WORK(M),IPIV(M),stat=error)
23+
24+
A = Ain
25+
26+
27+
if (error.ne.0)then
28+
print *,"error:not enough memory"
29+
stop
30+
end if
31+
32+
call ZGETRF(M,M,A,M,IPIV,info)
33+
if(info .eq. 0) then
34+
! write(*,*)"succeded"
35+
else
36+
write(*,*)"failed"
37+
end if
38+
39+
call ZGETRI(M,A,M,IPIV,WORK,M,info)
40+
if(info .eq. 0) then
41+
! write(*,*)"succeded"
42+
else
43+
write(*,*)"failed"
44+
end if
45+
46+
res = A
47+
48+
deallocate(A,IPIV,WORK)
49+
50+
51+
52+
end function findinv
53+
54+
end module linalg

‎fortran_implementation/main_recursive.f90

+42
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,42 @@
1+
program recursive_dos
2+
use solver_class
3+
use t2g_hamiltonian_class
4+
use get_command_line_input
5+
implicit none
6+
7+
type(InputParameters) :: params !Declare an instance of the inputvariables
8+
type(Hamiltonian) :: h1 ! Declare an instance of the Hamiltonian
9+
type(Solver) :: solver_instance ! Declare the solver instance
10+
11+
!We will use double precision throughout
12+
integer,parameter :: dp = selected_real_kind(15,307)
13+
!Variables we need from the command line
14+
integer :: length, numsteps, polarAngle
15+
16+
17+
!1. Grab some input from the command line
18+
params = InputParameters()
19+
call params%grabinput()
20+
21+
22+
!2. Construct the Hamiltonian
23+
h1 = Hamiltonian(spin_orbit=params%spinOrbit, mu=params%mu, t2=params%t2, &
24+
theta=params%theta, delta=params%delta )
25+
26+
!3. Now construct the solver instance
27+
solver_instance = Solver(Ham=h1, Nsize=params%length, num_steps=params%nsteps, &
28+
min_inv_field = params%minInvField, &
29+
max_inv_field = params%maxInvField)
30+
31+
!4. And now, run the solver!!
32+
call solver_instance%run()
33+
34+
35+
36+
37+
38+
39+
end program recursive_dos
40+
41+
42+

‎fortran_implementation/makefile

+34
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
#makefile for qo code, using mac Accelerate framework
2+
RM = rm -f
3+
4+
OBJS = t2g_hamiltonian_class.o linalg.o get_command_line_input.o solver_class.o main_recursive.o
5+
6+
COMP = ifort
7+
8+
#LDFLAGS=-framework Accelerate
9+
LDFLAGS = -L/opt/intel/compilers_and_libraries_2016.2.146/mac/mkl/lib/intel64/ -lmkl_lapack95_lp64 -lmkl_intel_lp64 -lmkl_intel_thread -qopenmp -lmkl_core -lpthread
10+
11+
qo: $(OBJS)
12+
$(COMP) $(FLAGS) $(LDFLAGS) -o qo $(OBJS)
13+
14+
main_recursive.o: main_recursive.f90 t2g_hamiltonian_class.mod solver_class.mod get_command_line_input.mod
15+
$(COMP) $(FLAGS) $(LDFLAGS) -c $<
16+
17+
18+
t2g_hamiltonian_class.mod t2g_hamiltonian_class.o: t2g_hamiltonian_class.f90
19+
$(COMP) $(FLAGS) $(LDFLAGS) -c $<
20+
21+
linalg.mod linalg.o: linalg.f90
22+
$(COMP) $(FLAGS) $(LDFLAGS) -c $<
23+
24+
get_command_line_input.mod get_command_line_input.o: get_command_line_input.f90
25+
$(COMP) $(FLAGS) $(LDFLAGS) -c $<
26+
27+
solver_class.mod solver_class.o: solver_class.f90 t2g_hamiltonian_class.mod linalg.mod
28+
$(COMP) $(FLAGS) $(LDFLAGS) -c $<
29+
30+
clean:
31+
$(RM) *.o *.mod Thubb2
32+
33+
debug: $(OBJS)
34+
$(COMP) $(OBJS) $(LDFLAGS) -g -o qo_debug

‎fortran_implementation/solver_class.f90

+187
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,187 @@
1+
module solver_class
2+
3+
use t2g_hamiltonian_class
4+
use linalg
5+
6+
!The Hamiltonian class - returns the diagonal and off diagonal
7+
!blocks of the Hamiltonian
8+
9+
implicit none
10+
private
11+
real :: pi = 3.1415926535897931d0 ! Class-wide private constant
12+
integer,parameter :: dp = selected_real_kind(15,307)
13+
integer :: dim_fix
14+
integer :: system_size
15+
complex(kind=8), dimension(:,:,:), allocatable :: dLeft,cLeft,dRight,cRight
16+
17+
!Construct the Solver class
18+
type, public :: Solver
19+
type(Hamiltonian) :: Ham
20+
integer :: num_steps = 1500
21+
integer :: Nsize = 8192
22+
integer :: dim = 6
23+
real(kind=8) :: min_inv_field = 30.d0
24+
real(kind=8) :: max_inv_field = 6000.d0
25+
logical :: verbose = .true.
26+
contains
27+
procedure :: run => run_iter
28+
! procedure :: dos => get_dos
29+
end type Solver
30+
31+
32+
33+
contains
34+
35+
subroutine initialize(this)
36+
class(Solver), intent(in) :: this
37+
38+
system_size = this%Nsize
39+
dim_fix = this%dim
40+
41+
allocate(dLeft(system_size,dim_fix,dim_fix),cLeft(system_size,dim_fix,dim_fix))
42+
allocate(dRight(system_size,dim_fix,dim_fix),cRight(system_size,dim_fix,dim_fix))
43+
44+
end subroutine initialize
45+
46+
47+
48+
subroutine run_iter(this)
49+
class(Solver), intent(in) :: this
50+
integer :: b
51+
real(kind=8) :: dos, dfield, bfield, invbfield
52+
53+
54+
55+
!system_size = this%Nsize
56+
do b=1,this%num_steps
57+
call initialize(this)
58+
dos = 0.d0
59+
dfield = (this%max_inv_field - this%min_inv_field)/this%num_steps
60+
bfield = 1.d0/(1.d0*this%min_inv_field + 1.d0*(b-1)*dfield)
61+
invbfield = 1.d0/bfield
62+
63+
64+
dos = get_dos(this,bfield)
65+
66+
if (this%verbose) then
67+
write(6,"(3F30.15,1I10)"), bfield, invbfield, dos, system_size
68+
endif
69+
70+
deallocate(dLeft,cLeft,dRight,cRight)
71+
72+
enddo
73+
74+
75+
76+
end subroutine run_iter
77+
78+
79+
80+
81+
82+
83+
function get_dos(this, bfield) result(dos2)
84+
!Gets a diagonal block of the Hamiltonian
85+
class(Solver), intent(in) :: this
86+
real(kind=8), intent(in) :: bfield
87+
real(kind=8) :: dos2
88+
complex(kind=8), dimension(dim_fix,dim_fix) :: block
89+
90+
dLeft(:,:,:) = dcmplx(0.d0,0.d0)
91+
cLeft(:,:,:) = dcmplx(0.d0,0.d0)
92+
dRight(:,:,:) = dcmplx(0.d0,0.d0)
93+
cRight(:,:,:) = dcmplx(0.d0,0.d0)
94+
95+
96+
call downsweep(this, bfield)
97+
call upsweep(this, bfield)
98+
dos2 = calc_dos(this, bfield)
99+
100+
end function get_dos
101+
102+
103+
104+
subroutine downsweep(this, bfield)
105+
class(Solver), intent(in) :: this
106+
real(kind=8), intent(in) :: bfield
107+
complex(kind=8), dimension(dim_fix,dim_fix) :: extra_term
108+
complex(kind=8), dimension(dim_fix,dim_fix) :: dinv
109+
integer :: i
110+
111+
do i=1,system_size-1
112+
if (i.eq.1) then
113+
dLeft(i,:,:) = this%Ham%diag(i,bfield,0.d0,0.d0)
114+
else
115+
extra_term = matmul(cLeft(i-1,:,:),this%Ham%rBlock(i-1,bfield,0.d0,0.d0))
116+
dLeft(i,:,:) = this%Ham%diag(i,bfield,0.d0,0.d0) + extra_term
117+
endif
118+
119+
dinv = findinv(dLeft(i,:,:), dim_fix)
120+
cLeft(i,:,:) = dcmplx(-1.d0,0.d0)*matmul(this%Ham%lBlock(i+1,bfield,0.d0,0.d0),dinv)
121+
enddo
122+
123+
i = system_size
124+
dLeft(i,:,:) = this%Ham%diag(i, bfield,0.d0,0.d0) &
125+
+ matmul(cLeft(i-1,:,:),this%Ham%rBlock(i-1,bfield,0.d0,0.d0))
126+
127+
end subroutine downsweep
128+
129+
130+
subroutine upsweep(this, bfield)
131+
class(Solver), intent(in) :: this
132+
real(kind=8), intent(in) :: bfield
133+
complex(kind=8), dimension(dim_fix,dim_fix) :: extra_term
134+
complex(kind=8), dimension(dim_fix,dim_fix) :: dinv
135+
integer :: i,j
136+
137+
do j=0,system_size-2
138+
i = system_size - j
139+
if (i.eq.system_size) then
140+
dRight(i,:,:) = this%Ham%diag(i,bfield,0.d0,0.d0)
141+
else
142+
extra_term = matmul(cRight(i+1,:,:),this%Ham%lBlock(i+1,bfield,0.d0,0.d0))
143+
dRight(i,:,:) = this%Ham%diag(i,bfield,0.d0,0.d0) + extra_term
144+
endif
145+
146+
dinv = findinv(dRight(i,:,:), dim_fix)
147+
cRight(i,:,:) = dcmplx(-1.d0,0.d0)*matmul(this%Ham%rBlock(i-1,bfield,0.d0,0.d0),dinv)
148+
enddo
149+
150+
i = 1
151+
dRight(i,:,:) = this%Ham%diag(i, bfield,0.d0,0.d0) &
152+
+ matmul(cRight(i+1,:,:),this%Ham%lBlock(i+1,bfield,0.d0,0.d0))
153+
154+
end subroutine upsweep
155+
156+
157+
158+
159+
function calc_dos(this, bfield) result(dos1)
160+
class(Solver), intent(in) :: this
161+
real(kind=8), intent(in) :: bfield
162+
real(kind=8) :: dos1
163+
complex(kind=8) :: trace
164+
165+
complex(kind=8), dimension(dim_fix,dim_fix) :: ginv, g
166+
integer :: i,j
167+
168+
dos1=0.d0
169+
do i=1,system_size
170+
ginv = - this%Ham%diag(i, bfield,0.d0,0.d0) + dLeft(i,:,:) + dRight(i,:,:)
171+
g = findinv(ginv, dim_fix)
172+
173+
trace = cmplx(0.d0,0.d0)
174+
do j=1,dim_fix
175+
trace = trace + g(j,j)
176+
enddo
177+
dos1 = dos1 - ( 1.d0/ (pi * system_size) ) * aimag(trace)
178+
!write(6,*) i,dos1
179+
enddo
180+
181+
end function calc_dos
182+
183+
184+
185+
end module solver_class
186+
187+

‎fortran_implementation/t2g_hamiltonian_class.f90

+143
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,143 @@
1+
module t2g_hamiltonian_class
2+
!The Hamiltonian class - returns the diagonal and off diagonal
3+
!blocks of the Hamiltonian
4+
5+
6+
implicit none
7+
private
8+
real :: pi = 3.1415926535897931d0 ! Class-wide private constant
9+
integer,parameter :: dp = selected_real_kind(15,307)
10+
real(kind=8), parameter :: ta = 1.d0
11+
integer, parameter :: dim_fix = 6
12+
13+
!Construct the Hamiltonian class
14+
type, public :: Hamiltonian
15+
real(kind=8) :: spin_orbit
16+
real(kind=8) :: mu
17+
real(kind=8) :: t2
18+
real(kind=8) :: theta
19+
real(kind=8) :: delta
20+
contains
21+
procedure :: dim => print_dim
22+
procedure :: diag => get_diagonal_block
23+
procedure :: rBlock => get_right_block
24+
procedure :: lBlock => get_left_block
25+
end type Hamiltonian
26+
27+
28+
29+
contains
30+
subroutine print_dim(this)
31+
class(Hamiltonian), intent(in) :: this
32+
print *, "Dimension of Hamiltonian is", this%spin_orbit
33+
end subroutine print_dim
34+
35+
36+
function get_diagonal_block(this, x, b, ky, kz) result(block)
37+
!Gets a diagonal block of the Hamiltonian
38+
class(Hamiltonian), intent(in) :: this
39+
integer, intent(in) :: x
40+
real(kind=8), intent(in) :: b,ky,kz
41+
42+
complex(kind=8), dimension(6,6) :: block
43+
44+
block(:,:) = (0.d0,0.d0)
45+
46+
47+
48+
block(1,1) = dcmplx(2.d0*ta*cos(2.d0*pi*b*x*cos(this%theta*pi/180.d0) - ky), this%delta)
49+
block(1,1) = block(1,1) + dcmplx(2.d0*ta*cos(-2.d0*pi*b*x*sin(this%theta*pi/180.d0) - kz), 0.d0)
50+
block(1,1) = block(1,1) + dcmplx(this%mu, 0.d0)
51+
block(4,4) = dcmplx(2.d0*ta*cos(2.d0*pi*b*x*cos(this%theta*pi/180.d0) - ky) &
52+
+ 2.d0*ta*cos(-2.d0*pi*b*x*sin(this%theta*pi/180.d0) - kz) &
53+
+ this%mu, this%delta)
54+
55+
56+
block(2,2) = dcmplx(2.d0*this%t2*cos(2.d0*pi*b*x*cos(this%theta*pi/180.d0) - ky) &
57+
+ 2.d0*ta*cos(-2.d0*pi*b*x*sin(this%theta*pi/180.d0) - kz) &
58+
+ this%mu, this%delta)
59+
block(5,5) = dcmplx(2.d0*this%t2*cos(2.d0*pi*b*x*cos(this%theta*pi/180.d0) - ky) &
60+
+ 2.d0*ta*cos(-2.d0*pi*b*x*sin(this%theta*pi/180.d0) - kz) &
61+
+ this%mu, this%delta)
62+
63+
64+
block(3,3) = dcmplx(2.d0*ta*cos(2.d0*pi*b*x*cos(this%theta*pi/180.d0) - ky) &
65+
+ 2.d0*this%t2*cos(-2.d0*pi*b*x*sin(this%theta*pi/180.d0) - kz) &
66+
+ this%mu, this%delta)
67+
block(6,6) = dcmplx(2.d0*ta*cos(2.d0*pi*b*x*cos(this%theta*pi/180.d0) - ky) &
68+
+ 2.d0*this%t2*cos(-2.d0*pi*b*x*sin(this%theta*pi/180.d0) - kz) &
69+
+ this%mu, this%delta)
70+
71+
!Now add the (complex) spin orbit terms
72+
block(1,2) = dcmplx(0.d0, this%spin_orbit)
73+
block(2,1) = dcmplx(0.d0, -1.d0*this%spin_orbit)
74+
75+
block(1,6) = dcmplx(-1.d0*this%spin_orbit,0.d0)
76+
block(6,1) = dcmplx(-1.d0*this%spin_orbit,0.d0)
77+
78+
block(2,6) = dcmplx(0.d0, this%spin_orbit)
79+
block(6,2) = dcmplx(0.d0, -1.d0*this%spin_orbit)
80+
81+
block(3,4) = dcmplx(this%spin_orbit,0.d0)
82+
block(4,3) = dcmplx(this%spin_orbit,0.d0)
83+
84+
block(3,5) = dcmplx(0.d0, -1.d0*this%spin_orbit)
85+
block(5,3) = dcmplx(0.d0, this%spin_orbit)
86+
87+
block(4,5) = dcmplx(0.d0,-1.d0*this%spin_orbit)
88+
block(5,4) = dcmplx(0.d0, this%spin_orbit)
89+
90+
! write(6,"(8F20.10)") x, b, block(1,1),block(2,2),block(3,3)
91+
92+
end function get_diagonal_block
93+
94+
95+
96+
function get_right_block(this, x, b, ky, kz) result(block)
97+
!Gets a diagonal block of the Hamiltonian
98+
class(Hamiltonian), intent(in) :: this
99+
integer, intent(in) :: x
100+
real(kind=8), intent(in) :: b,ky,kz
101+
102+
complex(kind=8), dimension(6,6) :: block
103+
104+
block(:,:) = (0.d0,0.d0)
105+
106+
!Now add the (complex) spin orbit terms
107+
block(1,1) = dcmplx(this%t2,0.d0)
108+
block(4,4) = dcmplx(this%t2,0.d0)
109+
110+
block(2,2) = dcmplx(ta,0.d0)
111+
block(5,5) = dcmplx(ta,0.d0)
112+
113+
block(3,3) = dcmplx(ta,0.d0)
114+
block(6,6) = dcmplx(ta,0.d0)
115+
116+
end function get_right_block
117+
118+
119+
function get_left_block(this, x, b, ky, kz) result(block)
120+
!Gets a diagonal block of the Hamiltonian
121+
class(Hamiltonian), intent(in) :: this
122+
integer, intent(in) :: x
123+
real(kind=8), intent(in) :: b,ky,kz
124+
125+
complex(kind=8), dimension(6,6) :: block
126+
127+
block(:,:) = (0.d0,0.d0)
128+
129+
!Now add the (complex) spin orbit terms
130+
block(1,1) = dcmplx(this%t2,0.d0)
131+
block(4,4) = dcmplx(this%t2,0.d0)
132+
133+
block(2,2) = dcmplx(ta,0.d0)
134+
block(5,5) = dcmplx(ta,0.d0)
135+
136+
block(3,3) = dcmplx(ta,0.d0)
137+
block(6,6) = dcmplx(ta,0.d0)
138+
139+
end function get_left_block
140+
141+
end module t2g_hamiltonian_class
142+
143+

0 commit comments

Comments
 (0)
Please sign in to comment.