TABLE OF CONTENTS
ABINIT/dfptnl_mv [ Functions ]
NAME
dfptnl_mv
FUNCTION
Compute the finite difference expression of the k-point derivative using the PEAD formulation of the third-order energy (see Nunes and Gonze PRB 63, 155107 (2001) Eq. 102) and the finite difference formula of Marzari and Vanderbilt (see Marzari and Vanderbilt, PRB 56, 12847 (1997), Appendix B)
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (MVeithen) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
cg(2,mpw*nspinor*mband*mkmem*nsppol) = array for planewave
coefficients of wavefunctions
cgindex(nkpt2,nsppol) = for each k-point, cgindex stores the location
of the WF in the cg array
cg1 = first-order wavefunction relative to the perturbations i1pert
cg3 = first-order wavefunction relative to the perturbations i3pert
dtfil <type(datafiles_type)>=variables related to files
dtset <type(dataset_type)>=all input variables for this dataset
gmet(3,3)=reciprocal space metric tensor in bohr**-2
i1pert,i3pert = type of perturbation that has to be computed
i1dir,i3dir=directions of the corresponding perturbations
kneigh(30,nkpt2) = index of the neighbours of each k-point
kg_neigh(30,nkpt2,3) = necessary to construct the vector joining a k-point
to its nearest neighbour in case of a single k-point,
a line of k-points or a plane of k-points.
See getshell.F90 for details
kptindex(2,nkpt3)= index of the k-points in the reduced BZ
related to a k-point in the full BZ
kpt3(3,nkpt3) = reduced coordinates of k-points in the full BZ
mband = maximum number of bands
mkmem = maximum number of k points which can fit in core memory
mkmem_max = maximal number of k-points on each processor (MPI //)
mk1mem = maximum number of k points for first-order WF
which can fit in core memory
mpi_enreg=MPI-parallelisation information
mpw = maximum number of planewaves in basis sphere (large number)
mvwtk(30,nkpt) = weights to compute the finite difference ddk
natom = number of atoms in unit cell
nkpt2 = number of k-points in the reduced part of the BZ
nkpt2 = nkpt/2 in case of time-reversal symmetry (kptopt = 2)
nkpt3 = number of k-points in the full BZ
nneigh = total number of neighbours required to evaluate the finite
difference formula
npwarr(nkpt) = array holding npw for each k point
nspinor = number of spinorial components of the wavefunctions
nsppol = number of channels for spin-polarization (1 or 2)
pwind(mpw,nneigh,mkmem) = array used to compute the overlap matrix smat
between k-points
OUTPUT
d3_berry(2,3) = Berry-phase part of the third-order energy
SIDE EFFECTS
mpi_enreg=MPI-parallelisation information
NOTES
For a given set of values of i1pert,i3pert,i1dir and i3dir, the routine computes the k-point derivatives for 12dir = 1,2,3
TODO
PARENTS
dfptnl_loop
CHILDREN
dzgedi,dzgefa,mpi_recv,mpi_send,status,wrtout,xmpi_sum
SOURCE
81 #if defined HAVE_CONFIG_H 82 #include "config.h" 83 #endif 84 85 #include "abi_common.h" 86 87 88 subroutine dfptnl_mv(cg,cgindex,cg1,cg3,dtset,dtfil,d3_berry,gmet,& 89 & i1pert,i3pert,i1dir,i3dir,kneigh,kg_neigh,kptindex,& 90 & kpt3,mband,mkmem,mkmem_max,mk1mem,mpi_enreg,& 91 & mpw,mvwtk,natom,nkpt2,nkpt3,nneigh,npwarr,nspinor,& 92 & nsppol,pwind) 93 94 use m_profiling_abi 95 96 use defs_basis 97 use defs_abitypes 98 use m_xmpi 99 #if defined HAVE_MPI2 100 use mpi 101 #endif 102 103 use m_abilasi, only : dzgedi, dzgefa 104 105 !This section has been created automatically by the script Abilint (TD). 106 !Do not modify the following lines by hand. 107 #undef ABI_FUNC 108 #define ABI_FUNC 'dfptnl_mv' 109 use interfaces_14_hidewrite 110 use interfaces_32_util 111 !End of the abilint section 112 113 implicit none 114 115 #if defined HAVE_MPI1 116 include 'mpif.h' 117 #endif 118 119 !Arguments ------------------------------------ 120 ! 121 !--- Arguments : integer scalars 122 integer, intent(in) :: i1dir,i1pert,i3dir,i3pert,mband,mk1mem 123 integer, intent(in) :: mkmem,mkmem_max,mpw,natom 124 integer, intent(in) :: nkpt2,nkpt3,nneigh,nspinor,nsppol 125 ! 126 !--- Arguments : integer arrays 127 integer, intent(in) :: cgindex(nkpt2,nsppol) 128 integer, intent(in) :: kneigh(30,nkpt2),kg_neigh(30,nkpt2,3),kptindex(2,nkpt3) 129 integer, intent(in) :: npwarr(nkpt2),pwind(mpw,nneigh,mkmem) 130 ! 131 !--- Arguments : real(dp) scalars 132 ! 133 !--- Arguments : real(dp) arrays 134 real(dp), intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol) 135 real(dp), intent(in) :: cg1(2,mpw*nspinor*mband*mk1mem*nsppol) 136 real(dp), intent(in) :: cg3(2,mpw*nspinor*mband*mk1mem*nsppol) 137 real(dp), intent(in) :: gmet(3,3),kpt3(3,nkpt3) 138 real(dp), intent(in) :: mvwtk(30,nkpt2) 139 real(dp), intent(out) :: d3_berry(2,3) 140 ! 141 !--- Arguments : structured datatypes 142 type(MPI_type), intent(in) :: mpi_enreg 143 type(datafiles_type), intent(in) :: dtfil 144 type(dataset_type), intent(in) :: dtset 145 146 !Local variables------------------------------- 147 ! 148 !---- Local variables : integer scalars 149 integer :: count,counter,count1,iband,icg 150 integer :: ierr,iexit,ii,ikpt,ikpt_loc,ikpt2 151 integer :: ikpt_rbz,ineigh,info,ipw,isppol,jband,jcg,jj,jkpt,job,jpw, jkpt2, jkpt_rbz 152 integer :: lband,lpband,nband_occ,npw_k,npw_k1,my_source,his_source,dest,tag 153 integer :: spaceComm 154 integer,parameter :: level=52 155 integer :: bdtot_index 156 ! 157 !---- Local variables : integer arrays 158 integer,allocatable :: ipvt(:) 159 integer, allocatable :: bd_index(:,:) 160 ! 161 !---- Local variables : real(dp) scalars 162 real(dp) :: dotnegi,dotnegr,dotposi,dotposr 163 ! real(dp) :: c1,c2 ! appear commented out below 164 ! 165 !---- Local variables : real(dp) arrays 166 real(dp) :: d3_aux(2,3),det(2,2),dk(3),dk_(3) 167 real(dp) :: z1(2),z2(2) 168 real(dp),allocatable :: buffer(:,:),cgq(:,:),cg1q(:,:),cg3q(:,:) 169 real(dp),allocatable :: qmat(:,:,:),s13mat(:,:,:),s1mat(:,:,:),s3mat(:,:,:) 170 real(dp),allocatable :: smat(:,:,:),zgwork(:,:) 171 ! 172 !---- Local variables : character variables 173 character(len=500) :: message 174 ! 175 !---- Local variables : structured datatypes 176 177 #if defined HAVE_MPI 178 integer :: status1(MPI_STATUS_SIZE) 179 !BEGIN TF_CHANGES 180 spaceComm=mpi_enreg%comm_cell 181 !END TF_CHANGES 182 #endif 183 184 185 ! *********************************************************************** 186 187 !DEBUG 188 !write(std_out,*)' dfptnl_mv : enter ' 189 !flush(6) 190 !stop 191 !ENDDEBUG 192 193 194 call status(0,dtfil%filstat,iexit,level,'enter ') 195 196 write(message,'(8a)') ch10,& 197 & ' dfptnl_mv : finite difference expression of the k-point derivative',ch10,& 198 & ' is performed using the PEAD formulation of ',& 199 & 'the third-order energy',ch10,& 200 & ' (see Nunes and Gonze PRB 63, 155107 (2001) Eq. 102)',ch10 201 !call wrtout(ab_out,message,'COLL') 202 call wrtout(std_out, message,'COLL') 203 204 205 !fab: I think that the following restriction must be eliminated: 206 !isppol = 1 207 208 ikpt_loc = 0 209 d3_aux(:,:) = 0_dp 210 211 ABI_ALLOCATE(s13mat,(2,mband,mband)) 212 ABI_ALLOCATE(smat,(2,mband,mband)) 213 ABI_ALLOCATE(s1mat,(2,mband,mband)) 214 ABI_ALLOCATE(qmat,(2,mband,mband)) 215 ABI_ALLOCATE(ipvt,(mband)) 216 ABI_ALLOCATE(s3mat,(2,mband,mband)) 217 ABI_ALLOCATE(zgwork,(2,mband)) 218 ABI_ALLOCATE(bd_index, (nkpt2, nsppol)) 219 220 bdtot_index = 0 221 do isppol = 1, nsppol 222 do ikpt_rbz = 1, nkpt2 223 bd_index(ikpt_rbz,isppol) = bdtot_index 224 bdtot_index = bdtot_index + dtset%nband(ikpt_rbz+nkpt2*(isppol-1)) 225 end do 226 end do 227 228 !fab: I think here I have to add the loop over spin 229 230 do isppol = 1, nsppol 231 232 ! Loop over k-points 233 ! COMMENT: Every processor has to make mkmem_max iterations 234 ! even if mkmem < mkemem_max. This is due to the fact 235 ! that it still has to communicate its wavefunctions 236 ! to other processors even if it has no more overlap 237 ! matrices to compute. 238 239 ikpt_loc = 0 ; ikpt = 0 240 241 do while (ikpt_loc < mkmem_max) 242 243 if (ikpt_loc < mkmem) ikpt = ikpt + 1 244 245 if (xmpi_paral == 1) then 246 ! if ((minval(abs(mpi_enreg%proc_distrb(ikpt,1:mband,1:dtset%nsppol) & 247 ! & - mpi_enreg%me)) /= 0).and.(ikpt_loc < mkmem)) cycle 248 if(ikpt>nkpt2)then 249 ikpt_loc=mkmem_max 250 cycle 251 end if 252 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,-1,mpi_enreg%me)) then 253 if(ikpt==nkpt2) ikpt_loc=mkmem_max 254 cycle 255 end if 256 end if 257 258 ikpt_loc = ikpt_loc + 1 259 npw_k = npwarr(ikpt) 260 counter = 100*ikpt 261 call status(counter,dtfil%filstat,iexit,level,'loop over k ') 262 263 ii = cgindex(ikpt,isppol) 264 265 ! Loop on the neighbours 266 267 do ineigh = 1,nneigh 268 269 s13mat(:,:,:) = zero 270 smat(:,:,:) = zero 271 s1mat(:,:,:) = zero 272 s3mat(:,:,:) = zero 273 qmat(:,:,:) = zero 274 275 ikpt2 = kneigh(ineigh,ikpt) 276 ikpt_rbz = kptindex(1,ikpt2) ! index of the k-point in the reduced BZ 277 jj = cgindex(ikpt_rbz,isppol) 278 ! previous fixed value for nband_k now called nband_occ: 279 !nband_occ = dtset%nband(ikpt_rbz+nkpt2*(isppol-1)) 280 ! TODO: check if all these bands are occupied in nsppol = 2 case 281 nband_occ = 0 282 do iband = 1, dtset%nband(ikpt_rbz+nkpt2*(isppol-1)) 283 if (dtset%occ_orig(bd_index(ikpt_rbz,isppol) + iband) > tol10) nband_occ = nband_occ + 1 284 end do 285 npw_k1 = npwarr(ikpt_rbz) 286 dk_(:) = kpt3(:,ikpt2) - dtset%kptns(:,ikpt) 287 dk(:) = dk_(:) - nint(dk_(:)) + real(kg_neigh(ineigh,ikpt,:),dp) 288 289 count = nspinor*mband*npw_k1 290 ABI_ALLOCATE(cgq,(2,count)) 291 ABI_ALLOCATE(cg1q,(2,count)) 292 ABI_ALLOCATE(cg3q,(2,count)) 293 294 #if defined HAVE_MPI 295 296 my_source = mpi_enreg%proc_distrb(ikpt_rbz,1,1) 297 298 ! do dest = 0, mpi_enreg%nproc-1 299 do dest = 0, maxval(mpi_enreg%proc_distrb(1:nkpt2,1:mband,1:dtset%nsppol)) 300 301 if ((dest==mpi_enreg%me).and.(ikpt_loc <= mkmem)) then 302 ! I am dest and have something to do 303 304 if (my_source == mpi_enreg%me) then 305 ! I am destination and source 306 jcg = cgindex(ikpt_rbz,isppol) 307 308 cgq(:,1:count) = cg(:,jcg+1:jcg+count) 309 cg1q(:,1:count) = cg1(:,jcg+1:jcg+count) 310 cg3q(:,1:count) = cg3(:,jcg+1:jcg+count) 311 312 else 313 ! I am the destination but not the source -> receive 314 315 tag = ikpt_rbz 316 317 ABI_ALLOCATE(buffer,(2,3*count)) 318 319 call MPI_RECV(buffer,2*3*count,MPI_DOUBLE_PRECISION,my_source,tag,spaceComm,status1,ierr) 320 321 cgq(:,1:count) = buffer(:,1:count) 322 cg1q(:,1:count) = buffer(:,count+1:2*count) 323 cg3q(:,1:count) = buffer(:,2*count+1:3*count) 324 ABI_DEALLOCATE(buffer) 325 326 end if 327 328 else if (ikpt_loc <= mpi_enreg%mkmem(dest)) then ! dest != me and the dest has a k-point to treat 329 330 jkpt=mpi_enreg%kpt_loc2ibz_sp(dest, ikpt_loc,1) 331 jkpt2 = kneigh(ineigh,jkpt) 332 jkpt_rbz = kptindex(1,jkpt2) ! index of the k-point in the reduced BZ 333 334 his_source = mpi_enreg%proc_distrb(jkpt_rbz,1,1) 335 336 if (his_source == mpi_enreg%me) then 337 338 jcg = cgindex(jkpt_rbz,isppol) 339 340 tag = jkpt_rbz 341 count1 = npwarr(jkpt_rbz)*mband*nspinor 342 ABI_ALLOCATE(buffer,(2,3*count1)) 343 buffer(:,1:count1) = cg(:,jcg+1:jcg+count1) 344 buffer(:,count1+1:2*count1) = cg1(:,jcg+1:jcg+count1) 345 buffer(:,2*count1+1:3*count1) = cg3(:,jcg+1:jcg+count1) 346 347 call MPI_SEND(buffer,2*3*count1,MPI_DOUBLE_PRECISION,dest,tag,spaceComm,ierr) 348 349 ABI_DEALLOCATE(buffer) 350 351 end if 352 353 end if 354 355 end do 356 ! 357 ! do jkpt = 1, nkpt2 358 ! 359 ! if ((jkpt == ikpt_rbz).and.(source /= mpi_enreg%me).and.& 360 ! & (ikpt_loc <= mkmem)) then 361 ! 362 ! tag = jkpt 363 ! 364 ! allocate(buffer(2,3*count)) 365 ! call MPI_RECV(buffer,2*3*count,MPI_DOUBLE_PRECISION,& 366 ! source,tag,spaceComm,status1,ierr) 367 ! 368 ! cgq(:,1:count) = buffer(:,1:count) 369 ! cg1q(:,1:count) = buffer(:,count+1:2*count) 370 ! cg3q(:,1:count) = buffer(:,2*count+1:3*count) 371 ! deallocate(buffer) 372 ! 373 ! end if 374 ! 375 ! ! ---------------------------------------------------------------------------- 376 ! ! --------------- Here: send the WF to all the cpus that need it ------------- 377 ! ! ---------------------------------------------------------------------------- 378 ! 379 ! do dest = 1, mpi_enreg%nproc 380 ! 381 ! if ((minval(abs(mpi_enreg%proc_distrb(jkpt,1:mband,1:dtset%nsppol) & 382 ! & - mpi_enreg%me)) == 0).and.& 383 ! & (mpi_enreg%kptdstrb(dest,ineigh,ikpt_loc) == jkpt)) then 384 ! 385 ! 386 ! 387 ! jcg = cgindex(jkpt,isppol) 388 ! 389 ! if (((dest-1) == mpi_enreg%me)) then 390 ! 391 ! cgq(:,1:count) = cg(:,jcg+1:jcg+count) 392 ! cg1q(:,1:count) = cg1(:,jcg+1:jcg+count) 393 ! cg3q(:,1:count) = cg3(:,jcg+1:jcg+count) 394 ! 395 ! else 396 ! 397 ! tag = jkpt 398 ! count1 = npwarr(jkpt)*mband*nspinor 399 ! allocate(buffer(2,3*count1)) 400 ! buffer(:,1:count1) = cg(:,jcg+1:jcg+count1) 401 ! buffer(:,count1+1:2*count1) = cg1(:,jcg+1:jcg+count1) 402 ! buffer(:,2*count1+1:3*count1) = cg3(:,jcg+1:jcg+count1) 403 ! 404 ! call MPI_SEND(buffer,2*3*count1,MPI_DOUBLE_PRECISION,(dest-1),tag,spaceComm,ierr) 405 ! 406 ! deallocate(buffer) 407 ! 408 ! end if 409 ! 410 ! end if 411 ! 412 ! end do ! loop over dest 413 ! 414 ! end do ! loop over jkpt 415 416 if (ikpt_loc > mkmem) then 417 ABI_DEALLOCATE(cgq) 418 ABI_DEALLOCATE(cg1q) 419 ABI_DEALLOCATE(cg3q) 420 cycle 421 end if 422 423 #else 424 ! no // over k-points 425 426 cgq(:,1:count) = cg(:,jj+1:jj+count) 427 cg1q(:,1:count) = cg1(:,jj+1:jj+count) 428 cg3q(:,1:count) = cg3(:,jj+1:jj+count) 429 430 #endif 431 432 ! Compute overlap matrices 433 434 if (kptindex(2,ikpt2) == 0) then ! no time-reversal symmetry 435 436 do ipw = 1, npw_k 437 438 jpw = pwind(ipw,ineigh,ikpt_loc) 439 if (jpw /= 0) then 440 441 do iband = 1, nband_occ 442 do jband = 1, nband_occ 443 444 icg = ii + (iband-1)*npw_k + ipw 445 jcg = (jband-1)*npw_k1 + jpw 446 447 smat(1,iband,jband) = smat(1,iband,jband) + & 448 & cg(1,icg)*cgq(1,jcg) + cg(2,icg)*cgq(2,jcg) 449 smat(2,iband,jband) = smat(2,iband,jband) + & 450 & cg(1,icg)*cgq(2,jcg) - cg(2,icg)*cgq(1,jcg) 451 452 s13mat(1,iband,jband) = s13mat(1,iband,jband) + & 453 & cg1(1,icg)*cg3q(1,jcg) + cg1(2,icg)*cg3q(2,jcg) 454 s13mat(2,iband,jband) = s13mat(2,iband,jband) + & 455 & cg1(1,icg)*cg3q(2,jcg) - cg1(2,icg)*cg3q(1,jcg) 456 457 s1mat(1,iband,jband) = s1mat(1,iband,jband) + & 458 & cg1(1,icg)*cgq(1,jcg) + cg1(2,icg)*cgq(2,jcg) + & 459 & cg(1,icg)*cg1q(1,jcg) + cg(2,icg)*cg1q(2,jcg) 460 s1mat(2,iband,jband) = s1mat(2,iband,jband) + & 461 & cg1(1,icg)*cgq(2,jcg) - cg1(2,icg)*cgq(1,jcg) + & 462 & cg(1,icg)*cg1q(2,jcg) - cg(2,icg)*cg1q(1,jcg) 463 464 s3mat(1,iband,jband) = s3mat(1,iband,jband) + & 465 & cg3(1,icg)*cgq(1,jcg) + cg3(2,icg)*cgq(2,jcg) + & 466 & cg(1,icg)*cg3q(1,jcg) + cg(2,icg)*cg3q(2,jcg) 467 s3mat(2,iband,jband) = s3mat(2,iband,jband) + & 468 & cg3(1,icg)*cgq(2,jcg) - cg3(2,icg)*cgq(1,jcg) + & 469 & cg(1,icg)*cg3q(2,jcg) - cg(2,icg)*cg3q(1,jcg) 470 471 end do 472 end do 473 474 end if 475 476 end do ! ipw 477 478 else ! use time-reversal symmetry 479 480 do ipw = 1,npw_k 481 482 jpw = pwind(ipw,ineigh,ikpt_loc) 483 if (jpw /= 0) then 484 485 do iband = 1, nband_occ 486 do jband = 1, nband_occ 487 488 icg = ii + (iband-1)*npw_k + ipw 489 jcg = (jband-1)*npw_k1 + jpw 490 491 smat(1,iband,jband) = smat(1,iband,jband) + & 492 & cg(1,icg)*cgq(1,jcg) - cg(2,icg)*cgq(2,jcg) 493 smat(2,iband,jband) = smat(2,iband,jband) - & 494 & cg(1,icg)*cgq(2,jcg) - cg(2,icg)*cgq(1,jcg) 495 496 s13mat(1,iband,jband) = s13mat(1,iband,jband) + & 497 & cg1(1,icg)*cg3q(1,jcg) - cg1(2,icg)*cg3q(2,jcg) 498 s13mat(2,iband,jband) = s13mat(2,iband,jband) - & 499 & cg1(1,icg)*cg3q(2,jcg) - cg1(2,icg)*cg3q(1,jcg) 500 501 s1mat(1,iband,jband) = s1mat(1,iband,jband) + & 502 & cg1(1,icg)*cgq(1,jcg) - cg1(2,icg)*cgq(2,jcg) + & 503 & cg(1,icg)*cg1q(1,jcg) - cg(2,icg)*cg1q(2,jcg) 504 s1mat(2,iband,jband) = s1mat(2,iband,jband) - & 505 & cg1(1,icg)*cgq(2,jcg) - cg1(2,icg)*cgq(1,jcg) - & 506 & cg(1,icg)*cg1q(2,jcg) - cg(2,icg)*cg1q(1,jcg) 507 508 s3mat(1,iband,jband) = s3mat(1,iband,jband) + & 509 & cg3(1,icg)*cgq(1,jcg) - cg3(2,icg)*cgq(2,jcg) + & 510 & cg(1,icg)*cg3q(1,jcg) - cg(2,icg)*cg3q(2,jcg) 511 s3mat(2,iband,jband) = s3mat(2,iband,jband) - & 512 & cg3(1,icg)*cgq(2,jcg) - cg3(2,icg)*cgq(1,jcg) - & 513 & cg(1,icg)*cg3q(2,jcg) - cg(2,icg)*cg3q(1,jcg) 514 515 end do 516 end do 517 518 end if 519 520 end do ! ipw 521 522 end if 523 524 ABI_DEALLOCATE(cgq) 525 ABI_DEALLOCATE(cg1q) 526 ABI_DEALLOCATE(cg3q) 527 528 ! Compute qmat, the inverse of smat 529 530 job = 1 ! compute inverse only 531 qmat(:,:,:) = smat(:,:,:) 532 533 call dzgefa(qmat,mband,nband_occ,ipvt,info) 534 call dzgedi(qmat,mband,nband_occ,ipvt,det,zgwork,job) 535 536 ! DEBUG 537 ! write(100,*) 538 ! write(100,*)'ikpt = ',ikpt,'ineigh = ',ineigh 539 ! do iband = 1,nband_occ 540 ! do jband = 1,nband_occ 541 ! c1 = 0_dp ; c2 = 0_dp 542 ! do lband = 1,nband_occ 543 ! c1 = c1 + smat(1,iband,lband)*qmat(1,lband,jband) - & 544 ! & smat(2,iband,lband)*qmat(2,lband,jband) 545 ! c2 = c2 + smat(1,iband,lband)*qmat(2,lband,jband) + & 546 ! & smat(2,iband,lband)*qmat(1,lband,jband) 547 ! end do 548 ! write(100,'(2(2x,i2),2(2x,f16.9))')iband,jband,& 549 ! & c1,c2 550 ! end do 551 ! end do 552 ! ENDDEBUG 553 554 555 556 ! Accumulate sum over bands 557 558 dotposr = 0_dp ; dotposi = 0_dp 559 dotnegr = 0_dp ; dotnegi = 0_dp 560 do iband = 1, nband_occ 561 do jband = 1, nband_occ 562 563 dotposr = dotposr + & 564 & s13mat(1,iband,jband)*qmat(1,jband,iband) - & 565 & s13mat(2,iband,jband)*qmat(2,jband,iband) 566 dotposi = dotposi + & 567 & s13mat(1,iband,jband)*qmat(2,jband,iband) + & 568 & s13mat(2,iband,jband)*qmat(1,jband,iband) 569 570 571 do lband = 1, nband_occ 572 do lpband= 1, nband_occ 573 574 z1(1) = s1mat(1,iband,jband)*qmat(1,jband,lband) - & 575 & s1mat(2,iband,jband)*qmat(2,jband,lband) 576 z1(2) = s1mat(1,iband,jband)*qmat(2,jband,lband) + & 577 & s1mat(2,iband,jband)*qmat(1,jband,lband) 578 579 z2(1) = s3mat(1,lband,lpband)*qmat(1,lpband,iband) - & 580 & s3mat(2,lband,lpband)*qmat(2,lpband,iband) 581 z2(2) = s3mat(1,lband,lpband)*qmat(2,lpband,iband) + & 582 & s3mat(2,lband,lpband)*qmat(1,lpband,iband) 583 584 dotnegr = dotnegr + & 585 & z1(1)*z2(1) - z1(2)*z2(2) 586 dotnegi = dotnegi + & 587 & z1(1)*z2(2) + z1(2)*z2(1) 588 589 end do ! lpband 590 end do ! lband 591 592 end do ! jband 593 end do ! iband 594 595 d3_aux(1,:) = d3_aux(1,:) + & 596 & dk(:)*mvwtk(ineigh,ikpt)*dtset%wtk(ikpt)*(2_dp*dotposr-dotnegr) 597 d3_aux(2,:) = d3_aux(2,:) + & 598 & dk(:)*mvwtk(ineigh,ikpt)*dtset%wtk(ikpt)*(2_dp*dotposi-dotnegi) 599 600 end do ! End loop over neighbours 601 602 603 end do ! End loop over k-points 604 605 end do ! fab: end loop over spin 606 607 608 609 610 call xmpi_sum(d3_aux,spaceComm,ierr) 611 612 613 ABI_DEALLOCATE(s13mat) 614 ABI_DEALLOCATE(smat) 615 ABI_DEALLOCATE(s1mat) 616 ABI_DEALLOCATE(qmat) 617 ABI_DEALLOCATE(ipvt) 618 ABI_DEALLOCATE(s3mat) 619 ABI_DEALLOCATE(zgwork) 620 ABI_DEALLOCATE(bd_index) 621 622 623 !fab: I think that in the following we have to make a distinction: 624 !for the spin unpolarized case we leave the PEAD expression as it is, while 625 !in the spin polarized case we have simply to divide by 2 626 !(see eq.19 di PRB 63,155107, eq. 7 di PRB 71,125107 and eq 13 di PRB 71, 125107... 627 !in this latter equation the 2 must be simply replaced by the sum over the spin components... 628 !and indeed we have inserted the loop over the spin, 629 !but there was a factor 2 already present in the routine due to spin degenracy that had to be removed) 630 631 632 if (nsppol==1) then 633 634 ! Take minus the imaginary part 635 636 d3_berry(1,:) = -1_dp*d3_aux(2,:) 637 d3_berry(2,:) = d3_aux(1,:) 638 639 d3_berry(2,:) = 0_dp 640 641 else 642 643 d3_berry(1,:) = -1_dp*d3_aux(2,:)/2._dp 644 d3_berry(2,:) = d3_aux(1,:)/2._dp 645 646 d3_berry(2,:) = 0_dp/2._dp 647 648 end if 649 650 651 652 653 !DEBUG 654 !write(100,*)'dfptnl_mv.f : d3_berry' 655 !write(100,*)'Perturbation',i1dir,i3dir 656 !write(100,*) 657 !write(100,*)'before transformation' 658 !write(100,*)'real part' 659 !write(100,'(3(2x,f20.9))')d3_berry(1,:) 660 !write(100,*) 661 !write(100,*)'imaginary part' 662 !write(100,'(3(2x,f20.9))')d3_berry(2,:) 663 !write(100,*) 664 !write(100,*)'after transformation' 665 !ENDDEBUG 666 667 !Compute the projection on the basis vectors of 668 !reciprocal space 669 670 d3_aux(:,:) = 0_dp 671 do ii = 1,3 672 do jj = 1,3 673 d3_aux(:,ii) = d3_aux(:,ii) + gmet(ii,jj)*d3_berry(:,jj) 674 end do 675 end do 676 d3_berry(:,:) = d3_aux(:,:) 677 678 !Write out the berryphase part of the third order energy 679 680 if (mpi_enreg%me == 0) then 681 682 write(message,'(a,a,a)')ch10,& 683 & ' Berryphase part of the third-order energy:',ch10 684 ! call wrtout(ab_out,message,'COLL') 685 call wrtout(std_out, message,'COLL') 686 687 if (i1pert < natom + 1) then 688 write(message,'(a,i3,a,i3)')& 689 & ' j1: Displacement of atom ',i1pert,& 690 & ' along direction ',i1dir 691 else if (i1pert == natom + 2) then 692 write(message,'(a,i3)')& 693 & ' j1: homogenous electric field along direction ',i1dir 694 end if 695 ! call wrtout(ab_out,message,'COLL') 696 call wrtout(std_out, message,'COLL') 697 698 write(message,'(a)')& 699 & ' j2: k-point derivative along direction i2dir ' 700 ! call wrtout(ab_out,message,'COLL') 701 call wrtout(std_out, message,'COLL') 702 703 if (i3pert < natom + 1) then 704 write(message,'(a,i3,a,i3,a)')& 705 & ' j3: Displacement of atom ',i3pert,& 706 & ' along direction ',i3dir,ch10 707 else if (i3pert == natom + 2) then 708 write(message,'(a,i3,a)')& 709 & ' j3: homogenous electric field along direction ',i3dir,ch10 710 end if 711 ! call wrtout(ab_out,message,'COLL') 712 call wrtout(std_out, message,'COLL') 713 714 ! write(ab_out,'(5x,a5,8x,a9,5x,a14)')'i2dir','real part','imaginary part' 715 write(std_out,'(5x,a5,8x,a9,5x,a14)')'i2dir','real part','imaginary part' 716 do ii = 1,3 717 write(std_out,'(7x,i1,3x,f16.9,3x,f16.9)')ii,& 718 & d3_berry(1,ii),d3_berry(2,ii) 719 write(std_out,'(7x,i1,3x,f16.9,3x,f16.9)')ii,& 720 & d3_berry(1,ii),d3_berry(2,ii) 721 end do 722 723 end if ! mpi_enreg%me == 0 724 725 !DEBUG 726 !write(100,*)'real part' 727 !write(100,'(3(2x,f20.9))')d3_berry(1,:) 728 !write(100,*) 729 !write(100,*)'imaginary part' 730 !write(100,'(3(2x,f20.9))')d3_berry(2,:) 731 !ENDDEBUG 732 733 734 735 !close(dtfil%unwff1) 736 !close(dtfil%unwff2) 737 738 call status(0,dtfil%filstat,iexit,level,'exit ') 739 740 !DEBUG 741 !write(std_out,*)' dfptnl_mv : exit ' 742 !ENDDEBUG 743 744 end subroutine dfptnl_mv