diff --git a/.gitignore b/.gitignore index bea26fa..2a27f18 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,5 @@ *.pyc *.egg-info* -docs/build/* \ No newline at end of file +docs/build/* +build/ +.vscode \ No newline at end of file diff --git a/README.md b/README.md index 66caf3f..7f0cdd9 100644 --- a/README.md +++ b/README.md @@ -3,7 +3,7 @@ # aiida-z2pack Official Z2pack plugin for AiiDA. -The repository as of aiida-z2pack>=2.0 is compatible with aiida-core>=1.0.0. +The repository as of aiida-z2pack>=2.0 is compatible with aiida-core>=1.0.0 (tested up to 1.6.5). For compatibility with older versions use aiida-z2pack==1.0. The plugin supports Quantum ESPRESSO only. @@ -12,16 +12,17 @@ The plugin supports Quantum ESPRESSO only. ### How do I get set up? ### The Z2pack plugin has the following dependencies: -* numpy~=1.17,<1.18 -* scipy>=1.4.1 -* scikit-learn>=0.22 +* numpy~=1.17 +* scipy~=1.4.1 +* scikit-learn~=0.22 * z2pack==2.1.1 -* aiida_quantumespresso==3.1.0 -* aiida_wannier90>=2.0.0 +* aiida-core~=1.5.0 +* aiida_quantumespresso~=3.1.0 +* aiida_wannier90~=2.0.0 Installing: * `pip install .` -* or `pip install .[dev]` to install the dependencies for developers (pre-commit, ...) +* or `pip install .[pre-commit,tests,docs]` to install the dependencies for developers (pre-commit, ...) ### Contribution guidelines ### diff --git a/modified_bands_QE6.x/bands.f90 b/modified_bands_QE6.x/bands.f90 new file mode 100644 index 0000000..e93e87f --- /dev/null +++ b/modified_bands_QE6.x/bands.f90 @@ -0,0 +1,664 @@ +! +! Copyright (C) 2001-2016 Quantum ESPRESSO group +! This file is distributed under the terms of the +! GNU General Public License. See the file `License' +! in the root directory of the present distribution, +! or http://www.gnu.org/copyleft/gpl.txt . +! +! +!----------------------------------------------------------------------- +PROGRAM do_bands + !----------------------------------------------------------------------- + ! + ! See files INPUT_BANDS.* in Doc/ directory for usage + ! + ! + USE io_files, ONLY : prefix, tmp_dir + USE mp_global, ONLY : mp_startup + USE mp_pools, ONLY : npool + USE control_flags, ONLY : gamma_only + USE environment, ONLY : environment_start, environment_end + USE wvfct, ONLY : nbnd + USE klist, ONLY : nkstot, two_fermi_energies + USE noncollin_module, ONLY : noncolin, i_cons + USE lsda_mod, ONLY : nspin + USE io_global, ONLY : ionode, ionode_id, stdout + USE mp, ONLY : mp_bcast + USE mp_images, ONLY : intra_image_comm + ! + IMPLICIT NONE + ! + CHARACTER(LEN=256), EXTERNAL :: trimcheck + ! + CHARACTER (len=256) :: filband, filp, outdir + LOGICAL :: lsigma(4), lsym, lp, no_overlap, plot_2d, parity + INTEGER :: spin_component, firstk, lastk + INTEGER :: ios + ! + NAMELIST / bands / outdir, prefix, filband, filp, spin_component, lsigma,& + lsym, parity, lp, filp, firstk, lastk, no_overlap, plot_2d + ! + ! initialise environment + ! +#if defined(__MPI) + CALL mp_startup ( ) +#endif + CALL environment_start ( 'BANDS' ) + ! + ! set default values for variables in namelist + ! + prefix = 'pwscf' + CALL get_environment_variable( 'ESPRESSO_TMPDIR', outdir ) + IF ( trim( outdir ) == ' ' ) outdir = './' + filband = 'bands.out' + lsym=.true. + parity=.false. + no_overlap=.true. + plot_2d=.false. + lsigma=.false. + lp=.false. + filp='p_avg.dat' + firstk=0 + lastk=10000000 + spin_component = 1 + ! + ios = 0 + ! + IF ( ionode ) THEN + ! + CALL input_from_file ( ) + ! + READ (5, bands, iostat = ios) + ! + lsigma(4)=.false. + tmp_dir = trimcheck (outdir) + ! + ENDIF + ! + ! + CALL mp_bcast( ios, ionode_id, intra_image_comm ) + IF (ios /= 0) CALL errore ('bands', 'reading bands namelist', abs(ios) ) + ! + ! ... Broadcast variables + ! + CALL mp_bcast( tmp_dir, ionode_id, intra_image_comm ) + CALL mp_bcast( prefix, ionode_id, intra_image_comm ) + CALL mp_bcast( filband, ionode_id, intra_image_comm ) + CALL mp_bcast( filp, ionode_id, intra_image_comm ) + CALL mp_bcast( spin_component, ionode_id, intra_image_comm ) + CALL mp_bcast( firstk, ionode_id, intra_image_comm ) + CALL mp_bcast( lastk, ionode_id, intra_image_comm ) + CALL mp_bcast( lp, ionode_id, intra_image_comm ) + CALL mp_bcast( lsym, ionode_id, intra_image_comm ) + CALL mp_bcast( lsigma, ionode_id, intra_image_comm ) + CALL mp_bcast( no_overlap, ionode_id, intra_image_comm ) + CALL mp_bcast( plot_2d, ionode_id, intra_image_comm ) + + IF (plot_2d) THEN + lsym=.false. + lp=.false. + no_overlap=.true. + ENDIF + IF (lsym) no_overlap=.true. + + IF ( npool > 1 .and..not.(lsym.or.no_overlap)) CALL errore('bands', & + 'pools not implemented',npool) + IF ( spin_component < 1 .OR. spin_component > 2 ) & + CALL errore('bands','incorrect spin_component',1) + ! + ! Now allocate space for pwscf variables, read and check them. + ! + CALL read_file() + ! + IF (gamma_only) CALL errore('bands','gamma_only case not implemented',1) + IF (two_fermi_energies.or.i_cons /= 0) & + CALL errore('bands',& + 'The bands code with constrained magnetization has not been tested',1) + IF ( ANY(lsigma) .AND. .NOT.noncolin ) & + CALL errore ('punch_band', 'lsigma requires noncollinear run', 1 ) + IF ( spin_component/=1 .and. nspin/=2 ) & + CALL errore('punch_bands','incorrect spin_component',1) + ! + CALL openfil_pp() + ! + IF (plot_2d) THEN + CALL punch_band_2d(filband,spin_component) + ELSE + CALL punch_band(filband,spin_component,lsigma,no_overlap) + IF (lsym) CALL sym_band(filband,spin_component,firstk,lastk) + IF (lp) CALL write_p_avg(filp,spin_component,firstk,lastk) + IF (parity) CALL parity_band(filband,spin_component,firstk,lastk) + END IF + ! + CALL environment_end ( 'BANDS' ) + ! + CALL stop_pp + STOP +END PROGRAM do_bands +! +!----------------------------------------------------------------------- +SUBROUTINE punch_band (filband, spin_component, lsigma, no_overlap) + !----------------------------------------------------------------------- + ! + ! This routine writes the band energies on a file. The routine orders + ! the eigenvalues using the overlap of the eigenvectors to give + ! an estimate crossing and anticrossing of the bands. This simplified + ! method works in many, but not in all the cases. + ! + ! + USE kinds, ONLY : dp + USE ions_base, ONLY : nat, ityp, ntyp => nsp + USE cell_base, ONLY : at + USE constants, ONLY : rytoev + USE gvect, ONLY : g, ngm + USE klist, ONLY : xk, nks, nkstot, ngk, igk_k + USE io_files, ONLY : iunpun, nwordwfc, iunwfc + USE wvfct, ONLY : nbnd, et, npwx + USE uspp, ONLY : nkb, vkb + USE uspp_param, ONLY : upf, nh, nhm + USE noncollin_module, ONLY : noncolin, npol + USE wavefunctions, ONLY : evc + USE io_global, ONLY : ionode, ionode_id, stdout + USE mp, ONLY : mp_bcast + USE mp_images, ONLY : intra_image_comm + USE becmod, ONLY : calbec, bec_type, allocate_bec_type, & + deallocate_bec_type, becp + + IMPLICIT NONE + CHARACTER (len=*) :: filband + INTEGER, INTENT(IN) :: spin_component + LOGICAL, INTENT(IN) :: lsigma(4), no_overlap + + ! becp : at current k-point + INTEGER :: ibnd, jbnd, i, ik, ig, ig1, ig2, ipol, npw, ngmax, jmax + INTEGER :: nks1tot, nks2tot, nks1, nks2 + INTEGER :: iunpun_sigma(4), ios(0:4), done(nbnd) + CHARACTER(len=256) :: nomefile + REAL(dp):: pscur, psmax, psr(nbnd) + COMPLEX(dp), ALLOCATABLE :: psi(:,:), spsi(:,:), ps(:,:) + INTEGER, ALLOCATABLE :: work(:), igg(:) + INTEGER, ALLOCATABLE :: closest_band(:,:)! index for band ordering + REAL(DP), ALLOCATABLE:: sigma_avg(:,:,:) ! expectation value of sigma + REAL(DP), ALLOCATABLE:: et_(:,:) ! reordered eigenvalues in eV + + + IF (filband == ' ') RETURN + + iunpun = 19 + ios(:) = 0 + ! + IF ( ionode ) THEN + ! + OPEN (unit = iunpun, file = filband, status = 'unknown', form = & + 'formatted', iostat = ios(0)) + REWIND (iunpun) + DO ipol=1,4 + IF (lsigma(ipol)) THEN + iunpun_sigma(ipol)=iunpun+ipol + WRITE(nomefile,'(".",i1)') ipol + OPEN (unit = iunpun_sigma(ipol), & + file = trim(filband)//trim(nomefile), & + status = 'unknown', form='formatted', iostat = ios(ipol)) + REWIND (iunpun_sigma(ipol)) + ENDIF + ENDDO + ! + ENDIF + ! + CALL mp_bcast( ios, ionode_id, intra_image_comm ) + IF ( ios(0) /= 0 ) & + CALL errore ('punch_band', 'Opening filband file', abs(ios(0)) ) + DO ipol=1,4 + IF ( ios(ipol) /= 0 ) & + CALL errore ('punch_band', 'Opening filband.N file ', ipol) + ENDDO + ! + CALL find_nks1nks2(1,nkstot,nks1tot,nks1,nks2tot,nks2,spin_component) + ! + ! index of largest G in plane-wave basis set across k-points + ! + ALLOCATE ( closest_band(nbnd,nkstot) ) + DO ik=1,nkstot + ! + ! default ordering: leave bands as they are + ! + DO ibnd=1,nbnd + closest_band(ibnd,ik) = ibnd + END DO + END DO + ! + IF ( noncolin ) ALLOCATE ( sigma_avg(4,nbnd,nkstot) ) + ALLOCATE ( psi(npwx*npol,nbnd), spsi(npwx*npol,nbnd), ps(nbnd,nbnd) ) + CALL allocate_bec_type(nkb, nbnd, becp) + ! + DO ik = nks1, nks2 + ! + ! read eigenfunctions + ! + CALL davcio (evc, 2*nwordwfc, iunwfc, ik, - 1) + ! + ! calculate becp = , needed to compute spsi = S|psi> + ! + npw = ngk(ik) + CALL init_us_2 (npw, igk_k(1,ik), xk(1,ik), vkb) + CALL calbec ( npw, vkb, evc, becp ) + ! + ! calculate average magnetization for noncolinear case + ! + IF (noncolin) & + CALL compute_sigma_avg(sigma_avg(1,1,ik),becp%nc,ik,lsigma) + ! + IF ( ik > nks1 .AND. .NOT. no_overlap ) THEN + ! + ! compute correspondence between k+G_i indices at current and previous k + ! + ngmax = MAXVAL ( igk_k(:,ik-1:ik) ) + ALLOCATE ( work(ngmax), igg(ngmax) ) + work(:) = 0 + DO ig1 = 1, ngk(ik-1) + work( igk_k(ig1,ik-1)) = ig1 + END DO + igg(:) = 0 + DO ig2 = 1, npw + ig1 = work( igk_k(ig2,ik)) + IF (ig1 > 0) igg(ig2) = ig1 + END DO + ! + ! compute overlap <\psi_k|S\psi_{k-1}> (without the Bloch factor) + ! psi(G) = \psi(k+G) (order of G-vectors same as for previous k) + ! + psi(:,:) = (0.0_dp,0.0_dp) + DO ig2 = 1, npw + IF ( igg(ig2) > 0 ) psi(igg(ig2),:) = evc(ig2,:) + END DO + IF ( noncolin) THEN + DO ig2 = 1, npw + IF ( igg(ig2) > 0 ) psi(npwx+igg(ig2),:) = evc(npwx+ig2,:) + END DO + END IF + DEALLOCATE (igg, work) + CALL calbec (ngk(ik-1), spsi, psi, ps ) + ! + ! ps(ibnd,jbnd) = + ! + ! assign bands on the basis of the relative overlap + ! simple cases first: large or very small overlap + ! + closest_band(:,ik) = -1 + done(:) =0 + !ndone = 0 + !nlost = 0 + DO ibnd=1,nbnd + ! + psr(:) = real(ps(ibnd,:))**2+aimag(ps(ibnd,:))**2 + psmax = MAXVAL( psr ) + ! + IF ( psmax > 0.75 ) THEN + ! simple case: large overlap with one specific band + closest_band(ibnd,ik) = MAXLOC( psr, 1 ) + ! record that this band at ik has been linked to a band at ik-1 + done( closest_band(ibnd,ik) ) = 1 + ! ndone = ndone + 1 + ! + ! ELSE IF ( psmax < 0.05 ) THEN + ! ! simple case: negligible overlap with all bands + ! closest_band(ibnd,ik) = 0 + ! nlost = nlost + 1 + ! + END IF + END DO + ! + ! assign remaining bands so as to maximise overlap + ! + DO ibnd=1,nbnd + ! + ! for unassigned bands ... + ! + IF ( closest_band(ibnd,ik) == -1 ) THEN + psmax = 0.0_dp + jmax = 0 + DO jbnd = 1, nbnd + ! + ! check if this band was already assigne ... + ! + IF ( done(jbnd) > 0 ) CYCLE + pscur = real(ps(ibnd,jbnd))**2+aimag(ps(ibnd,jbnd))**2 + IF ( pscur > psmax ) THEN + psmax = pscur + jmax = jbnd + END IF + END DO + closest_band(ibnd,ik) = jmax + done(jmax) = 1 + END IF + END DO + ENDIF + ! + IF ( ik < nks2 .AND. .NOT. no_overlap ) THEN + ! + ! compute S\psi_k to be used at next k-point + ! + CALL s_psi( npwx, npw, nbnd, evc, spsi ) + ! + END IF + ! + ENDDO + ! + IF (noncolin) CALL poolrecover(sigma_avg,4*nbnd,nkstot,nks) + ! + CALL deallocate_bec_type(becp) + DEALLOCATE ( psi, spsi, ps ) + ! + IF ( ionode ) THEN + ! + ! Re-order eigenvalues according to overlap + ! (ibnd=band index, jbnd=re-ordered index) + ! + ALLOCATE (et_(nbnd,nkstot)) + DO ibnd=1,nbnd + jbnd = ibnd + DO ik=nks1tot,nks2tot + et_(ibnd,ik) = et(jbnd,ik)* rytoev + jbnd = closest_band(jbnd,ik) + END DO + END DO + ! + CALL punch_plottable_bands ( filband, nks1tot, nks2tot, nkstot, nbnd, & + xk, et_ ) + ! + DO ik=nks1tot,nks2tot + IF (ik == nks1) THEN + WRITE (iunpun, '(" &plot nbnd=",i4,", nks=",i6," /")') & + nbnd, nks2tot-nks1tot+1 + DO ipol=1,4 + IF (lsigma(ipol)) WRITE(iunpun_sigma(ipol), & + '(" &plot nbnd=",i4,", nks=",i6," /")') & + nbnd, nks2tot-nks1tot+1 + ENDDO + ENDIF + WRITE (iunpun, '(10x,3f10.6)') xk(1,ik),xk(2,ik),xk(3,ik) + WRITE (iunpun, '(10f9.3)') (et_(ibnd, ik), ibnd = 1, nbnd) + DO ipol=1,4 + IF (lsigma(ipol)) THEN + WRITE (iunpun_sigma(ipol), '(10x,3f10.6)') & + xk(1,ik),xk(2,ik),xk(3,ik) + WRITE (iunpun_sigma(ipol), '(10f9.3)') & + (sigma_avg(ipol, ibnd, ik), ibnd = 1, nbnd) + ENDIF + ENDDO + ENDDO + ! + DEALLOCATE ( et_ ) + ! + ENDIF + ! + IF (noncolin) DEALLOCATE (sigma_avg) + ! + IF ( ionode ) THEN + CLOSE (iunpun) + WRITE ( stdout, & + '(5x,"Bands written to file ",A)') TRIM(filband) + DO ipol=1,4 + IF (lsigma(ipol)) CLOSE(iunpun_sigma(ipol)) + ENDDO + ENDIF + ! + RETURN + ! +END SUBROUTINE punch_band + +SUBROUTINE punch_band_2d(filband,spin_component) +! +! This routine opens a file for each band and writes on output +! kx, ky, energy, +! kx, ky, energy +! .., .., .. +! where kx and ky are proportional to the length +! of the vectors k_1 and k_2 specified in the input of the 2d plot. +! +! The k points are supposed to be in the form +! xk(i,j) = xk_0 + dkx *(i-1) + dky * (j-1) 1eps8.OR. & + ABS(xk_collect(2,j)-xkdum(2))>eps8.OR. & + ABS(xk_collect(3,j)-xkdum(3))>eps8) THEN + n2=j-1 + dkx(:)=xk_collect(:,j)-xk0(:) + EXIT loop_k + ENDIF + ENDDO loop_k + n1=nks_eff/n2 + IF (n1*n2 /= nks_eff) CALL errore('punch_band_2d',& + 'Problems with k points',1) + mdkx = sqrt( dkx(1)**2 + dkx(2)**2 + dkx(3)**2 ) + mdky = sqrt( dky(1)**2 + dky(2)**2 + dky(3)**2 ) +! +! write the output, a band per file +! + DO ibnd=1,nbnd + filename=TRIM(filband) // '.' // TRIM(int_to_char(ibnd)) + IF (ionode) & + open(unit=iuntmp,file=filename,status='unknown', err=100, iostat=ios) + CALL mp_bcast(ios,ionode_id, intra_image_comm) +100 CALL errore('punch_band_2d','Problem opening outputfile',ios) + ijk=0 + DO i1=1,n1 + DO i2=1,n2 + ijk=ijk+1 + IF (ionode) & + WRITE(iuntmp,'(3f16.6)') mdkx*(i1-1), mdky*(i2-1), & + et_collect(ibnd,ijk)*rytoev + ENDDO + ENDDO + IF (ionode) CLOSE(unit=iuntmp,status='KEEP') + ENDDO + + DEALLOCATE(xk_collect) + DEALLOCATE(et_collect) + + RETURN + END +! +!---------------------------------------------------------------------------- +SUBROUTINE punch_plottable_bands ( filband, nks1tot, nks2tot, nkstot, nbnd, & + xk, et ) + !--------------------------------------------------------------------------- + ! + USE kinds, ONLY : dp + IMPLICIT NONE + CHARACTER(LEN=*), INTENT(IN) :: filband + INTEGER, INTENT(IN) :: nks1tot, nks2tot, nkstot, nbnd + REAL(dp), INTENT(IN) :: xk(3,nkstot), et(nbnd,nkstot) + ! + INTEGER, PARAMETER :: max_lines = 100, stdout=6, iunpun0=18 + INTEGER:: ios, i, n, nlines, npoints(max_lines), point(max_lines) + LOGICAL :: high_symmetry(nkstot), opnd + REAL(dp):: k1(3), k2(3), kx(nkstot), ps, dxmod, dxmod_save + ! + ! + IF ( nks1tot < 1 .OR. nks2tot > nkstot .OR. nkstot < 1 .OR. nbnd < 1 ) THEN + CALL infomsg('punch_plottable_bands','incorrect input data, exiting') + RETURN + END IF + ios = 0 + OPEN (unit = iunpun0, file = TRIM(filband)//'.gnu', status = 'unknown',& + form = 'formatted', iostat = ios) + IF ( ios /= 0 ) THEN + WRITE ( stdout, & + '(/,5x,"Error opening plottable file ",A)') TRIM(filband)//'.gnu' + RETURN + END IF + ! + ! Find high-symmetry points (poor man's algorithm) + ! + DO n=nks1tot,nks2tot + IF (n==nks1tot .OR. n==nks2tot) THEN + high_symmetry(n) = .true. + ELSE + k1(:) = xk(:,n) - xk(:,n-1) + k2(:) = xk(:,n+1) - xk(:,n) + IF ( k1(1)*k1(1) + k1(2)*k1(2) + k1(3)*k1(3) < 1.0d-8 .OR. & + k2(1)*k2(1) + k2(2)*k2(2) + k2(3)*k2(3) < 1.0d-8 ) THEN + CALL infomsg('punch_plottable_bands','two consecutive same k, exiting') + RETURN + END IF + ps = ( k1(1)*k2(1) + k1(2)*k2(2) + k1(3)*k2(3) ) / & + sqrt( k1(1)*k1(1) + k1(2)*k1(2) + k1(3)*k1(3) ) / & + sqrt( k2(1)*k2(1) + k2(2)*k2(2) + k2(3)*k2(3) ) + high_symmetry(n) = (ABS(ps-1.d0) >1.0d-4) + ! + ! The gamma point is a high symmetry point + ! + IF (xk(1,n)**2+xk(2,n)**2+xk(3,n)**2 < 1.0d-8) high_symmetry(n)=.true. + ! + ! save the typical length of dk + ! + IF (n==nks1tot+1) dxmod_save = sqrt( k1(1)**2 + k1(2)**2 + k1(3)**2) + ENDIF + ENDDO + + kx(nks1tot) = 0.0_dp + DO n=nks1tot+1,nks2tot + dxmod=sqrt ( (xk(1,n)-xk(1,n-1))**2 + & + (xk(2,n)-xk(2,n-1))**2 + & + (xk(3,n)-xk(3,n-1))**2 ) + IF (dxmod > 5*dxmod_save) THEN + ! + ! A big jump in dxmod is a sign that the point xk(:,n) and xk(:,n-1) + ! are quite distant and belong to two different lines. We put them on + ! the same point in the graph + ! + kx(n)=kx(n-1) + ELSEIF (dxmod > 1.d-4) THEN + ! + ! This is the usual case. The two points xk(:,n) and xk(:,n-1) are in + ! the same path. + ! + kx(n) = kx(n-1) + dxmod + dxmod_save = dxmod + ELSE + ! + ! This is the case in which dxmod is almost zero. The two points + ! coincide in the graph, but we do not save dxmod. + ! + kx(n) = kx(n-1) + dxmod + ! + ENDIF + ENDDO + ! + ! Now we compute how many paths there are: nlines + ! The first point of this path: point(iline) + ! How many points are in each path: npoints(iline) + ! + DO n=nks1tot,nks2tot + IF (high_symmetry(n)) THEN + IF (n==nks1tot) THEN + ! + ! first point. Initialize the number of lines, and the number of point + ! and say that this line start at the first point + ! + nlines=1 + npoints(1)=1 + point(1)=1 + ELSEIF (n==nks2tot) THEN + ! + ! Last point. Here we save the last point of this line, but + ! do not increase the number of lines + ! + npoints(nlines) = npoints(nlines)+1 + point(nlines+1)=n + ELSE + ! + ! Middle line. The current line has one more points, and there is + ! a new line that has to be initialized. It has one point and its + ! first point is the current k. + ! + npoints(nlines) = npoints(nlines)+1 + nlines=nlines+1 + IF (nlines>max_lines) THEN + CALL infomsg('punch_plottable_bands','too many lines, exiting') + RETURN + END IF + npoints(nlines) = 1 + point(nlines)=n + ENDIF + ! + WRITE( stdout,'(5x,"high-symmetry point: ",3f7.4,& + &" x coordinate",f9.4)') (xk(i,n),i=1,3), kx(n) + ELSE + ! + ! This k is not an high symmetry line so we just increase the number + ! of points of this line. + ! + npoints(nlines) = npoints(nlines)+1 + ENDIF + ENDDO + ! + DO i=1,nbnd + WRITE (iunpun0,'(2f10.4)') (kx(n), et(i,n),n=nks1tot,nks2tot) + WRITE (iunpun0,*) + ENDDO + ! + WRITE ( stdout, & + '(/,5x,"Plottable bands (eV) written to file ",A)') TRIM(filband)//'.gnu' + CLOSE(unit=iunpun0, STATUS='KEEP') + + RETURN + ! +END SUBROUTINE punch_plottable_bands + diff --git a/modified_bands_QE6.x/sym_band.f90 b/modified_bands_QE6.x/sym_band.f90 new file mode 100644 index 0000000..4ef6d77 --- /dev/null +++ b/modified_bands_QE6.x/sym_band.f90 @@ -0,0 +1,1738 @@ +! +! Copyright (C) 2006-2007 Quantum ESPRESSO group +! This file is distributed under the terms of the +! GNU General Public License. See the file `License' +! in the root directory of the present distribution, +! or http://www.gnu.org/copyleft/gpl.txt . +! +! +!----------------------------------------------------------------------- +SUBROUTINE sym_band(filband, spin_component, firstk, lastk) + !----------------------------------------------------------------------- + ! + USE kinds, ONLY : DP + USE ions_base, ONLY : nat, ityp, ntyp => nsp + USE cell_base, ONLY : at, bg, ibrav + USE constants, ONLY : rytoev + USE fft_base, ONLY : dfftp + USE gvect, ONLY : ngm, g + USE lsda_mod, ONLY : nspin + USE wvfct, ONLY : et, nbnd, npwx + USE klist, ONLY : xk, nks, nkstot, ngk, igk_k + USE io_files, ONLY : nwordwfc, iunwfc + USE symm_base, ONLY : s, nsym, ft, t_rev, invs, sname + USE rap_point_group, ONLY : code_group, nclass, nelem, elem, which_irr, & + char_mat, name_rap, name_class, gname, ir_ram + USE rap_point_group_so, ONLY : nrap, nelem_so, elem_so, has_e, & + which_irr_so, char_mat_so, name_rap_so, & + name_class_so, d_spin, name_class_so1 + USE rap_point_group_is, ONLY : nsym_is, sr_is, ft_is, gname_is, & + sname_is, code_group_is + USE uspp, ONLY : nkb, vkb + USE spin_orb, ONLY : domag + USE noncollin_module, ONLY : noncolin + USE wavefunctions, ONLY : evc + USE io_global, ONLY : ionode, ionode_id, stdout + USE mp, ONLY : mp_bcast + USE mp_images, ONLY : intra_image_comm + ! + IMPLICIT NONE + ! + INTEGER :: ik, i, j, irot, iclass, ig, ibnd + INTEGER :: npw, spin_component, nks1, nks2, firstk, lastk + INTEGER :: nks1tot, nks2tot + INTEGER :: iunout, igroup, irap, dim_rap, ios + INTEGER :: sk(3,3,48), gk(3,48), sk_is(3,3,48), & + gk_is(3,48), invs_is(48), t_revk(48), invsk(48), nsymk, isym, ipol, jpol + REAL(dp) :: ftk(3,48) + LOGICAL :: is_complex, is_complex_so, is_symmorphic, search_sym + LOGICAL, ALLOCATABLE :: high_symmetry(:) + REAL(DP), PARAMETER :: accuracy=1.d-4 + COMPLEX(DP) :: d_spink(2,2,48), d_spin_is(2,2,48), zdotc + COMPLEX(DP),ALLOCATABLE :: times(:,:,:) + REAL(DP) :: dxk(3), dkmod, dkmod_save, modk1, modk2, k1(3), k2(3), ps + INTEGER, ALLOCATABLE :: rap_et(:,:), code_group_k(:) + INTEGER, ALLOCATABLE :: ngroup(:), istart(:,:) + CHARACTER(len=11) :: group_name + CHARACTER(len=45) :: snamek(48) + CHARACTER (len=256) :: filband, namefile + ! + IF (spin_component/=1.and.nspin/=2) & + CALL errore('sym_band','incorrect spin_component',1) + IF (spin_component<1.or.spin_component>2) & + CALL errore('sym_band','incorrect lsda spin_component',1) + + ALLOCATE(rap_et(nbnd,nkstot)) + ALLOCATE(code_group_k(nkstot)) + ALLOCATE(times(nbnd,24,nkstot)) + ALLOCATE(ngroup(nkstot)) + ALLOCATE(istart(nbnd+1,nkstot)) + ALLOCATE(high_symmetry(nkstot)) + + code_group_k=0 + rap_et=-1 + times=(0.0_DP,0.0_DP) + CALL find_nks1nks2(firstk,lastk,nks1tot,nks1,nks2tot,nks2,spin_component) + + ios=0 + IF ( ionode ) THEN + iunout=58 + namefile=trim(filband)//".rap" + OPEN (unit = iunout, file = namefile, status = 'unknown', form = & + 'formatted', iostat = ios) + REWIND (iunout) + ENDIF + + CALL mp_bcast ( ios, ionode_id, intra_image_comm ) + IF ( ios /= 0) CALL errore ('sym_band', 'Opening filband file', abs (ios) ) + + DO ik = nks1, nks2 + ! + npw = ngk(ik) + CALL init_us_2 (npw, igk_k(1,ik), xk (1, ik), vkb) + ! + ! read eigenfunctions + ! + CALL davcio (evc, 2*nwordwfc, iunwfc, ik, - 1) + ! + ! Find the small group of k + ! + CALL smallgk (xk(1,ik), at, bg, s, ft, t_rev, sname, nsym, sk, & + ftk, gk, t_revk, invsk, snamek, nsymk) + ! + ! character of the irreducible representations + ! + CALL find_info_group(nsymk,sk,t_revk,ftk,d_spink,gk,snamek,& + sk_is,d_spin_is,gk_is,invs_is,is_symmorphic,search_sym) + code_group_k(ik)=code_group + ! + IF (.not.search_sym) THEN + rap_et(:,ik)=-1 + GOTO 100 + ENDIF + ! + ! Find the symmetry of each state + ! + IF (noncolin) THEN + IF (domag) THEN + CALL find_band_sym_so(ik,evc,et(1,ik),nsym_is, & + sk_is,ft_is,d_spin_is,gk_is,invs_is,& + rap_et(1,ik),times(1,1,ik), & + ngroup(ik),istart(1,ik),accuracy) + ELSE + CALL find_band_sym_so(ik,evc,et(1,ik),nsymk,sk,ftk,d_spink,& + gk,invsk,rap_et(1,ik),times(1,1,ik),ngroup(ik),& + istart(1,ik),accuracy) + ENDIF + ELSE + CALL find_band_sym (ik,evc, et(1,ik), nsymk, sk, ftk, gk, invsk, & + rap_et(1,ik), times(1,1,ik), ngroup(ik),& + istart(1,ik),accuracy) + ENDIF + +100 CONTINUE + ENDDO + +#ifdef __MPI + ! + ! Only the symmetry of a set of k points is calculated by this + ! processor with pool. Here we collect the results into ionode + ! + CALL ipoolrecover(code_group_k,1,nkstot,nks) + CALL ipoolrecover(rap_et,nbnd,nkstot,nks) + CALL poolrecover(times,2*24*nbnd,nkstot,nks) + CALL ipoolrecover(ngroup,1,nkstot,nks) + CALL ipoolrecover(istart,nbnd+1,nkstot,nks) +#endif + IF (ionode) THEN + high_symmetry = .FALSE. + DO ik=nks1tot,nks2tot + IF ( ik==nks1tot .OR. ik==nks2tot ) THEN + high_symmetry(ik) = .TRUE. + ELSE + k1(:) = xk(:,ik) - xk(:,ik-1) + k2(:) = xk(:,ik+1) - xk(:,ik) + modk1=sqrt( k1(1)*k1(1) + k1(2)*k1(2) + k1(3)*k1(3) ) + modk2=sqrt( k2(1)*k2(1) + k2(2)*k2(2) + k2(3)*k2(3) ) + IF (modk1 <1.d-6 .OR. modk2 < 1.d-6) CYCLE + ps = ( k1(1)*k2(1) + k1(2)*k2(2) + k1(3)*k2(3) ) / & + modk1 / modk2 + high_symmetry(ik) = (ABS(ps-1.d0) >1.0d-4) +! +! The gamma point is a high symmetry point +! + IF (xk(1,ik)**2+xk(2,ik)**2+xk(3,ik)**2 < 1.0d-9) & + high_symmetry(ik)=.TRUE. + END IF + END DO +! + DO ik=nks1tot, nks2tot + CALL smallgk (xk(1,ik), at, bg, s, ft, t_rev, sname, & + nsym, sk, ftk, gk, t_revk, invsk, snamek, nsymk) + CALL find_info_group(nsymk,sk,t_revk,ftk,d_spink,gk,snamek,& + sk_is,d_spin_is,gk_is,invs_is,is_symmorphic,search_sym) + IF (code_group_k(ik) /= code_group) & + CALL errore('sym_band','problem with code_group',1) + WRITE(stdout, '(/,1x,74("*"))') + WRITE(stdout, '(/,20x,"xk=(",2(f10.5,","),f10.5," )")') & + xk(1,ik), xk(2,ik), xk(3,ik) + IF (.not.search_sym) THEN + WRITE(stdout,'(/,5x,"zone border point and non-symmorphic group ")') + WRITE(stdout,'(5x,"symmetry decomposition not available")') + WRITE(stdout, '(/,1x,74("*"))') + ENDIF + IF (ik == nks1tot) THEN + WRITE (iunout, '(" &plot_rap nbnd_rap=",i4,", nks_rap=",i4," /")') & + nbnd, nks2tot-nks1tot+1 + IF (search_sym) CALL write_group_info(.true.) + dxk(:) = xk(:,nks1tot+1) - xk(:,nks1tot) + dkmod_save = sqrt( dxk(1)**2 + dxk(2)**2 + dxk(3)**2 ) + ELSE + IF (code_group_k(ik)/=code_group_k(ik-1).and.search_sym) & + CALL write_group_info(.true.) +! +! When the symmetry changes the point must be considered a high +! symmetry point. If the previous point was also high_symmetry, there +! are two possibilities. The two points are distant and in this case +! both of them must be considered high symmetry. If they are close only +! the first point is a high symmetry point. First compute the distance +! + dxk(:) = xk(:,ik) - xk(:,ik-1) + dkmod= sqrt( dxk(1)**2 + dxk(2)**2 + dxk(3)**2 ) + IF (dkmod < 1.D-6) THEN + ! + ! In this case is_high_sym does not change because the point + ! is the same + high_symmetry(ik)=high_symmetry(ik-1) + ! + ELSE IF (dkmod < 5.0_DP * dkmod_save) THEN +! +! In this case the two points are considered close +! + IF ( .NOT. high_symmetry(ik-1) ) & + high_symmetry(ik)= ((code_group_k(ik)/=code_group_k(ik-1)) & + .OR. high_symmetry(ik) ) + IF (dkmod > 1.d-3) dkmod_save=dkmod + ELSE +! +! Points are distant. They are all high symmetry +! + high_symmetry(ik) = .TRUE. + ENDIF + ENDIF + WRITE (iunout, '(10x,3f10.6,l5)') xk(1,ik),xk(2,ik),xk(3,ik), & + high_symmetry(ik) + WRITE (iunout, '(10i8)') (rap_et(ibnd,ik), ibnd=1,nbnd) + IF (.not.search_sym) CYCLE + IF (noncolin) THEN + IF (domag) THEN + WRITE(stdout,'(/,5x,"Band symmetry, ",a11," [",a11, & + & "] magnetic double point group,")') gname, gname_is + WRITE(stdout,'(5x,"using ",a11,/)') gname_is + ELSE + WRITE(stdout,'(/,5x,"Band symmetry, ",a11,& + & " double point group:",/)') gname + ENDIF + ELSE + WRITE(stdout,'(/,5x,"Band symmetry, ",a11," point group:",/)') & + group_name(code_group_k(ik)) + ENDIF + + DO igroup=1,ngroup(ik) + dim_rap=istart(igroup+1,ik)-istart(igroup,ik) + DO irap=1,nclass + IF (noncolin) THEN + IF ((abs(nint(dble(times(igroup,irap,ik)))- & + dble(times(igroup,irap,ik))) > accuracy).or. & + (abs(aimag(times(igroup,irap,ik))) > accuracy) ) THEN + WRITE(stdout,'(5x,"e(",i3," -",i3,") = ",f12.5,2x,& + &"eV",3x,i3,3x, "--> ?")') & + istart(igroup,ik), istart(igroup+1,ik)-1, & + et(istart(igroup,ik),ik)*rytoev, dim_rap + exit + ELSEIF (abs(times(igroup,irap,ik)) > accuracy) THEN + IF (abs(nint(dble(times(igroup,irap,ik))-1.d0)) < & + accuracy) THEN + WRITE(stdout,'(5x, "e(",i3," -",i3,") = ",& + &f12.5,2x,"eV",3x,i3,3x,"--> ",a15)') & + istart(igroup,ik), istart(igroup+1,ik)-1, & + et(istart(igroup,ik),ik)*rytoev, & + dim_rap, name_rap_so(irap) + ELSE + WRITE(stdout,'(5x,"e(",i3," -",i3,") = ",& + &f12.5,2x,"eV",3x,i3,3x,"--> ",i3," ",a15)') & + istart(igroup,ik), istart(igroup+1,ik)-1, & + et(istart(igroup,ik),ik)*rytoev, dim_rap, & + nint(dble(times(igroup,irap,ik))), name_rap_so(irap) + ENDIF + ENDIF + ELSE + IF ((abs(nint(dble(times(igroup,irap,ik)))- & + dble(times(igroup,irap,ik))) > accuracy).or. & + (abs(aimag(times(igroup,irap,ik))) > accuracy) ) THEN + WRITE(stdout,'(5x,"e(",i3," -",i3,") = ",f12.5,2x,& + &"eV",3x,i3,3x, "--> ?")') & + istart(igroup,ik), istart(igroup+1,ik)-1, & + et(istart(igroup,ik),ik)*rytoev, dim_rap + exit + ELSEIF (abs(times(igroup,irap,ik)) > accuracy) THEN + IF (abs(nint(dble(times(igroup,irap,ik))-1.d0)) < & + accuracy) THEN + WRITE(stdout,'(5x, "e(",i3," -",i3,") = ",& + &f12.5,2x,"eV",3x,i3,3x,"--> ",a15)') & + istart(igroup,ik), istart(igroup+1,ik)-1, & + et(istart(igroup,ik),ik)*rytoev, & + dim_rap, name_rap(irap) + ELSE + WRITE(stdout,'(5x,"e(",i3," -",i3,") = ",& + &f12.5,2x,"eV",3x,i3,3x,"--> ",i3," ",a15)') & + istart(igroup,ik), istart(igroup+1,ik)-1, & + et(istart(igroup,ik),ik)*rytoev, dim_rap, & + nint(dble(times(igroup,irap,ik))), name_rap(irap) + ENDIF + ENDIF + ENDIF + ENDDO + ENDDO + WRITE( stdout, '(/,1x,74("*"))') + ENDDO + CLOSE(iunout) + ENDIF + + ! + DEALLOCATE(times) + DEALLOCATE(code_group_k) + DEALLOCATE(rap_et) + DEALLOCATE(ngroup) + DEALLOCATE(istart) + DEALLOCATE(high_symmetry) + ! + RETURN +END SUBROUTINE sym_band +! +SUBROUTINE parity_band(filband, spin_component, firstk, lastk) + !----------------------------------------------------------------------- + ! + USE kinds, ONLY : DP + USE ions_base, ONLY : nat, ityp, ntyp => nsp + USE cell_base, ONLY : at, bg, ibrav + USE constants, ONLY : rytoev + USE fft_base, ONLY : dfftp + USE gvect, ONLY : ngm, g + USE lsda_mod, ONLY : nspin + USE wvfct, ONLY : et, nbnd, npwx + USE klist, ONLY : xk, nks, nkstot, ngk, igk_k + USE io_files, ONLY : nwordwfc, iunwfc + USE symm_base, ONLY : s, nsym, ft, t_rev, invs, sname + USE rap_point_group, ONLY : code_group, nclass, nelem, elem, which_irr, & + char_mat, name_rap, name_class, gname, ir_ram + USE rap_point_group_so, ONLY : nrap, nelem_so, elem_so, has_e, & + which_irr_so, char_mat_so, name_rap_so, & + name_class_so, d_spin, name_class_so1 + USE rap_point_group_is, ONLY : nsym_is, sr_is, ft_is, gname_is, & + sname_is, code_group_is + USE uspp, ONLY : nkb, vkb + USE spin_orb, ONLY : domag + USE noncollin_module, ONLY : noncolin + USE wavefunctions, ONLY : evc + USE io_global, ONLY : ionode, ionode_id, stdout + USE mp, ONLY : mp_bcast + USE mp_images, ONLY : intra_image_comm + ! + ! + IMPLICIT NONE + ! + INTEGER :: npw, ik, i, j, irot, iclass, ig, ibnd + INTEGER :: spin_component, nks1, nks2, firstk, lastk + INTEGER :: nks1tot, nks2tot + INTEGER :: iunout, igroup, irap, dim_rap, ios + INTEGER :: sk(3,3,48), gk(3,48), sk_is(3,3,48), & + gk_is(3,48), invs_is(48), t_revk(48), invsk(48), nsymk, isym, ipol, jpol + REAL(dp) :: ftk(3,48) + LOGICAL :: is_complex, is_complex_so, is_symmorphic, search_sym + LOGICAL :: invsymk + LOGICAL, ALLOCATABLE :: high_symmetry(:) + REAL(DP), PARAMETER :: accuracy=1.d-4 + COMPLEX(DP) :: d_spink(2,2,48), d_spin_is(2,2,48), zdotc + COMPLEX(DP),ALLOCATABLE :: times(:,:,:) + REAL(DP) :: dxk(3), dkmod, dkmod_save, modk1, modk2, k1(3), k2(3), ps + INTEGER, ALLOCATABLE :: rap_par(:,:), code_group_k(:) + INTEGER, ALLOCATABLE :: ngroup(:), istart(:,:) + CHARACTER(len=11) :: group_name + CHARACTER(len=45) :: snamek(48) + CHARACTER (len=256) :: filband, namefile + ! + IF (spin_component/=1.and.nspin/=2) & + CALL errore('sym_band','incorrect spin_component',1) + IF (spin_component<1.or.spin_component>2) & + CALL errore('sym_band','incorrect lsda spin_component',1) + + ALLOCATE(rap_par(nbnd,nkstot)) + ALLOCATE(code_group_k(nkstot)) + ALLOCATE(times(nbnd,24,nkstot)) + ALLOCATE(ngroup(nkstot)) + ALLOCATE(istart(nbnd+1,nkstot)) + ALLOCATE(high_symmetry(nkstot)) + + + + IF ( ionode) THEN + write(stdout,*) "gnne" + ENDIF + + code_group_k=0 + rap_par=0 + times=(0.0_DP,0.0_DP) + CALL find_nks1nks2(firstk,lastk,nks1tot,nks1,nks2tot,nks2,spin_component) + + ios=0 + IF ( ionode ) THEN + iunout=59 + namefile=trim(filband)//".par" + OPEN (unit = iunout, file = namefile, status = 'unknown', form = & + 'formatted', iostat = ios) + REWIND (iunout) + ENDIF + + + CALL mp_bcast ( ios, ionode_id, intra_image_comm ) + IF ( ios /= 0) CALL errore ('sym_band', 'Opening filband file', abs (ios) ) + + DO ik = nks1, nks2 + ! + ! prepare the indices of this k point + ! + ! CALL gk_sort (xk (1, ik), ngm, g, ecutwfc / tpiba2, npw, & + ! igk, g2kin) + ! + npw = ngk(ik) + CALL init_us_2 (npw, igk_k(1,ik), xk (1, ik), vkb) + ! + ! read eigenfunctions + ! + CALL davcio (evc, 2*nwordwfc, iunwfc, ik, - 1) + ! + ! Find the small group of k + ! + CALL smallgk (xk(1,ik), at, bg, s, ft, t_rev, sname, nsym, sk, & + ftk, gk, t_revk, invsk, snamek, nsymk) + + + invsymk = all ( sk(:,:,nsymk/2+1) == -sk(:,:,1) ) + + + CALL find_info_group(nsymk,sk,t_revk,ftk,d_spink,gk,snamek,& + sk_is,d_spin_is,gk_is,invs_is,is_symmorphic,search_sym) + code_group_k(ik)=code_group + + IF (invsymk) THEN + CALL find_parity_sym_so(ik, evc,et(1,ik),nsymk,sk,ftk,d_spink,& + gk,invsk,rap_par(1,ik),times(1,1,ik),ngroup(ik),& + istart(1,ik),accuracy) + ENDIF + + + ENDDO + +#ifdef __MPI + ! + ! Only the symmetry of a set of k points is calculated by this + ! processor with pool. Here we collect the results into ionode + ! + CALL ipoolrecover(code_group_k,1,nkstot,nks) + CALL ipoolrecover(rap_par,nbnd,nkstot,nks) + CALL poolrecover(times,2*24*nbnd,nkstot,nks) + CALL ipoolrecover(ngroup,1,nkstot,nks) + CALL ipoolrecover(istart,nbnd+1,nkstot,nks) +#endif + + IF (ionode) THEN + high_symmetry = .FALSE. + + DO ik=nks1tot,nks2tot + IF ( ik==nks1tot .OR. ik==nks2tot ) THEN + high_symmetry(ik) = .TRUE. + ELSE + k1(:) = xk(:,ik) - xk(:,ik-1) + k2(:) = xk(:,ik+1) - xk(:,ik) + modk1=sqrt( k1(1)*k1(1) + k1(2)*k1(2) + k1(3)*k1(3) ) + modk2=sqrt( k2(1)*k2(1) + k2(2)*k2(2) + k2(3)*k2(3) ) + IF (modk1 <1.d-6 .OR. modk2 < 1.d-6) CYCLE + ps = ( k1(1)*k2(1) + k1(2)*k2(2) + k1(3)*k2(3) ) / & + modk1 / modk2 + high_symmetry(ik) = (ABS(ps-1.d0) >1.0d-4) +! +! The gamma point is a high symmetry point +! + IF (xk(1,ik)**2+xk(2,ik)**2+xk(3,ik)**2 < 1.0d-9) & + high_symmetry(ik)=.TRUE. + END IF + END DO +! + DO ik=nks1tot, nks2tot + CALL smallgk (xk(1,ik), at, bg, s, ft, t_rev, sname, & + nsym, sk, ftk, gk, t_revk, invsk, snamek, nsymk) + CALL find_info_group(nsymk,sk,t_revk,ftk,d_spink,gk,snamek,& + sk_is,d_spin_is,gk_is, & + invs_is,is_symmorphic,search_sym) + IF (code_group_k(ik) /= code_group) & + CALL errore('sym_band','problem with code_group',1) + WRITE(stdout, '(/,1x,74("*"))') + WRITE(stdout, '(/,20x,"xk=(",2(f10.5,","),f10.5," )")') & + xk(1,ik), xk(2,ik), xk(3,ik) + + IF (ik == nks1tot) THEN + WRITE (iunout, '(" &plot_par nbnd_par=",i4,", nks_par=",i4," /")') & + nbnd, nks2tot-nks1tot+1 + IF (search_sym) CALL write_group_info(.true.) + dxk(:) = xk(:,nks1tot+1) - xk(:,nks1tot) + dkmod_save = sqrt( dxk(1)**2 + dxk(2)**2 + dxk(3)**2 ) + ELSE + IF (code_group_k(ik)/=code_group_k(ik-1).and.search_sym) & + CALL write_group_info(.true.) + + + dxk(:) = xk(:,ik) - xk(:,ik-1) + dkmod= sqrt( dxk(1)**2 + dxk(2)**2 + dxk(3)**2 ) + IF (dkmod < 1.D-6) THEN + ! + ! In this case is_high_sym does not change because the point + ! is the same + high_symmetry(ik)=high_symmetry(ik-1) + ! + ELSE IF (dkmod < 5.0_DP * dkmod_save) THEN +! +! In this case the two points are considered close +! + IF ( .NOT. high_symmetry(ik-1) ) & + high_symmetry(ik)= ((code_group_k(ik)/=code_group_k(ik-1)) & + .OR. high_symmetry(ik) ) + IF (dkmod > 1.d-3) dkmod_save=dkmod + ELSE +! +! Points are distant. They are all high symmetry +! + high_symmetry(ik) = .TRUE. + ENDIF + ENDIF + !HERE TO CHANGE + WRITE (iunout, '(10x,3f10.6,l5)') xk(1,ik),xk(2,ik),xk(3,ik), & + high_symmetry(ik) + WRITE (iunout, '(10i8)') (rap_par(ibnd,ik), ibnd=1,nbnd) + + + + ENDDO + CLOSE(iunout) + ENDIF + + ! + DEALLOCATE(times) + DEALLOCATE(code_group_k) + DEALLOCATE(rap_par) + DEALLOCATE(ngroup) + DEALLOCATE(istart) + DEALLOCATE(high_symmetry) + ! + RETURN +END SUBROUTINE parity_band +! +SUBROUTINE find_parity_sym_so (ik,evc,et,nsym,s,ft,d_spin,gk, & + invs,rap_par,times,ngroup,istart,accuracy) + + ! + ! This subroutine finds the irreducible representations of the + ! double group which give the transformation properties of the + ! spinor wavefunctions evc. + ! Presently it does NOT work at zone border if the space group of + ! the crystal has fractionary translations (non-symmorphic space groups). + ! + ! + USE io_global, ONLY : stdout + USE kinds, ONLY : DP + USE constants, ONLY : rytoev + USE fft_base, ONLY : dfftp + USE rap_point_group, ONLY : code_group, nclass, gname + USE rap_point_group_so, ONLY : nrap, nelem_so, elem_so, has_e, which_irr_so, & + char_mat_so, name_rap_so, name_class_so, & + name_class_so1 + USE rap_point_group_is, ONLY : gname_is + USE gvect, ONLY : ngm + USE wvfct, ONLY : nbnd, npwx + USE klist, ONLY : ngk + USE uspp, ONLY : vkb, nkb, okvan + USE noncollin_module, ONLY : npol + USE becmod, ONLY : bec_type, becp, calbec, allocate_bec_type, deallocate_bec_type + USE mp_bands, ONLY : intra_bgrp_comm + USE mp, ONLY : mp_sum + + IMPLICIT NONE + + REAL(DP), INTENT(in) :: accuracy + + INTEGER, INTENT(in) :: ik + + INTEGER :: ftau(3) + INTEGER :: & + nsym, & + ngroup, & + istart(nbnd+1), & + rap_par(nbnd), & + gk(3,48), & + invs(48), & + s(3,3,48) + + REAL(DP) :: & + ft(3,48), & + et(nbnd) + + COMPLEX(DP) :: & + times(nbnd,24), & + d_spin(2,2,48), & + evc(npwx*npol, nbnd) + + REAL(DP), PARAMETER :: eps=1.d-5 + + INTEGER :: & + ibnd, & + igroup, & + dim_rap, & ! counters + irot, & + irap, & + shift, & + iclass, & + iclassI, & + na, i, j, ig, ipol, jpol, jrap, dimen, npw + + COMPLEX(DP) :: zdotc ! moltiplication factors + + ! INTEGER, ALLOCATABLE :: parity(:) + REAL(DP), ALLOCATABLE :: w1(:) ! list of energy eigenvalues in eV + COMPLEX(DP), ALLOCATABLE :: evcr(:,:), & ! the rotated of each wave function + trace(:,:) ! the trace of the symmetry matrix + ! within a given group + ! + ! Divide the bands on the basis of the band degeneracy. + ! + !ALLOCATE(parity(nbnd)) + ALLOCATE(w1(nbnd)) + ALLOCATE(evcr(npwx*npol,nbnd)) + ALLOCATE(trace(48,nbnd)) + IF (okvan) CALL allocate_bec_type ( nkb, nbnd, becp ) + + rap_par=0 + w1=et*rytoev + ! + ! divide the energies in groups of degenerate eigenvalues. Two eigenvalues + ! are assumed to be degenerate if their difference is less than 0.0001 eV. + ! + ngroup=1 + istart(ngroup)=1 + DO ibnd=2,nbnd + IF (abs(w1(ibnd)-w1(ibnd-1)) > 0.0001d0) THEN + ngroup=ngroup+1 + istart(ngroup)=ibnd + ENDIF + ENDDO + istart(ngroup+1)=nbnd+1 + + ! + ! Find the character of one symmetry operation per class + ! + trace=(0.d0,0.d0) + + iclassI = -1 + DO iclass=1,nclass + irot=elem_so(1,iclass) + + IF ( all(s(:,:,irot) == -s(:,:,1)) ) THEN + IF ( has_e(1,iclass) == 1 ) THEN + + ftau(1) = NINT(ft(1,invs(irot))*dfftp%nr1) + ftau(2) = NINT(ft(2,invs(irot))*dfftp%nr2) + ftau(3) = NINT(ft(3,invs(irot))*dfftp%nr3) + CALL rotate_all_psi_so(ik,evc,evcr,s(1,1,invs(irot)), & + ftau,d_spin(1,1,irot),has_e(1,iclass),gk(1,invs(irot))) + ! + ! and apply S in the US case. + ! + npw = ngk(ik) + IF ( okvan ) THEN + CALL calbec( npw, vkb, evcr, becp ) + CALL s_psi( npwx, npw, nbnd, evcr, evcr ) + ENDIF + + ! CALL rotate_all_psi_so(evc,evcr,s(1,1,irot), & + ! ftau(1,irot),d_spin(1,1,irot),has_e(1,iclass),gk(1,irot)) + + iclassI = iclass + ! + ! and apply S in the US case. + ! + ! Compute the trace of the representation for each group of bands + ! + + DO ibnd =1,nbnd + trace(iclass,ibnd)=trace(iclass,ibnd) + & + zdotc(2*npwx,evc(1,ibnd),1,evcr(1,ibnd),1) + ENDDO + ! DO igroup=1,ngroup + ! dim_rap=istart(igroup+1)-istart(igroup) + ! DO i=1,dim_rap + ! ibnd=istart(igroup)+i-1 + ! trace(iclass,igroup)=trace(iclass,igroup) + & + ! zdotc(2*npwx,evc(1,ibnd),1,evcr(1,ibnd),1) + ! ENDDO + ! ! write(6,*) igroup, iclass, dim_rap, trace(iclass,igroup) + ! ENDDO + + ENDIF + ENDIF + ENDDO + + CALL mp_sum(trace,intra_bgrp_comm) + + ! WRITE (stdout, '(/,1x,i8,74("*"))') ik + ! WRITE (stdout, '(10x,10f9.5)') (REAL(trace(iclassI,ibnd)), ibnd=1,nbnd) + + DO ibnd =1,nbnd + IF (iclassI == -1) THEN + rap_par(ibnd) = 0 + ELSE + IF (( ABS( 1 - ABS(REAL(trace(iclassI,ibnd))))) < 1.0d-5 ) THEN + rap_par(ibnd) = NINT(REAL(trace(iclassI,ibnd))) + ELSE + rap_par(ibnd) = 0 + ENDIF + ENDIF + ENDDO + + ! DO igroup=1,ngroup + ! dim_rap=istart(igroup+1)-istart(igroup) + ! DO i=1,dim_rap + ! ibnd=istart(igroup)+i-1 + ! IF (( ABS( 1 - ABS(REAL(trace(iclassI,igroup))))) < 1.0d-5 ) THEN + ! rap_par(ibnd) = NINT(REAL(trace(iclassI,igroup))) + ! ELSE + ! rap_par(ibnd) = 0 + ! ENDIF + ! ENDDO + ! ! write(6,*) igroup, iclass, dim_rap, trace(iclass,igroup) + ! ENDDO + + + !DEALLOCATE(parity) + DEALLOCATE(trace) + DEALLOCATE(w1) + DEALLOCATE(evcr) + IF (okvan) CALL deallocate_bec_type ( becp ) + RETURN +END SUBROUTINE find_parity_sym_so +! +SUBROUTINE find_band_sym (ik,evc,et,nsym,s,ft,gk,invs,rap_et,times,ngroup,& + istart,accuracy) + ! + ! This subroutine finds the irreducible representations which give + ! the transformation properties of the wavefunctions. + ! Presently it does NOT work at zone border if the space group of + ! the crystal has fractionary translations (non-symmorphic space groups). + ! + ! + USE io_global, ONLY : stdout + USE kinds, ONLY : DP + USE constants, ONLY : rytoev + USE rap_point_group, ONLY : code_group, nclass, nelem, elem, which_irr, & + char_mat, name_rap, name_class, gname + USE gvect, ONLY : ngm + USE wvfct, ONLY : nbnd, npwx + USE klist, ONLY : ngk, igk_k + USE uspp, ONLY : vkb, nkb, okvan + USE becmod, ONLY : bec_type, becp, calbec, & + allocate_bec_type, deallocate_bec_type + USE fft_base, ONLY : dfftp + USE fft_interfaces, ONLY : invfft + USE mp_bands, ONLY : intra_bgrp_comm + USE mp, ONLY : mp_sum + + IMPLICIT NONE + + INTEGER, INTENT(in) :: ik + REAL(DP), INTENT(in) :: accuracy + + INTEGER :: & + nsym, & + rap_et(nbnd), & + gk(3,48), & + s(3,3,48), & + invs(48), & + ngroup, & ! number of different frequencies groups + istart(nbnd+1) + + REAL(DP) :: & + ft(3,48), & + et(nbnd) + + COMPLEX(DP) :: & + times(nbnd,24), & + evc(npwx, nbnd) + + REAL(DP), PARAMETER :: eps=1.d-5 + + INTEGER :: & + ibnd, & + igroup, & + dim_rap, & + irot, & + irap, & + iclass, & + shift, & + na, i, j, ig, dimen, nrxx, npw + INTEGER :: ftau(3) + + COMPLEX(DP) :: zdotc + + REAL(DP), ALLOCATABLE :: w1(:) + COMPLEX(DP), ALLOCATABLE :: evcr(:,:), trace(:,:), psic(:,:) + ! + ! Divide the bands on the basis of the band degeneracy. + ! + nrxx=dfftp%nnr + ALLOCATE(w1(nbnd)) + ALLOCATE(evcr(npwx,nbnd)) + ALLOCATE(psic(nrxx,nbnd)) + ALLOCATE(trace(48,nbnd)) + IF (okvan) CALL allocate_bec_type ( nkb, nbnd, becp ) + + rap_et=-1 + w1=et*rytoev + + ngroup=1 + istart(ngroup)=1 + DO ibnd=2,nbnd + IF (abs(w1(ibnd)-w1(ibnd-1)) > 0.0001d0) THEN + ngroup=ngroup+1 + istart(ngroup)=ibnd + ENDIF + ENDDO + istart(ngroup+1)=nbnd+1 +! +! bring all the bands in real space +! + npw = ngk(ik) + psic=(0.0_DP,0.0_DP) + DO ibnd=1,nbnd + psic(dfftp%nl(igk_k(1:npw,ik)),ibnd) = evc(1:npw,ibnd) + CALL invfft ('Rho', psic(:,ibnd), dfftp) + ENDDO + ! + ! Find the character of one symmetry operation per class + ! + DO iclass=1,nclass + irot=elem(1,iclass) + ! + ! Rotate all the bands together. + ! NB: rotate_psi assume that s is in the small group of k. It does not + ! rotate the k point. + ! + ! + IF (irot==1) THEN + evcr=evc + ELSE + ftau(1) = NINT(ft(1,invs(irot))*dfftp%nr1) + ftau(2) = NINT(ft(2,invs(irot))*dfftp%nr2) + ftau(3) = NINT(ft(3,invs(irot))*dfftp%nr3) + CALL rotate_all_psi(ik,psic,evcr,s(1,1,invs(irot)), & + ftau,gk(1,invs(irot))) + ENDIF + ! + ! and apply S if necessary + ! + IF ( okvan ) THEN + CALL calbec( npw, vkb, evcr, becp ) + CALL s_psi( npwx, npw, nbnd, evcr, evcr ) + ENDIF + ! + ! Compute the trace of the representation for each group of bands + ! + DO igroup=1,ngroup + dim_rap=istart(igroup+1)-istart(igroup) + trace(iclass,igroup)=(0.d0,0.d0) + DO i=1,dim_rap + ibnd=istart(igroup)+i-1 + trace(iclass,igroup)=trace(iclass,igroup) + & + zdotc(npw,evc(1,ibnd),1,evcr(1,ibnd),1) + ENDDO + ! write(6,*) igroup, iclass, trace(iclass,igroup) + ENDDO + ENDDO + ! + CALL mp_sum( trace, intra_bgrp_comm ) + + !DO iclass=1,nclass + ! write(6,'(i5,3(2f11.8,1x))') iclass,trace(iclass,4),trace(iclass,5), & + ! trace(iclass,6) + !ENDDO + + ! + ! And now use the character table to identify the symmetry representation + ! of each group of bands + ! + !WRITE(stdout,'(/,5x,"Band symmetry, ",a11," point group:",/)') gname + + DO igroup=1,ngroup + dim_rap=istart(igroup+1)-istart(igroup) + shift=0 + DO irap=1,nclass + times(igroup,irap)=(0.d0,0.d0) + DO iclass=1,nclass + times(igroup,irap)=times(igroup,irap) & + +trace(iclass,igroup)*CONJG(char_mat(irap,which_irr(iclass)))& + *nelem(iclass) + ENDDO + times(igroup,irap)=times(igroup,irap)/nsym + IF ((abs(nint(dble(times(igroup,irap)))-dble(times(igroup,irap))) & + > accuracy).or. (abs(aimag(times(igroup,irap))) > eps) ) THEN + ! WRITE(stdout,'(5x,"e(",i3," -",i3,") = ",f12.5,2x,"eV",3x,i3,3x,& + ! & "--> ?")') & + ! istart(igroup), istart(igroup+1)-1, w1(istart(igroup)), dim_rap + ibnd=istart(igroup) + IF (rap_et(ibnd)==-1) THEN + DO i=1,dim_rap + ibnd=istart(igroup)+i-1 + rap_et(ibnd)=0 + ENDDO + ENDIF + GOTO 300 + ELSEIF (abs(times(igroup,irap)) > accuracy) THEN + ibnd=istart(igroup)+shift + dimen=nint(dble(char_mat(irap,1))) + IF (rap_et(ibnd)==-1) THEN + DO i=1,dimen*nint(dble(times(igroup,irap))) + ibnd=istart(igroup)+shift+i-1 + rap_et(ibnd)=irap + ENDDO + shift=shift+dimen*nint(dble(times(igroup,irap))) + ENDIF + ! IF (ABS(NINT(DBLE(times))-1.d0) < 1.d-4) THEN + ! WRITE(stdout,'(5x, "e(",i3," -",i3,") = ",f12.5,2x,"eV",3x,i3,& + ! & 3x,"--> ",a15)') & + ! istart(igroup), istart(igroup+1)-1, w1(istart(igroup)), & + ! dim_rap, name_rap(irap) + ! ELSE + ! WRITE(stdout,'(5x,"e(",i3," -",i3,") = ",f12.5,2x,"eV",3x,i3,3x,& + ! & "--> ",i3," ",a15)') & + ! istart(igroup), istart(igroup+1)-1, & + ! w1(istart(igroup)), dim_rap, NINT(DBLE(times)), name_rap(irap) + ! END IF + ENDIF + ENDDO +300 CONTINUE + ENDDO + !WRITE( stdout, '(/,1x,74("*"))') + + DEALLOCATE(trace) + DEALLOCATE(w1) + DEALLOCATE(evcr) + DEALLOCATE(psic) + IF (okvan) CALL deallocate_bec_type (becp) + + RETURN +END SUBROUTINE find_band_sym + + +SUBROUTINE rotate_all_psi(ik,psic,evcr,s,ftau,gk) + + USE kinds, ONLY : DP + USE constants, ONLY : tpi + USE gvect, ONLY : ngm + USE wvfct, ONLY : nbnd, npwx + USE klist, ONLY : ngk, igk_k + USE fft_base, ONLY : dfftp + USE scatter_mod, ONLY : cgather_sym_many, cscatter_sym_many + USE fft_interfaces, ONLY : fwfft, invfft + USE mp_bands, ONLY : intra_bgrp_comm + USE mp, ONLY : mp_sum + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: ik + INTEGER :: s(3,3), ftau(3), gk(3) + + COMPLEX(DP), ALLOCATABLE :: psir(:) + COMPLEX(DP) :: psic(dfftp%nnr,nbnd), evcr(npwx,nbnd) + COMPLEX(DP) :: phase + REAL(DP) :: arg + INTEGER :: i, j, k, ri, rj, rk, ir, rir, ipol, ibnd + INTEGER :: nr1, nr2, nr3, nr1x, nr2x, nr3x, nrxx, npw + LOGICAL :: zone_border + INTEGER :: start_band, last_band, my_nbnd_proc + INTEGER :: start_band_proc(dfftp%nproc), nbnd_proc(dfftp%nproc) + +#if defined (__MPI) + ! + COMPLEX (DP), ALLOCATABLE :: psir_collect(:) + COMPLEX (DP), ALLOCATABLE :: psic_collect(:,:) + ! +#endif + ! + nr1=dfftp%nr1 + nr2=dfftp%nr2 + nr3=dfftp%nr3 + nr1x=dfftp%nr1x + nr2x=dfftp%nr2x + nr3x=dfftp%nr3x + nrxx=dfftp%nnr + npw = ngk(ik) + ! + ALLOCATE(psir(nrxx)) + ! + zone_border=(gk(1)/=0.OR.gk(2)/=0.OR.gk(3)/=0) + ! + evcr= (0.0_DP, 0.0_DP) + ! +#if defined (__MPI) + + call divide (intra_bgrp_comm, nbnd, start_band, last_band) + start_band_proc=0 + start_band_proc(dfftp%mype+1)=start_band + nbnd_proc=0 + my_nbnd_proc=last_band-start_band+1 + nbnd_proc(dfftp%mype+1)=my_nbnd_proc + CALL mp_sum(start_band_proc, intra_bgrp_comm) + CALL mp_sum(nbnd_proc, intra_bgrp_comm) + ! + ALLOCATE (psic_collect(nr1x*nr2x*nr3x, my_nbnd_proc)) + ALLOCATE (psir_collect(nr1x*nr2x*nr3x)) + ! + CALL cgather_sym_many( dfftp, psic, psic_collect, nbnd, nbnd_proc, start_band_proc) + ! + DO ibnd = 1, my_nbnd_proc + psir_collect=(0.d0,0.d0) + ! + IF (zone_border) THEN + DO k = 1, nr3 + DO j = 1, nr2 + DO i = 1, nr1 + CALL ruotaijk (s, ftau, i, j, k, nr1, nr2, nr3, ri, rj, rk ) + ir=i+(j-1)*nr1x+(k-1)*nr1x*nr2x + rir=ri+(rj-1)*nr1x+(rk-1)*nr1x*nr2x + arg=tpi*( (gk(1)*(i-1))/dble(nr1)+(gk(2)*(j-1))/dble(nr2)+ & + (gk(3)*(k-1))/dble(nr3) ) + phase=cmplx(cos(arg),sin(arg),kind=DP) + psir_collect(ir)=psic_collect(rir,ibnd)*phase + ENDDO + ENDDO + ENDDO + ELSE + DO k = 1, nr3 + DO j = 1, nr2 + DO i = 1, nr1 + CALL ruotaijk (s, ftau, i, j, k, nr1, nr2, nr3, ri, rj, rk ) + ir=i+(j-1)*nr1x+(k-1)*nr1x*nr2x + rir=ri+(rj-1)*nr1x+(rk-1)*nr1x*nr2x + psir_collect(ir)=psic_collect(rir, ibnd) + ENDDO + ENDDO + ENDDO + ENDIF + psic_collect(:,ibnd)=psir_collect(:) + ENDDO + ! + DO ibnd=1, nbnd + CALL cscatter_sym_many( dfftp, psic_collect, psir, ibnd, nbnd, & + nbnd_proc, start_band_proc ) + ! + CALL fwfft ('Rho', psir, dfftp) + ! + evcr(1:npw,ibnd) = psir(dfftp%nl(igk_k(1:npw,ik))) + END DO + DEALLOCATE (psic_collect) + DEALLOCATE (psir_collect) + ! +#else + psir=(0.d0,0.d0) + DO ibnd=1,nbnd + IF (zone_border) THEN + DO k = 1, nr3 + DO j = 1, nr2 + DO i = 1, nr1 + CALL ruotaijk (s, ftau, i, j, k, nr1, nr2, nr3, ri, rj, rk ) + ir=i+(j-1)*nr1x+(k-1)*nr1x*nr2x + rir=ri+(rj-1)*nr1x+(rk-1)*nr1x*nr2x + arg=tpi*( (gk(1)*(i-1))/dble(nr1)+(gk(2)*(j-1))/dble(nr2)+ & + (gk(3)*(k-1))/dble(nr3) ) + phase=cmplx(cos(arg),sin(arg),kind=DP) + psir(ir)=psic(rir,ibnd)*phase + ENDDO + ENDDO + ENDDO + ELSE + DO k = 1, nr3 + DO j = 1, nr2 + DO i = 1, nr1 + CALL ruotaijk (s, ftau, i, j, k, nr1, nr2, nr3, ri, rj, rk ) + ir=i+(j-1)*nr1x+(k-1)*nr1x*nr2x + rir=ri+(rj-1)*nr1x+(rk-1)*nr1x*nr2x + psir(ir)=psic(rir,ibnd) + ENDDO + ENDDO + ENDDO + ENDIF + CALL fwfft ('Rho', psir, dfftp) + ! + evcr(1:npw,ibnd) = psir(dfftp%nl(igk_k(1:npw,ik))) + ENDDO + ! +#endif + ! + DEALLOCATE(psir) + ! + RETURN +END SUBROUTINE rotate_all_psi + +SUBROUTINE find_band_sym_so (ik,evc,et,nsym,s,ft,d_spin,gk, & + invs,rap_et,times,ngroup,istart,accuracy) + + ! + ! This subroutine finds the irreducible representations of the + ! double group which give the transformation properties of the + ! spinor wavefunctions evc. + ! Presently it does NOT work at zone border if the space group of + ! the crystal has fractionary translations (non-symmorphic space groups). + ! + ! + USE io_global, ONLY : stdout + USE kinds, ONLY : DP + USE constants, ONLY : rytoev + USE fft_base, ONLY : dfftp + USE rap_point_group, ONLY : code_group, nclass, gname + USE rap_point_group_so, ONLY : nrap, nelem_so, elem_so, has_e, which_irr_so, & + char_mat_so, name_rap_so, name_class_so, & + name_class_so1 + USE rap_point_group_is, ONLY : gname_is + USE gvect, ONLY : ngm + USE wvfct, ONLY : nbnd, npwx + USE klist, ONLY : ngk + USE uspp, ONLY : vkb, nkb, okvan + USE noncollin_module, ONLY : npol + USE becmod, ONLY : bec_type, becp, calbec, allocate_bec_type, deallocate_bec_type + USE mp_bands, ONLY : intra_bgrp_comm + USE mp, ONLY : mp_sum + + IMPLICIT NONE + + INTEGER, INTENT(in) :: ik + REAL(DP), INTENT(in) :: accuracy + + INTEGER :: & + nsym, & + ngroup, & + istart(nbnd+1), & + rap_et(nbnd), & + gk(3,48), & + invs(48), & + s(3,3,48) + + REAL(DP) :: & + ft(3,48), & + et(nbnd) + + COMPLEX(DP) :: & + times(nbnd,24), & + d_spin(2,2,48), & + evc(npwx*npol, nbnd) + + REAL(DP), PARAMETER :: eps=1.d-5 + + INTEGER :: ftau(3) + INTEGER :: & + ibnd, & + igroup, & + dim_rap, & ! counters + irot, & + irap, & + shift, & + iclass, & + na, i, j, ig, ipol, jpol, jrap, dimen, npw + + COMPLEX(DP) :: zdotc ! moltiplication factors + + REAL(DP), ALLOCATABLE :: w1(:) ! list of energy eigenvalues in eV + COMPLEX(DP), ALLOCATABLE :: evcr(:,:), & ! the rotated of each wave function + trace(:,:) ! the trace of the symmetry matrix + ! within a given group + ! + ! Divide the bands on the basis of the band degeneracy. + ! + ALLOCATE(w1(nbnd)) + ALLOCATE(evcr(npwx*npol,nbnd)) + ALLOCATE(trace(48,nbnd)) + IF (okvan) CALL allocate_bec_type ( nkb, nbnd, becp ) + + rap_et=-1 + w1=et*rytoev + ! + ! divide the energies in groups of degenerate eigenvalues. Two eigenvalues + ! are assumed to be degenerate if their difference is less than 0.0001 eV. + ! + ngroup=1 + istart(ngroup)=1 + DO ibnd=2,nbnd + IF (abs(w1(ibnd)-w1(ibnd-1)) > 0.0001d0) THEN + ngroup=ngroup+1 + istart(ngroup)=ibnd + ENDIF + ENDDO + istart(ngroup+1)=nbnd+1 + + trace=(0.d0,0.d0) + DO iclass=1,nclass + irot=elem_so(1,iclass) + ! + ! Rotate all the bands together. + ! NB: rotate_psi assumes that s is in the small group of k. It does not + ! rotate the k point. + ! + ftau(1) = NINT(ft(1,invs(irot))*dfftp%nr1) + ftau(2) = NINT(ft(2,invs(irot))*dfftp%nr2) + ftau(3) = NINT(ft(3,invs(irot))*dfftp%nr3) + CALL rotate_all_psi_so(ik,evc,evcr,s(1,1,invs(irot)), & + ftau,d_spin(1,1,irot),has_e(1,iclass),gk(1,invs(irot))) + ! + ! and apply S in the US case. + ! + npw = ngk(ik) + IF ( okvan ) THEN + CALL calbec( npw, vkb, evcr, becp ) + CALL s_psi( npwx, npw, nbnd, evcr, evcr ) + ENDIF + ! + ! Compute the trace of the representation for each group of bands + ! + DO igroup=1,ngroup + dim_rap=istart(igroup+1)-istart(igroup) + DO i=1,dim_rap + ibnd=istart(igroup)+i-1 + trace(iclass,igroup)=trace(iclass,igroup) + & + zdotc(2*npwx,evc(1,ibnd),1,evcr(1,ibnd),1) + ENDDO + ! write(6,*) igroup, iclass, dim_rap, trace(iclass,igroup) + ENDDO + ENDDO + ! + CALL mp_sum(trace,intra_bgrp_comm) + ! +! DO iclass=1,nclass +! write(6,'(i5,3(2f11.8,1x))') iclass,trace(iclass,1),trace(iclass,2), & +! trace(iclass,3) +! ENDDO + ! + ! And now use the character table to identify the symmetry representation + ! of each group of bands + ! + !IF (domag) THEN + ! WRITE(stdout,'(/,5x,"Band symmetry, ",a11," [",a11, & + ! & "] magnetic double point group,")') gname, gname_is + ! WRITE(stdout,'(5x,"using ",a11,/)') gname_is + !ELSE + ! WRITE(stdout,'(/,5x,"Band symmetry, ",a11," double point group:",/)') gname + !ENDIF + + DO igroup=1,ngroup + dim_rap=istart(igroup+1)-istart(igroup) + shift=0 + DO irap=1,nrap + times(igroup,irap)=(0.d0,0.d0) + DO iclass=1,nclass + times(igroup,irap)=times(igroup,irap) & + +trace(iclass,igroup)*CONJG(char_mat_so(irap, & + which_irr_so(iclass)))*DBLE(nelem_so(iclass)) + ENDDO + times(igroup,irap)=times(igroup,irap)/2/nsym + + IF ((abs(nint(dble(times(igroup,irap)))-dble(times(igroup,irap)))& + > accuracy).or. (abs(aimag(times(igroup,irap))) > accuracy) ) THEN + ! WRITE(stdout,'(5x,"e(",i3," -",i3,") = ",f12.5,2x,"eV",3x,i3,3x,& + ! & "--> ?")') & + ! istart(igroup), istart(igroup+1)-1, w1(istart(igroup)), dim_rap + ibnd=istart(igroup) + IF (rap_et(ibnd)==-1) THEN + DO i=1,dim_rap + ibnd=istart(igroup)+i-1 + rap_et(ibnd)=0 + ENDDO + ENDIF + GOTO 300 + ENDIF + IF (abs(times(igroup,irap)) > accuracy) THEN + dimen=nint(dble(char_mat_so(irap,1))) + ibnd=istart(igroup) + shift + IF (rap_et(ibnd)==-1) THEN + DO i=1,dimen*nint(dble(times(igroup,irap))) + ibnd=istart(igroup)+shift+i-1 + rap_et(ibnd)=irap + ENDDO + shift=shift+dimen*nint(dble(times(igroup,irap))) + ENDIF + ! IF (ABS(NINT(DBLE(times))-1.d0) < 1.d-4) THEN + ! WRITE(stdout,'(5x, "e(",i3," -",i3,") = ",f12.5,2x,"eV",3x,i3,3x,& + ! & "--> ",a15)') & + ! istart(igroup), istart(igroup+1)-1, w1(istart(igroup)), & + ! dim_rap, name_rap_so(irap) + ! ELSE + ! WRITE(stdout,'(5x,"e(",i3," -",i3,") = ",f12.5,2x,"eV",3x,i3,& + ! & 3x,"--> ",i3," ",a15)') & + ! istart(igroup), istart(igroup+1)-1, & + ! w1(istart(igroup)), dim_rap, NINT(DBLE(times)), name_rap_so(irap) + ! END IF + ENDIF + ENDDO +300 CONTINUE + ENDDO + !WRITE( stdout, '(/,1x,74("*"))') + + DEALLOCATE(trace) + DEALLOCATE(w1) + DEALLOCATE(evcr) + IF (okvan) CALL deallocate_bec_type ( becp ) + RETURN +END SUBROUTINE find_band_sym_so + +SUBROUTINE rotate_all_psi_so(ik,evc_nc,evcr,s,ftau,d_spin,has_e,gk) + ! + ! This subroutine rotates a spinor wavefunction according to the symmetry + ! s. d_spin contains the 2x2 rotation matrix in the spin space. + ! has_e=-1 means that also a 360 degrees rotation is applied in spin space. + ! + USE kinds, ONLY : DP + USE constants, ONLY : tpi + USE fft_base, ONLY : dfftp + USE scatter_mod, ONLY : cgather_sym_many, cscatter_sym_many + USE fft_interfaces, ONLY : fwfft, invfft + USE gvect, ONLY : ngm + USE wvfct, ONLY : nbnd, npwx + USE klist, ONLY : ngk, igk_k + USE noncollin_module, ONLY : npol + USE mp_bands, ONLY : intra_bgrp_comm + USE mp, ONLY : mp_sum + + IMPLICIT NONE + + INTEGER, INTENT(in) :: ik + INTEGER :: s(3,3), ftau(3), gk(3), has_e + COMPLEX(DP) :: evc_nc(npwx,2,nbnd), evcr(npwx,2,nbnd), d_spin(2,2) + + COMPLEX(DP), ALLOCATABLE :: psic(:,:), psir(:), evcr_save(:,:,:) + COMPLEX(DP) :: phase + REAL(DP) :: arg + INTEGER :: i, j, k, ri, rj, rk, ir, rir, ipol, jpol, ibnd + INTEGER :: nr1, nr2, nr3, nr1x, nr2x, nr3x, nrxx, npw + LOGICAL :: zone_border + INTEGER :: start_band, last_band, my_nbnd_proc + INTEGER :: start_band_proc(dfftp%nproc), nbnd_proc(dfftp%nproc) + ! +#if defined (__MPI) + ! + COMPLEX (DP), ALLOCATABLE :: psir_collect(:) + COMPLEX (DP), ALLOCATABLE :: psic_collect(:,:) + ! +#endif + + nr1=dfftp%nr1 + nr2=dfftp%nr2 + nr3=dfftp%nr3 + nr1x=dfftp%nr1x + nr2x=dfftp%nr2x + nr3x=dfftp%nr3x + nrxx=dfftp%nnr + +#if defined (__MPI) + call divide (intra_bgrp_comm, nbnd, start_band, last_band) + start_band_proc=0 + start_band_proc(dfftp%mype+1)=start_band + nbnd_proc=0 + my_nbnd_proc=last_band-start_band+1 + nbnd_proc(dfftp%mype+1)=my_nbnd_proc + CALL mp_sum(start_band_proc, intra_bgrp_comm) + CALL mp_sum(nbnd_proc, intra_bgrp_comm) + ALLOCATE (psic_collect(nr1x*nr2x*nr3x,my_nbnd_proc)) + ALLOCATE (psir_collect(nr1x*nr2x*nr3x)) +#endif + ! + ALLOCATE(psic(nrxx,nbnd)) + ALLOCATE(psir(nrxx)) + ALLOCATE(evcr_save(npwx,npol,nbnd)) + ! + zone_border=(gk(1)/=0.or.gk(2)/=0.or.gk(3)/=0) + ! + npw = ngk(ik) + evcr_save=(0.0_DP,0.0_DP) + DO ipol=1,npol + ! + psic = ( 0.D0, 0.D0 ) + psir = ( 0.D0, 0.D0 ) + ! + DO ibnd=1,nbnd + psic(dfftp%nl(igk_k(1:npw,ik)),ibnd) = evc_nc(1:npw,ipol,ibnd) + CALL invfft ('Rho', psic(:,ibnd), dfftp) + ENDDO + ! +#if defined (__MPI) + ! + ! + CALL cgather_sym_many( dfftp, psic, psic_collect, nbnd, nbnd_proc, & + start_band_proc ) + ! + psir_collect=(0.d0,0.d0) + DO ibnd=1,my_nbnd_proc + ! + IF (zone_border) THEN + DO k = 1, nr3 + DO j = 1, nr2 + DO i = 1, nr1 + CALL ruotaijk (s, ftau, i, j, k, nr1, nr2, nr3, ri, rj, rk ) + ir=i+(j-1)*nr1x+(k-1)*nr1x*nr2x + rir=ri+(rj-1)*nr1x+(rk-1)*nr1x*nr2x + arg=tpi*( (gk(1)*(i-1))/dble(nr1)+(gk(2)*(j-1))/dble(nr2)+ & + (gk(3)*(k-1))/dble(nr3) ) + phase=cmplx(cos(arg),sin(arg),kind=DP) + psir_collect(ir)=psic_collect(rir,ibnd)*phase + ENDDO + ENDDO + ENDDO + ELSE + DO k = 1, nr3 + DO j = 1, nr2 + DO i = 1, nr1 + CALL ruotaijk (s, ftau, i, j, k, nr1, nr2, nr3, ri, rj, rk ) + ir=i+(j-1)*nr1x+(k-1)*nr1x*nr2x + rir=ri+(rj-1)*nr1x+(rk-1)*nr1x*nr2x + psir_collect(ir)=psic_collect(rir,ibnd) + ENDDO + ENDDO + ENDDO + ENDIF + psic_collect(:,ibnd) = psir_collect(:) + ENDDO + DO ibnd=1,nbnd + ! + CALL cscatter_sym_many(dfftp, psic_collect, psir, ibnd, nbnd, nbnd_proc, & + start_band_proc) + CALL fwfft ('Rho', psir, dfftp) + ! + evcr_save(1:npw,ipol,ibnd) = psir(dfftp%nl(igk_k(1:npw,ik))) + ENDDO + ! +#else + DO ibnd=1,nbnd + IF (zone_border) THEN + DO k = 1, nr3 + DO j = 1, nr2 + DO i = 1, nr1 + CALL ruotaijk (s, ftau, i, j, k, nr1, nr2, nr3, ri, rj, rk ) + ir=i+(j-1)*nr1x+(k-1)*nr1x*nr2x + rir=ri+(rj-1)*nr1x+(rk-1)*nr1x*nr2x + arg=tpi*( (gk(1)*(i-1))/dble(nr1)+(gk(2)*(j-1))/dble(nr2)+ & + (gk(3)*(k-1))/dble(nr3) ) + phase=cmplx(cos(arg),sin(arg),kind=DP) + psir(ir)=psic(rir,ibnd)*phase + ENDDO + ENDDO + ENDDO + ELSE + DO k = 1, nr3 + DO j = 1, nr2 + DO i = 1, nr1 + CALL ruotaijk (s, ftau, i, j, k, nr1, nr2, nr3, ri, rj, rk ) + ir=i+(j-1)*nr1x+(k-1)*nr1x*nr2x + rir=ri+(rj-1)*nr1x+(rk-1)*nr1x*nr2x + psir(ir)=psic(rir,ibnd) + ENDDO + ENDDO + ENDDO + ENDIF + CALL fwfft ('Rho', psir(:), dfftp) + ! + evcr_save(1:npw,ipol,ibnd) = psir(dfftp%nl(igk_k(1:npw,ik))) + ENDDO + ! +#endif + ! + ! + ENDDO + + evcr=(0.d0,0.d0) + DO ibnd=1,nbnd + DO ipol=1,npol + DO jpol=1,npol + evcr(:,ipol,ibnd)=evcr(:,ipol,ibnd)+ & + d_spin(ipol,jpol)*evcr_save(:,jpol,ibnd) + ENDDO + ENDDO + ENDDO + IF (has_e==-1) evcr=-evcr + ! + DEALLOCATE(evcr_save) + DEALLOCATE(psic) + DEALLOCATE(psir) +#if defined (__MPI) + DEALLOCATE (psic_collect) + DEALLOCATE (psir_collect) +#endif + RETURN +END SUBROUTINE rotate_all_psi_so + +SUBROUTINE find_nks1nks2(firstk,lastk,nks1tot,nks1,nks2tot,nks2,spin_component) + ! + ! This routine selects the first (nks1) and last (nks2) k point calculated + ! by the current pool. + ! + USE lsda_mod, ONLY : nspin + USE klist, ONLY : nks, nkstot + USE mp_pools, ONLY : my_pool_id, npool, kunit + + IMPLICIT NONE + INTEGER, INTENT(out) :: nks1tot,nks1,nks2tot,nks2 + INTEGER, INTENT(in) :: firstk, lastk, spin_component + INTEGER :: nbase, rest + + IF (nspin==1.or.nspin==4) THEN + nks1tot=max(1,firstk) + nks2tot=min(nkstot, lastk) + ELSEIF (nspin==2) THEN + IF (spin_component == 1) THEN + nks1tot=max(1,firstk) + nks2tot=min(nkstot/2,lastk) + ELSEIF (spin_component==2) THEN + nks1tot=nkstot/2 + max(1,firstk) + nks2tot=nkstot/2 + min(nkstot/2,lastk) + ENDIF + ENDIF + IF (nks1tot>nks2tot) CALL errore('findnks1nks2','wrong nks1tot or nks2tot',1) + +#ifdef __MPI + nks = kunit * ( nkstot / kunit / npool ) + rest = ( nkstot - nks * npool ) / kunit + IF ( ( my_pool_id + 1 ) <= rest ) nks = nks + kunit + ! + ! ... calculates nbase = the position in the list of the first point that + ! ... belong to this npool - 1 + ! + nbase = nks * my_pool_id + IF ( ( my_pool_id + 1 ) > rest ) nbase = nbase + rest * kunit + + nks1=max(1,nks1tot-nbase) + IF (nks1>nks) nks1=nks+1 + nks2=min(nks,nks2tot-nbase) + IF (nks2<1) nks2=nks1-1 +#else + nks1=nks1tot + nks2=nks2tot +#endif + +END SUBROUTINE find_nks1nks2 + +SUBROUTINE find_info_group(nsym,s,t_rev,ft,d_spink,gk,sname, & + s_is,d_spin_is,gk_is, invs_is,is_symmorphic,search_sym) + ! + ! This routine receives as input a point group and sets the corresponding + ! variables for the description of the classes and of the irreducible + ! representations. It sets also the group name and code. + ! In the magnetic case it selects the invariat subgroup. + ! + USE kinds, ONLY : DP + USE cell_base, ONLY : at, bg + USE fft_base, ONLY : dfftp + USE noncollin_module, ONLY : noncolin + USE spin_orb, ONLY : domag + USE rap_point_group, ONLY : code_group, nclass, nelem, elem, which_irr, & + char_mat, name_rap, name_class, gname, ir_ram + USE rap_point_group_so, ONLY : nrap, nelem_so, elem_so, has_e, & + which_irr_so, char_mat_so, name_rap_so, & + name_class_so, d_spin, name_class_so1 + USE rap_point_group_is, ONLY : nsym_is, sr_is, ft_is, gname_is, & + sname_is, code_group_is + + IMPLICIT NONE + + INTEGER, INTENT(in) :: nsym, & ! dimension of the group + s(3,3,48), & ! rotation matrices + t_rev(48), & ! if time reversal is need + gk(3,48) + REAL(dp), INTENT(IN) :: ft(3,48)! fractionary translation + + INTEGER, INTENT(out) :: s_is(3,3,48), & ! rotation matrices + gk_is(3,48), invs_is(48) + + COMPLEX(DP),INTENT(out) :: d_spink(2,2,48), & ! rotation in spin space + d_spin_is(2,2,48) ! rotation in spin space + + LOGICAL, INTENT(out) :: is_symmorphic, & ! true if the gruop is symmorphic + search_sym ! true if gk + + CHARACTER(len=45), INTENT(in) :: sname(48) + + REAL(DP) :: sr(3,3,48) + INTEGER :: isym, jsym, ss(3,3) + LOGICAL :: found + + + is_symmorphic=.true. + search_sym=.true. + + DO isym=1,nsym + is_symmorphic=( is_symmorphic.and.(ft(1,isym)==0.d0).and. & + (ft(2,isym)==0.d0).and. & + (ft(3,isym)==0.d0) ) + ENDDO + + IF (.NOT.is_symmorphic) THEN + + DO isym=1,nsym + DO jsym=1,nsym + search_sym=search_sym.AND.(ABS(gk(1,isym)*ft(1,jsym)+ & + gk(2,isym)*ft(2,jsym)+gk(3,isym)*ft(3,jsym))<1.D-8) + ENDDO + ENDDO + ENDIF + ! + ! Set the group name, divide it in classes and set the irreducible + ! representations + ! + nsym_is=0 + DO isym=1,nsym + CALL s_axis_to_cart (s(1,1,isym), sr(1,1,isym), at, bg) + IF (noncolin) THEN + ! + ! In the noncollinear magnetic case finds the invariant subgroup of the point + ! group of k. Presently we use only this subgroup to classify the levels. + ! + IF (domag) THEN + IF (t_rev(isym)==0) THEN + nsym_is=nsym_is+1 + CALL s_axis_to_cart (s(1,1,isym), sr_is(1,1,nsym_is), at, bg) + CALL find_u(sr_is(1,1,nsym_is),d_spin_is(1,1,nsym_is)) + s_is(:,:,nsym_is)=s(:,:,isym) + gk_is(:,nsym_is)=gk(:,isym) + ft_is(:,nsym_is)=ft(:,isym) + sname_is(nsym_is)=sname(isym) + ENDIF + ELSE + CALL find_u(sr(1,1,isym),d_spink(1,1,isym)) + ENDIF + ENDIF + ENDDO + + IF (noncolin.AND.domag) THEN +! +! find the inverse of each element +! + DO isym = 1, nsym_is + found = .FALSE. + DO jsym = 1, nsym_is + ! + ss = MATMUL (s_is(:,:,jsym),s_is(:,:,isym)) + ! s(:,:,1) is the identity + IF ( ALL ( s_is(:,:,1) == ss(:,:) ) ) THEN + invs_is (isym) = jsym + found = .TRUE. + ENDIF + END DO + IF ( .NOT.found) CALL errore ('find_info_group', ' Not a group', 1) + ENDDO +! +! Recheck if we can compute the representations. The group is now smaller +! + is_symmorphic=.TRUE. + search_sym=.TRUE. + + DO isym=1,nsym_is + is_symmorphic=( is_symmorphic.AND.(ft_is(1,isym)==0.d0).AND. & + (ft_is(2,isym)==0.d0).AND. & + (ft_is(3,isym)==0.d0) ) + ENDDO + IF (.NOT.is_symmorphic) THEN + DO isym=1,nsym_is + DO jsym=1,nsym_is + search_sym=search_sym.AND.(ABS(gk_is(1,isym)*ft_is(1,jsym)+ & + gk_is(2,isym)*ft_is(2,jsym)+ & + gk_is(3,isym)*ft_is(3,jsym))<1.D-8) + ENDDO + ENDDO + ENDIF + END IF + ! + ! Set the group name, divide it in classes and set the irreducible + ! + CALL find_group(nsym,sr,gname,code_group) + IF (noncolin) THEN + IF (domag) THEN + CALL find_group(nsym_is,sr_is,gname_is,code_group_is) + CALL set_irr_rap_so(code_group_is,nclass,nrap,char_mat_so, & + name_rap_so,name_class_so,name_class_so1) + CALL divide_class_so(code_group_is,nsym_is,sr_is,d_spin_is,& + has_e,nclass,nelem_so,elem_so,which_irr_so) + ELSE + CALL set_irr_rap_so(code_group,nclass,nrap,char_mat_so, & + name_rap_so,name_class_so,name_class_so1) + CALL divide_class_so(code_group,nsym,sr,d_spink, & + has_e,nclass,nelem_so,elem_so,which_irr_so) + ENDIF + ELSE + CALL set_irr_rap(code_group,nclass,char_mat,name_rap,name_class,ir_ram) + CALL divide_class(code_group,nsym,sr,nclass,nelem,elem,which_irr) + ENDIF + + RETURN +END SUBROUTINE find_info_group +! +! Copyright (C) 2001 PWSCF group +! This file is distributed under the terms of the +! GNU General Public License. See the file `License' +! in the root directory of the present distribution, +! or http://www.gnu.org/copyleft/gpl.txt . +! +! +!---------------------------------------------------------------------- +SUBROUTINE s_axis_to_cart (s, sr, at, bg) + !---------------------------------------------------------------------- + ! + ! This routine transform a symmetry matrix expressed in the + ! basis of the crystal axis in the cartesian basis. + ! + ! last revised 3 may 1995 by A. Dal Corso + ! + ! + USE kinds + IMPLICIT NONE + ! + ! first the input parameters + ! + INTEGER :: s (3, 3) + ! input: matrix in crystal axis + real(DP) :: sr (3, 3), at (3, 3), bg (3, 3) + ! output: matrix in cartesian axis + ! input: direct lattice vectors + ! input: reciprocal lattice vectors + ! + ! here the local variable + ! + + INTEGER :: apol, bpol, kpol, lpol + ! + ! counters on polarizations + ! + DO apol = 1, 3 + DO bpol = 1, 3 + sr (apol, bpol) = 0.d0 + DO kpol = 1, 3 + DO lpol = 1, 3 + sr (apol, bpol) = sr (apol, bpol) + at (apol, kpol) * & + dble ( s (lpol, kpol) ) * bg (bpol, lpol) + ENDDO + ENDDO + ENDDO + ENDDO + RETURN +END SUBROUTINE s_axis_to_cart +