TABLE OF CONTENTS
ABINIT/fxphas [ Functions ]
NAME
fxphas
FUNCTION
Fix phase of all bands. Keep normalization but maximize real part (minimize imag part). Also fix the sign of real part by setting the first non-zero element to be positive.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, MT) 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,mcg)= contains the wavefunction |c> coefficients. gsc(2,mgsc)= if useoverlap==1, contains the S|c> coefficients, where S is an overlap matrix. icg=shift to be applied on the location of data in the array cg igsc=shift to be applied on the location of data in the array gsc istwfk=input option parameter that describes the storage of wfs (set to 1 if usual complex vectors) mcg=size of second dimension of cg mgsc=size of second dimension of gsc mpi_enreg=information about MPI parallelization nband_k=number of bands npw_k=number of planewaves useoverlap=describe the overlap of wavefunctions: 0: no overlap (S=Identi0,ty_matrix) 1: wavefunctions are overlapping
OUTPUT
cg(2,mcg)=same array with altered phase. gsc(2,mgsc)= same array with altered phase.
NOTES
When the sign of the real part was fixed (modif v3.1.3g.6), the test Tv3#5 , dataset 5, behaved differently than previously. This should be cleared up.
PARENTS
vtowfk
CHILDREN
timab,xmpi_sum
SOURCE
53 #if defined HAVE_CONFIG_H 54 #include "config.h" 55 #endif 56 57 #include "abi_common.h" 58 59 subroutine fxphas(cg,gsc,icg,igsc,istwfk,mcg,mgsc,mpi_enreg,nband_k,npw_k,useoverlap) 60 61 use defs_basis 62 use defs_abitypes 63 use m_errors 64 use m_profiling_abi 65 use m_xmpi 66 67 !This section has been created automatically by the script Abilint (TD). 68 !Do not modify the following lines by hand. 69 #undef ABI_FUNC 70 #define ABI_FUNC 'fxphas' 71 use interfaces_18_timing 72 !End of the abilint section 73 74 implicit none 75 76 !Arguments ------------------------------------ 77 !scalars 78 integer,intent(in) :: icg,igsc,istwfk,mcg,mgsc,nband_k,npw_k,useoverlap 79 type(MPI_type),intent(in) :: mpi_enreg 80 !arrays 81 real(dp),intent(inout) :: cg(2,mcg),gsc(2,mgsc*useoverlap) 82 83 !Local variables------------------------------- 84 !scalars 85 integer :: iband,ierr,ii,indx 86 real(dp) :: cim,cre,gscim,gscre,quotient,root1,root2,saa,sab,sbb,theta 87 real(dp) :: thppi,xx,yy 88 character(len=500) :: message 89 !arrays 90 real(dp) :: buffer2(nband_k,2),buffer3(nband_k,3),tsec(2) 91 real(dp),allocatable :: cimb(:),creb(:),saab(:),sabb(:),sbbb(:) !,sarr(:,:) 92 93 ! ************************************************************************* 94 95 !The general case, where a complex phase indeterminacy is present 96 if(istwfk==1)then 97 98 ABI_ALLOCATE(cimb,(nband_k)) 99 ABI_ALLOCATE(creb,(nband_k)) 100 ABI_ALLOCATE(saab,(nband_k)) 101 ABI_ALLOCATE(sabb,(nband_k)) 102 ABI_ALLOCATE(sbbb,(nband_k)) 103 cimb(:)=zero ; creb(:)=zero 104 105 ! Loop over bands 106 ! TODO: MG store saa arrays in sarr(3,nband_k) to reduce false sharing. 107 !$OMP PARALLEL DO DEFAULT(PRIVATE) SHARED(nband_k,icg,npw_k,cg,saab,sbbb,sabb) 108 do iband=1,nband_k 109 indx=icg+(iband-1)*npw_k 110 111 ! Compute several sums over Re, Im parts of c 112 saa=zero; sbb=zero; sab=zero 113 do ii=1+indx,npw_k+indx 114 saa=saa+cg(1,ii)*cg(1,ii) 115 sbb=sbb+cg(2,ii)*cg(2,ii) 116 sab=sab+cg(1,ii)*cg(2,ii) 117 end do 118 saab(iband)=saa 119 sbbb(iband)=sbb 120 sabb(iband)=sab 121 end do 122 123 ! XG030513 : MPIWF : should transmit saab,sbbb,sabb from the procs 124 ! of the WF group to the master processor of the WF group 125 if (mpi_enreg%paral_kgb == 1) then 126 buffer3(:,1)=saab(:) 127 buffer3(:,2)=sbbb(:) 128 buffer3(:,3)=sabb(:) 129 call timab(48,1,tsec) 130 call xmpi_sum(buffer3,mpi_enreg%comm_fft,ierr) 131 if (mpi_enreg%paral_spinor==1) then 132 call xmpi_sum(buffer3,mpi_enreg%comm_spinor,ierr) 133 end if 134 call timab(48,2,tsec) 135 saab(:)=buffer3(:,1) 136 sbbb(:)=buffer3(:,2) 137 sabb(:)=buffer3(:,3) 138 end if 139 140 ! XG030513 : MPIWF this loop should only be executed by the master of the WF group 141 142 if (mpi_enreg%paral_kgb==0.or.mpi_enreg%me_fft==0) then 143 do iband=1,nband_k 144 indx=icg+(iband-1)*npw_k 145 146 saa=saab(iband) 147 sbb=sbbb(iband) 148 sab=sabb(iband) 149 150 ! Get phase angle theta 151 if (sbb+saa>tol8)then 152 if(abs(sbb-saa)>tol8*(sbb+saa) .or. 2*abs(sab)>tol8*(sbb+saa))then 153 if (abs(sbb-saa)>tol8*abs(sab)) then 154 quotient=sab/(sbb-saa) 155 theta=0.5_dp*atan(2.0_dp*quotient) 156 else 157 ! Taylor expansion of the atan in terms of inverse of its argument. Correct up to 1/x2, included. 158 theta=0.25_dp*(pi-(sbb-saa)/sab) 159 end if 160 ! Check roots to get theta for max Re part 161 root1=cos(theta)**2*saa+sin(theta)**2*sbb-2.0_dp*cos(theta)*sin(theta)*sab 162 thppi=theta+0.5_dp*pi 163 root2=cos(thppi)**2*saa+sin(thppi)**2*sbb-2.0_dp*cos(thppi)*sin(thppi)*sab 164 if (root2>root1) theta=thppi 165 else 166 ! The real part vector and the imaginary part vector are orthogonal, and of same norm. Strong indeterminacy. 167 ! Will determine the first non-zero coefficient, and fix its phase 168 ! Hypothesis : there is at least one non-zero element on the master node ... 169 do ii=1+indx,npw_k+indx 170 cre=cg(1,ii) 171 cim=cg(2,ii) 172 if(cre**2+cim**2>tol8**2*(saa+sbb))then 173 if(cre**2>tol8**2**cim**2)then 174 theta=atan(cim/cre) 175 else 176 ! Taylor expansion of the atan in terms of inverse of its argument. Correct up to 1/x2, included. 177 theta=pi/2-cre/cim 178 end if 179 exit 180 end if 181 end do 182 end if 183 else 184 write(message,'(a,i0,5a)')& 185 & 'The eigenvector with band ',iband,' has zero norm.',ch10,& 186 & 'This usually happens when the number of bands (nband) is comparable to the number of planewaves (mpw)',ch10,& 187 & 'Action: Check the parameters of the calculation. If nband ~ mpw, then decrease nband or, alternatively, increase ecut' 188 MSG_ERROR(message) 189 end if 190 191 xx=cos(theta) 192 yy=sin(theta) 193 194 ! Here, set the first non-zero element to be positive 195 ! Comment the next nine lines to recover the behaviour of pre v3.1.3g 196 ! Hypothesis : there is at least one non-zero element on the master node ... 197 do ii=1+indx,npw_k+indx 198 cre=cg(1,ii) 199 cim=cg(2,ii) 200 cre=xx*cre-yy*cim 201 if(abs(cre)>tol8)exit 202 end do 203 if(cre<zero)then 204 xx=-xx ; yy=-yy 205 end if 206 207 creb(iband)=xx 208 cimb(iband)=yy 209 210 end do 211 end if 212 213 ! XG030513 : MPIWF : should transmit creb(:),cimb(:) of the master 214 ! processor of the WF group to the others procs of the WF group 215 if (mpi_enreg%paral_kgb == 1) then 216 call timab(48,1,tsec) 217 buffer2(:,1)=creb(:) 218 buffer2(:,2)=cimb(:) 219 call xmpi_sum(buffer2,mpi_enreg%comm_fft,ierr) 220 if (mpi_enreg%paral_spinor==1) then 221 call xmpi_sum(buffer2,mpi_enreg%comm_spinor,ierr) 222 end if 223 creb(:)=buffer2(:,1) 224 cimb(:)=buffer2(:,2) 225 end if 226 227 ! MG TODO: Scaling can be done with zscal 228 !$OMP PARALLEL DO PRIVATE(indx,xx,yy,cre,cim,gscre,gscim) 229 do iband=1,nband_k 230 indx=icg+(iband-1)*npw_k 231 232 xx=creb(iband) 233 yy=cimb(iband) 234 ! Alter phase of array |cg> 235 do ii=1+indx,npw_k+indx 236 cre=cg(1,ii) 237 cim=cg(2,ii) 238 cg(1,ii)=xx*cre-yy*cim 239 cg(2,ii)=xx*cim+yy*cre 240 end do 241 242 ! Alter phase of array S|cg> 243 if (useoverlap==1) then 244 indx=igsc+(iband-1)*npw_k 245 do ii=1+indx,npw_k+indx 246 gscre=gsc(1,ii) 247 gscim=gsc(2,ii) 248 gsc(1,ii)=xx*gscre-yy*gscim 249 gsc(2,ii)=xx*gscim+yy*gscre 250 end do 251 end if 252 end do ! iband 253 254 ABI_DEALLOCATE(cimb) 255 ABI_DEALLOCATE(creb) 256 ABI_DEALLOCATE(saab) 257 ABI_DEALLOCATE(sabb) 258 ABI_DEALLOCATE(sbbb) 259 260 else ! if istwfk/=1. Storages that take into account the time-reversal symmetry : the freedom is only a sign freedom 261 262 ABI_ALLOCATE(creb,(nband_k)) 263 creb(:)=zero 264 ! XG030513 : MPIWF : this loop should be done only by the master processor of the WF group 265 266 if (mpi_enreg%paral_kgb==0.or.mpi_enreg%me_fft==0) then 267 268 ! Loop over bands 269 do iband=1,nband_k 270 271 indx=icg+(iband-1)*npw_k 272 273 ! Here, set the first non-zero real element to be positive 274 do ii=1+indx,npw_k+indx 275 cre=cg(1,ii) 276 if(abs(cre)>tol8)exit 277 end do 278 creb(iband)=cre 279 280 end do ! iband 281 282 end if 283 ! XG030513 : MPIWF : should transmit cre(:) of the master processor of the WF group to the others 284 if (mpi_enreg%paral_kgb == 1) then 285 call timab(48,1,tsec) 286 call xmpi_sum(creb,mpi_enreg%comm_fft,ierr) 287 if (mpi_enreg%paral_spinor==1) then 288 call xmpi_sum(creb,mpi_enreg%comm_spinor,ierr) 289 end if 290 call timab(48,2,tsec) 291 end if 292 293 do iband=1,nband_k 294 cre=creb(iband) 295 if(cre<zero)then 296 indx=icg+(iband-1)*npw_k 297 do ii=1+indx,npw_k+indx 298 cg(1,ii)=-cg(1,ii) 299 cg(2,ii)=-cg(2,ii) 300 end do 301 if(useoverlap==1)then 302 indx=igsc+(iband-1)*npw_k 303 do ii=1+indx,npw_k+indx 304 gsc(1,ii)=-gsc(1,ii) 305 gsc(2,ii)=-gsc(2,ii) 306 end do 307 end if 308 end if 309 end do ! iband 310 311 ABI_DEALLOCATE(creb) 312 end if ! istwfk 313 314 end subroutine fxphas