Fortran结合MPI并行计算自旋极化率
在之前的博客Julia的MPI并行计算极化率(重复Bilayer Two-Orbital Model of La$_3$Ni$_2$O$_7$ under Pressure)利用Julia计算了一个模型的极化率，但是在撒点数量增加的时候计算速度还是堪忧。这里就返祖利用Fortran语言再结合MPI并行来重写code。根据测试发现计算速度显著提升。
{:.info}
前言
虽然使用Julia在计算的时候速度已经相当快了，但是如果涉及到计算极化率这种需要在布里渊区撒点很多才能精确计算的量，Julia在计算的时候需要的时间还是有点久。况且想要速度更快，就需要仔细的去对代码进行优化，想想还是挺复杂的。最后考虑直接返祖，就用Fortran(公式翻译器)来计算极化率。
关于公式具体的内容可以参考Julia的MPI并行计算极化率(重复Bilayer Two-Orbital Model of La$_3$Ni$_2$O$_7$ under Pressure)这篇Blog，或者去查看原文Bilayer Two-Orbital Model of La$_3$Ni$_2$O$_7$ under Pressure，下面直接上代码。
代码
这里在写的时候偷懒了，将哈密顿量设置为全局变量了，而且对角化厄米矩阵的函数并没有进行封装，调用之后就是直接对角化哈密顿量。实际上正确且安全的写法就是通过子过程返回哈密顿量，并将其传递给对角化函数计算本征值和本征矢量，这样才能让程序具有通用性。不过事情我是知道的，但这个程序很简单，就先不做这样做了，后续写大程序的时候就会规范了。
- 计算耗时
这里撒点数量为 256 * 256，调用了64个核，计算时间如下下面就是完整的代码了1
2
3======== Job starts at 2024-04-15 15:25:16 on n26 ========
Start Fortran code
======== Job ends at 2024-04-15 15:25:37 on n26 ========计算结果1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242module param
implicit none
integer kn,hn
parameter(kn = 64,hn = 4)
real,parameter::pi = 3.1415926535897
real,parameter::omega = 0.00
complex,parameter::im = (0.,1.) !Imagine unit
complex Ham(hn,hn),Umat(hn,hn),ones(hn,hn) ! Hamiltonian and interaction Matrix
real t0,t1x,t1z,t2x,t2z,t3xz,t4xz,tvx,tvz,ex,ez ! 哈密顿量参数
real U0,J0 ! 相互作用参数
real valmesh(2,-kn:kn-1,-kn:kn-1,hn)
complex vecmesh(2,-kn:kn-1,-kn:kn-1,hn,hn)
complex chi0(-kn:kn-1,-kn:kn-1,hn,hn),chi(-kn:kn-1,-kn:kn-1,hn,hn)
! LAPACK PACKAGE PARAM
integer::lda = hn
integer,parameter::lwmax = 2*hn+hn**2
real,allocatable::val(:)
complex,allocatable::work(:)
real,allocatable::rwork(:)
integer,allocatable::iwork(:)
integer lwork ! at least 2*hn+N**2
integer lrwork ! at least 1 + 5*hn +2*hn**2
integer liwork ! at least 3 +5*hn
integer info
end module param
!================================================================================================================================================================================================
program main
use param
use mpi
!------------------------------------
complex temp1(hn,hn),temp2(hn,hn),temp3
integer iky,ikx
real qx,qy,local_rechi(-kn:kn-1,-kn:kn-1,4),rechi(-kn:kn-1,-kn:kn-1,4)
!---------------------------------
integer numcore,indcore,ierr
character(len=20)::filename,char1,char2
!---------------------------------
external::cheevd
allocate(val(hn))
allocate(work(lwmax))
allocate(rwork(1+5*hn+2*hn**2))
allocate(iwork(3+5*hn))
!---------------------------------
call MPI_INIT(ierr) ! 初始化进程
call MPI_COMM_RANK(MPI_COMM_WORLD, indcore, ierr) ! 得到本进程在通信空间中的rank值,即在组中的逻辑编号(该 rank值为0到p-1间的整数,相当于进程的ID。)
call MPI_COMM_SIZE(MPI_COMM_WORLD, numcore, ierr) !获得进程个数 size, 这里用变量p保存
call MPI_Barrier(MPI_COMM_WORLD,ierr)
nki = floor(indcore * (2.0 * kn)/numcore) - kn
nkf = floor((indcore + 1) * (2.0 * kn)/numcore) - kn - 1
! 遍历BZ求极化率
do iky = nki,nkf
qy = pi * iky/kn
do ikx = -kn,kn - 1
qx = pi * ikx/kn
call chi0cal(qx,qy,chi0(iky,ikx,:,:)) ! 得到裸极化率
call inv(ones - matmul(chi0(iky,ikx,:,:),Umat),temp2) ! 矩阵求逆
chi(iky,ikx,:,:) = matmul(temp2,chi0(iky,ikx,:,:))
temp3 = sum(chi(iky,ikx,:,:))
local_rechi(iky,ikx,1) = qx
local_rechi(iky,ikx,2) = qy
local_rechi(iky,ikx,3) = real(temp3)
local_rechi(iky,ikx,4) = aimag(temp3)
end do
end do
call MPI_Barrier(MPI_COMM_WORLD,ierr)
! call MPI_Gather(local_rechi, 2 * kn, MPI_COMPLEX, rechi, 2 * kn, MPI_COMPLEX, 0, MPI_COMM_WORLD,ierr)
call MPI_Reduce(local_rechi, rechi, (2 * kn)**2 * 4, MPI_REAL, MPI_SUM, 0, MPI_COMM_WORLD,ierr)
if (indcore .eq. 0) then
char1 = "fortran-chi-"
write(char2,"(I3.3)")kn
filename = trim(char1)//trim(char2)
char1 = ".dat"
filename = trim(filename)//trim(char1)
open(12,file = filename)
do iky = -kn,kn - 1
do ikx = -kn,kn - 1
write(12,"(4F5.3)")rechi(iky,ikx,1),rechi(iky,ikx,2),abs(rechi(iky,ikx,3)),abs(rechi(iky,ikx,4))
end do
end do
close(12)
end if
call MPI_Finalize(ierr)
stop
end program main
!================================================================================================================================================================================================
subroutine matset(kx,ky)
! 矩阵赋值
use param
real kx,ky
integer k0
t0 = 1.0
t1x = -0.483 * t0
t1z = -0.110 * t0
t2x = 0.069 * t0
t2z = -0.017 * t0
t3xz = 0.239 * t0
t4xz = -0.034 * t0
tvx = 0.005 * t0
tvz = -0.635 * t0
ex = 0.776 * t0
ez = 0.409 * t0
Ham = 0.0
Ham(1, 1) = 2 * t1x * (cos(kx) + cos(ky)) + 4 * t2x * cos(kx) * cos(ky) + ex
Ham(2, 2) = 2 * t1z * (cos(kx) + cos(ky)) + 4 * t2z * cos(kx) * cos(ky) + ez
Ham(1, 2) = 2 * t3xz * (cos(kx) - cos(ky))
Ham(2, 1) = 2 * t3xz * (cos(kx) - cos(ky))
Ham(3, 3) = 2 * t1x * (cos(kx) + cos(ky)) + 4 * t2x * cos(kx) * cos(ky) + ex
Ham(4, 4) = 2 * t1z * (cos(kx) + cos(ky)) + 4 * t2z * cos(kx) * cos(ky) + ez
Ham(3, 4) = 2 * t3xz * (cos(kx) - cos(ky))
Ham(4, 3) = 2 * t3xz * (cos(kx) - cos(ky))
Ham(1, 3) = tvx
Ham(1, 4) = 2 * t4xz * (cos(kx) - cos(ky))
Ham(2, 3) = 2 * t4xz * (cos(kx) - cos(ky))
Ham(2, 4) = tvz
Ham(3, 1) = tvx
Ham(4, 1) = 2 * t4xz * (cos(kx) - cos(ky))
Ham(3, 2) = 2 * t4xz * (cos(kx) - cos(ky))
Ham(4, 2) = tvz
!---------------------------------------------------------------------
! 相互作用矩阵赋值
U0 = 3.0
J0 = 0.4
Umat(1,1) = U0
Umat(2,2) = U0
Umat(3,3) = U0
Umat(4,4) = U0
Umat(1,2) = J0/2
Umat(2,1) = J0/2
Umat(3,4) = J0/2
Umat(4,3) = J0/2
!---------------------------------------------------------------------
! 单位矩阵
do k0 = 1,hn
ones(k0,k0) = 1
end do
return
end subroutine
!================================================================================================================================================================================================
subroutine chi0cal(qx,qy,re1)
! 计算极化率 返回到re1
use param
integer ikx,iky,l1,l2,e1,e2
real qx,qy,kx,ky
complex re1(hn,hn)
do iky = -kn,kn - 1
ky = pi * iky/kn
do ikx = -kn,kn - 1
kx = pi * ikx/kn
! k
call matset(kx,ky)
call eigSol()
valmesh(1,iky,ikx,:) = val(:)
vecmesh(1,iky,ikx,:,:) = Ham(:,:)
! k + q
call matset(kx + qx,ky + qy)
call eigSol()
valmesh(2,iky,ikx,:) = val(:)
vecmesh(2,iky,ikx,:,:) = Ham(:,:)
! 计算极化率
do l1 = 1,hn ! orbit ondex
do l2 = 1,hn
do e1 = 1,hn ! band index
do e2 = 1,hn
re1(l1,l2) = re1(l1,l2) + (fermi(valmesh(1,iky,ikx,e1)) - fermi(valmesh(2,iky,ikx,e2)))/(im * (omega + 0.0001) + valmesh(1,iky,ikx,e1) - valmesh(2,iky,ikx,e2))&
* conjg(vecmesh(2,iky,ikx,l1,e2)) * vecmesh(2,iky,ikx,l2,e2) * conjg(vecmesh(1,iky,ikx,l2,e1)) * vecmesh(1,iky,ikx,l1,e1)
end do
end do
end do
end do
end do
end do
re1 = re1/(2 * kn)**2
return
end subroutine
!================================================================================================================================================================================================
function fermi(ek)
! 费米分布函数
implicit none
real fermi,ek,kbt
kbt = 0.001
fermi = 1/(exp(ek/kbt) + 1)
return
end
!================================================================================================================================================================================================
function equivkpq(i0)
! 找到BZ中k+q等价与k的索引
use param
integer equivkpq,i0
if (i0 <= kn/2 .and. i0 > -kn/2) equivkpq = i0
if (i0 > kn/2) equivkpq = i0 - kn
if (i0 <= -kn/2) equivkpq = i0 + kn
end
!================================================================================================================================================================================================
subroutine eigSol()
! 对角化得到本征值w和本征矢量Ham
use param
integer m
lwork = -1
liwork = -1
lrwork = -1
call cheevd('V','U',hn,Ham,lda,val,work,lwork,rwork,lrwork,iwork,liwork,info)
lwork = min(2 * hn + hn**2, int( work( 1 ) ) )
lrwork = min(1 + 5 * hn + 2 * hn**2, int( rwork( 1 ) ) )
liwork = min(3 + 5 * hn, iwork( 1 ) )
call cheevd('V','U',hn,Ham,lda,val,work,lwork,rwork,lrwork,iwork,liwork,info)
if( info .GT. 0 ) then
open(11,file = "mes.dat",status = "unknown")
write(11,*)'The algorithm failed to compute eigenvalues.'
close(11)
end if
! open(12,file="eigval.dat",status="uknnown")
! do m = 1,N
! write(12,*)m,val(m)
! end do
! close(12)
return
end subroutine eigSol
!================================================================================================================================================================================================
subroutine inv(matin,matout)
! 矩阵求逆
use param
complex,intent(in) :: matin(hn,hn)
complex:: matout(size(matin,1),size(matin,2))
real:: work2(size(matin,1)) ! work2 array for LAPACK
integer::ipiv(size(matin,1)) ! pivot indices
! Store matin in matout to prevent it from being overwritten by LAPACK
matout = matin
! SGETRF computes an LU factorization of a general M - by - N matrix A
! using partial pivoting with row interchanges .
call CGETRF(hn,hn,matout,hn,ipiv,info)
if (info.ne.0) stop 'Matrix is numerically singular!'
! SGETRI computes the inverse of a matrix using the LU factorization
! computed by SGETRF.
call CGETRI(hn,matout,hn,ipiv,work2,hn,info)
if (info.ne.0) stop 'Matrix inversion failed!'
return
end subroutine inv
!================================================================================================================================================================================================
绘图程序
1 | def plotchi(numk): |
公众号
相关内容均会在公众号进行同步，若对该Blog感兴趣，欢迎关注微信公众号。
{:.info}
|
yxliphy@gmail.com |
本博客所有文章除特别声明外，均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来源 Yu-Xuan's Blog！
wechat
alipay
相关推荐
2021-06-26
快速格林函数方法计算Chern绝缘体边界态
这篇博客整理一下如何利用格林函数方法来计算Chern绝缘体不同边界上的边界态.{:.info}在Chern Insulator边界态及Chern数计算这篇博客中提供了计算Chern绝缘体边界态的程序,但是因为那个方法中是在一个cylinder上进行计算的,所以会存在两个边界,从而也就会在能谱中看到有两个边界态,这在有时候的研究中是不太方便的,这里就像通过格林函数的方法,计算一个半无限大的系统,因为只存在一个边界,所以对于Chern绝缘体来说此时就可以得到只有一个边界态的能谱图,而且还可以分别计算左右两端的边界态,可以发现其对应的费米速度是相反的. 计算公式这里使用的就是Highly convergent schemes for the calculation of bulk and surface Green functions这篇文章中的计算方法,不过需要对写程序中一些具体内容进行一下说明. 当将一个动量空间中的哈密顿量沿某一个方向取开边界,另一个方向取周期边界的时候,对应的哈密顿量矩阵为 H=\left[\begin{array}{ccccc}H_{00}&H_{01}&0&0&0\\ H_{10}&H_{11}&H_{12}&0&0\\0&H_{21}&H_{22}&H_{23}&0\\ 0&0&H_{32}&H_{33}&\cdots\\ 0&0&0&\cdots&\cdots \end{array}\right]因为哈密顿量是厄米的,所以就会有$H_{01}=H_{i,i+1}=H^\dagger_{i+1,i},H_{00}=H_{ii}=H_{i+1.i+1}$. 想要得到格林函数 (\omega-H)G=I可以通过下面的迭代方程进行 \begin{equation}\begin{aligned}\alpha_i&=\alpha_{i-1}(\omega-\epsilon_{i-1})^{-1}\alpha_{i-1}\\ \beta_i&=\beta_{i-1}(\omega-\epsilon_{i-1})^{-1}\beta_{i-1}\\...
2024-04-19
Fortran结合MPI并行
该Bolg整理在Fortran用MPI编写并行程序，用在科学计算中速度提升肉眼可见，{:.info} 代码解析先给一个完整的代码，是用来计算自旋极化率的，实际上在Fortran结合MPI并行计算自旋极化率中已经出现过了，这里就是解读一下其中的并行部分是如何用MPI写的。 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242module param implicit none integer kn,hn parameter(kn = 64,hn = 4) real,parameter::pi = 3.1415926535897 real,parameter::omega = 0.00 complex,parameter::im = (0.,1.) !Imagine unit complex Ham(hn,hn),Umat(hn,hn),ones(hn,hn)...
2019-01-01
p-wave 超导体Vortex中的Majorana zero mode
整理研究生接触的第一个科研工作 研究生开始时，导师给的第一个任务就是重复一篇关于p-wave超导Vortex中Majorana零能模的文章。从刚开始一无所知的小白，到完整复现文章内容还是学到了很多东西，最后将自己学习中的想法和代码整理成一个项目，如果可以帮到别人最好了。整理的过程中也是自己对所学内容的进一步梳理，可以更好的理解文章内容和计算方法。{:.success} BdG equation solve在超导的问题中，通常会涉及到BdG方程的求解。本代码是用来重复Vortex pinning by the point potential in topological superconductors:A scheme for braiding Majorana bound states论文中关于Majorana费米子在拓扑超导中存在的一些问题。 BdG equation哈密顿量本征值问题求解的关键问题在于哈密顿量的对角化，而BdG方程的关键点在于引入Bogoliubov变换，引 进了准粒子算符来对非对角化的哈密顿量进行对角化处理，该准粒子算符具有粒子空穴对称性。 Important关于如何将论文中的模型一步步的转换成Fortran程序，在这里提供一份尚未完善的帮助手册，读者可以借助这份手册同时对照提供的程序，可以帮助读者更好的明白程序中每一步的含义，以及如何去重复该论文的结果，同时也可以在该程序的基础上进行进一步的数值计算。 文件地址为Manual 程序介绍 整个项目中有多个文件夹，按英文序号进行命名，记录了作者一步一步对代码的编写，同时里面也包括了对论文中公式到矩阵形式进行转换的想法。 final_version中有两个文件夹Linux 和Windows，分别是程序的Linux版本和windows版本 每个文件夹中都有程序源代码，已经通过了检验，同时也附上了计算得到的最终数据。 同时作者还提供了一个由Mathematic编写的简单的绘图程序，可以用来及时的展示计算得到的结果（visual.nb） 在运行Visual.nb文件时，只需将该文件和数据结果放在同一个文件夹之中即可 final_version中同时有origin的绘图文件pic，用来比较清晰的展示计算结果 ...
评论
