TABLE OF CONTENTS


ABINIT/fxphas [ Functions ]

[ Top ] [ 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