使用Fortran计算Hall电导
整理用Fortran计算某一条能带的Hall电导，这里使用的QWZ模型。{:.info} 前言这里使用的是最简单的QWZ模型，其哈密顿量为 H = (m_0 - t_x \cos(k_x) - t_y \cos(k_y))\sigma_z + a_x \sin(k_x) \sigma_x + a_y * \sin(k_y) * sigma_y程序实现上目前是分别计算每一条能带的Hall电导，而不是计算费米面以下电子贡献的电导。当化学势落在能隙中的时候，这里计算的Hall电导就是单独一条能带贡献，如果化学势与能带有交点，那么Hall电导的计算应该是考虑所有费米面以下的贡献，这里先不考虑这件事情。 代码串行版123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185module code_param implicit none integer, parameter :: dp = kind(1.0d0) ! 双精度 real(dp),parameter::pi = acos(-1.0_dp) complex(dp),parameter::im = (0.,1.) ! Imagine unit integer Nk,numk_bz ...
Fortran踩坑记录(持续更新中)
最近在用Fortran写程序的时候，在意想不到的地方踩了坑，这篇Bolg整理一下。{:.info} 前言用Fortran写程序主打一个方便，只要把公式翻译成代码就行了，但是在码代码的时候还是遇到了一些意想不到的问题，这里整理一下，方便自己查阅闭坑。 费米分布函数在把费米分布转成代码为123456789function fermi(ek) ! 费米分布函数 integer,parameter::dp = kind(1.0) real(dp) fermi,ek,kbt kbt = 0.01 fermi = 1.0/(exp(ek/kbt) + 1) ! fermi = ek returnend 但这里有$e^x$指数函数，如果$x$足够大的话，如果是在单精度下面(dp = kind(1.0))，此时分布函数就会溢出。不过就算使用双精度(kind(1.0d0))，还是无法完全避免数据溢出这个问题，所以最好的方式就是给定一个阈值，判断在哪些范围内可以使用费米分布函数，其他范围就是0或者1，取决于该量子态是占据还是非占据12345678910function fermi(ek) ! 万万不能直接用费米分布函数，会存在浮点溢出 integer,parameter::dp = kind(1.0) real(dp) fermi,ek,kbt kbt = 0.001 if(ek/kbt>-40 .and. ek/kbt<40) fermi = 1.0/(exp(ek/kbt) + 1) if(ek/kbt <-40) fermi = 1.0 if(ek/kbt >40) fermi = 0.0 returnend 除了上面的截断方式，还可以在能量大于零和小于零的情形下，对费米分布函数做一下变形处理，在数值上避免发散 f(E) = \begin{cases} \frac{e^{-E/(k_B T)}}{1 + e^{-E/(k_B T)}}, & E > 0 \\ \frac{1}{1 + e^{E/(k_B T)}}, & E \leq 0 \end{cases}123456789101112real(dp) function...
正方晶格超导序参量自洽
最近在学习计算超导体系的超流权重，其中要涉及到在确定相互作用$U$下自洽出序参量，虽然之前也会，为了保证所有结果的正确性，这里重新将一些过程记录一下，方便后面查阅理解。{:.info} 超导序参量自洽 代码Fortran Version123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265!==============================================================================================================================================================! 计算自由电子的超流权重! H(k) = t(cos...
平带中的紧致局域态(compact localized states in flat band)
平带最近受到很大的关注，最近正好在研究量子几何，里面涉及到平带相关的内容，这篇笔记整理了一下在学习平带时相关的紧致局域态到底是个什么东西。{:.info} Main 参考文献 Designing flat-band tight-binding models with tunable multifold band touching points= General construction of flat bands with and without band crossings based on wave function singularity 公众号相关内容均会在公众号进行同步，若对该Blog感兴趣，欢迎关注微信公众号。{:.info} Email yxliphy@gmail.com
Fortran常用函数的封装与收集
这里整理一下在用Fortran编写物理程序的时候经常用的几个函数封装，目前只是简单的几个矩阵对角化函数，后续有新的函数也会继续增加。{:.info} 前言最近因为遇到了计算瓶颈，需要用Fortran来重写代码，其中就遇到了对矩阵的一些操作，虽然lapack已经提供了很好的函数，但是在实际程序编写中还是有些麻烦，这里就对目前我用到的函数进行了一些封装，方便调用。 一般方矩阵对角化1234567891011121314151617181920212223242526subroutine diagonalize_general_matrix(matdim,matin,matout,mateigval) use param,only:dp,im implicit none integer matdim, LDA, LDVL, LDVR, LWMAX, INFO, LWORK complex(dp) mateigval(matdim),matin(matdim,matdim),matout(matdim,matdim) real,allocatable::VL(:,:),VR(:,:),WR(:),WI(:),WORK(:) LDA = matdim LDVL = matdim LDVR = matdim LWMAX = 2 * matdim + matdim**2 allocate(VL( LDVL, matdim), VR(LDVL,matdim), WR(matdim), WI(matdim), WORK(LWMAX)) LWORK = -1 matout = matin CALL SGEEV( 'V', 'V', matdim, matout, LDA, WR, WI, VL, LDVL, VR, LDVR, WORK, LWORK, INFO ) LWORK = MIN( LWMAX, INT( WORK( 1 ) ) ) CALL SGEEV( 'V', 'V', matdim, matout, LDA, WR, WI, VL, LDVL, VR, LDVR, WORK, LWORK,...
