1! 2! Copyright (C) 2004 Vanderbilt's group at Rutgers University, NJ 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8! April 2012, A. Dal Corso: parallelization for gdir /= 3 imported 9! from c_phase_field.f90 10! May 2012, A. Dal Corso: Noncollinear/spin-orbit case allowed (experimental). 11! January 2019 Ronald E Cohen, fixed mods and averages on same branch 12! 13!##############################################################################! 14!# #! 15!# #! 16!# This is the main one of a set of Fortran 90 files designed to compute #! 17!# the electrical polarization in a crystaline solid. #! 18!# #! 19!# #! 20!# AUTHORS #! 21!# ~~~~~~~ #! 22!# This set of subprograms is based on code written in an early Fortran #! 23!# 77 version of PWSCF by Alessio Filippetti. These were later ported #! 24!# into another version by Lixin He. Oswaldo Dieguez, in collaboration #! 25!# with Lixin He and Jeff Neaton, ported these routines into Fortran 90 #! 26!# version 1.2.1 of PWSCF. He, Dieguez, and Neaton were working at the #! 27!# time in David Vanderbilt's group at Rutgers, The State University of #! 28!# New Jersey, USA. #! 29!# #! 30!# #! 31!# LIST OF FILES #! 32!# ~~~~~~~~~~~~~ #! 33!# The complete list of files added to the PWSCF distribution is: #! 34!# * ../PW/bp_calc_btq.f90 #! 35!# * ../PW/bp_c_phase.f90 #! 36!# * ../PW/bp_qvan3.f90 #! 37!# * ../PW/bp_strings.f90 #! 38!# #! 39!# The PWSCF files that needed (minor) modifications were: #! 40!# * ../PW/electrons.f90 #! 41!# * ../PW/input.f90 #! 42!# * ../PW/pwcom.f90 #! 43!# * ../PW/setup.f90 #! 44!# #! 45!# Present in the original version and later removed: #! 46!# * bp_ylm_q.f bp_dbess.f bp_radin.f bp_bess.f #! 47!# #! 48!# BRIEF SUMMARY OF THE METHODOLOGY #! 49!# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #! 50!# The spontaneous polarization has two contibutions, electronic #! 51!# and ionic. With these additional routines, PWSCF will output both. #! 52!# #! 53!# The ionic contribution is relatively trivial to compute, requiring #! 54!# knowledge only of the atomic positions and core charges. The new #! 55!# subroutines focus mainly on evaluating the electronic contribution, #! 56!# computed as a Berry phase, i.e., a global phase property that can #! 57!# be computed from inner products of Bloch states at neighboring #! 58!# points in k-space. #! 59!# #! 60!# The standard procedure would be for the user to first perform a #! 61!# self-consistent (sc) calculation to obtain a converged charge density. #! 62!# With well-converged sc charge density, the user would then run one #! 63!# or more non-self consistent (or "band structure") calculations, #! 64!# using the same main code, but with a flag to ask for the polarization. #! 65!# Each such run would calculate the projection of the polarization #! 66!# onto one of the three primitive reciprocal lattice vectors. In #! 67!# cases of high symmetry (e.g. a tetragonal ferroelectric phase), one #! 68!# such run would suffice. In the general case of low symmetry, the #! 69!# user would have to submit up to three jobs to compute the three #! 70!# components of polarization, and would have to obtain the total #! 71!# polarization "by hand" by summing these contributions. #! 72!# #! 73!# Accurate calculation of the electronic or "Berry-phase" polarization #! 74!# requires overlaps between wavefunctions along fairly dense lines (or #! 75!# "strings") in k-space in the direction of the primitive G-vector for #! 76!# which one is calculating the projection of the polarization. The #! 77!# code would use a higher-density k-mesh in this direction, and a #! 78!# standard-density mesh in the two other directions. See below for #! 79!# details. #! 80!# #! 81!# #! 82!# FUNCTIONALITY/COMPATIBILITY #! 83!# ~~~~~~~~~~~~~~~~~~~~~~~~~~~ #! 84!# * Berry phases for a given G-vector. #! 85!# #! 86!# * Contribution to the polarization (in relevant units) for a given #! 87!# G-vector. #! 88!# #! 89!# * Spin-polarized systems supported. #! 90!# #! 91!# * Ultrasoft and norm-conserving pseudopotentials supported. #! 92!# #! 93!# * The value of the "polarization quantum" and the ionic contribution #! 94!# to the polarization are reported. #! 95!# #! 96!# #! 97!# NEW INPUT PARAMETERS #! 98!# ~~~~~~~~~~~~~~~~~~~~ #! 99!# * lberry (.TRUE. or .FALSE.) #! 100!# Tells PWSCF that a Berry phase calcultion is desired. #! 101!# #! 102!# * gdir (1, 2, or 3) #! 103!# Specifies the direction of the k-point strings in reciprocal space. #! 104!# '1' refers to the first reciprocal lattice vector, '2' to the #! 105!# second, and '3' to the third. #! 106!# #! 107!# * nppstr (integer) #! 108!# Specifies the number of k-points to be calculated along each #! 109!# symmetry-reduced string. #! 110!# #! 111!# #! 112!# EXPLANATION OF K-POINT MESH #! 113!# ~~~~~~~~~~~~~~~~~~~~~~~~~~~ #! 114!# If gdir=1, the program takes the standard input specification of the #! 115!# k-point mesh (nk1 x nk2 x nk3) and stops if the k-points in dimension #! 116!# 1 are not equally spaced or if its number is not equal to nppstr, #! 117!# working with a mesh of dimensions (nppstr x nk2 x nk3). That is, for #! 118!# each point of the (nk2 x nk3) two-dimensional mesh, it works with a #! 119!# string of nppstr k-points extending in the third direction. Symmetry #! 120!# will be used to reduce the number of strings (and assign them weights) #! 121!# if possible. Of course, if gdir=2 or 3, the variables nk2 or nk3 will #! 122!# be overridden instead, and the strings constructed in those #! 123!# directions, respectively. #! 124!# #! 125!# #! 126!# BIBLIOGRAPHY #! 127!# ~~~~~~~~~~~~ #! 128!# The theory behind this implementation is described in: #! 129!# [1] R D King-Smith and D Vanderbilt, "Theory of polarization of #! 130!# crystaline solids", Phys Rev B 47, 1651 (1993). #! 131!# #! 132!# Other relevant sources of information are: #! 133!# [2] D Vanderbilt and R D King-Smith, "Electronic polarization in the #! 134!# ultrasoft pseudopotential formalism", internal report (1998), #! 135!# [3] D Vanderbilt, "Berry phase theory of proper piezoelectric #! 136!# response", J Phys Chem Solids 61, 147 (2000). #! 137!# #! 138!# #! 139!# dieguez@physics.rutgers.edu #! 140!# 09 June 2003 #! 141!# #! 142!# #! 143!##############################################################################! 144!# Updated by Ronald Cohen, Carnegie Institution, 2019 145!# To make sure average over strings is done so that each straing is 146!# on the same branch cut. Also computations are kept as phases as long 147!# as possible before converting to polarization lattice 148 149 150!======================================================================! 151 152!--------------------------------------------------------------------- 153SUBROUTINE c_phase 154 !---------------------------------------------------------------------- 155 !! Geometric phase calculation along a strip of \(\text{nppstr}\) k-points 156 !! averaged over a 2D grid of \(\text{nkort}\) orthogonal k-points. 157 ! 158 USE kinds, ONLY : DP 159 USE io_global, ONLY : stdout 160 USE io_files, ONLY : iunwfc, nwordwfc 161 USE buffers, ONLY : get_buffer 162 USE ions_base, ONLY : nat, ntyp => nsp, ityp, tau, zv, atm 163 USE cell_base, ONLY : at, alat, tpiba, omega 164 USE constants, ONLY : pi, tpi 165 USE gvect, ONLY : ngm, g, gcutm, ngm_g, ig_l2g 166 USE fft_base, ONLY : dfftp 167 USE uspp, ONLY : nkb, vkb, okvan 168 USE uspp_param, ONLY : upf, lmaxq, nbetam, nh, nhm 169 USE lsda_mod, ONLY : nspin 170 USE klist, ONLY : nelec, degauss, nks, xk, wk, igk_k, ngk 171 USE wvfct, ONLY : npwx, nbnd, wg 172 USE wavefunctions, ONLY : evc 173 USE bp, ONLY : gdir, nppstr, mapgm_global, pdl_tot 174 USE becmod, ONLY : calbec, bec_type, allocate_bec_type, & 175 deallocate_bec_type 176 USE noncollin_module, ONLY : noncolin, npol, nspin_lsda 177 USE spin_orb, ONLY : lspinorb 178 USE mp_bands, ONLY : intra_bgrp_comm, nproc_bgrp 179 USE mp, ONLY : mp_sum 180 USE qes_libs_module, ONLY : qes_reset 181 USE qexsd_init, ONLY : qexsd_init_berryPhaseOutput, qexsd_bp_obj 182! --- Avoid implicit definitions --- 183 IMPLICIT NONE 184 185! --- Internal definitions --- 186 INTEGER :: i 187 INTEGER :: igk1(npwx) 188 INTEGER :: igk0(npwx) 189 INTEGER :: ig 190 INTEGER :: ind1 191 INTEGER :: info 192 INTEGER :: is 193 INTEGER :: istring 194 INTEGER :: iv 195 INTEGER :: ivpt(nbnd) 196 INTEGER :: j 197 INTEGER :: jkb 198 INTEGER :: jkb_bp 199 INTEGER :: jkb1 200 INTEGER :: job 201 INTEGER :: jv 202 INTEGER :: kindex 203 INTEGER :: kort 204 INTEGER :: kpar 205 INTEGER :: kpoint 206 INTEGER :: kstart 207 INTEGER :: mb 208 INTEGER :: mk1 209 INTEGER :: mk2 210 INTEGER :: mk3 211 INTEGER , ALLOCATABLE :: mod_elec(:) 212 INTEGER , ALLOCATABLE :: ln(:,:,:) 213 INTEGER :: mod_elec_dw 214 INTEGER :: mod_elec_tot 215 INTEGER :: mod_elec_up 216 INTEGER :: mod_ion(nat) 217 INTEGER :: mod_ion_tot 218 INTEGER :: mod_tot 219 INTEGER :: n1 220 INTEGER :: n2 221 INTEGER :: n3 222 INTEGER :: na 223 INTEGER :: nb 224 INTEGER :: ng 225 INTEGER :: nhjkb 226 INTEGER :: nhjkbm 227 INTEGER :: nkbtona(nkb) 228 INTEGER :: nkbtonh(nkb) 229 INTEGER :: nkort 230 INTEGER :: np 231 INTEGER :: npw1 232 INTEGER :: npw0 233 INTEGER :: nstring 234 INTEGER :: nbnd_occ 235 INTEGER :: nt 236 INTEGER, ALLOCATABLE :: map_g(:) 237 LOGICAL :: lodd 238 LOGICAL :: l_para 239 LOGICAL, ALLOCATABLE :: l_cal(:) ! flag for occupied/empty states 240 REAL(DP) :: t1,t !!REC 241 REAL(DP) :: dk(3) 242 REAL(DP) :: dkmod 243 REAL(DP) :: el_loc 244 REAL(DP) :: eps 245 REAL(DP) :: fac 246 REAL(DP) :: gpar(3) 247 REAL(DP) :: gtr(3) 248 REAL(DP) :: gvec 249 REAL(DP), ALLOCATABLE :: pdl_elec(:) 250 REAL(DP), ALLOCATABLE :: phik(:) 251 REAL(DP) :: phik_ave 252 REAL(DP) :: qrad_dk(nbetam,nbetam,lmaxq,ntyp) 253 REAL(DP) :: weight 254 REAL(DP) :: upol(3) 255 REAL(DP) :: pdl_elec_dw 256 REAL(DP) :: pdl_elec_tot 257 REAL(DP) :: pdl_elec_up 258 REAL(DP) :: pdl_ion(nat) 259 REAL(DP) :: pdl_ion_tot 260 REAL(DP) :: phidw 261 REAL(DP) :: phiup 262 REAL(DP) :: rmod 263 REAL(DP), ALLOCATABLE :: wstring(:) 264 REAL(DP) :: ylm_dk(lmaxq*lmaxq) 265 REAL(DP) :: zeta_mod 266 COMPLEX(DP), ALLOCATABLE :: aux(:) 267 COMPLEX(DP), ALLOCATABLE :: aux_g(:) 268 COMPLEX(DP), ALLOCATABLE :: aux0(:) 269 TYPE (bec_type) :: becp0 270 TYPE (bec_type) :: becp_bp 271 COMPLEX(DP) :: cave 272 COMPLEX(DP) , ALLOCATABLE :: cphik(:) 273 COMPLEX(DP) :: det 274 COMPLEX(DP) :: dtheta 275 COMPLEX(DP) :: mat(nbnd,nbnd) 276 COMPLEX(DP) :: pref 277 COMPLEX(DP), ALLOCATABLE :: psi(:,:) 278 COMPLEX(DP), ALLOCATABLE :: q_dk_so(:,:,:,:) 279 COMPLEX(DP) :: q_dk(nhm,nhm,ntyp) 280 COMPLEX(DP) :: struc(nat) 281 COMPLEX(DP) :: theta0 282 COMPLEX(DP) :: zeta 283 284! ------------------------------------------------------------------------- ! 285! INITIALIZATIONS 286! ------------------------------------------------------------------------- ! 287 ALLOCATE (psi(npwx*npol,nbnd)) 288 ALLOCATE (aux(ngm*npol)) 289 ALLOCATE (aux0(ngm*npol)) 290 IF (okvan) THEN 291 CALL allocate_bec_type ( nkb, nbnd, becp0 ) 292 CALL allocate_bec_type ( nkb, nbnd, becp_bp ) 293 IF (lspinorb) ALLOCATE(q_dk_so(nhm,nhm,4,ntyp)) 294 END IF 295 296 l_para= (nproc_bgrp > 1 .AND. gdir /= 3) 297 IF (l_para) THEN 298 ALLOCATE ( aux_g(ngm_g*npol) ) 299 ELSE 300 ALLOCATE ( map_g(ngm) ) 301 ENDIF 302 303! --- Write header --- 304 WRITE( stdout,"(/,/,/,15X,50('='))") 305 WRITE( stdout,"(28X,'POLARIZATION CALCULATION')") 306 WRITE( stdout,"(25X,'!!! NOT THOROUGHLY TESTED !!!')") 307 WRITE( stdout,"(15X,50('-'),/)") 308 309! --- Check that we are working with an insulator with no empty bands --- 310 IF ( degauss > 0.0_dp ) CALL errore('c_phase', & 311 'Polarization only for insulators',1) 312 313! --- Define a small number --- 314 eps=1.0E-6_dp 315 316! --- Recalculate FFT correspondence (see ggen.f90) --- 317 ALLOCATE (ln (-dfftp%nr1:dfftp%nr1, -dfftp%nr2:dfftp%nr2, -dfftp%nr3:dfftp%nr3) ) 318 DO ng=1,ngm 319 mk1=nint(g(1,ng)*at(1,1)+g(2,ng)*at(2,1)+g(3,ng)*at(3,1)) 320 mk2=nint(g(1,ng)*at(1,2)+g(2,ng)*at(2,2)+g(3,ng)*at(3,2)) 321 mk3=nint(g(1,ng)*at(1,3)+g(2,ng)*at(2,3)+g(3,ng)*at(3,3)) 322 ln(mk1,mk2,mk3) = ng 323 END DO 324 325 if(okvan) then 326! --- Initialize arrays --- 327 jkb_bp=0 328 DO nt=1,ntyp 329 DO na=1,nat 330 IF (ityp(na).eq.nt) THEN 331 DO i=1, nh(nt) 332 jkb_bp=jkb_bp+1 333 nkbtona(jkb_bp) = na 334 nkbtonh(jkb_bp) = i 335 END DO 336 END IF 337 END DO 338 END DO 339 endif 340! --- Get the number of strings --- 341 nstring=nks/nppstr 342 nkort=nstring/nspin_lsda 343 344! --- Allocate memory for arrays --- 345 ALLOCATE(phik(nstring)) 346 ALLOCATE(cphik(nstring)) 347 ALLOCATE(wstring(nstring)) 348 ALLOCATE(pdl_elec(nstring)) 349 ALLOCATE(mod_elec(nstring)) 350 351! ------------------------------------------------------------------------- ! 352! electronic polarization: set values for k-points strings ! 353! ------------------------------------------------------------------------- ! 354 355! --- Find vector along strings --- 356 gpar(1)=xk(1,nppstr)-xk(1,1) 357 gpar(2)=xk(2,nppstr)-xk(2,1) 358 gpar(3)=xk(3,nppstr)-xk(3,1) 359 gvec=dsqrt(gpar(1)**2+gpar(2)**2+gpar(3)**2)*tpiba 360 361! --- Find vector between consecutive points in strings --- 362 dk(1)=xk(1,2)-xk(1,1) 363 dk(2)=xk(2,2)-xk(2,1) 364 dk(3)=xk(3,2)-xk(3,1) 365 dkmod=SQRT(dk(1)**2+dk(2)**2+dk(3)**2)*tpiba 366 IF (ABS(dkmod-gvec/(nppstr-1)) > eps) & 367 CALL errore('c_phase','Wrong k-strings?',1) 368 369! --- Check that k-points form strings --- 370 DO i=1,nspin_lsda*nkort 371 DO j=2,nppstr 372 kindex=j+(i-1)*nppstr 373 IF (ABS(xk(1,kindex)-xk(1,kindex-1)-dk(1)) > eps) & 374 CALL errore('c_phase','Wrong k-strings?',1) 375 IF (ABS(xk(2,kindex)-xk(2,kindex-1)-dk(2)) > eps) & 376 CALL errore('c_phase','Wrong k-strings?',1) 377 IF (ABS(xk(3,kindex)-xk(3,kindex-1)-dk(3)) > eps) & 378 CALL errore('c_phase','Wrong k-strings?',1) 379 IF (ABS(wk(kindex)-wk(kindex-1)) > eps) & 380 CALL errore('c_phase','Wrong k-strings weights?',1) 381 END DO 382 END DO 383 384! ------------------------------------------------------------------------- ! 385! electronic polarization: weight strings ! 386! ------------------------------------------------------------------------- ! 387 388! --- Calculate string weights, normalizing to 1 (no spin or noncollinear) 389! or 1+1 (spin) --- 390 DO is=1,nspin_lsda 391 weight=0.0_dp 392 DO kort=1,nkort 393 istring=kort+(is-1)*nkort 394 wstring(istring)=wk(nppstr*istring) 395 weight=weight+wstring(istring) 396 END DO 397 DO kort=1,nkort 398 istring=kort+(is-1)*nkort 399 wstring(istring)=wstring(istring)/weight 400 END DO 401 END DO 402 403! ------------------------------------------------------------------------- ! 404! electronic polarization: structure factor ! 405! ------------------------------------------------------------------------- ! 406 407! --- Calculate structure factor e^{-i dk*R} --- 408 DO na=1,nat 409 fac=(dk(1)*tau(1,na)+dk(2)*tau(2,na)+dk(3)*tau(3,na))*tpi 410 struc(na)=CMPLX(cos(fac),-sin(fac),kind=DP) 411 END DO 412! ------------------------------------------------------------------------- ! 413! electronic polarization: form factor ! 414! ------------------------------------------------------------------------- ! 415 if(okvan) then 416! --- Calculate Bessel transform of Q_ij(|r|) at dk [Q_ij^L(|r|)] --- 417 CALL calc_btq(dkmod,qrad_dk,0) 418! --- Calculate the q-space real spherical harmonics at dk [Y_LM] --- 419 dkmod=dk(1)**2+dk(2)**2+dk(3)**2 420 CALL ylmr2(lmaxq*lmaxq, 1, dk, dkmod, ylm_dk) 421! --- Form factor: 4 pi sum_LM c_ij^LM Y_LM(Omega) Q_ij^L(|r|) --- 422 q_dk = (0.d0, 0.d0) 423 DO np =1, ntyp 424 if( upf(np)%tvanp ) then 425 DO iv = 1, nh(np) 426 DO jv = iv, nh(np) 427 call qvan3(iv,jv,np,pref,ylm_dk,qrad_dk) 428 q_dk(iv,jv,np) = omega*pref 429 q_dk(jv,iv,np) = omega*pref 430 ENDDO 431 ENDDO 432 endif 433 ENDDO 434 IF (lspinorb) CALL transform_qq_so(q_dk,q_dk_so) 435 endif 436 437! ------------------------------------------------------------------------- ! 438! electronic polarization: strings phases ! 439! ------------------------------------------------------------------------- ! 440 441 el_loc=0.d0 442 kpoint=0 443 ALLOCATE ( l_cal(nbnd) ) 444 CALL weights() 445 446! --- Start loop over spin --- 447 DO is=1,nspin_lsda 448 449 ! l_cal(n) = .true./.false. if n-th state is occupied/empty 450 nbnd_occ=0 451 DO nb = 1, nbnd 452 l_cal(nb) = (wg(nb,1+nks*(is-1)/2) > eps) 453 IF (l_cal(nb)) nbnd_occ = nbnd_occ + 1 454 END DO 455 456! --- Start loop over orthogonal k-points --- 457 DO kort=1,nkort 458 459! --- Index for this string --- 460 istring=kort+(is-1)*nkort 461 462! --- Initialize expectation value of the phase operator --- 463 zeta=(1.d0,0.d0) 464 zeta_mod = 1.d0 465 466! --- Start loop over parallel k-points --- 467 DO kpar = 1,nppstr 468 469! --- Set index of k-point --- 470 kpoint = kpoint + 1 471 472! --- Calculate dot products between wavefunctions and betas --- 473 IF (kpar /= 1) THEN 474 475! --- Dot wavefunctions and betas for PREVIOUS k-point --- 476 npw0 = ngk(kpoint-1) 477 igk0(:) = igk_k(:,kpoint-1) 478 CALL get_buffer (psi,nwordwfc,iunwfc,kpoint-1) 479 if (okvan) then 480 CALL init_us_2 (npw0,igk0,xk(1,kpoint-1),vkb) 481 CALL calbec (npw0, vkb, psi, becp0) 482 endif 483! --- Dot wavefunctions and betas for CURRENT k-point --- 484 IF (kpar /= nppstr) THEN 485 npw1 = ngk(kpoint) 486 igk1(:) = igk_k(:,kpoint) 487 CALL get_buffer(evc,nwordwfc,iunwfc,kpoint) 488 if (okvan) then 489 CALL init_us_2 (npw1,igk1,xk(1,kpoint),vkb) 490 CALL calbec (npw1, vkb, evc, becp_bp) 491 endif 492 ELSE 493 kstart = kpoint-nppstr+1 494 npw1 = ngk(kstart) 495 igk1(:) = igk_k(:,kstart) 496 CALL get_buffer(evc,nwordwfc,iunwfc,kstart) 497 if (okvan) then 498 CALL init_us_2 (npw1,igk1,xk(1,kstart),vkb) 499 CALL calbec(npw1, vkb, evc, becp_bp) 500 endif 501 ENDIF 502 503 IF (kpar == nppstr .AND. .NOT. l_para) THEN 504 map_g(:) = 0 505!$omp parallel & 506!$omp shared ( npw1, gpar, ln, gcutm, g, eps, map_g, at ) & 507!$omp private ( ig, gtr, n1, n2, n3, ng ) 508!$omp do 509 DO ig=1,npw1 510! --- If k'=k+G_o, the relation psi_k+G_o (G-G_o) --- 511! --- = psi_k(G) is used, gpar=G_o, gtr = G-G_o --- 512 513 gtr(1)=g(1,igk1(ig)) - gpar(1) 514 gtr(2)=g(2,igk1(ig)) - gpar(2) 515 gtr(3)=g(3,igk1(ig)) - gpar(3) 516! --- Find crystal coordinates of gtr, n1,n2,n3 --- 517! --- and the position ng in the ngm array --- 518 519 IF (gtr(1)**2+gtr(2)**2+gtr(3)**2 <= gcutm) THEN 520 n1=NINT(gtr(1)*at(1,1)+gtr(2)*at(2,1) & 521 +gtr(3)*at(3,1)) 522 n2=NINT(gtr(1)*at(1,2)+gtr(2)*at(2,2) & 523 +gtr(3)*at(3,2)) 524 n3=NINT(gtr(1)*at(1,3)+gtr(2)*at(2,3) & 525 +gtr(3)*at(3,3)) 526 ng=ln(n1,n2,n3) 527 528 IF ( (ABS(g(1,ng)-gtr(1)) > eps) .OR. & 529 (ABS(g(2,ng)-gtr(2)) > eps) .OR. & 530 (ABS(g(3,ng)-gtr(3)) > eps) ) THEN 531 WRITE(6,*) ' error: translated G=', & 532 gtr(1),gtr(2),gtr(3), & 533 & ' with crystal coordinates',n1,n2,n3, & 534 & ' corresponds to ng=',ng,' but G(ng)=', & 535 & g(1,ng),g(2,ng),g(3,ng) 536 WRITE(6,*) ' probably because G_par is NOT', & 537 & ' a reciprocal lattice vector ' 538 WRITE(6,*) ' Possible choices as smallest ', & 539 ' G_par:' 540 DO i=1,50 541 WRITE(6,*) ' i=',i,' G=', & 542 g(1,i),g(2,i),g(3,i) 543 ENDDO 544 CALL errore('c_phase','wrong g',1) 545 ENDIF 546 ELSE 547 WRITE(6,*) ' |gtr| > gcutm for gtr=', & 548 gtr(1),gtr(2),gtr(3) 549 CALL errore('c_phase','wrong gtr',1) 550 END IF 551 map_g(ig)=ng 552 END DO 553!$omp end do 554!$omp end parallel 555 END IF 556 557! --- Matrix elements calculation --- 558 559 mat(:,:) = (0.d0, 0.d0) 560 DO mb=1,nbnd 561 IF ( .NOT. l_cal(mb) ) THEN 562 mat(mb,mb)=(1.d0, 0.d0) 563 ELSE 564 aux(:) = (0.d0, 0.d0) 565 IF (kpar /= nppstr) THEN 566 DO ig=1,npw1 567 aux(igk1(ig))=evc(ig,mb) 568 ENDDO 569 IF (noncolin) THEN 570 DO ig=1,npw1 571 aux(igk1(ig)+ngm)=evc(ig+npwx,mb) 572 ENDDO 573 ENDIF 574 ELSEIF (.NOT. l_para) THEN 575 DO ig=1,npw1 576 aux(map_g(ig))=evc(ig,mb) 577 ENDDO 578 IF (noncolin) THEN 579 DO ig=1,npw1 580 aux(map_g(ig)+ngm)=evc(ig+npwx,mb) 581 ENDDO 582 ENDIF 583 ELSE 584! 585! In this case this processor might not have the G-G_0 586! 587 aux_g=(0.d0,0.d0) 588 DO ig=1,npw1 589 aux_g(mapgm_global(ig_l2g(igk1(ig)),gdir)) & 590 =evc(ig,mb) 591 ENDDO 592 IF (noncolin) THEN 593 DO ig=1,npw1 594 aux_g(mapgm_global(ig_l2g(igk1(ig)),gdir) & 595 + ngm_g) =evc(ig+npwx,mb) 596 ENDDO 597 ENDIF 598 CALL mp_sum(aux_g(:), intra_bgrp_comm ) 599 DO ig=1,ngm 600 aux(ig) = aux_g(ig_l2g(ig)) 601 ENDDO 602 IF (noncolin) THEN 603 DO ig=1,ngm 604 aux(ig+ngm) = aux_g(ig_l2g(ig)+ngm_g) 605 ENDDO 606 ENDIF 607 ENDIF 608! 609 DO nb=1,nbnd 610 IF ( l_cal(nb) ) THEN 611 aux0(:)= (0.d0, 0.d0) 612 DO ig=1,npw0 613 aux0(igk0(ig))=psi(ig,nb) 614 END DO 615 IF (noncolin) THEN 616 DO ig=1,npw0 617 aux0(igk0(ig)+ngm)=psi(ig+npwx,nb) 618 END DO 619 ENDIF 620 mat(nb,mb) = dot_product(aux0,aux) 621 END IF 622 END DO 623 END IF 624 END DO 625 ! 626 call mp_sum( mat, intra_bgrp_comm ) 627 ! 628 629 DO nb=1,nbnd 630!$omp parallel & 631!$omp shared ( nbnd, l_cal, nb, okvan, nkb, nkbtonh, ityp, nh, nkbtona ) & 632!$omp shared ( noncolin, lspinorb, becp0, q_dk_so, struc, mat ) & 633!$omp private ( mb, pref, jkb, nhjkb, na, np, nhjkbm, jkb1, j ) 634!$omp do 635 DO mb=1,nbnd 636! --- Calculate the augmented part: ij=KB projectors, --- 637! --- R=atom index: SUM_{ijR} q(ijR) <u_nk|beta_iR> --- 638! --- <beta_jR|u_mk'> e^i(k-k')*R = --- 639! --- also <u_nk|beta_iR>=<psi_nk|beta_iR> = becp^* --- 640 IF ( l_cal(nb) .AND. l_cal(mb) ) THEN 641 if (okvan) then 642 pref = (0.d0,0.d0) 643 DO jkb=1,nkb 644 nhjkb = nkbtonh(jkb) 645 na = nkbtona(jkb) 646 np = ityp(na) 647 nhjkbm = nh(np) 648 jkb1 = jkb - nhjkb 649 DO j = 1,nhjkbm 650 IF (noncolin) THEN 651 IF (lspinorb) THEN 652 pref = pref+(CONJG(becp0%nc(jkb,1,nb))* & 653 becp_bp%nc(jkb1+j,1,mb) & 654 *q_dk_so(nhjkb,j,1,np) & 655 +CONJG(becp0%nc(jkb,1,nb))* & 656 becp_bp%nc(jkb1+j,2,mb) & 657 *q_dk_so(nhjkb,j,2,np) & 658 +CONJG(becp0%nc(jkb,2,nb))* & 659 becp_bp%nc(jkb1+j,1,mb) & 660 *q_dk_so(nhjkb,j,3,np) & 661 +CONJG(becp0%nc(jkb,2,nb))* & 662 becp_bp%nc(jkb1+j,2,mb) & 663 *q_dk_so(nhjkb,j,4,np))*struc(na) 664 ELSE 665 pref = pref+(CONJG(becp0%nc(jkb,1,nb))* & 666 becp_bp%nc(jkb1+j,1,mb) + & 667 CONJG(becp0%nc(jkb,2,nb))* & 668 becp_bp%nc(jkb1+j,2,mb)) & 669 *q_dk(nhjkb,j,np)*struc(na) 670 END IF 671 ELSE 672 pref = pref+CONJG(becp0%k(jkb,nb))* & 673 becp_bp%k(jkb1+j,mb) & 674 *q_dk(nhjkb,j,np)*struc(na) 675 END IF 676 ENDDO 677 ENDDO 678 mat(nb,mb) = mat(nb,mb) + pref 679 endif 680 endif 681 ENDDO 682!$omp end do 683!$omp end parallel 684 ENDDO 685 686! --- Calculate matrix determinant --- 687 CALL ZGETRF (nbnd,nbnd,mat,nbnd,ivpt,info) 688 CALL errore('c_phase','error in factorization',abs(info)) 689 det=(1.d0,0.d0) 690 do nb=1,nbnd 691 det = det*mat(nb,nb) 692 if(nb.ne.ivpt(nb)) det=-det 693 enddo 694! --- Multiply by the already calculated determinants --- 695 zeta=zeta*det 696 697! --- End of dot products between wavefunctions and betas --- 698 ENDIF 699 700! --- End loop over parallel k-points --- 701 END DO 702 703! --- Calculate the phase for this string --- 704 phik(istring)=AIMAG(LOG(zeta)) 705 cphik(istring)=COS(phik(istring))*(1.0_dp,0.0_dp) & 706 +SIN(phik(istring))*(0.0_dp,1.0_dp) 707 708! --- Calculate the localization for current kort --- 709 zeta_mod= DBLE(CONJG(zeta)*zeta) 710!REC if zeta_mod=0 then angle is zero! 711 if(zeta_mod.le.eps)then 712 phik(istring)=0d0 713 cphik(istring)=cmplx(1d0,0d0) 714 endif 715 716! --- End loop over orthogonal k-points --- 717 END DO !kort 718 719! --- End loop over spin --- 720 END DO 721 DEALLOCATE ( l_cal ) 722 723! ------------------------------------------------------------------------- ! 724! electronic polarization: phase average ! 725! ------------------------------------------------------------------------- ! 726 727! --- Start loop over spins --- 728 DO is=1,nspin_lsda 729 730! --- Initialize average of phases as complex numbers --- 731 cave=(0.0_dp,0.0_dp) 732 phik_ave=(0.0_dp,0.0_dp) 733 734! --- Start loop over strings with same spin --- 735 DO kort=1,nkort 736 737! --- Calculate string index --- 738 istring=kort+(is-1)*nkort 739 740! --- Average phases as complex numbers --- 741 cave=cave+wstring(istring)*cphik(istring) 742 743! --- End loop over strings with same spin --- 744 END DO 745 746! --- Get the angle corresponding to the complex numbers average --- 747 theta0=atan2(AIMAG(cave), DBLE(cave)) 748! --- Put the phases in an around theta0 --- 749 DO kort=1,nkort 750 istring=kort+(is-1)*nkort 751 cphik(istring)=cphik(istring)/cave 752 dtheta=atan2(AIMAG(cphik(istring)), DBLE(cphik(istring))) 753 phik(istring)=theta0+dtheta 754 end do 755!REC First you need to multiply phase by two if only summed over 1 set of bands for non-spin-polarized case NO--summed below 756! if(nspin_lsda.eq.1)phik(1:istring)=2d0*phik(1:istring) 757!REC Second you need to take mod so phase is -Pi to Pi 758 DO kort=1,nkort 759 istring=kort+(is-1)*nkort 760 phik(istring)=phik(istring)-tpi*nint(phik(istring)/tpi) 761 enddo 762!REC Third you need to fix jumps before you take average 763 t1=phik(1)/tpi 764 DO kort=1,nkort 765 istring=kort+(is-1)*nkort 766 t=phik(istring)/tpi 767 if(abs(t+1d0-t1).lt.abs(t-t1))phik(istring)=phik(istring)+tpi 768 if(abs(t-1d0-t1).lt.abs(t-t1))phik(istring)=phik(istring)-tpi 769 pdl_elec(istring) = phik(istring)/tpi 770 phik_ave=phik_ave+wstring(istring)*phik(istring) 771 enddo 772! --- Assign this angle to the corresponding spin phase average --- 773 IF (nspin == 1) THEN 774 phiup=phik_ave !theta0+dtheta 775 phidw=phik_ave !theta0+dtheta 776 pdl_elec = pdl_elec*2._DP 777 ELSE IF (nspin == 2) THEN 778 IF (is == 1) THEN 779 phiup=phik_ave !theta0+dtheta 780 ELSE IF (is == 2) THEN 781 phidw=phik_ave !theta0+dtheta 782 END IF 783 ELSE IF (nspin==4 ) THEN 784 phiup=phik_ave 785 phidw=0.0_DP 786 END IF 787! --- End loop over spins 788 END DO 789 790! ------------------------------------------------------------------------- ! 791! ------------------------------------------------------------------------- ! 792 pdl_elec_up=phiup/tpi 793 pdl_elec_dw=phidw/tpi 794 pdl_elec_tot=pdl_elec_up+pdl_elec_dw 795! you need to do mod again! 796 pdl_elec_tot=pdl_elec_tot-nint(pdl_elec_tot) 797 pdl_elec_up=pdl_elec_up-nint(pdl_elec_up) 798 pdl_elec_dw=pdl_elec_dw-nint(pdl_elec_dw) 799 800! ------------------------------------------------------------------------- ! 801! ionic polarization ! 802! ------------------------------------------------------------------------- ! 803 804! --- Look for ions with odd number of charges --- 805 mod_ion=2 806 lodd=.FALSE. 807 DO na=1,nat 808 IF (MOD(NINT(zv(ityp(na))),2) == 1) THEN 809 mod_ion(na)=1 810 lodd=.TRUE. 811 END IF 812 END DO 813 814! --- Calculate ionic polarization phase for every ion --- 815 pdl_ion=0.0_dp 816 DO na=1,nat 817 DO i=1,3 818 pdl_ion(na)=pdl_ion(na)+zv(ityp(na))*tau(i,na)*gpar(i) 819 ENDDO 820 IF (mod_ion(na) == 1) THEN 821 pdl_ion(na)=pdl_ion(na)-1.0_dp*nint(pdl_ion(na)/1.0_dp) 822 ELSE IF (mod_ion(na) == 2) THEN 823 pdl_ion(na)=pdl_ion(na)-2.0_dp*nint(pdl_ion(na)/2.0_dp) 824 END IF 825 ENDDO 826! --- Add up the phases modulo 2 iff the ionic charges are even numbers --- 827 828 !REC You don't need a correction for jumps for the ionic part 829 ! This doesn't do anything since there is not an average but a sum 830 pdl_ion_tot=SUM(pdl_ion(1:nat)) 831 IF (lodd) THEN 832 pdl_ion_tot=pdl_ion_tot-1.0d0*nint(pdl_ion_tot/1.0d0) 833 mod_ion_tot=1 834 else 835 pdl_ion_tot=pdl_ion_tot-2.d0*nint(pdl_ion_tot/2.d0) 836 mod_ion_tot=2 837 endif 838 839! ------------------------------------------------------------------------- ! 840! total polarization ! 841! ------------------------------------------------------------------------- ! 842 843! --- Add electronic and ionic contributions to total phase --- 844 pdl_tot=pdl_elec_tot+pdl_ion_tot 845 IF ((.NOT.lodd).AND.(nspin == 1)) THEN 846 mod_tot=2 847 ELSE 848 mod_tot=1 849 END IF 850! ------------------------------------------------------------------------- ! 851! write output information ! 852! ------------------------------------------------------------------------- ! 853 854! --- Information about the k-points string used --- 855 WRITE( stdout,"(/,21X,'K-POINTS STRINGS USED IN CALCULATIONS')") 856 WRITE( Stdout,"(21X,37('~'),/)") 857 WRITE( stdout,"(7X,'G-vector along string (2 pi/a):',3F9.5)") & 858 gpar(1),gpar(2),gpar(3) 859 WRITE( stdout,"(7X,'Modulus of the vector (1/bohr):',F9.5)") & 860 gvec 861 WRITE( stdout,"(7X,'Number of k-points per string:',I4)") nppstr 862 WRITE( stdout,"(7X,'Number of different strings :',I4)") nkort 863 864! --- Information about ionic polarization phases --- 865 WRITE( stdout,"(2/,31X,'IONIC POLARIZATION')") 866 WRITE( stdout,"(31X,18('~'),/)") 867! WRITE( stdout,"(8X,'Note: (mod 1) means that the phases (angles ranging from' & 868! & /,8X,'-pi to pi) have been mapped to the interval [-1/2,+1/2) by'& 869! & /,8X,'dividing by 2*pi; (mod 2) refers to the interval [-1,+1)',& 870! & /)") 871 WRITE( stdout,"(2X,76('='))") 872 WRITE( stdout,"(4X,'Ion',4X,'Species',4X,'Charge',14X, & 873 & 'Position',16X,'Phase')") 874 WRITE( stdout,"(2X,76('-'))") 875 DO na=1,nat 876 WRITE( stdout,"(3X,I3,8X,A2,F12.3,5X,3F8.4,F12.5,' (mod ',I1,')')") & 877 & na,atm(ityp(na)),zv(ityp(na)), & 878 & tau(1,na),tau(2,na),tau(3,na),pdl_ion(na),mod_ion(na) 879 END DO 880 WRITE( stdout,"(2X,76('-'))") 881 WRITE( stdout,"(47X,'IONIC PHASE: ',F9.5,' (mod ',I1,')')") pdl_ion_tot,mod_ion_tot 882 WRITE( stdout,"(2X,76('='))") 883 884! --- Information about electronic polarization phases --- 885 WRITE( stdout,"(2/,28X,'ELECTRONIC POLARIZATION')") 886 WRITE( stdout,"(28X,23('~'),/)") 887! WRITE( stdout,"(8X,'Note: (mod 1) means that the phases (angles ranging from' & 888! & /,8X,'-pi to pi) have been mapped to the interval [-1/2,+1/2) by',& 889! & /,8X,'dividing by 2*pi; (mod 2) refers to the interval [-1,+1)',& 890! & /)") 891 WRITE( stdout,"(2X,76('='))") 892 IF(nspin_lsda.gt.1)THEN 893 WRITE( stdout,"(3X,'Spin',4X,'String',5X,'Weight',6X, & 894 & 'First k-point in string',9X,'Phase')") 895 WRITE( stdout,"(2X,76('-'))") 896 DO istring=1,nstring/nspin_lsda 897 ind1=1+(istring-1)*nppstr 898 WRITE( stdout,"(3X,' up ',3X,I5,F14.6,4X,3(F8.4),F12.5)") & 899 & istring,wstring(istring), & 900 & xk(1,ind1),xk(2,ind1),xk(3,ind1),pdl_elec(istring) 901 END DO 902 WRITE( stdout,"(2X,76('-'))") 903 DO istring=nstring/2+1,nstring 904 ind1=1+(istring-1)*nppstr 905 WRITE( stdout,"(3X,'down',3X,I4,F15.6,4X,3(F8.4),F12.5)") & 906 & istring,wstring(istring), xk(1,ind1),xk(2,ind1),xk(3,ind1), & 907 & pdl_elec(istring) 908 END DO 909 WRITE( stdout,"(2X,76('-'))") 910 WRITE( stdout,"(40X,'Average phase (up): ',F9.5)") & 911 pdl_elec_up 912 WRITE( stdout,"(38X,'Average phase (down): ',F9.5)")& 913 pdl_elec_dw 914 ENDIF 915! --- Treat unpolarized/polarized spin cases --- 916 IF (nspin_lsda == 1) THEN 917 WRITE( stdout,"(4X,'String',5X,'Weight',6X, & 918 & 'First k-point in string',9X,'Phase')") 919 WRITE( stdout,"(2X,76('-'))") 920 DO istring=1,nstring/nspin_lsda 921 ind1=1+(istring-1)*nppstr 922 WRITE( stdout,"(3X,I5,F14.6,4X,3(F8.4),F12.5)") & 923 & istring,wstring(istring), & 924 & xk(1,ind1),xk(2,ind1),xk(3,ind1),pdl_elec(istring) 925 END DO 926 WRITE( stdout,"(2X,76('-'))") 927 END IF 928 IF (noncolin) THEN 929 WRITE( stdout,"(42X,'Average phase : ',F9.5)") & 930 pdl_elec_up 931 ELSE 932 WRITE( stdout,"(42X,'ELECTRONIC PHASE: ',F9.5)") & 933 pdl_elec_tot 934 ENDIF 935 WRITE( stdout,"(2X,76('='))") 936 937! --- Information about total phase --- 938 WRITE( stdout,"(2/,31X,'SUMMARY OF PHASES')") 939 WRITE( stdout,"(31X,17('~'),/)") 940 WRITE( stdout,"(26X,'Ionic Phase:',F9.5)") & 941 pdl_ion_tot 942 WRITE( stdout,"(21X,'Electronic Phase:',F9.5)") & 943 pdl_elec_tot 944 WRITE( stdout,"(26X,'TOTAL PHASE:',F9.5, ' MOD_TOT:',i2)") & 945 pdl_tot,mod_tot 946 947! --- Information about the value of polarization --- 948 WRITE( stdout,"(2/,29X,'VALUES OF POLARIZATION')") 949 WRITE( stdout,"(29X,22('~'),/)") 950 WRITE( stdout,"( & 951 & 8X,'The calculation of phases done along the direction of vector ',I1, & 952 & /,8X,'of the reciprocal lattice gives the following contribution to', & 953 & /,8X,'the polarization vector (in different units, and being Omega', & 954 & /,8X,'the volume of the unit cell):')") & 955 gdir 956! --- Calculate direction of polarization and modulus of lattice vector --- 957 rmod=SQRT(at(1,gdir)*at(1,gdir)+at(2,gdir)*at(2,gdir) & 958 +at(3,gdir)*at(3,gdir)) 959 upol(:)=at(:,gdir)/rmod 960 rmod=alat*rmod 961! --- Give polarization in units of (e/Omega).bohr --- 962 fac=rmod 963 WRITE( stdout,"(/,11X,'P = ',F11.7,' (mod ',F11.7,') (e/Omega).bohr')") & 964 fac*pdl_tot,fac*DBLE(mod_tot) 965! --- Give polarization in units of e.bohr --- 966 fac=rmod/omega 967 WRITE( stdout,"(/,11X,'P = ',F11.7,' (mod ',F11.7,') e/bohr^2')") & 968 fac*pdl_tot,fac*DBLE(mod_tot) 969! --- Give polarization in SI units (C/m^2) --- 970 fac=(rmod/omega)*(1.60097E-19_dp/5.29177E-11_dp**2) 971 WRITE( stdout,"(/,11X,'P = ',F11.7,' (mod ',F11.7,') C/m^2')") & 972 fac*pdl_tot,fac*DBLE(mod_tot) 973! --- Write polarization direction --- 974 WRITE( stdout,"(/,8X,'The polarization direction is: (', & 975 & F8.5,' ,',F8.5,' ,',F8.5,' )')") upol(1),upol(2),upol(3) 976 977! --- End of information relative to polarization calculation --- 978 WRITE( stdout,"(/,/,15X,50('=')/,/)") 979 980 !------------------------------------------------------------------------------ 981! INITIALIZE QEXSD OUTPUT ELEMENT 982! Here we write all output information in a berry_phase_type variable to print 983! them in the XML output P.D. april 2016 984!------------------------------------------------------------------------------ 985 CALL qes_reset( qexsd_bp_obj ) 986 CALL qexsd_init_berryPhaseOutput(qexsd_bp_obj, gpar, gvec, nppstr, nkort, xk, pdl_ion, mod_ion, & 987 pdl_ion_tot, mod_ion_tot, nstring, pdl_elec , mod_elec, wstring, & 988 pdl_elec_up, mod_elec_up, pdl_elec_dw, mod_elec_dw, pdl_elec_tot,& 989 mod_elec_tot, pdl_tot, mod_tot, upol, rmod) 990! ------------------------------------------------------------------------- ! 991! finalization ! 992! ------------------------------------------------------------------------- ! 993 994! --- Free memory --- 995 DEALLOCATE(mod_elec) 996 DEALLOCATE(pdl_elec) 997 DEALLOCATE(wstring) 998 DEALLOCATE(cphik) 999 DEALLOCATE(phik) 1000 DEALLOCATE(ln) 1001 DEALLOCATE(aux) 1002 DEALLOCATE(aux0) 1003 DEALLOCATE(psi) 1004 IF (l_para) THEN 1005 DEALLOCATE ( aux_g ) 1006 ELSE 1007 DEALLOCATE ( map_g ) 1008 ENDIF 1009 IF (okvan) THEN 1010 CALL deallocate_bec_type ( becp0 ) 1011 CALL deallocate_bec_type ( becp_bp ) 1012 IF (lspinorb) DEALLOCATE(q_dk_so) 1013 END IF 1014 1015!------------------------------------------------------------------------------! 1016 1017END SUBROUTINE c_phase 1018 1019!==============================================================================! 1020