1
!{\src2tex{textfont=tt}}
7
!! This routine is called in scfcv.f to compute the derivative of
8
!! ground-state wavefunctions with respect to k (du/dk) by finite differencing
9
!! on neighbouring k points
10
!! Work for nsppol=1 or 2, but only accept nspinor=1,
13
!! Copyright (C) 2001-2007 ABINIT group (NSAI).
16
!! bdberry(4)=band limits for Berry phase contributions (or du/dk)
17
!! spin up and spin down (bdberry(3:4) is irrelevant when nsppol=1)
18
!! cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions
19
!! gprimd(3,3)=reciprocal space dimensional primitive translations
20
!! hdr <type(hdr_type)>=the header of wf, den and pot files
21
!! istwfk(nkpt_)=input option parameter that describes the storage of wfs
22
!! kberry(3,20)= different delta k for Berry phases(or du/dk),
23
!! in unit of kptrlatt only kberry(1:3,1:nberry) is relevant
24
!! kg(3,mpw*mkmem)=reduced planewave coordinates
25
!! kpt_(3,nkpt_)=reduced coordinates of k points generated by ABINIT,
26
!! kpt_ samples half the BZ if time-reversal symetrie is used
27
!! kptopt=2 when time-reversal symmetry is used
28
!! kptrlatt(3,3)=k-point lattice specification
29
!! mband=maximum number of bands
30
!! mgfft=maximum size of 1D FFTs
31
!! mkmem=number of k points which can fit in memory; set to 0 if use disk
32
!! mpi_enreg=informations about MPI parallelization
33
!! mpw=maximum dimensioned size of npw
34
!! natom=number of atoms in cell
35
!! nband(nkpt*nsppol)=number of bands at each k point, for each polarization
36
!! nberry=number of Berry phases(or du/dk) to be computed
37
!! nkpt=number of k points
38
!! npwarr(nkpt)=number of planewaves in basis at this k point
39
!! nspinor=number of spinorial components of the wavefunctions
40
!! nsppol=1 for unpolarized, 2 for spin-polarized
41
!! rprimd(3,3)=dimensional real space primitive translations (bohr)
42
!! unddk=unit number for ddk file
43
!! unkg=unit number for (k+G) basis sphere data
44
!! wffnow=struct info for wf disk file
47
!! (the ddk wavefunctions are written on disk)
52
!! Cleaning, checking for rules
53
!! Should allow for time-reversal symmetry (istwfk)
54
!! WARNING : the use of nspinor is completely erroneous
58
!! cmatrix(:,:,:)= overlap matrix of size maxband*maxband
59
!! cg_index(:,:,:)= unpacked cg index array for specific band,
60
!! k point and polarization.
61
!! det(2,2)= intermediate output of Lapack routine zgedi.f
62
!! dk(3)= step taken to the next k mesh point along the kberry direction
63
!! gpard(3)= dimensionalreciprocal lattice vector G along which the
64
!! polarization is computed
65
!! kg_kpt(:,:,:)= unpacked reduced planewave coordinates with subscript of
66
!! planewave and k point
67
!! kpt(3,nkpt)=reduced coordinates of k-point grid that samples the whole BZ
68
!! kpt_flag(nkpt)=kpt_flag(ikpt)=0 when the wf was generated by the ABINIT
70
!! kpt_flag(ikpt) gives the indices of the k-point
71
!! related to ikpt by time revers
72
!! kpt_mark(nkpt)= 0, if k point is unmarked; 1, if k point has been marked
73
!! maxband/minband= control the minimum and maximum band calculated in the
75
!! npw_k= npwarr(ikpt), number of planewaves in basis at this k point
76
!! shift_g_2(nkpt,nkpt)= .true. if the k point should be shifted by a G vector;
78
!! tr(2)=variable that changes k to -k
80
!! $c_g$ to $c_g^*$ when time-reversal symetrie is used
86
!! appdig,cgedi,cgefa,handle_ncerr,hdr_io,hdr_io_netcdf,hdr_skip
87
!! ini_wf_netcdf,leave_new,matr3inv,rwwf,waveformat,wffclose,wffopen
88
!! wrtout,xdefineoff,zgedi,zgefa
92
#if defined HAVE_CONFIG_H
96
subroutine uderiv(bdberry,cg,gprimd,hdr,istwfk,kberry,kg,kpt_,kptopt,kptrlatt,&
97
& mband,mgfft,mkmem,mpi_enreg,mpw,natom,nband,nberry,npwarr,nspinor,nsppol,nkpt_,&
98
& rprimd,unddk,unkg,wffnow,filnam)
102
#if defined HAVE_NETCDF
106
!This section has been created automatically by the script Abilint (TD). Do not modify these by hand.
107
#ifdef HAVE_FORTRAN_INTERFACES
108
use interfaces_01manage_mpi
109
use interfaces_11util
110
use interfaces_13io_mpi
111
use interfaces_13ionetcdf
112
use interfaces_15common, except_this_one => uderiv
113
use interfaces_lib01hidempi
117
!End of the abilint section
121
!Arguments ------------------------------------
123
integer,intent(in) :: kptopt,mband,mgfft,mkmem,mpw,natom,nberry,nkpt_,nspinor
124
integer,intent(in) :: nsppol,unddk,unkg
125
type(MPI_type),intent(inout) :: mpi_enreg
126
type(hdr_type),intent(inout) :: hdr
127
type(wffile_type),intent(inout) :: wffnow
129
integer,intent(in) :: bdberry(4),istwfk(nkpt_),kberry(3,20),kg(3,mpw*mkmem)
130
integer,intent(in) :: kptrlatt(3,3),nband(nkpt_*nsppol),npwarr(nkpt_)
131
real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol),gprimd(1:3,1:3)
132
real(dp),intent(in) :: kpt_(3,nkpt_),rprimd(3,3)
133
character(len=fnlen),intent(in) :: filnam(5)
135
!Local variables -------------------------
137
integer :: accesswff,band_in,cg_index_iband,cg_index_jband,fform,flag,flag1
138
integer :: formeig,iband,iberry,icg,idir,ierr,ifor,ii,iii,ikpt,ikpt2,ikpt_
139
integer :: index,index1,info,ipert,ipw,isgn,isp,isppol,jband,jj,jkpt,jkpt_,jpw
140
integer :: jsppol,lkpt,master,maxband,mcg_disk,me,minband,nband_diff,nband_k
141
integer :: ncerr,ncid_hdr,nkpt,npw_k,option,pertcase,rdwr,read_k,spaceComm
143
real(dp) :: gmod,twodk
144
character(len=500) :: message
145
character(len=fnlen) :: fiwf1o,wff2nm
146
type(wffile_type) :: wffddk
148
integer :: kpt_flag(2*nkpt_),kpt_mark(2*nkpt_)
149
integer,allocatable :: cg_index(:,:,:),ikpt_dk(:,:),ipvt(:),kg_dum(:,:)
150
integer,allocatable :: kg_jl(:,:,:),kg_kpt(:,:,:)
151
real(dp) :: det(2,2),diffk(3),diffk2(3),dk(3),gpard(3),klattice(3,3)
152
real(dp) :: kptrlattr(3,3),tr(2)
153
real(dp),allocatable :: cg_disk(:,:,:),cg_disk_1(:,:),cmatrix(:,:,:),dudk(:,:)
154
real(dp),allocatable :: eig_dum(:),eig_dum_2(:),kpt(:,:),occ_dum(:)
155
real(dp),allocatable :: occ_dum_2(:),phi(:,:,:),u_tilde(:,:,:,:),zgwork(:,:)
156
logical,allocatable :: shift_g_2(:,:)
158
! *************************************************************************
160
!write(6,*)' uderiv : enter '
164
write(6,*)' uderiv : does not yet work for nspinor=2'
168
if(maxval(istwfk(:))/=1)then
169
write(message, '(a,a,a,a,a,a)' )ch10,&
170
& ' berryphase : BUG -',ch10,&
171
& ' Sorry, this routine does not work yet with istwfk/=1.',ch10,&
172
& ' This should have been tested previously ...'
173
call wrtout(6,message,'COLL')
174
call leave_new('COLL')
179
allocate(kpt(3,nkpt))
181
else if (kptopt==2) then
183
allocate(kpt(3,nkpt))
186
kpt(:,ikpt)=kpt_(:,ikpt)
189
do ikpt = (nkpt/2+1),nkpt
192
if (((abs(kpt_(1,ikpt-nkpt/2)+kpt_(1,jkpt))<1.0d-8).or.&
193
& (abs(1-abs(kpt_(1,ikpt-nkpt/2)+kpt_(1,jkpt)))<1.0d-8))&
194
& .and.((abs(kpt_(2,ikpt-nkpt/2)+kpt_(2,jkpt))<1.0d-8).or.&
195
& (abs(1-abs(kpt_(2,ikpt-nkpt/2)+kpt_(2,jkpt)))<1.0d-8))&
196
& .and.((abs(kpt_(3,ikpt-nkpt/2)+kpt_(3,jkpt))<1.0d-8).or.&
197
& (abs(1-abs(kpt_(3,ikpt-nkpt/2)+kpt_(3,jkpt)))<1.0d-8))) then
204
kpt_flag(ikpt-index)=ikpt-nkpt/2
205
kpt(:,ikpt-index)=-kpt_(:,ikpt-nkpt/2)
212
!write(101,*) 'beginning write kpt'
214
!write(101,*) kpt(:,ikpt)
218
allocate(shift_g_2(nkpt,nkpt))
220
!Compute primitive vectors of the k point lattice
222
kptrlattr(:,:)=kptrlatt(:,:)
223
!Go to reciprocal space (in reduced coordinates)
224
call matr3inv(kptrlattr,klattice)
228
!**************************************************************************
229
! Determine the appended index for ddk 1WF files
232
if (kberry(idir,iberry) ==1) then
234
pertcase=idir+(ipert-1)*3
240
wff2nm=trim(filnam(4))//'_1WF'
242
call appdig(pertcase,wff2nm,fiwf1o)
244
#if defined HAVE_NETCDF
245
if(accesswff==2) then
246
! Create empty netCDF file
247
ncerr = nf90_create(path=wff2nm, cmode=NF90_CLOBBER, ncid=ncid_hdr)
248
call handle_ncerr(ncerr," create netcdf wavefunction file")
249
ncerr = nf90_close(ncid_hdr)
250
call handle_ncerr(ncerr," close netcdf wavefunction file")
251
else if(accesswff==3) then
252
write (std_out,*) "FIXME: ETSF I/O support in uderiv"
256
spaceComm=abinit_comm_serial ; me=0 ; master=0 ; accesswff=0
257
call WffOpen(accesswff,spaceComm,fiwf1o,ierr,wffddk,master,me,unddk)
260
if (wffddk%accesswff /= 2) then
261
call hdr_io(fform,hdr,rdwr,wffddk)
262
#if defined HAVE_NETCDF
263
else if (wffddk%accesswff == 2) then
264
call hdr_io_netcdf(fform,hdr,rdwr,wffddk)
266
call ini_wf_netcdf(mpw,wffddk%unwff,1)
267
else if (wffddk%accesswff==3) then
268
write (std_out,*) "FIXME: ETSF I/O support in uderiv"
272
! Define offsets, in case of MPI I/O
273
call xdefineOff(formeig,wffddk,mpi_enreg,nband,npwarr,nspinor,nsppol,nkpt_)
275
! *****************************************************************************
276
! Calculate dimensional recip lattice vector along which P is calculated
277
! dk = step to the nearest k point along that direction
278
! in reduced coordinates
280
dk(:)=dble(kberry(1,iberry))*klattice(:,1)+&
281
& dble(kberry(2,iberry))*klattice(:,2)+&
282
& dble(kberry(3,iberry))*klattice(:,3)
285
if (dk(idir)/=0) then
290
gpard(:)=dk(1)*gprimd(:,1)+dk(2)*gprimd(:,2)+dk(3)*gprimd(:,3)
291
gmod=sqrt(dot_product(gpard,gpard))
293
!******************************************************************************
294
! Select the k grid points along the direction to compute dudk
295
! dk = step to the nearest k point along that direction
297
! For each k point, find k_prim such that k_prim= k + dk mod(G)
298
! where G is a vector of the reciprocal lattice
299
allocate(ikpt_dk(2,nkpt))
300
ikpt_dk(1:2,1:nkpt)=0
301
shift_g_2(:,:)= .false.
305
diffk(:)=abs(kpt(:,ikpt2)-kpt(:,ikpt)-dk(:))
306
diffk2(:)=abs(kpt(:,ikpt2)-kpt(:,ikpt)+dk(:))
307
if (sum(abs(diffk(:)-nint(diffk(:))))<3*tol8)then
308
ikpt_dk(1,ikpt)=ikpt2
309
if(sum(diffk(:))>=3*tol8) shift_g_2(ikpt,ikpt2) = .true.
311
if (sum(abs(diffk2(:)-nint(diffk2(:))))<3*tol8)then
312
ikpt_dk(2,ikpt)=ikpt2
313
if(sum(diffk2(:))>=3*tol8) shift_g_2(ikpt,ikpt2) = .true.
318
write(message,'(a,a,a,3f9.5,a,a,3f9.5,a)')ch10,&
319
& ' Computing the derivative for reciprocal vector:',ch10,&
320
& dk(:),' (in reduced coordinates)',ch10,&
321
& gpard(1:3),' (in cartesian coordinates - atomic units)'
322
call wrtout(ab_out,message,'COLL')
323
call wrtout(6,message,'COLL')
326
write(message, '(a,i5,a,i5)')&
327
& ' From band number',bdberry(1),' to band number',bdberry(2)
329
write(message, '(a,i5,a,i5,a,a,a,i5,a,i5,a)')&
330
& ' From band number',bdberry(1),' to band number',bdberry(2),' for spin up,',&
332
& ' from band number',bdberry(3),' to band number',bdberry(4),' for spin down.'
334
call wrtout(ab_out,message,'COLL')
335
call wrtout(6,message,'COLL')
337
!*****************************************************************************
338
allocate(dudk(2,mpw*nspinor*mband*nsppol))
339
allocate(eig_dum_2((2*mband)**formeig*mband))
340
allocate(occ_dum_2((2*mband)**formeig*mband))
347
! Find the location of each wavefunction
349
allocate(cg_index(mband,nkpt_,nsppol))
353
nband_k=nband(ikpt+(isppol-1)*nkpt_)
356
cg_index(iband,ikpt,isppol)=(iband-1)*npw_k*nspinor+icg
358
icg=icg+npw_k*nspinor*nband(ikpt)
362
! Find the planewave vectors for each k point
363
! SHOULD BE REMOVED WHEN ANOTHER INDEXING TECHNIQUE WILL BE USED FOR kg
364
allocate(kg_kpt(3,mpw*nspinor,nkpt_))
369
do ipw=1,npw_k*nspinor
370
kg_kpt(1:3,ipw,ikpt)=kg(1:3,ipw+index1)
372
index1=index1+npw_k*nspinor
377
call hdr_skip(wffnow,ierr)
379
! Should Define offsets, in case of MPI I/O (Also, other hdr_skip calls !!)
380
! For the time being does not work with MPI I/O
381
! call xdefineOff(formeig,wffnow,mpi_enreg,nband,npwarr,nspinor,nsppol,nkpt_)
383
allocate(cg_disk(2,mpw*nspinor*mband,2))
385
allocate(kg_jl(3,mpw,2))
389
!*************************************************************************
393
minband=bdberry(2*isppol-1)
394
maxband=bdberry(2*isppol)
397
write(message,'(a,a,a,a,i5,a)')ch10,&
398
& ' berryphase : BUG - ',ch10,&
399
& ' The band limit minband=',minband,', is lower than 0.'
400
call wrtout(6,message,'COLL')
401
call leave_new('COLL')
404
write(message,'(a,a,a,a,i5,a)')ch10,&
405
& ' berryphase : BUG - ',ch10,&
406
& ' The band limit maxband=',maxband,', is lower than 0.'
407
call wrtout(6,message,'COLL')
408
call leave_new('COLL')
410
if(maxband<minband)then
411
write(message,'(4a,i5,a,i5)')ch10,&
412
& ' berryphase : BUG - ',ch10,&
413
& ' maxband=',maxband,', is lower than minband=',minband
414
call wrtout(6,message,'COLL')
415
call leave_new('COLL')
427
if (read_k == 0) then
428
if (kpt_flag(ikpt_)/=0) then
430
ikpt= kpt_flag(ikpt_)
433
if (kpt_flag(ikpt_)/=0) then
434
tr(-1*read_k+3) = -1.0_dp
435
ikpt= kpt_flag(ikpt_)
440
nband_k=nband(ikpt+(isppol-1)*nkpt_)
442
if(nband_k<maxband)then
443
write(message,'(4a,i5,a,i5)')ch10,&
444
& ' uderiv : BUG - ',ch10,&
445
& ' maxband=',maxband,', is larger than nband(i,isppol)=',nband_k
446
call wrtout(6,message,'COLL')
447
call leave_new('COLL')
452
allocate(u_tilde(2,npw_k,maxband,2))
453
u_tilde(1:2,1:npw_k,1:maxband,1:2)=0.0_dp
455
! ifor = 1,2 represents forward and backward neighbouring k points of ikpt
456
! respectively along dk direction
460
allocate(phi(2,mpw,mband),cmatrix(2,maxband,maxband))
461
phi(1:2,1:mpw,1:mband)=0.0_dp; cmatrix(1:2,1:maxband,1:maxband)=0.0_dp
464
jkpt_= ikpt_dk(ifor,ikpt_)
471
if (read_k == 0) then
472
if (kpt_flag(jkpt_)/=0) then
474
jkpt= kpt_flag(jkpt_)
477
if (kpt_flag(jkpt_)/=0) then
479
jkpt= kpt_flag(jkpt_)
486
! if read_k = 0,read first k point of string
487
! ******************************************
493
do while(index < ikpt)
502
read(unkg) kg_jl(1:3,1:npw_k,read_k)
505
allocate(eig_dum(mband),occ_dum(mband))
506
mcg_disk=mpw*nspinor*mband
507
allocate(cg_disk_1(2,mcg_disk),kg_dum(3,0))
509
call hdr_skip(wffnow,ierr)
512
if(isp==isppol .and. lkpt==ikpt)then
513
option=-2 ! will read cg only
515
option=-1 ! will skip
517
call rwwf(cg_disk_1,eig_dum,0,0,0,lkpt,isppol,kg_dum,mband,mcg_disk,nband_k, &
518
& nband_k,npw_k,nspinor,occ_dum,option,0,tim_rwwf,wffnow)
522
cg_disk(1:2,:,read_k)=cg_disk_1(1:2,:)
523
deallocate(cg_disk_1,eig_dum,kg_dum,occ_dum)
529
! read the next k point
530
! *********************
531
if (read_k /= 0) then
535
do while(index < jkpt)
544
read(unkg) kg_jl(1:3,1:npw_k,read_k)
547
allocate(eig_dum(mband),occ_dum(mband))
548
mcg_disk=mpw*nspinor*mband
549
allocate(cg_disk_1(2,mcg_disk))
551
call hdr_skip(wffnow,ierr)
554
if(isp==isppol .and. lkpt==jkpt)then
555
option=-2 ! will read cg only
557
option=-1 ! will skip
559
call rwwf(cg_disk_1,eig_dum,0,0,0,lkpt,isppol,kg_dum,mband,mcg_disk,nband_k, &
560
& nband_k,npw_k,nspinor,occ_dum,option,0,tim_rwwf,wffnow)
564
cg_disk(1:2,:,read_k)=cg_disk_1(1:2,:)
565
deallocate(cg_disk_1,eig_dum,kg_dum,occ_dum)
570
if (ifor==1) read_k = 2
575
call waveformat(cg,cg_disk,cg_index,phi,dk,ii,ikpt,&
576
& ikpt_,isgn,isppol,jj,jkpt,jkpt_,kg_kpt,kpt,kg_jl,maxband,mband,&
577
& minband,mkmem,mpw,nkpt,nkpt_,npwarr,nsppol,nspinor,shift_g_2,tr)
579
! Compute the overlap matrix <u_k|u_k+b>
583
do iband=minband,maxband
584
cg_index_iband= (iband-1)*npwarr(ikpt)
585
do jband=minband,maxband
586
do ipw=1,npwarr(ikpt)
587
cmatrix(1,iband,jband)=cmatrix(1,iband,jband)+&
588
& cg_disk(1,ipw+cg_index_iband,ii)*phi(1,ipw,jband)+&
589
& tr(ii)*cg_disk(2,ipw+cg_index_iband,ii)*tr(jj)*phi(2,ipw,jband)
591
cmatrix(2,iband,jband)=cmatrix(2,iband,jband)+&
592
& cg_disk(1,ipw+cg_index_iband,ii)*tr(jj)*phi(2,ipw,jband)-&
593
& tr(ii)*cg_disk(2,ipw+cg_index_iband,ii)*phi(1,ipw,jband)
600
do iband=minband,maxband
601
cg_index_iband=cg_index(iband,ikpt,isppol)
602
do jband=minband,maxband
603
do ipw=1,npwarr(ikpt)
604
cmatrix(1,iband,jband)=cmatrix(1,iband,jband)+&
605
& cg(1,ipw+cg_index_iband)*phi(1,ipw,jband)+&
606
& tr(ii)*cg(2,ipw+cg_index_iband)*tr(jj)*phi(2,ipw,jband)
608
cmatrix(2,iband,jband)=cmatrix(2,iband,jband)+&
609
& cg(1,ipw+cg_index_iband)*tr(jj)*phi(2,ipw,jband)-&
610
& tr(ii)*cg(2,ipw+cg_index_iband)*phi(1,ipw,jband)
616
! Compute the inverse of cmatrix(1:2,minband:maxband, minband:maxband)
618
band_in = maxband - minband + 1
619
allocate(ipvt(maxband))
620
allocate(zgwork(2,1:maxband))
622
! Last argument of zgedi means calculate inverse only
624
call cgefa(cmatrix(1,minband,minband),maxband, band_in,ipvt,info)
625
call cgedi(cmatrix(1,minband,minband),maxband, band_in,ipvt,det,zgwork,01)
627
call zgefa(cmatrix(1,minband,minband),maxband, band_in,ipvt,info)
628
call zgedi(cmatrix(1,minband,minband),maxband, band_in,ipvt,det,zgwork,01)
631
deallocate(zgwork,ipvt)
633
! Compute the product of Inverse overlap matrix with the wavefunction
635
do iband=minband,maxband
636
do ipw=1,npwarr(ikpt)
637
u_tilde(1,ipw,iband,ifor)= &
638
& dot_product(cmatrix(1,minband:maxband,iband),&
639
& phi(1,ipw,minband:maxband))-&
640
& dot_product(cmatrix(2,minband:maxband,iband),&
641
& tr(jj)*phi(2,ipw,minband:maxband))
642
u_tilde(2,ipw,iband,ifor)= &
643
& dot_product(cmatrix(1,minband:maxband,iband),&
644
& tr(jj)*phi(2,ipw,minband:maxband))+&
645
& dot_product(cmatrix(2,minband:maxband,iband),&
646
& phi(1,ipw,minband:maxband))
649
deallocate(cmatrix,phi)
653
! Compute dudk for ikpt
657
do iband=minband,maxband
659
icg=(iband-minband)*npw_k
661
dudk(1,1+icg:npw_k+icg)=(u_tilde(1,1:npw_k,iband,1)-&
662
& u_tilde(1,1:npw_k,iband,2))/twodk
664
dudk(2,1+icg:npw_k+icg)=(u_tilde(2,1:npw_k,iband,1)-&
665
& u_tilde(2,1:npw_k,iband,2))/twodk
670
mcg_disk=mpw*nspinor*mband
671
nband_diff=maxband-minband+1
672
call rwwf(dudk,eig_dum_2,formeig,0,0,ikpt,isppol,kg_kpt(:,:,ikpt),&
673
& mband,mcg_disk,nband_diff,nband_diff,&
674
& npw_k,nspinor,occ_dum_2,2,1,tim_rwwf,wffddk)
681
deallocate(eig_dum_2,occ_dum_2)
684
call WffClose(wffddk,ierr)
687
deallocate(cg_disk,kg_jl)
689
deallocate(kg_kpt,cg_index)
695
deallocate(shift_g_2,kpt)
697
write(6,*) 'uderiv: exit '
699
end subroutine uderiv
added
removed
src/15common/uderiv.F90
