Merge branch 'develop' into 'us_acc'

# Conflicts:
#   GWW/head/lanczos_k.f90
#   GWW/head/solve_head.f90
#   LR_Modules/lr_sm1_psi.f90
#   PP/src/projwfc.f90
This commit is contained in:
giannozz 2021-09-28 19:00:39 +00:00
commit 42674ef499
52 changed files with 972 additions and 243 deletions

View File

@ -72,7 +72,7 @@ build:cmake-nvhpc:
- cmake --version
- mkdir build
- cd build
- cmake -DBUILD_SHARED_LIBS=ON -DCMAKE_Fortran_COMPILER=mpif90 -DCMAKE_C_COMPILER=mpicc
- cmake -DBUILD_SHARED_LIBS=OFF -DCMAKE_Fortran_COMPILER=mpif90 -DCMAKE_C_COMPILER=mpicc
-DQE_ENABLE_CUDA=ON -DQE_ENABLE_OPENACC=ON .. && make
#build:centos:

View File

@ -1,6 +1,7 @@
New in development version:
* RMM-DIIS for CPU (S. Nisihara) and GPU (E. de Paoli, P. Delugas)
* DFT-D3: MPI parallelization and GPU acceleration with OPenACC
* projwfc.x can be used to compute the PDOS in a local basis (I. Timrov)
Fixed in development version:
* Some build problems occurring under special circumstances
@ -11,6 +12,9 @@ Fixed in development version:
* In DFT+U (lda_plus_u_kind = 0 and 1) the pw.x code was printing the squared
eigenvectors instead of simply eigenvectors. Now it prints the
eigenvectors (consistent with lda_plus_u_kind = 2).
* plotband.x wasn't correctly plotting the bands, under some not-so-special
circumstances
* CP with DFT+U could crash when writing the xml file (since v.6.6)
Incompatible changes in develoment version:
* Changes to Makefiles and to file "make.inc"

View File

@ -606,7 +606,6 @@
!
USE kinds, ONLY : DP
USE ions_base, ONLY : nat, ityp, ntyp => nsp
USE spin_orb, ONLY : lspinorb
USE cell_base, ONLY : tpiba2, omega, tpiba
USE gvect, ONLY : ngm, gg, g, eigts1, eigts2, eigts3, mill
USE scf, ONLY : v, vltot

View File

@ -376,8 +376,8 @@ subroutine do_polarization_lanczos(tf,options,ispin)
call read_compact_q_lanczos(cql, 1)!just for obtaining numpw,numt...poor man solution
allocate(cql_save(1,1,1))
endif
allocate(vtl_save(1,1,1))
allocate(ttl_save(1,1,1))
else
!put all matrices vtl and ttl in memory, distributed according to valence state
@ -581,9 +581,9 @@ subroutine do_polarization_lanczos(tf,options,ispin)
if(options%l_t_wannier) then
deallocate(cql_save)
call free_memory_compact_q_lanczos(cql)
else
deallocate(vtl_save,ttl_save)
endif
deallocate(vtl_save,ttl_save)
! deallocate(e_mat)
deallocate(af,occ)

View File

@ -36,7 +36,8 @@ subroutine bcast_ph_input ( )
USE io_global, ONLY : ionode_id
USE run_info, ONLY : title
USE wannier_gw, ONLY : l_head, omega_gauss, n_gauss, grid_type, nsteps_lanczos,&
&second_grid_n,second_grid_i,l_scissor,scissor,len_head_block_freq,len_head_block_wfc
&second_grid_n,second_grid_i,l_scissor,scissor,len_head_block_freq,&
l_easy,len_head_block_wfc
USE control_lr, ONLY : lgamma, lrpa
@ -106,6 +107,7 @@ subroutine bcast_ph_input ( )
call mp_bcast(scissor, ionode_id, world_comm)
call mp_bcast(len_head_block_freq, ionode_id, world_comm)
call mp_bcast(len_head_block_wfc, ionode_id,world_comm)
call mp_bcast(l_easy , ionode_id,world_comm)
#endif
return

View File

@ -79,7 +79,7 @@ SUBROUTINE phq_readin()
USE ahc, ONLY : elph_ahc, ahc_dir, ahc_nbnd, ahc_nbndskip, &
skip_upperfan
USE wannier_gw, ONLY : l_head, omega_gauss, n_gauss, grid_type, nsteps_lanczos,second_grid_n,second_grid_i,&
&l_scissor,scissor, len_head_block_freq, len_head_block_wfc
&l_scissor,scissor, len_head_block_freq, len_head_block_wfc, l_easy
!
IMPLICIT NONE
!
@ -128,7 +128,7 @@ SUBROUTINE phq_readin()
wpot_dir, ahc_dir, ahc_nbnd, ahc_nbndskip, &
skip_upperfan, &
l_head, omega_gauss, n_gauss, grid_type,nsteps_lanczos,l_scissor,scissor,&
second_grid_n,second_grid_i,len_head_block_wfc,len_head_block_freq
second_grid_n,second_grid_i,len_head_block_wfc,len_head_block_freq, l_easy
! tr2_ph : convergence threshold
! amass : atomic masses
@ -348,6 +348,8 @@ SUBROUTINE phq_readin()
d2ns_type = 'full'
len_head_block_freq=0
len_head_block_wfc=0
l_easy=.false.
!
! ... reading the namelist inputph
!

View File

@ -149,13 +149,19 @@ PROGRAM head
!
! ... cleanup of the variables for the next q point
!
write(stdout,*) 'DEBUG 1'
CALL clean_pw_ph(iq)
write(stdout,*) 'DEBUG 2'
!
END DO
CALL ph_writefile('init',0,0,ierr)
write(stdout,*) 'DEBUG 3'
CALL ph_writefile('init',0,0,ierr)
write(stdout,*) 'DEBUG 4'
CALL collect_grid_files()
write(stdout,*) 'DEBUG 5'
CALL destroy_status_run()
write(stdout,*) 'DEBUG 6'
!
IF (bands_computed) CALL print_clock_pw()
!

View File

