Commit 47c5de62 authored by giannozz's avatar giannozz

More USPP_related variables moved to Modules/uspp.f90

Note that lqx => lmaxq for consistency with other names
(those ending in x are static dimensioning)
Beware unexpected side effects (PG)


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@933 c92efa57-630b-4861-b058-cf58834340f0
parent 598adf55
......@@ -17,6 +17,7 @@ subroutine d0rhod2v (ipert, drhoscf)
USE io_global, ONLY : stdout
USE io_files, ONLY : iunigk
USE kinds, only : DP
USE uspp_param, ONLY: nh
use pwcom
USE wavefunctions_module, ONLY : evc
use phcom
......
......@@ -20,7 +20,7 @@ subroutine d3_summary
USE io_global, ONLY : stdout
USE kinds, only : DP
use pwcom
USE uspp_param, ONLY : rinner, nqlc, nqf, lll, nbeta, psd, iver
USE uspp_param, ONLY : rinner, nqlc, nqf, lll, nbeta, psd, iver, tvanp
USE atom, ONLY: mesh, numeric, xmin, dx, nlcc
USE control_flags, ONLY : iverbosity
use phcom
......
......@@ -16,6 +16,7 @@ subroutine d3vrho
#include "machine.h"
USE kinds, only : DP
use pwcom
USE uspp_param, ONLY: nh
USE wavefunctions_module, ONLY: evc
USE io_files, ONLY : iunigk
use phcom
......
......@@ -17,6 +17,7 @@ subroutine dqrhod2v (ipert, drhoscf)
#include "machine.h"
USE kinds, only : DP
use pwcom
USE uspp_param, ONLY: nh
USE wavefunctions_module, ONLY: evc
USE io_files, ONLY : iunigk
use phcom
......
......@@ -17,6 +17,7 @@ subroutine dvdpsi (nu_i, xq_, dvloc, vkb_, vkbq_, psi_, dvpsi_)
use pwcom
use phcom
use d3com
USE uspp_param, ONLY: nh
!
implicit none
integer :: nu_i
......
......@@ -19,6 +19,7 @@ subroutine incdrhoscf2 (drhoscf, weight, ik, dbecsum, mode, flag)
USE kinds, only : DP
use pwcom
USE wavefunctions_module, ONLY: evc
USE uspp_param, ONLY: nhm
use phcom
implicit none
......
......@@ -12,7 +12,7 @@ subroutine A_h(e,h,ah)
#include "machine.h"
USE kinds, only: DP
USE cell_base,ONLY : alat, omega, tpiba2
USE us, ONLY : vkb, nkb
USE uspp, ONLY : vkb, nkb
USE lsda_mod, ONLY : current_spin, nspin
USE wvfct, ONLY: nbnd, npwx, npw, g2kin, igk
USE wavefunctions_module, ONLY: evc, psic
......
......@@ -16,6 +16,7 @@ subroutine dvpsi_e(kpoint,ipol)
#include "machine.h"
USE kinds, only: DP
use pwcom
USE uspp_param, ONLY: nh
USE wavefunctions_module, ONLY: evc
USE rbecmod, ONLY: becp
use cgcom
......
......@@ -14,6 +14,7 @@ subroutine dvpsi_kb(kpoint,nu)
#include "machine.h"
USE kinds, only: DP
use pwcom
USE uspp_param, ONLY: nh, nhm
USE atom, ONLY: nlcc, mesh, dx, r, rab, rho_atc, numeric
USE wavefunctions_module, ONLY: evc, psic
use cgcom
......
......@@ -14,7 +14,7 @@ subroutine H_h(e,h,Ah)
USE kinds, only: DP
USE wvfct, ONLY: nbnd, npwx, npw, g2kin, igk
USE gvect, ONLY : gstart
USE us, ONLY : vkb, nkb
USE uspp, ONLY : vkb, nkb
USE lsda_mod, ONLY : current_spin
USE scf, ONLY : vrs
use rbecmod, only: becp
......
......@@ -15,6 +15,7 @@ subroutine rhod2vkb(dyn0)
#include "machine.h"
use pwcom
USE wavefunctions_module, ONLY: evc, psic
USE uspp_param, only: nh
use cgcom
!
implicit none
......
......@@ -1067,17 +1067,14 @@ MODULE input_parameters
! CARDS parameters
!=----------------------------------------------------------------------------=!
!
! Note: See file read_cards.f90 for the cards syntax and usage
! Note: See file read_cards.f90 for card syntax and usage
!
! ATOMIC_SPECIES
!
CHARACTER(LEN=4) :: atom_label(nsx) = 'XX' ! label of the atomic species being read
CHARACTER(LEN=80) :: atom_pfile(nsx) = 'YY' ! pseudopotential file name
REAL(dbl) :: atom_mass(nsx) = 0.0d0 ! atomic mass
INTEGER :: atom_ptyp(nsx) = 0 ! pseudopotential type
! unsorted atomic masses
REAL(dbl) :: atom_mass(nsx) = 0.0d0 ! atomic mass of the i-th atomic species
! in atomic mass units: 1 a.m.u. = 1822.9 a.u. = 1.6605 * 10^-27 kg
! atomic mass of the i-th atomic species
LOGICAL :: taspc = .FALSE.
!
......@@ -1090,7 +1087,7 @@ MODULE input_parameters
INTEGER :: na_inp(nsx) = 0 ! number of atom for each specie
LOGICAL :: tapos = .FALSE.
CHARACTER(LEN=80) :: atomic_positions = 'crystal'
! atomic_positions = 'bohr' | 'armstrong' | 'crystal' | 'alat'
! atomic_positions = 'bohr' | 'angstrong' | 'crystal' | 'alat'
! select the units for the atomic positions being read from stdin
!
......
......@@ -327,11 +327,7 @@ MODULE read_cards_module
DO is = 1, ntyp
!
CALL read_line( input_line )
IF ( prog == 'CP' ) THEN
READ( input_line, * ) lb_pos, atom_mass(is), psfile, atom_ptyp(is)
ELSE
READ( input_line, * ) lb_pos, atom_mass(is), psfile
END IF
READ( input_line, * ) lb_pos, atom_mass(is), psfile
atom_pfile(is) = TRIM( psfile )
lb_pos = ADJUSTL( lb_pos )
atom_label(is) = TRIM( lb_pos )
......
......@@ -27,6 +27,8 @@ MODULE uspp_param
rinner(lqmax,npsx) ! values of r_L
INTEGER :: &
nbeta(npsx), &! number of beta functions
nh(npsx), &! number of beta functions per atomic type
nhm, &! max number of different beta functions per atom
kkbeta(npsx), &! point where the beta are zero
nqf(npsx), &! number of coefficients for Q
nqlc(npsx), &! number of angular momenta in Q
......@@ -35,21 +37,25 @@ MODULE uspp_param
iver(3,npsx) ! version of the atomic code
INTEGER :: &
lmaxkb, &! max angular momentum
lqx ! max angular momentum + 1 for Q functions
lmaxq ! max angular momentum + 1 for Q functions
LOGICAL :: &
tvanp(npsx), &! if .TRUE. the atom is of Vanderbilt type
newpseudo(npsx) ! if .TRUE. multiple projectors are allowed
END MODULE uspp_param
!
!
MODULE uspp
!
! Ultrasoft PPs:
! - Clebsch-Gordan coefficients "ap", auxiliary variables "lpx", "lpl"
! - beta and q functions of the solid
!
USE kinds, ONLY: DP
USE parameters, ONLY: lmaxx, lqmax
IMPLICIT NONE
PRIVATE
SAVE
PUBLIC :: nlx, lpx, lpl, ap, aainit
PUBLIC :: nlx, lpx, lpl, ap, aainit, indv, nhtol, nhtolm, nkb, vkb, dvan, &
deeq, qq, becsum, nhtoj, deallocate_uspp, beta
INTEGER, PARAMETER :: &
nlx = (lmaxx+1)**2, &! maximum number of combined angular momentum
mx = 2*lqmax-1 ! maximum magnetic angular momentum of Q
......@@ -60,6 +66,25 @@ MODULE uspp
ap(lqmax*lqmax,nlx,nlx)
! Clebsch-Gordan coefficients for spherical harmonics
!
INTEGER :: nkb ! total number of beta functions, with struct.fact.
!
INTEGER, ALLOCATABLE ::&
indv(:,:), &! indes linking atomic beta's to beta's in the solid
nhtol(:,:), &! correspondence n <-> angular momentum l
nhtolm(:,:) ! correspondence n <-> combined lm index for (l,m)
!
COMPLEX(KIND=DP), ALLOCATABLE, TARGET :: &
vkb(:,:), &! all beta functions in reciprocal space
dvan(:,:,:,:), &! the D functions of the solid
deeq(:,:,:,:) ! the integral of V_eff and Q_{nm}
REAL(KIND=DP), ALLOCATABLE :: &
qq(:,:,:), &! the q functions in the solid
becsum(:,:,:), &! the sum of bec = <beta|psi>
nhtoj(:,:) ! correspondence n <-> total angular momentum
!
REAL(KIND=DP), ALLOCATABLE :: &
beta(:,:,:) ! beta functions for CP (without struct.factor)
!
CONTAINS
!
!-----------------------------------------------------------------------
......@@ -213,5 +238,15 @@ CONTAINS
return
end function compute_ap
subroutine deallocate_uspp()
IF( ALLOCATED( nhtol ) ) DEALLOCATE( nhtol )
IF( ALLOCATED( indv ) ) DEALLOCATE( indv )
IF( ALLOCATED( nhtolm ) ) DEALLOCATE( nhtolm )
IF( ALLOCATED( nhtoj ) ) DEALLOCATE( nhtoj )
IF( ALLOCATED( vkb ) ) DEALLOCATE( vkb )
IF( ALLOCATED( qq ) ) DEALLOCATE( qq )
IF( ALLOCATED( dvan ) ) DEALLOCATE( dvan )
end subroutine deallocate_uspp
end module uspp
......@@ -10,6 +10,7 @@ subroutine add_dkmds(kpoint, uact, jpol, dvkb)
use pwcom
USE kinds, only : DP
USE wavefunctions_module, ONLY : evc
USE uspp_param, only: nh
use phcom
implicit none
......
......@@ -8,6 +8,7 @@ subroutine add_for_charges (ik, uact)
use pwcom
USE kinds, only : DP
USE uspp_param, only: nh
use phcom
implicit none
!
......
subroutine adddvepsi_us(becp2,ipol,kpoint)
! This subdoutine adds to dvpsi the terms which depend on the augmentation
! charge. It assume that the variable dpqq, has been set.
!
subroutine adddvepsi_us(becp2,ipol,kpoint)
! This subdoutine adds to dvpsi the terms which depend on the augmentation
! charge. It assume that the variable dpqq, has been set.
!
#include "machine.h"
use pwcom
USE kinds, only : DP
use phcom
implicit none
integer, intent(in) :: ipol, kpoint
complex(kind=dp), intent(in) :: becp2(nkb,nbnd)
real(kind=dp) :: fact
complex(kind=dp), allocatable :: ps(:)
integer:: ijkb0, nt, na, ih, jh, ikb, jkb, ibnd
allocate (ps(nbnd))
ijkb0 = 0
do nt = 1, ntyp
do na = 1, nat
if (ityp(na).eq.nt) then
do ih = 1, nh (nt)
ikb = ijkb0 + ih
ps = (0.d0,0.d0)
do jh = 1, nh (nt)
jkb = ijkb0 + jh
fact=at(1,ipol)*dpqq(ih,jh,1,nt)+ &
at(2,ipol)*dpqq(ih,jh,2,nt)+ &
at(3,ipol)*dpqq(ih,jh,3,nt)
do ibnd=1, nbnd_occ(kpoint)
ps(ibnd) = ps(ibnd) &
+ becp2(jkb,ibnd)*(0.d0,1.d0)*qq(ih,jh,nt)+ &
becp1(jkb,ibnd,kpoint)*fact
enddo
enddo
do ibnd = 1, nbnd_occ (kpoint)
call ZAXPY(npw,ps(ibnd),vkb(1,ikb),1,dvpsi(1,ibnd),1)
enddo
enddo
ijkb0=ijkb0+nh(nt)
endif
enddo
enddo
if (jkb.ne.nkb) call errore ('adddvepsi_us', 'unexpected error', 1)
deallocate(ps)
return
end
use pwcom
USE kinds, only : DP
USE uspp_param, only: nh
use phcom
implicit none
integer, intent(in) :: ipol, kpoint
complex(kind=dp), intent(in) :: becp2(nkb,nbnd)
real(kind=dp) :: fact
complex(kind=dp), allocatable :: ps(:)
integer:: ijkb0, nt, na, ih, jh, ikb, jkb, ibnd
allocate (ps(nbnd))
ijkb0 = 0
do nt = 1, ntyp
do na = 1, nat
if (ityp(na).eq.nt) then
do ih = 1, nh (nt)
ikb = ijkb0 + ih
ps = (0.d0,0.d0)
do jh = 1, nh (nt)
jkb = ijkb0 + jh
fact=at(1,ipol)*dpqq(ih,jh,1,nt)+ &
at(2,ipol)*dpqq(ih,jh,2,nt)+ &
at(3,ipol)*dpqq(ih,jh,3,nt)
do ibnd=1, nbnd_occ(kpoint)
ps(ibnd) = ps(ibnd) &
+ becp2(jkb,ibnd)*(0.d0,1.d0)*qq(ih,jh,nt)+ &
becp1(jkb,ibnd,kpoint)*fact
enddo
enddo
do ibnd = 1, nbnd_occ (kpoint)
call ZAXPY(npw,ps(ibnd),vkb(1,ikb),1,dvpsi(1,ibnd),1)
enddo
enddo
ijkb0=ijkb0+nh(nt)
endif
enddo
enddo
if (jkb.ne.nkb) call errore ('adddvepsi_us', 'unexpected error', 1)
deallocate(ps)
return
end subroutine adddvepsi_us
......@@ -17,8 +17,10 @@ subroutine adddvscf (ipert, ik)
#include "machine.h"
USE kinds, ONLY : DP
USE uspp_param, ONLY : nh, tvanp
USE uspp, ONLY : vkb
! modules from pwcom
USE us, ONLY : okvan, tvanp, nh, vkb
USE us, ONLY : okvan
USE lsda_mod, ONLY : lsda, current_spin, isk
USE basis, ONLY : ntyp, nat, ityp
USE wvfct, ONLY : nbnd
......
......@@ -18,6 +18,7 @@ subroutine addusdbec (ik, wgt, psi, dbecsum)
use pwcom
USE kinds, only : DP
USE uspp_param, only: nh, tvanp, nhm
use phcom
implicit none
!
......
......@@ -21,7 +21,8 @@ subroutine addusddens (drhoscf, dbecsum, irr, mode0, npe, iflag)
USE wavefunctions_module, ONLY: psic
use phcom
USE kinds, only : DP
USE uspp_param, ONLY: lqx
USE uspp_param, ONLY: lmaxq, nh, nhm, tvanp
implicit none
!
! the dummy variables
......@@ -72,7 +73,7 @@ subroutine addusddens (drhoscf, dbecsum, irr, mode0, npe, iflag)
call start_clock ('addusddens')
allocate (aux( ngm , nspin , npertx))
allocate (sk ( ngm))
allocate (ylmk0(ngm , lqx * lqx))
allocate (ylmk0(ngm , lmaxq * lmaxq))
allocate (qgm( ngm))
allocate (qmod( ngm))
if (.not.lgamma) allocate (qpg( 3 , ngm))
......@@ -82,12 +83,12 @@ subroutine addusddens (drhoscf, dbecsum, irr, mode0, npe, iflag)
!
if (.not.lgamma) then
call setqmod (ngm, xq, g, qmod, qpg)
call ylmr2 (lqx * lqx, ngm, qpg, qmod, ylmk0)
call ylmr2 (lmaxq * lmaxq, ngm, qpg, qmod, ylmk0)
do ig = 1, ngm
qmod (ig) = sqrt (qmod (ig) )
enddo
else
call ylmr2 (lqx * lqx, ngm, g, gg, ylmk0)
call ylmr2 (lmaxq * lmaxq, ngm, g, gg, ylmk0)
do ig = 1, ngm
qmod (ig) = sqrt (gg (ig) )
enddo
......
......@@ -20,7 +20,7 @@ subroutine addusddense (drhoscf, dbecsum)
use pwcom
use phcom
USE kinds, only : DP
USE uspp_param, ONLY: lqx
USE uspp_param, ONLY: lmaxq, nh, nhm, tvanp
implicit none
!
! the dummy variables
......@@ -55,13 +55,13 @@ subroutine addusddense (drhoscf, dbecsum)
allocate (aux( ngm, nspin, 3))
allocate (sk ( ngm))
allocate (qg ( nrxx))
allocate (ylmk0(ngm , lqx * lqx))
allocate (ylmk0(ngm , lmaxq * lmaxq))
allocate (qgm (ngm))
allocate (qmod (ngm))
!
! And then we compute the additional charge in reciprocal space
!
call ylmr2 (lqx * lqx, ngm, g, gg, ylmk0)
call ylmr2 (lmaxq * lmaxq, ngm, g, gg, ylmk0)
do ig = 1, ngm
qmod (ig) = sqrt (gg (ig) )
enddo
......
......@@ -20,6 +20,7 @@ subroutine addusdynmat (dynwrk)
use pwcom
USE kinds, only : DP
USE uspp_param, only: tvanp, nh
use phcom
implicit none
......
......@@ -16,7 +16,7 @@ subroutine addusldos (ldos, becsum1)
#include "machine.h"
use pwcom
USE wavefunctions_module, ONLY: psic
USE uspp_param, ONLY: lqx
USE uspp_param, ONLY: lmaxq, tvanp, nh, nhm
implicit none
complex(kind=DP) :: ldos (nrxx, nspin)
! local density of states
......@@ -38,12 +38,12 @@ subroutine addusldos (ldos, becsum1)
! work space
allocate (aux ( ngm , nspin))
allocate (ylmk0(ngm , lqx * lqx))
allocate (ylmk0(ngm , lmaxq * lmaxq))
allocate (qgm ( ngm))
allocate (qmod( ngm))
aux (:,:) = (0.d0,0.d0)
call ylmr2 (lqx * lqx, ngm, g, gg, ylmk0)
call ylmr2 (lmaxq * lmaxq, ngm, g, gg, ylmk0)
do ig = 1, ngm
qmod (ig) = sqrt (gg (ig) )
enddo
......
......@@ -21,6 +21,7 @@ subroutine allocate_phq
use phcom
use el_phon
USE becmod, ONLY: becp
USE uspp_param, ONLY: nhm
implicit none
!
! allocate space for the quantities needed in the phonon program
......
......@@ -19,6 +19,7 @@ subroutine compute_alphasum
#include "machine.h"
use pwcom
USE kinds, only : DP
USE uspp_param, ONLY: nh, tvanp
use phcom
implicit none
......
......@@ -19,6 +19,7 @@ subroutine compute_becsum
#include "machine.h"
use pwcom
USE kinds, only : DP
USE uspp_param, ONLY: nh, tvanp
use phcom
implicit none
......
......@@ -20,6 +20,7 @@ subroutine compute_drhous (drhous, dbecsum, wgg, becq, alpq)
USE wavefunctions_module, ONLY: evc
USE io_files, ONLY: iunigk
USE kinds, only : DP
USE uspp_param, ONLY: nhm
use phcom
implicit none
!
......
......@@ -19,6 +19,7 @@ subroutine compute_nldyn (wdyn, wgg, becq, alpq)
use pwcom
USE kinds, only : DP
USE uspp_param, ONLY: nh
use phcom
implicit none
......
......@@ -7,10 +7,9 @@ subroutine compute_qdipol
USE constants, ONLY: fpi
USE atom, ONLY: r, rab
USE ions_base, ONLY: ntyp => nsp
USE us, only: nh, nhtol, nhtolm, indv, tvanp
USE uspp, only: nhtol, nhtolm, indv, nlx, ap
USE uspp_param, only: nbrx, nbeta, lll, kkbeta, qfunc, rinner, &
qfcoef, nqf
USE uspp, ONLY: nlx, ap
qfcoef, nqf, tvanp, nh
USE phus, ONLY: dpqq
implicit none
......
......@@ -20,6 +20,7 @@ subroutine drho
use pwcom
USE kinds, only : DP
USE uspp_param, only: nhm
use phcom
implicit none
......
......@@ -17,6 +17,7 @@ subroutine drhodvnl (ik, ikk, nper, nu_i0, wdyn, dbecq, dalpq)
#include "machine.h"
use pwcom
USE kinds, only : DP
USE uspp_param, only: nh
use phcom
implicit none
integer :: ik, ikk, nper, nu_i0
......
......@@ -21,7 +21,7 @@ subroutine dvanqq
use pwcom
USE kinds, only : DP
use phcom
USE uspp_param, ONLY: lqx
USE uspp_param, ONLY: lmaxq, nh, tvanp
implicit none
!
! And the local variables
......@@ -62,10 +62,10 @@ subroutine dvanqq
allocate (aux5( ngm ,nat, 3 ))
allocate (qmodg( ngm))
allocate (veff ( nrxx , nspin))
allocate (ylmk0( ngm , lqx * lqx))
allocate (ylmk0( ngm , lmaxq * lmaxq))
allocate (qgm ( ngm))
if (.not.lgamma) then
allocate (ylmkq(ngm , lqx * lqx))
allocate (ylmkq(ngm , lmaxq * lmaxq))
allocate (qpg (3, ngm))
allocate (qmod( ngm))
allocate (qgmq( ngm))
......@@ -75,13 +75,13 @@ subroutine dvanqq
!
! compute spherical harmonics
!
call ylmr2 (lqx * lqx, ngm, g, gg, ylmk0)
call ylmr2 (lmaxq * lmaxq, ngm, g, gg, ylmk0)
do ig = 1, ngm
qmodg (ig) = sqrt (gg (ig) )
enddo
if (.not.lgamma) then
call setqmod (ngm, xq, g, qmod, qpg)
call ylmr2 (lqx * lqx, ngm, qpg, qmod, ylmkq)
call ylmr2 (lmaxq * lmaxq, ngm, qpg, qmod, ylmkq)
do ig = 1, ngm
qmod (ig) = sqrt (qmod (ig) )
enddo
......
......@@ -10,6 +10,7 @@ subroutine dvkb3(kpoint,dvkb)
use pwcom
USE kinds, only : DP
USE wavefunctions_module, ONLY : evc
USE uspp_param, ONLY: nh
use phcom
implicit none
......
......@@ -23,6 +23,7 @@ subroutine dvpsi_e (kpoint, ipol)
USE wavefunctions_module, ONLY: evc
USE kinds, only : DP
USE becmod, ONLY: becp
USE uspp_param, ONLY: nh
use phcom
implicit none
!
......
......@@ -19,6 +19,7 @@ subroutine dvqpsi_us_only (ik, mode, uact)
use pwcom
USE kinds, only : DP
USE uspp_param, ONLY: nh
use phcom
implicit none
!
......
......@@ -19,6 +19,7 @@ subroutine dynmat_us
USE wavefunctions_module, ONLY: evc
USE io_files, ONLY: iunigk
USE kinds, only : DP
USE uspp_param, ONLY: nh
use phcom
#ifdef __PARA
use para
......
......@@ -19,6 +19,7 @@ subroutine incdrhoscf (drhoscf, weight, ik, dbecsum, mode)
use pwcom
USE wavefunctions_module, ONLY: evc
USE kinds, only : DP
USE uspp_param, ONLY: nhm
use phcom
implicit none
......
......@@ -18,6 +18,7 @@ subroutine incdrhous (drhoscf, weight, ik, dbecsum, evcr, wgg, becq, &
use pwcom
USE kinds, only : DP
USE uspp_param, ONLY: nhm, nh
use phcom
implicit none
......
......@@ -21,6 +21,7 @@ subroutine localdos (ldos, ldoss, dos_ef)
use pwcom
USE wavefunctions_module, ONLY: evc, psic
USE kinds, only : DP
USE uspp_param, ONLY: nh, nhm, tvanp
use phcom
USE io_files, ONLY: iunigk
implicit none
......
......@@ -19,7 +19,7 @@ subroutine newdq (dvscf, npe)
use pwcom
USE kinds, only : DP
use phcom
USE uspp_param, ONLY: lqx
USE uspp_param, ONLY: nh, nhm, tvanp, lmaxq
implicit none
!
! The dummy variables
......@@ -54,7 +54,7 @@ subroutine newdq (dvscf, npe)
allocate (aux1 (ngm))
allocate (aux2 (ngm , nspin))
allocate (veff (nrxx))
allocate (ylmk0(ngm , lqx * lqx))
allocate (ylmk0(ngm , lmaxq * lmaxq))
allocate (qgm (ngm))
allocate (qmod (ngm))
......@@ -64,12 +64,12 @@ subroutine newdq (dvscf, npe)
!
if (.not.lgamma) then
call setqmod (ngm, xq, g, qmod, qg)
call ylmr2 (lqx * lqx, ngm, qg, qmod, ylmk0)
call ylmr2 (lmaxq * lmaxq, ngm, qg, qmod, ylmk0)
do ig = 1, ngm
qmod (ig) = sqrt (qmod (ig) )