@ -29,6 +29,7 @@ subroutine lanczos_state_k(ik,nstates, nsteps,in_states,d,f,omat,dpsi_ipol, t_ou
USE klist, ONLY : xk,igk_k, ngk
USE noncollin_module, ONLY : noncolin, npol
USE uspp_init, ONLY : init_us_2
USE wvfct, ONLY : current_k
implicit none
@ -60,7 +61,9 @@ subroutine lanczos_state_k(ik,nstates, nsteps,in_states,d,f,omat,dpsi_ipol, t_ou
allocate(alpha(nstates),beta(nstates),gamma(nstates),n_1(nstates),delta(nstates))
allocate(c(nstates))
allocate(spsi(npwx,nstates))
current_k=ik
npw = ngk(ik)
t_out(:,:,:)=(0.d0,0.d0)

View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001-2013 Quantum ESPRESSO group
! Copyright (C) 2001-2021 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,
@ -27,7 +27,7 @@ subroutine solve_head
USE wannier_gw, ONLY : n_gauss, omega_gauss, grid_type,&
nsteps_lanczos,second_grid_n,second_grid_i,&
l_scissor,scissor,len_head_block_freq, &
len_head_block_wfc
len_head_block_wfc, l_easy
USE control_ph, ONLY : tr2_ph
USE gvect, ONLY : ngm, ngm_g, ig_l2g, gstart, g
USE gvecs, ONLY : doublegrid
@ -46,6 +46,7 @@ subroutine solve_head
use qpoint, ONLY : npwq, nksq
use control_lr, ONLY : nbnd_occ, lgamma
use uspp_init, ONLY : init_us_2
use SCF
implicit none
@ -93,6 +94,7 @@ subroutine solve_head
COMPLEX(kind=DP), ALLOCATABLE :: z_dl(:),z_d(:),z_du(:),z_b(:)
COMPLEX(kind=DP) :: csca, csca1
COMPLEX(kind=DP), ALLOCATABLE :: t_out(:,:,:), psi_tmp(:)
TYPE(scf_type) :: wing
INTEGER :: n
INTEGER :: npwx_g
@ -145,6 +147,9 @@ subroutine solve_head
do i=1,n_gauss
freqs(1+i)=omega_gauss*dble(i)/dble(n_gauss)
enddo
if(l_easy) then
freqs(1)=omega_gauss/dble(n_gauss)/4.d0
endif
else if(grid_type==4) then!equally spaced grid shifted of 1/2
freqs(1) = 0.d0
do i=1,n_gauss
@ -472,18 +477,19 @@ subroutine solve_head
#else
call syme (pola_charge(:,:,:,i))
#endif
call create_scf_type ( wing, .true. )
do ipol=1,3
CALL fwfft ('Rho', pola_charge(1:dfftp%nnr,1,ipol,i), dfftp)
tmp_g(:)=(0.d0,0.d0)
tmp_g(gstart:ngm)=pola_charge(dfftp%nl(gstart:ngm),1,ipol,i)
wing%of_g(1:ngm,1)=-4.d0*tmp_g(1:ngm)
call write_wing ( wing, nspin,ipol,i)
!loop on frequency
do ig=gstart,ngm
e_head_pol(ig,i,ipol)=-4.d0*tmp_g(ig)
enddo
enddo
call destroy_scf_type (wing )
enddo
@ -569,9 +575,74 @@ subroutine solve_head
call mp_barrier( world_comm )
write(stdout,*) 'ATT2'
write(stdout,*) 'THIS IS THE END'
call stop_clock ('solve_head')
return
end subroutine solve_head
SUBROUTINE write_wing ( rho, nspin,ipol,iw)
USE kinds, ONLY : DP
USE io_files, ONLY : create_directory
USE io_base, ONLY : write_rhog, read_rhog
!
USE paw_variables, ONLY : okpaw
USE ldaU, ONLY : lda_plus_u
USE noncollin_module, ONLY : noncolin
USE spin_orb, ONLY : domag
USE scf, ONLY : scf_type
!
USE cell_base, ONLY : bg, tpiba
USE gvect, ONLY : ig_l2g, mill
USE control_flags, ONLY : gamma_only
USE io_files, ONLY : seqopn, tmp_dir, prefix, postfix
USE io_global, ONLY : ionode, ionode_id, stdout
USE mp_pools, ONLY : my_pool_id
USE mp_bands, ONLY : my_bgrp_id, root_bgrp_id, &
root_bgrp, intra_bgrp_comm
USE mp_images, ONLY : intra_image_comm
USE mp, ONLY : mp_bcast
!
IMPLICIT NONE
TYPE(scf_type), INTENT(IN) :: rho
INTEGER, INTENT(IN) :: nspin
INTEGER, INTENT(IN) :: ipol!direction
INTEGER, INTENT(IN) :: iw !frequency
!
CHARACTER (LEN=256) :: dirname
LOGICAL :: lexist
INTEGER :: nspin_, iunocc, iunpaw, ierr
INTEGER, EXTERNAL :: find_free_unit
CHARACTER(5) :: nfile
CHARACTER :: npol
write(nfile,'(5i1)') &
& iw/10000,mod(iw,10000)/1000,mod(iw,1000)/100,mod(iw,100)/10,mod(iw,10)
write(npol,'(1i1)') ipol
dirname = TRIM(tmp_dir) // TRIM(prefix) // postfix
CALL create_directory( dirname )
! in the following case do not read or write polarization
IF ( noncolin .AND. .NOT.domag ) THEN
nspin_ = 1
ELSE
nspin_ = nspin
ENDIF
! Write G-space density
IF ( my_pool_id == 0 .AND. my_bgrp_id == root_bgrp_id ) &
CALL write_rhog( TRIM(dirname) // "wing_" // npol // "_" //nfile, &
root_bgrp, intra_bgrp_comm, &
bg(:,1)*tpiba, bg(:,2)*tpiba, bg(:,3)*tpiba, &
gamma_only, mill, ig_l2g, rho%of_g(:,1:nspin_) )
RETURN
END SUBROUTINE write_wing

View File

@ -881,7 +881,7 @@ MODULE convergence_gw
write(stdout,*) iw,dble(cmat(iw,ct%nf))
enddo
if(ionode.and. kk==ct%iband .and. l_verbose) close(iun)
sigmac=0.d0
do dw=-ct%nf+1,ct%nf
dww=dw+ct%nf
@ -1030,7 +1030,7 @@ MODULE convergence_gw
deallocate(vprods,vprods_im,wterms,wterms_im)
deallocate(gwmat_rr,gwmat_ri,gwmat_ir,gwmat_ii)
if(ionode) close(iun)
enddo
enddo
nr_counter=nr_counter+1
@ -1574,7 +1574,7 @@ MODULE convergence_gw
REAL(kind=DP) :: energy0, energy1, energy,energy_old,mod_force, mod_force_old
COMPLEX(kind=DP), ALLOCATABLE :: force(:),phi_old(:)
COMPLEX(kind=DP), ALLOCATABLE :: psi_g_phi(:,:)
REAL(kind=DP) :: alpha_old,par_a,par_b,par_c, alpha_new
REAL(kind=DP) :: alpha_old,par_a,par_b=1.d0,par_c=1.d0, alpha_new
LOGICAL :: l_updated
COMPLEX(kind=DP), ALLOCATABLE :: h_diag2(:,:), s_diag2(:,:)
INTEGER :: is

View File

@ -74,6 +74,19 @@
TYPE(exchange_cus) :: exx_cus
INTERFACE
SUBROUTINE energies_xc( lda, n, m, psi, e_xc, e_h,ispin, v_states )
USE kinds, ONLY : DP
USE fft_base, ONLY : dffts
USE lsda_mod, ONLY : nspin
INTEGER :: lda, n, m
COMPLEX(kind=DP) :: psi(lda,m)
REAL(kind=DP) :: e_xc(m), e_h(m)
INTEGER, INTENT(in) :: ispin !spin 1,2
REAL(kind=DP), OPTIONAL :: v_states(dffts%nnr,m, nspin)
END SUBROUTINE energies_xc
END INTERFACE
! interface
! subroutine fake_conduction_wannier(fcw_n,fcw_s,fcw_m,cut,s_cut)

View File

@ -14,7 +14,6 @@ subroutine product_basis
USE uspp_param, ONLY : upf, nh
USE noncollin_module, ONLY: npol, noncolin
USE mp_world, ONLY : world_comm
USE spin_orb, ONLY: lspinorb
USE ions_base, ONLY : nat, nsp, ityp
USE io_global, ONLY : stdout, ionode
USE input_simple
@ -290,7 +289,6 @@ SUBROUTINE optimal_gram_schmidt_nc(num_in,wfcs,thres,num_out)
USE noncollin_module, ONLY: npol, noncolin
USE input_simple, ONLY : npw_max,vkb_max
USE noncollin_module, ONLY: npol, noncolin
USE spin_orb, ONLY: lspinorb
USE ions_base, ONLY : nat, nsp, ityp

View File

@ -9,7 +9,6 @@ subroutine v_product
USE klist, ONLY : nks,ngk,xk
USE noncollin_module, ONLY: npol, noncolin
USE mp_world, ONLY : world_comm,mpime
USE spin_orb, ONLY: lspinorb
USE io_global, ONLY : stdout, ionode
USE input_simple
USE gvect, ONLY : ngm, gstart,gg, g

View File

@ -12,7 +12,6 @@ subroutine wfc_basis
USE uspp_param, ONLY : upf, nh
USE noncollin_module, ONLY: npol, noncolin
USE mp_world, ONLY : world_comm
USE spin_orb, ONLY: lspinorb
USE ions_base, ONLY : nat, nsp, ityp
USE io_global, ONLY : stdout, ionode
USE input_simple
@ -541,7 +540,6 @@ SUBROUTINE optimal_gram_schmidt_z(num_in,wfcs,ithres,thres,num_out)
USE uspp, ONLY : nkb, vkb, becsum, nhtol, nhtoj, indv, okvan
USE uspp_param, ONLY : upf, nh
USE noncollin_module, ONLY: npol, noncolin
USE spin_orb, ONLY: lspinorb
USE ions_base, ONLY : nat, nsp, ityp

View File

@ -46,22 +46,19 @@ subroutine mix_potential (ndim, vout, vin, alphamix, dr2, tr2, &
!
! Here the local variables
!
! max number of iterations used in mixing: n_iter < maxter must hold.
integer, parameter :: maxter = 8
!
integer :: iunmix, n, i, j, iwork (maxter), info, iter_used, &
ipos, inext, ndimtot
! work space containing info from previous iterations:
! must be kept in memory and saved between calls if file_extension=' '
real(DP), allocatable, save :: df (:,:), dv (:,:)
!
real(DP), allocatable :: vinsave (:)
real(DP) :: beta (maxter, maxter), gamma, work (maxter), norm
logical :: saveonfile, exst
integer :: iunmix, n, i, j, info, iter_used, ipos, inext, ndimtot
real(DP) :: gamma, norm
real(DP), allocatable, save :: df (:,:), dv (:,:)
!! work space containing info from previous iterations:
!! must be kept in memory and saved between calls if file_extension=' '
!
real(DP) :: w(maxter) = 1.d0
real(DP) :: w0 = 0.01d0
!! adjustable parameters as suggested in the original paper
integer, allocatable :: iwork(:)
real(DP), allocatable :: vinsave(:), beta(:, :), work(:), w(:)
!
real(DP) :: w0 = 0.01_dp
real(DP) :: w1 = 1.0_dp
!! adjustable parameters as suggested in the original paper. w(:) is set to w1.
!
INTEGER, EXTERNAL :: find_free_unit
REAL(DP), EXTERNAL :: ddot, dnrm2
@ -69,7 +66,6 @@ subroutine mix_potential (ndim, vout, vin, alphamix, dr2, tr2, &
call start_clock ('mix_pot')
!
if (iter < 1) call errore ('mix_potential', 'iter must be positive', 1)
if (n_iter > maxter) call errore ('mix_potential', 'n_iter too big', 1)
if (ndim < 1) call errore ('mix_potential', 'ndim must be positive', 3)
!
saveonfile = file_extension /= ' '
@ -161,29 +157,41 @@ subroutine mix_potential (ndim, vout, vin, alphamix, dr2, tr2, &
call DCOPY (ndim, vin, 1, vinsave, 1)
endif
!
do i = 1, iter_used
do j = i + 1, iter_used
beta (i, j) = w (i) * w (j) * ddot (ndim, df (1, j), 1, df (1, i), 1)
call mp_sum ( beta (i, j), intra_bgrp_comm )
if ( iter_used > 0 ) then
!
ALLOCATE(beta(iter_used, iter_used))
ALLOCATE(w(iter_used))
ALLOCATE(work(iter_used))
ALLOCATE(iwork(iter_used))
w(:) = w1
!
beta = 0.0_dp
do i = 1, iter_used
do j = i + 1, iter_used
beta (i, j) = w (i) * w (j) * ddot (ndim, df (1, j), 1, df (1, i), 1)
call mp_sum ( beta (i, j), intra_bgrp_comm )
enddo
beta (i, i) = w0**2 + w (i) **2
enddo
beta (i, i) = w0**2 + w (i) **2
enddo
!
call DSYTRF ('U', iter_used, beta, maxter, iwork, work, maxter, info)
call errore ('broyden', 'factorization', info)
call DSYTRI ('U', iter_used, beta, maxter, iwork, work, info)
call errore ('broyden', 'DSYTRI', info)
!
do i = 1, iter_used
do j = i + 1, iter_used
beta (j, i) = beta (i, j)
!
call DSYTRF ('U', iter_used, beta, iter_used, iwork, work, iter_used, info)
call errore ('broyden', 'factorization', info)
call DSYTRI ('U', iter_used, beta, iter_used, iwork, work, info)
call errore ('broyden', 'DSYTRI', info)
deallocate ( iwork )
!
do i = 1, iter_used
do j = i + 1, iter_used
beta (j, i) = beta (i, j)
enddo
enddo
enddo
!
do i = 1, iter_used
work (i) = ddot (ndim, df (1, i), 1, vout, 1)
enddo
call mp_sum ( work(1:iter_used), intra_bgrp_comm )
!
do i = 1, iter_used
work (i) = ddot (ndim, df (1, i), 1, vout, 1)
enddo
call mp_sum ( work, intra_bgrp_comm )
!
end if
!
do n = 1, ndim
vin (n) = vin (n) + alphamix * vout (n)
@ -200,6 +208,12 @@ subroutine mix_potential (ndim, vout, vin, alphamix, dr2, tr2, &
enddo
enddo
!
IF ( iter_used > 0 ) THEN
deallocate(beta)
deallocate(w)
deallocate(work)
ENDIF
!
if (saveonfile) then
close (iunmix, status='keep')
deallocate(dv)

View File

@ -46,7 +46,7 @@ USE kinds, ONLY : DP
USE ions_base, ONLY : nat, ityp
USE uspp_param, ONLY : nh, nhm
USE noncollin_module, ONLY : npol, nspin_mag
USE spin_orb, ONLY : fcoef, domag
USE spin_orb, ONLY : fcoef
USE lrus, ONLY : intq_nc
!
IMPLICIT NONE

View File

@ -453,11 +453,13 @@ CONTAINS
IF (PRESENT(starting_ns)) CALL init_starting_ns(starting_ns_ , label)
IF (PRESENT(Hub_ns)) CALL init_Hubbard_ns(Hubbard_ns_ , label)
IF (PRESENT(Hub_ns_nc)) CALL init_Hubbard_ns(Hubbard_ns_nc_ , label)
IF (ANY(is_hubbard_back) .AND. PRESENT(hubb_l_back)) &
IF (PRESENT(is_hubbard_back)) THEN
IF (ANY(is_hubbard_back) .AND. PRESENT(hubb_l_back)) &
CALL init_Hubbard_back(is_hubbard_back, Hub_back_, hubb_l_back, backall, hubb_l1_back)
IF (ANY(is_hubbard_back) .AND. .NOT. PRESENT (hubb_l_back)) &
IF (ANY(is_hubbard_back) .AND. .NOT. PRESENT (hubb_l_back)) &
CALL errore('qexsd_init_dft:',&
'internal error background is set to true but hubb_l_back is not present',1)
END IF
!
CALL qes_init (obj, "dftU", lda_plus_u_kind, U_, J0_, alpha_, beta_, J_, starting_ns_, Hubbard_ns_, &
U_projection_type, Hub_back_, U_back_, alpha_back_, Hubbard_ns_nc_)

View File

@ -43,8 +43,6 @@ function wgauss (x, n)
! counter on 2n
real(DP), parameter :: maxarg = 200.d0
! maximum value for the argument of the exponential
real(DP), parameter :: c = 0.7071067811865475_DP
! c = sqrt(1/2)
! Fermi-Dirac smearing
if (n.eq. - 99) then
@ -67,9 +65,13 @@ function wgauss (x, n)
return
endif
! Methfessel-Paxton
!gauss_freq(x) = 0.5_DP * ERFC( - x * c)
wgauss = 0.5_DP * ERFC( - x * sqrt (2.0d0) * c)
! Methfessel-Paxton and plain gaussian cases
arg = -x
IF (arg .LT. sqrt(maxarg)) THEN
wgauss = 0.5_DP * ERFC( arg)
ELSE
wgauss = 0._DP
END IF
if (n.eq.0) return
hd = 0.d0
arg = min (maxarg, x**2)

View File

@ -87,7 +87,10 @@ input_description -distribution {Quantum ESPRESSO} -package PHonon -program ph.x
var nmix_ph -type INTEGER {
default { 4 }
info { Number of iterations used in potential mixing. }
info {
Number of iterations used in potential mixing. Using a larger value (8~20)
can significantly speed up convergence, at the cost of using more memory.
}
}
var verbosity -type CHARACTER {

View File

@ -30,7 +30,6 @@ SUBROUTINE dynmat_us()
USE cell_base, ONLY : omega, tpiba2
USE uspp_param, ONLY : nh, nhm
USE noncollin_module, ONLY : noncolin, npol
USE spin_orb, ONLY : lspinorb
USE becmod, ONLY : calbec, bec_type, allocate_bec_type, &
deallocate_bec_type, beccopy
USE modes, ONLY : u

View File

@ -65,7 +65,7 @@ MODULE io_dyn_mat
!
! ... open XML descriptor
!
iunout = xml_openfile (TRIM( fildyn ) // '.xml' )
iunout = xml_open_file (TRIM( fildyn ) // '.xml' )
!
ENDIF
CALL mp_bcast( iunout, ionode_id, intra_image_comm )
@ -260,7 +260,7 @@ MODULE io_dyn_mat
!
! Open XML descriptor
!
IF (ionode) iunout = xml_openfile( TRIM(fildyn) // '.xml')
IF (ionode) iunout = xml_open_file( TRIM(fildyn) // '.xml')
!
CALL mp_bcast(iunout, ionode_id, intra_image_comm)
IF ( iunout == -1 ) &

View File

@ -833,7 +833,7 @@ MODULE ph_restart
!! Error code
INTEGER :: iun
!
iun = xml_openfile (filename)
iun = xml_open_file (filename)
IF ( iun == -1 ) then
ierr = 1
return
@ -1072,7 +1072,7 @@ MODULE ph_restart
filename1=TRIM(filename) // TRIM(int_to_char(irr)) // '.xml'
INQUIRE(FILE=TRIM(filename1), EXIST=exst)
IF (.NOT.exst) CYCLE
iunout = xml_openfile( filename1 )
iunout = xml_open_file( filename1 )
IF (iunout == -1 ) THEN
ierr = 1
GOTO 100
@ -1097,7 +1097,7 @@ MODULE ph_restart
filename1=TRIM(filename) // TRIM(int_to_char(irr)) // '.xml'
INQUIRE(FILE=TRIM(filename1), EXIST=exst)
IF (.NOT.exst) CYCLE
iunout = xml_openfile( filename1 )
iunout = xml_open_file( filename1 )
IF (iunout == -1 ) THEN
ierr = 1
GOTO 100
@ -1359,7 +1359,7 @@ MODULE ph_restart
IF (.NOT.exst) GOTO 100
ENDIF
iunpun = xml_openfile( filename )
iunpun = xml_open_file( filename )
!
exst = (iunpun /= -1)
IF (.NOT.exst) GOTO 100

View File

@ -404,8 +404,6 @@ SUBROUTINE phq_readin()
ENDDO
IF (niter_ph.LT.1.OR.niter_ph.GT.maxter) CALL errore ('phq_readin', &
' Wrong niter_ph ', 1)
IF (nmix_ph.LT.1.OR.nmix_ph.GT.5) CALL errore ('phq_readin', ' Wrong &
&nmix_ph ', 1)
IF (iverbosity.NE.0.AND.iverbosity.NE.1) CALL errore ('phq_readin', &
&' Wrong iverbosity ', 1)
IF (fildyn.EQ.' ') CALL errore ('phq_readin', ' Wrong fildyn ', 1)

View File

@ -73,6 +73,20 @@ input_description -distribution {Quantum Espresso} -package PWscf -program projw
}
}
var diag_basis -type LOGICAL {
default { .false. }
info {
if @b .false. the projections of Kohn-Sham states are
done on the orthogonalized atomic orbitals
in the global XYZ coordinate frame.
if @b .true. the projections of Kohn-Sham states are
done on the orthogonalized atomic orbitals
that are rotated to the basis in which the
atomic occupation matrix is diagonal
(i.e. local XYZ coordinate frame).
}
}
var pawproj -type LOGICAL {
default { .false. }
info {

View File

@ -26,7 +26,7 @@ SUBROUTINE atomic_wfc_nc_proj (ik, wfcatom)
USE wvfct, ONLY : npwx, nbnd
USE uspp_param, ONLY : upf, nwfcm
USE noncollin_module, ONLY : noncolin, npol, angle1, angle2
USE spin_orb, ONLY : lspinorb, rot_ylm, fcoef, lmaxx
USE spin_orb, ONLY : lspinorb, rot_ylm, lmaxx
!
IMPLICIT NONE
!

View File

@ -42,7 +42,7 @@ PROGRAM plotband
real :: emin = 1.e10, emax =-1.e10, etic, eref, deltaE, Ef
INTEGER, PARAMETER :: max_lines=99
INTEGER, PARAMETER :: max_lines=999
real :: mine, dxmod, dxmod_save
INTEGER :: point(max_lines+1), nrap(max_lines)
INTEGER :: ilines, irap, ibnd, ipoint, jnow
@ -75,15 +75,15 @@ PROGRAM plotband
ENDIF
filename1=trim(filename)//".rap"
!!! replace with inquire statement?
exist_rap=.true.
OPEN(unit=21,file=filename1,form='formatted',status='old',err=100,iostat=ios)
100 IF (ios /= 0) THEN
exist_rap=.false.
ENDIF
!!!
!
inquire(file=filename1, exist=exist_rap)
IF (exist_rap) THEN
OPEN(unit=21,file=filename1,form='formatted',status='old',iostat=ios)
IF ( ios /= 0 ) THEN
WRITE(6,'("error opening file with representations")')
exist_rap=.false.
GO TO 100
ENDIF
READ (21, plot_rap, iostat=ios)
IF (nks_rap/=nks.or.nbnd_rap/=nbnd.or.ios/=0) THEN
WRITE(6,'("file with representations not compatible with bands")')
@ -91,6 +91,8 @@ PROGRAM plotband
ENDIF
ENDIF
!
100 CONTINUE
!
ALLOCATE (e(nbnd,nks))
ALLOCATE (k(3,nks), e_in(nks), kx(nks), npoints(nks), high_symmetry(nks))
ALLOCATE (is_in_range(nbnd))
@ -105,7 +107,6 @@ PROGRAM plotband
ENDIF
!!!
filename2=trim(filename)//".proj"
exist_proj=.false.
INQUIRE(file=filename2,exist=exist_proj)
IF (exist_proj) THEN
OPEN(UNIT=22, FILE=filename2, FORM='formatted', STATUS='old', IOSTAT=ios)
@ -136,7 +137,6 @@ PROGRAM plotband
ENDIF
!!!
high_symmetry=.false.
DO n=1,nks
@ -225,7 +225,7 @@ PROGRAM plotband
dxmod=sqrt ( (k(1,n)-k(1,n-1))**2 + &
(k(2,n)-k(2,n-1))**2 + &
(k(3,n)-k(3,n-1))**2 )
IF (dxmod > 5*dxmod_save) THEN
IF (dxmod >10*dxmod_save) THEN
!
! A big jump in dxmod is a sign that the point k(:,n) and k(:,n-1)
! are quite distant and belong to two different lines. We put them on
@ -238,7 +238,6 @@ PROGRAM plotband
! 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

View File

@ -36,6 +36,7 @@ PROGRAM do_projwfc
USE cell_base, ONLY : at, bg
USE start_k, ONLY : k1, k2, k3, nk1, nk2, nk3
USE lsda_mod, ONLY : lsda
USE control_flags, ONLY : gamma_only
!
IMPLICIT NONE
!
@ -46,7 +47,7 @@ PROGRAM do_projwfc
REAL (DP) :: Emin, Emax, DeltaE, degauss1, ef_0
INTEGER :: nks2, ngauss1, ios
LOGICAL :: lwrite_overlaps, lbinary_data, needwf = .TRUE.
LOGICAL :: lsym, kresolveddos, tdosinboxes, plotboxes, pawproj
LOGICAL :: lsym, kresolveddos, tdosinboxes, plotboxes, pawproj, diag_basis
INTEGER, PARAMETER :: N_MAX_BOXES = 999
INTEGER :: n_proj_boxes, irmin(3,N_MAX_BOXES), irmax(3,N_MAX_BOXES)
LOGICAL :: lgww !if .true. use GW QP energies from file bands.dat
@ -54,7 +55,7 @@ PROGRAM do_projwfc
NAMELIST / projwfc / outdir, prefix, ngauss, degauss, lsym, &
Emin, Emax, DeltaE, filpdos, filproj, filowdin, lgww, &
kresolveddos, tdosinboxes, n_proj_boxes, irmin, irmax, plotboxes, &
lwrite_overlaps, lbinary_data, pawproj, lforcet, ef_0
lwrite_overlaps, lbinary_data, pawproj, lforcet, ef_0, diag_basis
!
! initialise environment
!
@ -75,6 +76,7 @@ PROGRAM do_projwfc
DeltaE = 0.01d0
ngauss = 0
lsym = .true.
diag_basis = .false.
degauss= 0.d0
lgww = .false.
pawproj= .false.
@ -117,6 +119,7 @@ PROGRAM do_projwfc
CALL mp_bcast( degauss1, ionode_id, intra_image_comm )
CALL mp_bcast( DeltaE, ionode_id, intra_image_comm )
CALL mp_bcast( lsym, ionode_id, intra_image_comm )
CALL mp_bcast( diag_basis,ionode_id, intra_image_comm )
CALL mp_bcast( Emin, ionode_id, intra_image_comm )
CALL mp_bcast( Emax, ionode_id, intra_image_comm )
CALL mp_bcast( lwrite_overlaps, ionode_id, intra_image_comm )
@ -148,6 +151,12 @@ PROGRAM do_projwfc
IF ( lforcet .AND. tdosinboxes ) CALL errore ('projwfc','incompatible options',3)
IF ( lforcet .AND. lsym ) CALL errore ('projwfc','incompatible options',4)
!
IF (diag_basis) THEN
IF ( pawproj ) CALL errore ('projwfc','diag_basis=.true. is not available for pawproj=.true.',1)
IF ( noncolin ) CALL errore ('projwfc','diag_basis=.true. is not implemented for noncolin=.true.',1)
IF ( gamma_only ) CALL errore ('projwfc','diag_basis=.true. is not implemented for gamma_only',1)
ENDIF
!
! Tetrahedron method
!
IF ( ltetra .AND. degauss1==0.d0 ) THEN
@ -210,7 +219,7 @@ PROGRAM do_projwfc
ELSE IF ( pawproj ) THEN
CALL projwave_paw ( )
ELSE
CALL projwave(filproj, filowdin, lsym, lwrite_overlaps)
CALL projwave(filproj, filowdin, lsym, diag_basis, lwrite_overlaps)
IF ( lforcet ) CALL force_theorem ( ef_0, filproj )
ENDIF
!
@ -269,27 +278,42 @@ SUBROUTINE get_et_from_gww ( nbnd, et )
ENDIF
END SUBROUTINE get_et_from_gww
!
SUBROUTINE print_lowdin ( unit, nat, lmax_wfc, nspin, charges, charges_lm )
SUBROUTINE print_lowdin ( unit, nat, lmax_wfc, nspin, diag_basis, charges, charges_lm )
!
USE kinds, ONLY : dp
USE kinds, ONLY : DP
USE io_global, ONLY : stdout, ionode
USE klist, ONLY: nelec
USE klist, ONLY : nelec
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: unit, nat, lmax_wfc, nspin
LOGICAL, INTENT(IN) :: diag_basis
REAL(DP), INTENT(in) :: charges (nat, 0:lmax_wfc, nspin )
REAL(DP), INTENT(in), OPTIONAL :: charges_lm (nat, 0:lmax_wfc, 1:2*lmax_wfc+1, nspin )
!
INTEGER :: is, l, m, na
REAL(DP) :: totcharge(2), psum
CHARACTER (len=1) :: l_label(0:3)=(/'s','p','d','f'/)
CHARACTER (len=7) :: lm_label(1:7,1:3)=reshape( (/ &
CHARACTER (len=7) :: lm_label(1:7,1:3)
CHARACTER (len=7) :: lm_label_global_frame(1:7,1:3)=reshape( (/ &
'z ','x ','y ',' ',' ',' ',' ', &
'z2 ','xz ','yz ','x2-y2 ','xy ',' ',' ', &
'z3 ','xz2 ','yz2 ','zx2-zy2','xyz ','x3-3xy2','3yx2-y3' /), (/7,3/) )
! TODO: think of a better way how to automatically label states in the diagonalized basis
! (i.e. eg, t2g, etc.)
CHARACTER (len=7) :: lm_label_diag(1:7,1:3)=reshape( (/ &
'1 ','2 ','3 ',' ',' ',' ',' ', &
'1 ','2 ','3 ','4 ','5 ',' ',' ', &
'1 ','2 ','3 ','4 ','5 ','6 ','7 ' /), (/7,3/) )
!
IF ( ionode ) THEN
!
IF (diag_basis) THEN
lm_label = lm_label_diag
ELSE
lm_label = lm_label_global_frame
ENDIF
!
WRITE( unit, '(/"Lowdin Charges: "/)')
!
DO na = 1, nat
@ -302,7 +326,7 @@ SUBROUTINE print_lowdin ( unit, nat, lmax_wfc, nspin, charges, charges_lm )
IF (l /= 0 .AND. present(charges_lm)) THEN
DO m = 1, 2*l+1
WRITE( unit,'(A1,A,"=",F8.4,", ")',advance='no') &
l_label(l), trim(lm_label(m,l)), charges_lm(na,l,m,1)
l_label(l), trim(lm_label(m,l)), charges_lm(na,l,m,1)
ENDDO
ENDIF
WRITE(unit,*)
@ -697,7 +721,7 @@ SUBROUTINE sym_proj_nc ( proj0, proj_out )
!
END SUBROUTINE sym_proj_nc
!-----------------------------------------------------------------------
SUBROUTINE print_proj ( lmax_wfc, proj, lowdin_unit )
SUBROUTINE print_proj ( lmax_wfc, proj, lowdin_unit, diag_basis )
!-----------------------------------------------------------------------
!
USE kinds, ONLY : DP
@ -715,6 +739,7 @@ SUBROUTINE print_proj ( lmax_wfc, proj, lowdin_unit )
IMPLICIT NONE
INTEGER, INTENT(in) :: lmax_wfc, lowdin_unit
REAL(DP), INTENT(IN) :: proj(natomwfc,nbnd,nkstot)
LOGICAL, INTENT(IN) :: diag_basis
!
INTEGER :: nspin0, nwfc, ibnd, i, j, ik, na, l, m
INTEGER, ALLOCATABLE :: idx(:)
@ -725,9 +750,10 @@ SUBROUTINE print_proj ( lmax_wfc, proj, lowdin_unit )
CHARACTER (len=1) :: plus
!
INTERFACE
SUBROUTINE print_lowdin ( unit, nat, lmax_wfc, nspin, charges, charges_lm )
SUBROUTINE print_lowdin ( unit, nat, lmax_wfc, nspin, diag_basis, charges, charges_lm )
IMPORT :: DP
INTEGER, INTENT(IN) :: unit, nat, lmax_wfc, nspin
LOGICAL, INTENT(IN) :: diag_basis
REAL(DP), INTENT(in) :: charges (nat, 0:lmax_wfc, nspin )
REAL(DP), INTENT(in), OPTIONAL :: charges_lm (nat, 0:lmax_wfc, 1:2*lmax_wfc+1, nspin )
END SUBROUTINE print_lowdin
@ -848,10 +874,10 @@ SUBROUTINE print_proj ( lmax_wfc, proj, lowdin_unit )
ENDDO
!
IF ( nspin /= 4 ) THEN
CALL print_lowdin ( lowdin_unit, nat, lmax_wfc, nspin, charges, charges_lm )
CALL print_lowdin ( lowdin_unit, nat, lmax_wfc, nspin, diag_basis, charges, charges_lm )
DEALLOCATE (charges_lm)
ELSE
CALL print_lowdin ( lowdin_unit, nat, lmax_wfc, nspin0, charges )
CALL print_lowdin ( lowdin_unit, nat, lmax_wfc, nspin0, diag_basis, charges )
END IF
DEALLOCATE (charges)
!
@ -1093,7 +1119,7 @@ END FUNCTION compute_mj
! projwave with distributed matrixes
!
!-----------------------------------------------------------------------
SUBROUTINE projwave( filproj, filowdin, lsym, lwrite_ovp )
SUBROUTINE projwave( filproj, filowdin, lsym, diag_basis, lwrite_ovp )
!-----------------------------------------------------------------------
!
USE kinds, ONLY : DP
@ -1108,16 +1134,17 @@ SUBROUTINE projwave( filproj, filowdin, lsym, lwrite_ovp )
USE uspp, ONLY : nkb, vkb
USE becmod, ONLY : bec_type, becp, calbec, allocate_bec_type, deallocate_bec_type
USE io_files, ONLY : prefix, restart_dir, tmp_dir
USE control_flags, ONLY : gamma_only, use_para_diag
USE control_flags, ONLY : gamma_only, use_para_diag, io_level
USE pw_restart_new,ONLY : read_collected_wfc
USE wavefunctions, ONLY : evc
!
USE projections, ONLY: nlmchi, fill_nlmchi, proj, proj_aux, ovps_aux
!
USE io_files, ONLY: nd_nmbr
USE io_files, ONLY: nd_nmbr, nwordatwfc
USE mp, ONLY: mp_bcast
USE mp_pools, ONLY: me_pool, root_pool, intra_pool_comm
USE uspp_init, ONLY : init_us_2
USE buffers, ONLY : open_buffer, save_buffer, get_buffer, close_buffer
!
IMPLICIT NONE
!
@ -1125,6 +1152,7 @@ SUBROUTINE projwave( filproj, filowdin, lsym, lwrite_ovp )
!
CHARACTER (len=*), INTENT(IN) :: filproj, filowdin
LOGICAL, INTENT(IN) :: lsym
LOGICAL, INTENT(IN) :: diag_basis
LOGICAL, INTENT(INOUT) :: lwrite_ovp
!
LOGICAL :: ionode_pool
@ -1139,9 +1167,9 @@ SUBROUTINE projwave( filproj, filowdin, lsym, lwrite_ovp )
REAL (DP), ALLOCATABLE ::roverlap_d(:,:)
!
INTEGER :: nksinit, nkslast
LOGICAL :: freeswfcatom
LOGICAL :: freeswfcatom, exst, lrotated
!
INTEGER :: iunaux, lowdin_unit
INTEGER :: iunaux, lowdin_unit, iuwfc
INTEGER, EXTERNAL :: find_free_unit
CHARACTER(len=256) :: auxname
!
@ -1170,6 +1198,11 @@ SUBROUTINE projwave( filproj, filowdin, lsym, lwrite_ovp )
CALL fill_nlmchi ( natomwfc, lmax_wfc )
!
ALLOCATE( proj (natomwfc, nbnd, nkstot) )
IF (diag_basis) THEN
iuwfc = find_free_unit()
nwordatwfc = npwx*natomwfc*npol
CALL open_buffer( iuwfc, 'wfcrot', nwordatwfc, io_level, exst )
ENDIF
!
IF (.not. ALLOCATED(swfcatom)) THEN
ALLOCATE(swfcatom (npwx*npol , natomwfc ) )
@ -1340,6 +1373,12 @@ SUBROUTINE projwave( filproj, filowdin, lsym, lwrite_ovp )
ENDIF
DEALLOCATE( overlap_d )
!
! Save O^{-1/2} \hat S | phi_j> for a given k
IF (diag_basis) THEN
CALL save_buffer (wfcatom, nwordatwfc, iuwfc, ik)
GOTO 100
ENDIF
!
! make the projection <psi_i| O^{-1/2} \hat S | phi_j>,
! symmetrize the projections if required
!
@ -1374,8 +1413,37 @@ SUBROUTINE projwave( filproj, filowdin, lsym, lwrite_ovp )
DEALLOCATE (proj0)
!
ENDIF
! on k-points
ENDDO
!
100 CONTINUE
!
ENDDO ! ik
!
! Compute the projections in a local frame if requested
!
IF (diag_basis) THEN
CALL rotate_basis (iuwfc, lrotated)
ALLOCATE (proj0(natomwfc,nbnd))
DO ik = 1, nks
npw = ngk(ik)
! Read the KSwavefunction evc at this k
CALL read_collected_wfc (restart_dir(), ik, evc)
! Read the rotated orbital wfcatom at this k
CALL get_buffer (wfcatom, nwordatwfc, iuwfc, ik)
! Calculate proj0 = <wfcatom|evc> at this k
CALL calbec (npw, wfcatom, evc, proj0)
IF (ionode_pool) WRITE( iunaux ) proj0
! Symmetrization (lsym=.true.) must not be used here
! if the rotation of orbitals was done because otherwise
! the results will be wrong
IF (lsym .AND. .NOT.lrotated) THEN
CALL sym_proj_k (proj0, proj(:,:,ik))
ELSE
proj(:,:,ik)=abs(proj0(:,:))**2
ENDIF
ENDDO
DEALLOCATE (proj0)
CALL close_buffer (iuwfc, 'DELETE')
ENDIF
!
DEALLOCATE (e)
DEALLOCATE (wfcatom)
@ -1402,7 +1470,7 @@ SUBROUTINE projwave( filproj, filowdin, lsym, lwrite_ovp )
lowdin_unit = stdout
END IF
!
CALL print_proj( lmax_wfc, proj, lowdin_unit )
CALL print_proj( lmax_wfc, proj, lowdin_unit, diag_basis )
!
IF (TRIM(filowdin) /= ' ') CLOSE( unit=lowdin_unit )
!
@ -1715,3 +1783,364 @@ CONTAINS
!
END SUBROUTINE projwave
!
SUBROUTINE rotate_basis (iuwfc, lrotated)
!
!! This routine rotates the ortho-atomic orbitals to the basis
!! where the occupation matrix is diagonal. This is useful e.g.
!! for determining the eg and t2g states in chemical elements
!! containing the d-type electrons.
!!
!! Input: wfcatom = O^{-1/2} \hat S | phi_j>
!! Output: wfcatom = T O^{-1/2} \hat S | phi_j>
!!
!! here: T is the transformation matrix composed of the eigenvectors
!! of the occupation matrix;
!! O is the orbital overlap matrix;
!! S is the ultrasoft/PAW operator;
!! |phi_j> is the j-th atomic orbital.
!!
! Written by I. Timrov (September 2021).
!
USE kinds, ONLY : DP
USE basis, ONLY : natomwfc
USE wvfct, ONLY : nbnd, wg, npwx
USE symm_base, ONLY : nsym, irt, d1, d2, d3
USE projections, ONLY : nlmchi
USE ions_base, ONLY : nat
USE klist, ONLY : nks, ngk
USE mp, ONLY : mp_sum
USE mp_pools, ONLY : inter_pool_comm
USE lsda_mod, ONLY : lsda, current_spin, nspin, isk
USE io_global, ONLY : stdout
USE pw_restart_new, ONLY : read_collected_wfc
USE wavefunctions, ONLY : evc
USE io_files, ONLY : restart_dir, nwordatwfc
USE becmod, ONLY : calbec
USE buffers, ONLY : save_buffer, get_buffer
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: iuwfc
! unit from which to read wfcatom
LOGICAL, INTENT(OUT) :: lrotated
!
INTEGER :: ibnd, nwfc, na, nb, l, m, m1, m2, m0, m00, &
n, ik, is, isym, ldim, npw, ig
COMPLEX(DP), ALLOCATABLE :: proj(:,:), proj_aux(:,:,:,:), f(:,:) , vet(:,:), &
wfcatom_aux(:,:,:,:), wfcatomk_rot(:), &
transformation_matrix(:,:,:,:,:), wfcatom(:,:)
REAL(DP), ALLOCATABLE :: ns(:,:,:,:,:), nr(:,:,:,:,:), lambda(:)
INTEGER, ALLOCATABLE :: n_max(:), orbital_quantum_number(:,:)
LOGICAL, ALLOCATABLE :: diagonalize(:,:)
INTEGER, PARAMETER :: nmax = 10, & ! max number of shells per atom
mmax = 7 ! max number of m's per shell
REAL(DP) :: psum
!
lrotated = .TRUE.
!
ALLOCATE (n_max(nat))
ALLOCATE (orbital_quantum_number(nat,nmax))
ALLOCATE (proj_aux(nat,nmax,mmax,nbnd))
ALLOCATE (nr(nat,nmax,mmax,mmax,nspin))
ALLOCATE (proj(natomwfc,nbnd))
ALLOCATE (wfcatom(npwx,natomwfc))
ALLOCATE (diagonalize(nat,nmax))
nr = 0.0d0
!
DO ik = 1, nks
!
IF (lsda) current_spin = isk(ik)
!
npw = ngk(ik)
!
! Read the Kohn-Sham wavefunctions evc at this k
CALL read_collected_wfc (restart_dir(), ik, evc)
!
! Read the orbitals wfcatom at this k
CALL get_buffer (wfcatom, nwordatwfc, iuwfc, ik)
!
! Calculate proj = <wfcatom|evc> at this k
CALL calbec (npw, wfcatom, evc, proj)
!
! Rewrite the arrays proj and wfcatomk in a different way
! (in terms of n,l,m)
!
n_max(:) = 0
orbital_quantum_number = -1
proj_aux = (0.0d0, 0.0d0)
DO nwfc = 1, natomwfc
na= nlmchi(nwfc)%na
n = nlmchi(nwfc)%n
l = nlmchi(nwfc)%l
m = nlmchi(nwfc)%m
n_max(na) = MAX(n,n_max(na))
orbital_quantum_number(na,n) = l
DO ibnd = 1, nbnd
proj_aux(na,n,m,ibnd) = proj(nwfc,ibnd)
ENDDO
ENDDO
!
! Compute the occupation matrix
!
DO na = 1, nat
! n are the shells of atom na
DO n = 1, n_max(na)
l = orbital_quantum_number(na,n)
ldim = 2*l+1
DO m1 = 1, ldim
DO m2 = 1, ldim
DO ibnd = 1, nbnd
nr(na,n,m1,m2,current_spin) = nr(na,n,m1,m2,current_spin) + &
wg(ibnd,ik) * &
DBLE(proj_aux(na,n,m2,ibnd) * &
CONJG(proj_aux(na,n,m1,ibnd)))
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
!
ENDDO
!
CALL mp_sum( nr, inter_pool_comm )
!
DEALLOCATE (proj)
DEALLOCATE (proj_aux)
ALLOCATE (ns(nat,nmax,mmax,mmax,nspin))
ns = 0.0d0
!
! Symmetrization
!
CALL d_matrix (d1, d2, d3)
DO na = 1, nat
DO n = 1, n_max(na)
l = orbital_quantum_number(na,n)
ldim = 2*l+1
DO is = 1, nspin
DO m1 = 1, ldim
DO m2 = 1, ldim
DO isym = 1, nsym
nb = irt (isym, na)
DO m0 = 1, ldim
DO m00 = 1, ldim
IF (l == 0) THEN
ns(na,n,m1,m2,is) = ns(na,n,m1,m2,is) + &
nr(nb,n,m0,m00,is) / nsym
ELSEIF (l == 1) THEN
ns(na,n,m1,m2,is) = ns(na,n,m1,m2,is) + &
d1(m0 ,m1,isym) * nr(nb,n,m0,m00,is) * &
d1(m00,m2,isym) / nsym
ELSEIF (l == 2) THEN
ns(na,n,m1,m2,is) = ns(na,n,m1,m2,is) + &
d2(m0 ,m1,isym) * nr(nb,n,m0,m00,is) * &
d2(m00,m2,isym) / nsym
ELSEIF (l == 3) THEN
ns(na,n,m1,m2,is) = ns(na,n,m1,m2,is) + &
d3(m0 ,m1,isym) * nr(nb,n,m0,m00,is) * &
d3(m00,m2,isym) / nsym
ELSE
CALL errore( 'rotate_basis', &
'angular momentum not implemented', &
ABS(l) )
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
!
DEALLOCATE (nr)
!
! Now we make the matrix ns strictly Hermitean
!
DO na = 1, nat
DO n = 1, n_max(na)
l = orbital_quantum_number(na,n)
ldim = 2*l+1
DO is = 1, nspin
DO m1 = 1, ldim
DO m2 = m1, ldim
psum = ABS( ns(na,n,m1,m2,is) - ns(na,n,m2,m1,is) )
IF (psum > 1.d-10) THEN
WRITE( stdout, * ) na, n, m1, m2, is
WRITE( stdout, * ) ns(na,n,m1,m2,is)
WRITE( stdout, * ) ns(na,n,m2,m1,is)
CALL errore( 'rotate_basis', 'non Hermitean matrix', 1 )
ELSE
ns(na,n,m1,m2,is) = 0.5d0 * (ns(na,n,m1,m2,is) + &
ns(na,n,m2,m1,is) )
ns(na,n,m2,m1,is) = ns(na,n,m1,m2,is)
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
!
! Before diagonalizng the occupation matrix we need to check
! whether there are non-zero off-diagonal matrix elements. If the
! matrix is already diagonal, then we do not diagonalize it.
!
diagonalize(:,:) = .FALSE.
DO na = 1, nat
DO n = 1, n_max(na)
l = orbital_quantum_number(na,n)
IF (l>0) THEN
ldim = 2*l+1
DO is = 1, nspin
DO m1 = 1, ldim-1
DO m2 = m1+1, ldim
! If any off-diagonal element is non-zero then proceed
IF (ABS(ns(na,n,m1,m2,is)) > 1.d-3) THEN
diagonalize(na,n) = .TRUE.
GO TO 10
ENDIF
ENDDO
ENDDO
ENDDO
ENDIF
10 CONTINUE
ENDDO
ENDDO
!
! Diagonalization of the occupation matrix ns
!
ALLOCATE (transformation_matrix(nat,nmax,mmax,mmax,nspin))
DO na = 1, nat
WRITE( stdout,'(/5x,21("-")," ATOM ",i4,1x,22("-"))') na
DO n = 1, n_max(na)
l = orbital_quantum_number(na,n)
ldim = 2*l+1
WRITE( stdout,'(/5x,"Orbital quantum number l = ", i1)') l
ALLOCATE (f(ldim,ldim), vet(ldim,ldim), lambda(ldim))
DO is = 1, nspin
IF (nspin /= 1) WRITE( stdout,'(5x,"SPIN ",i2)') is
DO m1 = 1, ldim
DO m2 = 1, ldim
f(m1,m2) = ns(na,n,m1,m2,is)
ENDDO
ENDDO
WRITE( stdout,'(5x,"occupation matrix ns (before diag.):")')
DO m1 = 1, ldim
WRITE( stdout,'(5x,7f7.3)') ( DBLE(ns(na,n,m1,m2,is)), m2=1, ldim )
ENDDO
IF (diagonalize(na,n)) THEN
! Diagonalize ns
CALL cdiagh( ldim, f, ldim, lambda, vet )
DO m1 = 1, ldim
DO m2 = 1, ldim
transformation_matrix(na,n,m1,m2,is) = vet(m1,m2)
ENDDO
ENDDO
ELSE
! Do not diagonalize ns
DO m1 = 1, ldim
lambda(m1) = f(m1,m1)
ENDDO
transformation_matrix(na,n,:,:,is) = (0.0d0, 0.0d0)
DO m1 = 1, ldim
transformation_matrix(na,n,m1,m1,is) = (1.0d0, 0.d0)
ENDDO
ENDIF
WRITE( stdout,'(5x,"eigenvalues:")')
WRITE( stdout,'(5x,7f7.3)') (lambda(m1), m1=1, ldim)
WRITE( stdout,'(5x,"eigenvectors (columns):")')
DO m1 = 1, ldim
WRITE( stdout,'(5x,7f7.3)') ( DBLE(transformation_matrix(na,n,m1,m2,is)), m2=1, ldim )
ENDDO
ENDDO
DEALLOCATE (f, vet, lambda)
ENDDO
ENDDO
!
IF (.NOT. ANY(diagonalize(:,:))) THEN
WRITE( stdout,'(/5x,"All occupation matrices are already diagonal. No rotation will be performed!")')
lrotated = .FALSE.
GO TO 11
ENDIF
!
WRITE( stdout,'(/5x,"Rotating the orbitals to the diagonal representation of ns...")')
!
! Rotate the orbitals using the transformation matrix composed of
! the eigenvectors of the occupation matrix
!
ALLOCATE (wfcatom_aux(npwx,nat,nmax,mmax))
DO ik = 1, nks
!
npw = ngk(ik)
!
IF (lsda) current_spin = isk(ik)
!
! Read the orbitals wfcatom at this k
CALL get_buffer (wfcatom, nwordatwfc, iuwfc, ik)
!
! Rewrite the array wfcatom in a different way (in terms of na,n,m)
!
DO nwfc = 1, natomwfc
na= nlmchi(nwfc)%na
n = nlmchi(nwfc)%n
m = nlmchi(nwfc)%m
DO ig = 1, npw
wfcatom_aux(ig,na,n,m) = wfcatom(ig,nwfc)
ENDDO
ENDDO
!
nwfc = 0
wfcatom(:,:) = (0.0d0, 0.0d0)
DO na = 1, nat
DO n = 1, n_max(na)
l = orbital_quantum_number(na,n)
ldim = 2*l+1
!
! Rotate the orbitals
!
ALLOCATE (wfcatomk_rot(ldim))
DO ig = 1, npw
wfcatomk_rot(:) = (0.0d0, 0.0d0)
DO m1 = 1, ldim
DO m2 = 1, ldim
wfcatomk_rot(m1) = wfcatomk_rot(m1) + &
transformation_matrix(na,n,m2,m1,current_spin) * &
wfcatom_aux(ig,na,n,m2)
ENDDO
ENDDO
DO m = 1, ldim
wfcatom_aux(ig,na,n,m) = wfcatomk_rot(m)
ENDDO
ENDDO
DEALLOCATE (wfcatomk_rot)
!
! Copy wfcatom_aux to the original array wfcatom
!
DO m = 1, ldim
nwfc = nwfc + 1
DO ig = 1, npw
wfcatom(ig,nwfc) = wfcatom_aux(ig,na,n,m)
ENDDO
ENDDO
!
ENDDO
ENDDO
!
! Write the rotated orbitals to file
CALL save_buffer (wfcatom, nwordatwfc, iuwfc, ik)
!
ENDDO
!
DEALLOCATE (wfcatom_aux)
!
11 CONTINUE
!
DEALLOCATE (transformation_matrix)
DEALLOCATE (n_max)
DEALLOCATE (orbital_quantum_number)
DEALLOCATE (ns)
DEALLOCATE (wfcatom)
DEALLOCATE (diagonalize)
!
RETURN
!
END SUBROUTINE rotate_basis

View File

@ -3583,7 +3583,6 @@ SUBROUTINE write_kih (kih_file_name, vxc_hybrid_file_name, diag_nmin, diag_nmax,
USE ions_base, ONLY : nat, ityp, ntyp => nsp
USE wvfct, ONLY: npwx, current_k !FZ:
USE uspp_param, ONLY: upf, nh !FZ:
USE spin_orb, ONLY: lspinorb ,domag !FZ: added domag
USE becmod, ONLY : bec_type, becp, calbec, &
allocate_bec_type, deallocate_bec_type
USE fft_base, ONLY : dffts !FZ: test

View File

@ -13,7 +13,7 @@ CONTAINS
nspin, nelec, ef, xk, wk, et, projs, ovps )
!-----------------------------------------------------------------------
USE kinds, ONLY : dp
USE xmltools, ONLY : xml_openfile, xml_closefile, xmlr_readtag,&
USE xmltools, ONLY : xml_open_file, xml_closefile, xmlr_readtag,&
xmlr_opentag, xmlr_closetag, get_attr
!
IMPLICIT NONE
@ -38,7 +38,7 @@ CONTAINS
CALL infomsg ('read_xml_proj', 'xml data file not found')
END IF
!
iun = xml_openfile (filename)
iun = xml_open_file (filename)
IF ( iun == -1 ) THEN
ierr = 2
CALL infomsg ('read_xml_proj', 'xml data file not readable')

View File

@ -29,7 +29,6 @@ SUBROUTINE sum_band_kin(kin_r)
USE uspp_param, ONLY : upf, nh, nhm
USE wavefunctions, ONLY : evc, psic, psic_nc
USE noncollin_module, ONLY : noncolin, npol, nspin_mag
USE spin_orb, ONLY : lspinorb, fcoef
USE wvfct, ONLY : nbnd, npwx, wg, et
USE mp_pools, ONLY : inter_pool_comm
USE mp_bands, ONLY : inter_bgrp_comm, intra_bgrp_comm, nbgrp

View File

@ -16,7 +16,7 @@ SUBROUTINE write_xml_proj (filename, projs, lwrite_ovp, ovps )
USE lsda_mod, ONLY : nspin
USE ener, ONLY : ef
USE wvfct, ONLY : et, nbnd
USE xmltools, ONLY : xml_openfile, xml_closefile, xmlw_writetag,&
USE xmltools, ONLY : xml_open_file, xml_closefile, xmlw_writetag,&
xmlw_opentag, xmlw_closetag, add_attr
IMPLICIT NONE
@ -29,7 +29,7 @@ SUBROUTINE write_xml_proj (filename, projs, lwrite_ovp, ovps )
INTEGER :: ik, ik_eff, is, nwfc, ibnd, nspin_lsda, num_k_points
!
!
iunpun = xml_openfile ( TRIM( restart_dir() )//TRIM(filename) )
iunpun = xml_open_file ( TRIM( restart_dir() )//TRIM(filename) )
IF ( iunpun == -1 ) RETURN
!
nspin_lsda = 1

View File

@ -22,8 +22,7 @@ SUBROUTINE atomic_wfc( ik, wfcatom )
USE wvfct, ONLY : npwx
USE uspp_param, ONLY : upf, nwfcm
USE noncollin_module, ONLY : noncolin, npol, angle1, angle2
USE spin_orb, ONLY : lspinorb, rot_ylm, fcoef, lmaxx, domag, &
starting_spin_angle
USE spin_orb, ONLY : rot_ylm, lmaxx, domag, starting_spin_angle
!
implicit none
INTEGER, INTENT(IN) :: ik

View File

@ -22,8 +22,7 @@ SUBROUTINE atomic_wfc_gpu( ik, wfcatom_d )
USE wvfct, ONLY : npwx
USE uspp_param, ONLY : upf, nwfcm
USE noncollin_module, ONLY : noncolin, npol, angle1, angle2
USE spin_orb, ONLY : lspinorb, rot_ylm, fcoef, lmaxx, domag, &
starting_spin_angle
USE spin_orb, ONLY : rot_ylm, lmaxx, domag, starting_spin_angle
!
IMPLICIT NONE
!

View File

@ -158,7 +158,7 @@ SUBROUTINE force_ew( alat, nat, ntyp, ityp, zv, at, bg, tau, omega, &
IF ( mykey > 0 ) GOTO 100
rmax = 5.d0 / (SQRT(alpha) * alat)
!
! with this choice terms up to ZiZj*erfc(5) are counted (erfc(5)=2x10^-1
! with this choice terms up to ZiZj*erfc(5) are counted (erfc(5)=2x10^-12)
!
!$omp parallel do default(none) shared(na_s, na_e, nat, at, bg, ityp, zv, alpha, forceion, alat, rmax, tau)&
!$omp &private(nb, dtau, nrm, r, r2, rr, fact)

View File

@ -24,7 +24,6 @@ SUBROUTINE force_us( forcenl )
USE symme, ONLY : symvector
USE wavefunctions, ONLY : evc
USE noncollin_module, ONLY : npol, noncolin
USE spin_orb, ONLY : lspinorb
USE io_files, ONLY : iunwfc, nwordwfc
USE buffers, ONLY : get_buffer
USE becmod, ONLY : calbec, becp, bec_type, allocate_bec_type, &

View File

@ -25,7 +25,6 @@ SUBROUTINE force_us_gpu( forcenl )
USE wavefunctions, ONLY : evc
USE wavefunctions_gpum, ONLY : evc_d, using_evc, using_evc_d
USE noncollin_module, ONLY : npol, noncolin
USE spin_orb, ONLY : lspinorb
USE io_files, ONLY : iunwfc, nwordwfc
USE buffers, ONLY : get_buffer
USE becmod, ONLY : calbec, becp, bec_type, allocate_bec_type, &

View File

@ -37,7 +37,6 @@ SUBROUTINE newq_gpu(vr,deeq_d,skip_vltot)
USE control_flags, ONLY : gamma_only
USE wavefunctions, ONLY : psic
USE wavefunctions_gpum, ONLY : psic_d
USE spin_orb, ONLY : lspinorb, domag
USE noncollin_module, ONLY : nspin_mag
USE mp_bands, ONLY : intra_bgrp_comm
USE mp_pools, ONLY : inter_pool_comm

View File

@ -1682,7 +1682,6 @@ MODULE paw_onecenter
!! potential in the spherical basis. It receives as input the charge
!! density and its variation.
!
USE spin_orb, ONLY : domag
USE noncollin_module, ONLY : nspin_mag
USE lsda_mod, ONLY : nspin
USE atom, ONLY : g => rgrid

View File

@ -300,8 +300,6 @@ SUBROUTINE atomic_wfc_nc_updown( ik, wfcatom )
USE wvfct, ONLY : npwx, nbnd
USE uspp_param, ONLY : upf, nwfcm
USE noncollin_module, ONLY : noncolin, npol, angle1, angle2
USE spin_orb, ONLY : lspinorb, rot_ylm, fcoef, domag, &
starting_spin_angle
!
IMPLICIT NONE
!

View File

@ -33,7 +33,7 @@ SUBROUTINE sum_band()
USE uspp_param, ONLY : nh, nhm
USE wavefunctions, ONLY : evc, psic, psic_nc
USE noncollin_module, ONLY : noncolin, npol, nspin_mag
USE spin_orb, ONLY : lspinorb, domag, fcoef
USE spin_orb, ONLY : domag
USE wvfct, ONLY : nbnd, npwx, wg, et, btype
USE mp_pools, ONLY : inter_pool_comm
USE mp_bands, ONLY : inter_bgrp_comm, intra_bgrp_comm, nbgrp

View File

@ -36,7 +36,7 @@ SUBROUTINE sum_band_gpu()
USE uspp_param, ONLY : upf, nh, nhm
USE wavefunctions, ONLY : evc, psic
USE noncollin_module, ONLY : noncolin, npol, nspin_mag
USE spin_orb, ONLY : lspinorb, domag, fcoef
USE spin_orb, ONLY : domag
USE wvfct, ONLY : nbnd, npwx, wg, et, btype
USE mp_pools, ONLY : inter_pool_comm
USE mp_bands, ONLY : inter_bgrp_comm, intra_bgrp_comm, nbgrp

View File

@ -131,7 +131,7 @@ SUBROUTINE write_ns
!
CALL cdiagh( ldim, f, ldim, lambda, vet )
!
IF (nspin /= 1) WRITE( stdout,'(5x"SPIN ",i2)') is
IF (nspin /= 1) WRITE( stdout,'(5x,"SPIN ",i2)') is
WRITE( stdout,'(5x,"eigenvalues:")')
WRITE( stdout,'(5x,7f7.3)') (lambda(m1), m1=1, ldim)
!
@ -191,7 +191,7 @@ SUBROUTINE write_ns
!
CALL cdiagh(ldim, f, ldim, lambda, vet)
!
IF (nspin /= 1) WRITE( stdout,'(5x"SPIN ",i2)') is
IF (nspin /= 1) WRITE( stdout,'(5x,"SPIN ",i2)') is
WRITE( stdout,'(5x,"eigenvalues:")')
WRITE( stdout,'(5x,7f7.3)') (lambda(m1), m1=1, ldim)
!
@ -584,7 +584,7 @@ SUBROUTINE write_nsg
!
CALL cdiagh(ldim1, f, ldim1, lambda, vet)
!
IF (nspin /= 1) WRITE( stdout,'(5x"SPIN ",i2)') is
IF (nspin /= 1) WRITE( stdout,'(5x,"SPIN ",i2)') is
WRITE( stdout,'(5x,"eigenvalues:")')
WRITE( stdout,'(5x,7f7.3)') (lambda(m1), m1=1, ldim1)
!

View File

@ -28,7 +28,6 @@ SUBROUTINE lr_apply_liouvillian_eels ( evc1, evc1_new, interaction )
USE noncollin_module, ONLY : noncolin, npol, nspin_mag
USE uspp, ONLY : okvan
USE mp_bands, ONLY : ntask_groups, me_bgrp
USE spin_orb, ONLY : domag
USE buffers, ONLY : get_buffer
USE qpoint, ONLY : ikks, ikqs, nksq
USE eqv, ONLY : evq, dpsi, dvpsi

View File

@ -35,7 +35,7 @@ SUBROUTINE dgcxc_unpol( length, r_in, s2_in, vrrx, vsrx, vssx, vrrc, vsrc, vssc
!! This routine computes the derivative of the exchange and correlation
!! potentials of GGA family.
!
USE qe_drivers_gga, ONLY: gcxc
USE qe_drivers_gga, ONLY: gcxc, gcxc_beef
!
IMPLICIT NONE
!
@ -92,6 +92,9 @@ SUBROUTINE dgcxc_unpol( length, r_in, s2_in, vrrx, vsrx, vssx, vrrc, vsrc, vssc
!
CALL gcxc( length*4, raux, s2aux, sx, sc, v1x, v2x, v1c, v2c )
!
IF ( igcx==43 .OR. igcc==14 ) CALL gcxc_beef( length*4, raux, s2aux, &
sx, sc, v1x, v2x, v1c, v2c )
!
! ... to avoid NaN in the next operations
WHERE( r_in<=small .OR. s2_in<=small )
dr = 1._DP ; ds = 1._DP ; s = 1._DP
@ -126,7 +129,7 @@ SUBROUTINE dgcxc_spin( length, r_in, g_in, vrrx, vrsx, vssx, vrrc, vrsc, &
!! This routine computes the derivative of the exchange and correlation
!! potentials in the spin-polarized case.
!
USE qe_drivers_gga, ONLY: gcx_spin, gcc_spin
USE qe_drivers_gga, ONLY: gcx_spin, gcc_spin, gcx_spin_beef, gcc_spin_beef
!
IMPLICIT NONE
!
@ -236,7 +239,8 @@ SUBROUTINE dgcxc_spin( length, r_in, g_in, vrrx, vrsx, vssx, vrrc, vrsc, &
raux(i8:f8,:) = r ; s2aux(i8:f8,:) = (s-dsdw)**2
!
!
CALL gcx_spin( length*8, raux, s2aux, sx, v1x, v2x )
IF ( igcx/=43 ) CALL gcx_spin( length*8, raux, s2aux, sx, v1x, v2x )
IF ( igcx==43 ) CALL gcx_spin_beef( length*8, raux, s2aux, sx, v1x, v2x )
!
! ... up
vrrx(:,1) = 0.5_DP * (v1x(i1:f1,1) - v1x(i2:f2,1)) / drup(:,1)
@ -307,7 +311,8 @@ SUBROUTINE dgcxc_spin( length, r_in, g_in, vrrx, vrsx, vssx, vrrc, vrsc, &
rtaux(i5:f5) = rt ; s2taux(i5:f5) = s2t ; zetaux(i5:f5) = zeta+dz
rtaux(i6:f6) = rt ; s2taux(i6:f6) = s2t ; zetaux(i6:f6) = zeta-dz
!
CALL gcc_spin( length*6, rtaux, zetaux, s2taux, sc, v1c, v2c )
IF ( igcc/=14 ) CALL gcc_spin( length*6, rtaux, zetaux, s2taux, sc, v1c, v2c )
IF ( igcc==14 ) CALL gcc_spin_beef( length*6, rtaux, zetaux, s2taux, sc, v1c, v2c )
!
vrrc(:,1) = 0.5_DP * (v1c(i1:f1,1) - v1c(i2:f2,1)) / dr * null_v(:,1)
vrrc(:,2) = 0.5_DP * (v1c(i1:f1,2) - v1c(i2:f2,2)) / dr * null_v(:,1)

View File

@ -25,14 +25,15 @@ MODULE qe_drivers_gga
!
PRIVATE
!
PUBLIC :: gcxc, gcx_spin, gcc_spin, gcc_spin_more
PUBLIC :: gcxc, gcx_spin, gcc_spin, gcc_spin_more, gcxc_beef, gcx_spin_beef, &
gcc_spin_beef
!
!
CONTAINS
!
!-----------------------------------------------------------------------
SUBROUTINE gcxc( length, rho_in, grho_in, sx_out, sc_out, v1x_out, &
v2x_out, v1c_out, v2c_out )
v2x_out, v1c_out, v2c_out )
!---------------------------------------------------------------------
!! Gradient corrections for exchange and correlation - Hartree a.u.
!! See comments at the beginning of module for implemented cases
@ -81,10 +82,6 @@ SUBROUTINE gcxc( length, rho_in, grho_in, sx_out, sc_out, v1x_out, &
!
!
#if defined(_OPENACC)
!
IF (igcx==43 .OR. igcc==14) CALL xclib_error( 'gcxc', 'BEEF not available with&
& OpenACC enabled', 1 )
!
!$acc data copyin(rho_in,grho_in), copyout(sx_out,sc_out,v1x_out,v2x_out,v1c_out,v2c_out)
!$acc parallel loop
#endif
@ -301,13 +298,6 @@ SUBROUTINE gcxc( length, rho_in, grho_in, sx_out, sc_out, v1x_out, &
v2x = (1.0_DP - exx_fraction) * v2x
ENDIF
!
#if !defined(_OPENACC)
CASE( 43 ) ! 'beefx' --- BEEF unavailable with OpenACC-- Will be enabled soon
! last parameter = 0 means do not add LDA (=Slater) exchange
! (espresso) will add it itself
CALL beefx(rho, grho, sx, v1x, v2x, 0)
#endif
!
CASE( 44 ) ! 'RPBE'
!
CALL pbex( rho, grho, 8, sx, v1x, v2x )
@ -379,13 +369,6 @@ SUBROUTINE gcxc( length, rho_in, grho_in, sx_out, sc_out, v1x_out, &
v2c = 0.871_DP * v2c
ENDIF
!
#if !defined(_OPENACC)
CASE( 14 ) ! BEEF unavailable with OpenACC-- Will be enabled soon
! last parameter 0 means: do not add lda contributions
! espresso will do that itself
call beeflocalcorr(rho, grho, sc, v1c, v2c, 0)
#endif
!
CASE DEFAULT
!
sc = 0.0_DP
@ -399,14 +382,12 @@ SUBROUTINE gcxc( length, rho_in, grho_in, sx_out, sc_out, v1x_out, &
v2x_out(ir) = v2x ; v2c_out(ir) = v2c
!
ENDDO
#if defined(_OPENMP) && !defined(_OPENACC)
#if defined(_OPENACC)
!$acc end data
#else
!$omp end do
!$omp end parallel
#endif
#if defined(_OPENACC)
!$acc end data
#endif
!
!
RETURN
!
@ -462,14 +443,9 @@ SUBROUTINE gcx_spin( length, rho_in, grho2_in, sx_tot, v1x_out, v2x_out )
sx_tot = 0.0_DP
!
#if defined(_OPENACC)
!
IF (igcx==43) CALL xclib_error( 'gcx_spin', 'BEEF not available with&
& OpenACC enabled', 1 )
!
!$acc data copyin(rho_in, grho2_in), copyout(sx_tot, v1x_out, v2x_out)
!$acc parallel loop
#endif
#if defined(_OPENMP) && !defined(_OPENACC)
#else
!$omp parallel if(ntids==1) default(none) &
!$omp private( rho_up, rho_dw, grho2_up, grho2_dw, rnull_up, rnull_dw, &
!$omp sx_up, sx_dw, sxsr_up, sxsr_dw, v1xsr_up, v1xsr_dw, &
@ -480,7 +456,6 @@ SUBROUTINE gcx_spin( length, rho_in, grho2_in, sx_tot, v1x_out, v2x_out )
!$omp gau_parameter )
!$omp do
#endif
!
DO ir = 1, length
!
rho_up = rho_in(ir,1)
@ -846,24 +821,11 @@ SUBROUTINE gcx_spin( length, rho_in, grho2_in, sx_tot, v1x_out, v2x_out )
v2x_dw = (1.0_DP - exx_fraction) * v2x_dw
ENDIF
!
#if !defined(_OPENACC)
CASE( 43 ) ! 'beefx' --- BEEF unavailable with OpenACC-- Will be enabled soon
!
rho_up = 2.0_DP * rho_up ; rho_dw = 2.0_DP * rho_dw
grho2_up = 4.0_DP * grho2_up ; grho2_dw = 4.0_DP * grho2_dw
! case igcx == 5 (HCTH) and 6 (OPTX) not implemented
! case igcx == 7 (meta-GGA) must be treated in a separate call to another
! routine: needs kinetic energy density in addition to rho and grad rho
!
CALL beefx(rho_up, grho2_up, sx_up, v1x_up, v2x_up, 0)
CALL beefx(rho_dw, grho2_dw, sx_dw, v1x_dw, v2x_dw, 0)
!
sx_tot(ir) = 0.5_DP * (sx_up*rnull_up + sx_dw*rnull_dw)
v2x_up = 2.0_DP * v2x_up
v2x_dw = 2.0_DP * v2x_dw
#endif
!
! case igcx == 5 (HCTH) and 6 (OPTX) not implemented
! case igcx == 7 (meta-GGA) must be treated in a separate call to another
! routine: needs kinetic energy density in addition to rho and grad rho
!
CASE DEFAULT
!
sx_tot(ir) = 0.0_DP
@ -878,12 +840,11 @@ SUBROUTINE gcx_spin( length, rho_in, grho2_in, sx_tot, v1x_out, v2x_out )
v2x_out(ir,2) = v2x_dw * rnull_dw
!
ENDDO
#if defined(_OPENMP) && !defined(_OPENACC)
!$omp end do
!$omp end parallel
#endif
#if defined(_OPENACC)
!$acc end data
#else
!$omp end do
!$omp end parallel
#endif
!
!
@ -933,13 +894,9 @@ SUBROUTINE gcc_spin( length, rho_in, zeta_io, grho_in, sc_out, v1c_out, v2c_out
#endif
!
#if defined(_OPENACC)
IF (igcc==14) CALL xclib_error( 'gcc_spin', 'BEEF not available with&
& OpenACC enabled', 1 )
!
!$acc data copyin(rho_in, grho_in), copyout(sc_out, v1c_out, v2c_out), copy(zeta_io)
!$acc parallel loop
#endif
#if defined(_OPENMP) && !defined(_OPENACC)
#else
!$omp parallel if(ntids==1) default(none) &
!$omp private( rho, zeta, grho, sc, v1c_up, v1c_dw, v2c ) &
!$omp shared( igcc, sc_out, v1c_out, v2c_out, &
@ -987,12 +944,6 @@ SUBROUTINE gcc_spin( length, rho_in, zeta_io, grho_in, sc_out, v1c_out, v2c_out
!
CALL pbec_spin( rho, zeta, grho, 2, sc, v1c_up, v1c_dw, v2c )
!
#if !defined(_OPENACC)
CASE( 14 ) !*****BEEF unavailable with OpenACC-- Will be enabled soon
!
call beeflocalcorrspin(rho, zeta, grho, sc, v1c_up, v1c_dw, v2c, 0)
#endif
!
CASE DEFAULT
!
sc = 0.0_DP
@ -1008,12 +959,11 @@ SUBROUTINE gcc_spin( length, rho_in, zeta_io, grho_in, sc_out, v1c_out, v2c_out
v2c_out(ir) = v2c
!
ENDDO
#if defined(_OPENMP) && !defined(_OPENACC)
!$omp end do
!$omp end parallel
#endif
#if defined(_OPENACC)
!$acc end data
#else
!$omp end do
!$omp end parallel
#endif
!
RETURN
@ -1076,8 +1026,7 @@ SUBROUTINE gcc_spin_more( length, rho_in, grho_in, grho_ud_in, &
#if defined(_OPENACC)
!$acc data copyin(rho_in, grho_in, grho_ud_in), copyout(sc, v1c, v2c, v2c_ud)
!$acc parallel loop
#endif
#if defined(_OPENMP) && !defined(_OPENACC)
#else
!$omp parallel if(ntids==1) default(none) &
!$omp private( rho_up, rho_dw, grho_up, grho_dw, grho_ud ) &
!$omp shared( length, rho_in, grho_in, grho_ud_in, &
@ -1139,17 +1088,228 @@ SUBROUTINE gcc_spin_more( length, rho_in, grho_in, grho_ud_in, &
END SELECT
!
ENDDO
#if defined(_OPENMP) && !defined(_OPENACC)
#if defined(_OPENACC)
!$acc end data
#else
!$omp end do
!$omp end parallel
#endif
#if defined(_OPENACC)
!$acc end data
#endif
!
RETURN
!
END SUBROUTINE gcc_spin_more
!
!
! ========> BEEF GGA DRIVERS <========================
!
!------------------------------------------------------------------------
SUBROUTINE gcxc_beef( length, rho_in, grho_in, sx_out, sc_out, v1x_out, &
v2x_out, v1c_out, v2c_out )
!---------------------------------------------------------------------
!! Driver for BEEF gga xc terms. Unpolarized case.
!
USE beef_interface, ONLY: beefx, beeflocalcorr
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: length
REAL(DP), INTENT(IN), DIMENSION(length) :: rho_in
REAL(DP), INTENT(IN), DIMENSION(length) :: grho_in
REAL(DP), INTENT(OUT), DIMENSION(length) :: sx_out, sc_out
REAL(DP), INTENT(OUT), DIMENSION(length) :: v1x_out, v2x_out
REAL(DP), INTENT(OUT), DIMENSION(length) :: v1c_out, v2c_out
!
! ... local variables
!
INTEGER :: ir
REAL(DP) :: rho, grho
REAL(DP) :: sx, v1x, v2x
REAL(DP) :: sc, v1c, v2c
!
IF (igcx == 43 .OR. igcc == 14) THEN
!
DO ir = 1, length
grho = grho_in(ir)
IF ( rho_in(ir) <= rho_threshold_gga .OR. grho <= grho_threshold_gga ) THEN
sx_out(ir) = 0.0_DP ; sc_out(ir) = 0.0_DP
v1x_out(ir) = 0.0_DP ; v1c_out(ir) = 0.0_DP
v2x_out(ir) = 0.0_DP ; v2c_out(ir) = 0.0_DP
CYCLE
ENDIF
!
rho = ABS(rho_in(ir))
!
IF ( igcx == 43 ) THEN
! last parameter = 0 means do not add LDA (=Slater) exchange
! (espresso) will add it itself
CALL beefx( rho, grho, sx, v1x, v2x, 0 )
sx_out(ir) = sx
v1x_out(ir) = v1x
v2x_out(ir) = v2x
ENDIF
!
IF ( igcc == 14 ) THEN
! last parameter 0 means: do not add lda contributions
! espresso will do that itself
CALL beeflocalcorr( rho, grho, sc, v1c, v2c, 0 )
sc_out(ir) = sc
v1c_out(ir) = v1c
v2c_out(ir) = v2c
ENDIF
ENDDO
!
ENDIF
!
RETURN
!
END SUBROUTINE gcxc_beef
!
!-----------------------------------------------------------------------------
SUBROUTINE gcx_spin_beef( length, rho_in, grho2_in, sx_tot, v1x_out, v2x_out )
!--------------------------------------------------------------------------
!! Driver for BEEF gga exchange term. Polarized case.
!
USE beef_interface, ONLY: beefx
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: length
REAL(DP), INTENT(IN), DIMENSION(length,2) :: rho_in
REAL(DP), INTENT(IN), DIMENSION(length,2) :: grho2_in
REAL(DP), INTENT(OUT), DIMENSION(length) :: sx_tot
REAL(DP), INTENT(OUT), DIMENSION(length,2) :: v1x_out
REAL(DP), INTENT(OUT), DIMENSION(length,2) :: v2x_out
!
! ... local variables
!
INTEGER :: ir, is, iflag
REAL(DP) :: rho_up, rho_dw, grho2_up, grho2_dw
REAL(DP) :: v1x_up, v1x_dw, v2x_up, v2x_dw
REAL(DP) :: sx_up, sx_dw, rnull_up, rnull_dw
REAL(DP) :: sxsr_up, sxsr_dw
REAL(DP) :: v1xsr_up, v1xsr_dw, v2xsr_up, v2xsr_dw
!
REAL(DP), PARAMETER :: small=1.D-10
REAL(DP), PARAMETER :: rho_trash=0.5_DP, grho2_trash=0.2_DP
! ... Temporary workaround for BEEF, which does not support GPU.
IF ( igcx == 43 ) THEN
!
DO ir = 1, length
!
rho_up = rho_in(ir,1)
rho_dw = rho_in(ir,2)
grho2_up = grho2_in(ir,1)
grho2_dw = grho2_in(ir,2)
rnull_up = 1.0_DP
rnull_dw = 1.0_DP
!
IF ( rho_up+rho_dw <= small ) THEN
sx_tot(ir) = 0.0_DP
v1x_out(ir,1) = 0.0_DP
v2x_out(ir,1) = 0.0_DP
v1x_out(ir,2) = 0.0_DP
v2x_out(ir,2) = 0.0_DP
CYCLE
ELSE
IF ( rho_up<=small .OR. SQRT(ABS(grho2_up))<=small ) THEN
rho_up = rho_trash
grho2_up = grho2_trash
rnull_up = 0.0_DP
ENDIF
IF ( rho_dw<=small .OR. SQRT(ABS(grho2_dw))<=small ) THEN
rho_dw = rho_trash
grho2_dw = grho2_trash
rnull_dw = 0.0_DP
ENDIF
ENDIF
!
rho_up = 2.0_DP * rho_up ; rho_dw = 2.0_DP * rho_dw
grho2_up = 4.0_DP * grho2_up ; grho2_dw = 4.0_DP * grho2_dw
!
CALL beefx(rho_up, grho2_up, sx_up, v1x_up, v2x_up, 0)
CALL beefx(rho_dw, grho2_dw, sx_dw, v1x_dw, v2x_dw, 0)
!
sx_tot(ir) = 0.5_DP * (sx_up*rnull_up + sx_dw*rnull_dw)
v2x_up = 2.0_DP * v2x_up
v2x_dw = 2.0_DP * v2x_dw
!
v1x_out(ir,1) = v1x_up * rnull_up
v1x_out(ir,2) = v1x_dw * rnull_dw
v2x_out(ir,1) = v2x_up * rnull_up
v2x_out(ir,2) = v2x_dw * rnull_dw
!
ENDDO
!
ENDIF
!
RETURN
!
END SUBROUTINE gcx_spin_beef
!
!
!--------------------------------------------------------------------------------
SUBROUTINE gcc_spin_beef( length, rho_in, zeta_io, grho_in, sc_out, v1c_out, v2c_out )
!-------------------------------------------------------------------------------
!! BEEF Gradient correction.
!
USE beef_interface, ONLY: beeflocalcorrspin
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: length
REAL(DP), INTENT(IN), DIMENSION(length) :: rho_in
REAL(DP), INTENT(INOUT), DIMENSION(length) :: zeta_io
REAL(DP), INTENT(IN), DIMENSION(length) :: grho_in
REAL(DP), INTENT(OUT), DIMENSION(length) :: sc_out
REAL(DP), INTENT(OUT), DIMENSION(length,2) :: v1c_out
REAL(DP), INTENT(OUT), DIMENSION(length) :: v2c_out
!
! ... local variables
!
INTEGER :: ir
REAL(DP) :: rho, zeta, grho
REAL(DP) :: sc, v1c_up, v1c_dw, v2c
!
IF ( igcc == 14 ) THEN
!
DO ir = 1, length
!
rho = rho_in(ir)
grho = grho_in(ir)
IF ( ABS(zeta_io(ir))<=1.0_DP ) zeta_io(ir) = SIGN( MIN(ABS(zeta_io(ir)), &
(1.0_DP-rho_threshold_gga)), zeta_io(ir) )
zeta = zeta_io(ir)
!
IF ( ABS(zeta)>1.0_DP .OR. rho<=rho_threshold_gga .OR. &
SQRT(ABS(grho))<=rho_threshold_gga ) THEN
sc_out(ir) = 0.0_DP
v1c_out(ir,1) = 0.0_DP ; v2c_out(ir) = 0.0_DP
v1c_out(ir,2) = 0.0_DP
CYCLE
ENDIF
rho = rho_in(ir)
grho = grho_in(ir)
zeta = zeta_io(ir)
!
IF ( ABS(zeta)>1.0_DP .OR. rho<=rho_threshold_gga .OR. &
SQRT(ABS(grho))<=rho_threshold_gga ) CYCLE
!
CALL beeflocalcorrspin( rho, zeta, grho, sc, v1c_up, v1c_dw, v2c, 0 )
!
sc_out(ir) = sc
v1c_out(ir,1) = v1c_up
v1c_out(ir,2) = v1c_dw
v2c_out(ir) = v2c
!
ENDDO
!
ENDIF
!
RETURN
!
END SUBROUTINE gcc_spin_beef
!
!
END MODULE qe_drivers_gga

View File

@ -86,7 +86,6 @@ SUBROUTINE xc_gcx_( length, ns, rho, grho, ex, ec, v1x, v2x, v1c, v2c, v2c_ud )
ec = 0.0_DP ; v1c = 0.0_DP ; v2c = 0.0_DP
IF ( PRESENT(v2c_ud) ) v2c_ud = 0.0_DP
!
!
#if defined(__LIBXC)
!
fkind_x = -1
@ -136,7 +135,10 @@ SUBROUTINE xc_gcx_( length, ns, rho, grho, ex, ec, v1x, v2x, v1c, v2c, v2c_ud )
!
IF ( ns==1 .AND. ANY(.NOT.is_libxc(3:4)) ) THEN
!
CALL gcxc( length, ABS(rho(:,1)), sigma, ex, ec, v1x(:,1), v2x(:,1), v1c(:,1), v2c(:,1) )
CALL gcxc( length, ABS(rho(:,1)), sigma, ex, ec, v1x(:,1), v2x(:,1), v1c(:,1), v2c(:,1) )
!
IF ( igcx==43 .OR. igcc==14 ) CALL gcxc_beef( length, ABS(rho(:,1)), grho2(:,1), ex, ec, &
v1x(:,1), v2x(:,1), v1c(:,1), v2c(:,1) )
!
DO k = 1, length
sgn(1) = SIGN(1._DP, rho(k,1))
@ -215,7 +217,8 @@ SUBROUTINE xc_gcx_( length, ns, rho, grho, ex, ec, v1x, v2x, v1c, v2c, v2c_ud )
( grho(2,:,1) + grho(2,:,2) )**2 + &
( grho(3,:,1) + grho(3,:,2) )**2
!
CALL gcc_spin( length, rh, zeta, grho2(:,1), ec, v1c, v2c(:,1) )
IF ( igcc/=14 ) CALL gcc_spin( length, rh, zeta, grho2(:,1), ec, v1c, v2c(:,1) )
IF ( igcc==14 ) CALL gcc_spin_beef( length, rh, zeta, grho2(:,1), ec, v1c, v2c(:,1) )
!
v2c(:,2) = v2c(:,1)
IF ( PRESENT(v2c_ud) ) v2c_ud(:) = v2c(:,1)
@ -263,7 +266,8 @@ SUBROUTINE xc_gcx_( length, ns, rho, grho, ex, ec, v1x, v2x, v1c, v2c, v2c_ud )
grho2(:,is) = grho(1,:,is)**2 + grho(2,:,is)**2 + grho(3,:,is)**2
ENDDO
!
CALL gcx_spin( length, rho, grho2, ex, v1x, v2x )
IF ( igcx/=43 ) CALL gcx_spin( length, rho, grho2, ex, v1x, v2x )
IF ( igcx==43 ) CALL gcx_spin_beef( length, rho, grho2, ex, v1x, v2x )
!
ENDIF
!
@ -293,6 +297,9 @@ SUBROUTINE xc_gcx_( length, ns, rho, grho, ex, ec, v1x, v2x, v1c, v2c, v2c_ud )
!
CALL gcxc( length, ABS(rho(:,1)), grho2(:,1), ex, ec, v1x(:,1), v2x(:,1), v1c(:,1), v2c(:,1) )
!
IF ( igcx==43 .OR. igcc==14 ) CALL gcxc_beef( length, ABS(rho(:,1)), grho2(:,1), ex, ec, &
v1x(:,1), v2x(:,1), v1c(:,1), v2c(:,1) )
!
DO k = 1, length
sgn(1) = SIGN(1._DP, rho(k,1))
ex(k) = ex(k) * sgn(1)
@ -305,7 +312,8 @@ SUBROUTINE xc_gcx_( length, ns, rho, grho, ex, ec, v1x, v2x, v1c, v2c, v2c_ud )
grho2(:,is) = grho(1,:,is)**2 + grho(2,:,is)**2 + grho(3,:,is)**2
ENDDO
!
CALL gcx_spin( length, rho, grho2, ex, v1x, v2x )
IF ( igcx/=43 ) CALL gcx_spin( length, rho, grho2, ex, v1x, v2x )
IF ( igcx==43 ) CALL gcx_spin_beef( length, rho, grho2, ex, v1x, v2x )
!
IF (igcc==3 .OR. igcc==7 .OR. igcc==13 ) THEN
!
@ -333,7 +341,9 @@ SUBROUTINE xc_gcx_( length, ns, rho, grho, ex, ec, v1x, v2x, v1c, v2c, v2c_ud )
( grho(2,:,1) + grho(2,:,2) )**2 + &
( grho(3,:,1) + grho(3,:,2) )**2
!
CALL gcc_spin( length, rh, zeta, grho2(:,1), ec, v1c, v2c(:,1) )
IF ( igcc/=14 ) CALL gcc_spin( length, rh, zeta, grho2(:,1), ec, v1c, v2c(:,1) )
IF ( igcc==14 ) CALL gcc_spin_beef( length, rh, zeta, grho2(:,1), ec, v1c, v2c(:,1) )
!
v2c(:,2) = v2c(:,1)
IF ( PRESENT(v2c_ud) ) v2c_ud(:) = v2c(:,1)

View File

@ -2,23 +2,27 @@
include ../make.inc
QEMODS = libupf.a
MODFLAGS= $(MOD_FLAG)../UtilXlib $(MOD_FLAG)../external/devxlib/src
# list of modules
# list of modules:
# OBJS_NODEP do not depend upon any external modules
# OBJS_DEP depend upon modules in UtilXlib/
# OBJS_GPU depend upon modules in devxlib
OBJS= \
atom.o \
atomic_number.o \
dqvan2.o \
dylmr2.o \
gth.o \
OBJS_DEP= \
init_us_0.o \
init_us_b0.o \
init_us_1.o \
init_tab_atwfc.o \
init_tab_beta.o \
init_tab_qrad.o \
init_tab_qrad.o
OBJS_NODEP= \
atom.o \
atomic_number.o \
dqvan2.o \
dylmr2.o \
gth.o \
init_us_2_base.o \
interp_atwfc.o \
paw_variables.o \
@ -56,7 +60,7 @@ xmltools.o \
ylmr2.o
# GPU versions of routines
OBJS += \
OBJS_GPU= \
dylmr2_gpu.o \
init_us_2_base_gpu.o \
interp_atwfc_gpu.o \
@ -67,17 +71,21 @@ OBJS += \
all : libupf.a virtual_v2.x upfconv.x casino2upf.x
libupf.a: $(MODULES) $(OBJS)
libupf.a: $(OBJS_DEP) $(OBJS_NODEP) $(OBJS_GPU)
$(AR) $(ARFLAGS) $@ $?
$(RANLIB) $@
virtual_v2.x : virtual_v2.o $(QEMODS)
$(LD) $(LDFLAGS) -o $@ virtual_v2.o $(QEMODS) $(LAPACK_LIBS) $(BLAS_LIBS)
# The following utilities are linked with $(OBJS_NODEP) instead of libupf.a
# because the latter has external dependencies that produce "missing symbols"
# errors for some linkers that try to resolve symbols even if never used
casino2upf.x : casino2upf.o casino_pp.o $(QEMODS)
$(LD) $(LDFLAGS) -o $@ casino2upf.o casino_pp.o $(QEMODS) $(LAPACK_LIBS) $(BLAS_LIBS)
upfconv.x : upfconv.o casino_pp.o $(QEMODS)
$(LD) $(LDFLAGS) -o $@ upfconv.o casino_pp.o $(QEMODS) $(LAPACK_LIBS) $(BLAS_LIBS)
virtual_v2.x : virtual_v2.o $(OBJS_NODEP)
$(LD) $(LDFLAGS) -o $@ virtual_v2.o $(OBJS_NODEP) $(LAPACK_LIBS) $(BLAS_LIBS)
casino2upf.x : casino2upf.o casino_pp.o $(OBJS_NODEP)
$(LD) $(LDFLAGS) -o $@ casino2upf.o casino_pp.o $(OBJS_NODEP) $(LAPACK_LIBS) $(BLAS_LIBS)
upfconv.x : upfconv.o casino_pp.o $(OBJS_NODEP)
$(LD) $(LDFLAGS) -o $@ upfconv.o casino_pp.o $(OBJS_NODEP) $(LAPACK_LIBS) $(BLAS_LIBS)
clean :
- /bin/rm -f *.o *.a *.d *.i *~ *_tmp.f90 *.mod *.L

View File

@ -19,7 +19,7 @@ http://www.quantum-espresso.org/pseudopotentials/unified-pseudopotential-format
The xml schema for the newer UPF definition can be found here:
http://www.quantum-espresso.org/ns/qes/qe_pp-1.0.xsd
In addition to the `libupf.a` library, two executable utilities are produced:
In addition to the `libupf.a` library, executable utilities are produced:
- `upfconv.x`, converting pseudopotentials in other formats into UPF:
see `upfconv.x -h` for more
@ -27,6 +27,8 @@ In addition to the `libupf.a` library, two executable utilities are produced:
- `virtual_v2.x`, courtesy Jingyang Wang (jw598@cornell.edu), generates
an averaged pseudopotential suitable for Virtual Crystal Approximation
- `casino2upf.x`, courtesy Mike Towler (see below)
A python script `fixfile.py` is also present, to remove undesired `&`
characters from UPF files that hinder their parsing by xml tools.

View File

@ -40,7 +40,7 @@ CONTAINS
!! ierr=0 : xml schema, ierr=-2: UPF v.2
!! ierr=-81: error reading PP file
!
iun = xml_openfile ( filename )
iun = xml_open_file ( filename )
IF ( iun == -1 ) CALL upf_error('read_upf', 'cannot open file',1)
call xmlr_opentag ( 'qe_pp:pseudo', IERR = ierr )
if ( ierr == 0 ) then

View File

@ -523,7 +523,7 @@ subroutine read_pseudo_pswfc (upf, iunps)
character (len=75) :: dummy
integer :: nb, ir
ALLOCATE( upf%chi( upf%mesh, MAX( upf%nwfc, 1 ) ) )
ALLOCATE( upf%chi( upf%mesh, upf%nwfc ) )
upf%chi = 0.0_DP
do nb = 1, upf%nwfc
read (iunps, *, err=100, end=100) dummy !Wavefunction labels

View File

@ -54,7 +54,7 @@ CONTAINS
v2 = .true.
END SELECT
!
iun = xml_openfile ( filename )
iun = xml_open_file ( filename )
IF ( iun == -1 ) CALL upf_error('write_upf', 'cannot open file',1)
!
! The starting line, header and info sections differ a lot

View File

@ -60,7 +60,7 @@ MODULE xmltools
!
PRIVATE
! general subroutines
PUBLIC :: xml_openfile, xml_closefile
PUBLIC :: xml_open_file, xml_closefile
! subroutines for writing
PUBLIC :: add_attr
PUBLIC :: xmlw_writetag, xmlw_opentag, xmlw_closetag
@ -262,7 +262,7 @@ CONTAINS
!
END SUBROUTINE add_c_attr
!
FUNCTION xml_openfile ( filexml ) RESULT (iun)
FUNCTION xml_open_file ( filexml ) RESULT (iun)
!
! returns on output the opened unit number if opened successfully
! returns -1 otherwise
@ -289,7 +289,7 @@ CONTAINS
print "('file ',a,' opened with unit ',i5)",trim(filexml),iun
#endif
!
END FUNCTION xml_openfile
END FUNCTION xml_open_file
!
SUBROUTINE xml_closefile ( )
!