scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

band_interpolation input with namelists and cards

This commit is contained in:
Ivan Carnimeo 2022-01-20 16:41:44 +01:00
parent 73235b943c
commit fd8699db75
8 changed files with 259 additions and 178 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Modules/read_cards.f90
View File

@ -26,7 +26,7 @@ MODULE read_cards_module
!
PRIVATE
!
PUBLIC :: read_cards
PUBLIC :: read_cards, card_kpoints
!
! ... end of module-scope declarations
!

2
PP/examples/exx_interpolated_bands_example/input_fourier
View File

@ -1,7 +1,7 @@
Method
fourier
User-given star functions
1
1 .true.
0.5 0.5 0.5
Interpolation parameters
NMax 7

2
PP/examples/exx_interpolated_bands_example/input_fourier-diff
View File

@ -1,7 +1,7 @@
Method
fourier-diff
User-given star functions
1
1 .true.
0.5 0.5 0.5
Interpolation parameters
NMax 7

122
PP/src/band_interpolation.f90
View File

@ -11,31 +11,135 @@
!
!----------------------------------------------------------------------------
program band_interpolation
use globalmod
use idwmod
use fouriermod
use fourierdiffmod
use globalmod, ONLY : card_user_stars, card_roughness, read_xml_input, t, build_kpath, at, bg, &
ek, eq, Nb, Nq, NSym, q, Op, print_bands, deallocate_global
use idwmod, ONLY : idw, p_metric, scale_sphere
use fouriermod, ONLY : fourier, miller_max, check_periodicity, RoughC, RoughN, NUser
use fourierdiffmod, ONLY : fourierdiff
use input_parameters, ONLY : k_points , xk, nkstot
use parser, ONLY : read_line
USE mp_global, ONLY : mp_startup
USE io_global, ONLY : stdout
USE read_cards_module, ONLY : card_kpoints
implicit none
!
Call read_input ()
integer, parameter :: iunit = 5
integer :: ios
integer :: i, iii, ik
CHARACTER(len=256) :: input_line
LOGICAL :: tend
CHARACTER(len=80) :: card
CHARACTER(len=1), EXTERNAL :: capital
!
CHARACTER(len=80) :: method = ' '
!! a string describing the method used for interpolation
!
NAMELIST / interpolation / method, miller_max, check_periodicity, p_metric, scale_sphere
!
!
write(*,*) 'PROGRAM: band_interpolation '
!
#if defined(__MPI)
CALL mp_startup ( )
#endif
!
! Set defaults
!
method = 'fourier-diff'
check_periodicity = .false.
p_metric = 2
scale_sphere = 4.0d0
miller_max = 6
RoughN = 1
allocate( RoughC(RoughN) )
RoughC(1) = 1.0d0
NUser = 0
k_points = 'none'
!
! Read input file
!
ios = 0
READ( iunit, interpolation, iostat = ios )
!
100 CALL read_line( input_line, end_of_file=tend )
!
IF( tend ) GOTO 120
IF( input_line == ' ' .OR. input_line(1:1) == '#' .OR. &
input_line(1:1) == '!' ) GOTO 100
!
READ (input_line, *) card
!
DO i = 1, len_trim( input_line )
input_line( i : i ) = capital( input_line( i : i ) )
ENDDO
!
IF ( trim(card) == 'ROUGHNESS' ) THEN
!
CALL card_roughness( input_line )
!
ELSEIF ( trim(card) == 'USER_STARS' ) THEN
!
CALL card_user_stars( input_line )
!
ELSEIF ( trim(card) == 'K_POINTS' ) THEN
!
CALL card_kpoints( input_line )
!
ELSE
!
WRITE( *,'(A)') 'Warning: card '//trim(input_line)//' ignored'
!
ENDIF
!
GOTO 100
!
120 CONTINUE
!
if(k_points.ne.'tpiba') then
write(stdout,'(A,A)') 'k_points = ', k_points
Call errore('band_interpolation ' , ' K_POINTS card must be specified with tpiba_b ', 1)
end if
!
if(nkstot.le.0) then
write(stdout,'(A,I5)') 'nkstot = ', nkstot
Call errore('band_interpolation ' , ' wrong number of k-points ', 1)
end if
!
Write(stdout,'(A,A)') 'Interpolation method: ', method
if( TRIM(method).ne.'idw'.and.TRIM(method).ne.'idw-sphere'&
.and.TRIM(method).ne.'fourier'.and.TRIM(method).ne.'fourier-diff' ) &
Call errore('band_interpolation', 'Wrong interpolation method ', 1)
!
! Build the abscissa values for band plotting
!
allocate( t(nkstot) )
t = 0.0d0
do ik= 2, nkstot
t(ik) = t(ik-1) + sqrt(dot_product(xk(:,ik)-xk(:,ik-1),xk(:,ik)-xk(:,ik-1)))
write(*,'(I5,f12.6)') ik, t(ik)
end do
!
! Read xml file
!
Call read_xml_input ()
!
ek = 0.0d0
!
if(TRIM(method).eq.'idw') then
!
Call idw (1, Nb, Nq, q, eq, Nk, k, ek, at, bg)
Call idw (1, Nb, Nq, q, eq, nkstot, xk, ek, at, bg)
!
elseif(TRIM(method).eq.'idw-sphere') then
!
Call idw (2, Nb, Nq, q, eq, Nk, k, ek, at, bg)
Call idw (2, Nb, Nq, q, eq, nkstot, xk, ek, at, bg)
!
elseif(TRIM(method).eq.'fourier') then
!
Call fourier (Nb, Nq, q, eq, Nk, k, ek, Nsym, at, bg, Op)
Call fourier (Nb, Nq, q, eq, nkstot, xk, ek, Nsym, at, bg, Op)
!
elseif(TRIM(method).eq.'fourier-diff') then
!
Call fourierdiff (Nb, Nq, q, eq, Nk, k, ek, Nsym, at, bg, Op)
Call fourierdiff (Nb, Nq, q, eq, nkstot, xk, ek, Nsym, at, bg, Op)
!
else
!

9
PP/src/fourierdiffmod.f90
View File

@ -57,6 +57,7 @@ implicit none
!
write(*,'(A)') '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,'(A)') 'Fourier difference interpolation method'
if(check_periodicity) write(*,*) 'Checking Star functions periodicity (WARNING: time consuming)'
!
Na = Nq - 1 ! dimension of the linear system
!
@ -71,8 +72,10 @@ implicit none
write(*,'(A)') 'Creating Star functions...'
Call find_stars(NSym, Op, at, .true.)
!
write(*,*) 'Checking Star functions periodicity...'
Call check_stars(Nq, q, NSym, Op, bg)
if(check_periodicity) then
write(*,*) 'Checking Star functions periodicity...'
Call check_stars(Nq, q, NSym, Op, bg)
end if
!
! fStarsOnQ = [S_m(q_i)-S_m(q_Nq)] / sqrt(rho_m)
write(*,*) 'Computing fStarsOnQ...'
@ -129,7 +132,7 @@ implicit none
ek(:,:) = dble(ek_c(:,:))
!
deallocate( matX, matB, matA, matC, matC1, ek_c )
deallocate( C )
deallocate( RoughC )
deallocate( VecStars )
deallocate( fStarsOnQ )
deallocate( fStarsOnK )

47
PP/src/fouriermod.f90
View File

@ -17,10 +17,12 @@ save
real(dp), parameter :: eps = 0.000010d0, Zero = 0.0d0, One = 1.0d0, Two = 2.0d0, Four = 4.0d0
real(dp), parameter :: Pi = Four*atan(One)
!
logical :: check_periodicity = .false.
!
! the largest Miller index used to generate all the lattice vectors inside an outer shell
integer :: NMax
integer :: NC
real(dp), allocatable :: C(:)
integer :: miller_max
integer :: RoughN
real(dp), allocatable :: RoughC(:)
!
integer :: NStars ! total number of Star functions
real(dp), allocatable :: VecStars(:,:) ! symmetry inequivalent lattice vectors generating the Star functions (one per Star)
@ -31,7 +33,7 @@ CONTAINS
!----------------------------------------------------------------------------
subroutine allocate_fourier( )
implicit none
allocate ( C(NC) )
allocate ( RoughC(RoughN) )
return
end subroutine
!----------------------------------------------------------------------------
@ -64,12 +66,15 @@ implicit none
real(dp), allocatable :: rmatJ(:,:)
!
write(*,'(A)') '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,'(A)') 'Fourier interpolation method'
write(*,'(A)') 'Fourier interpolation method'
if(check_periodicity) write(*,*) 'Checking Star functions periodicity (WARNING: time consuming)'
write(*,*) 'Creating Star functions...'
Call find_stars(NSym, Op, at)
!
write(*,*) 'Checking Star functions periodicity...'
Call check_stars(Nq, q, NSym, Op, bg)
if(check_periodicity) then
write(*,*) 'Checking Star functions periodicity...'
Call check_stars(Nq, q, NSym, Op, bg)
end if
!
! fStarsOnQ = S_m(q) / sqrt(rho_m)
write(*,*) 'Computing fStarsOnQ...'
@ -133,7 +138,7 @@ implicit none
!
ek(:,:) = dble(ek_c(:,:))
!
deallocate( C )
deallocate( RoughC )
deallocate( ek_c, eq_c )
deallocate( VecStars )
deallocate( matQQ )
@ -156,7 +161,7 @@ implicit none
! local variables
logical :: Skip000
!
! arrays containing the lattice vectors inside the outer shell defined by NMax
! arrays containing the lattice vectors inside the outer shell defined by miller_max
integer :: NAll, NPrint
real(dp), allocatable :: VecAll(:,:) ! all lattice vectors
real(dp), allocatable :: VecInq(:,:) ! symmetry inequivalent lattice vectors
@ -174,7 +179,7 @@ implicit none
Skip000 = .false.
endif
!
NAll = (2 * NMax + 1 )**3 ! from -Nmax to Nmax is (2*Nmax + 1), for 3 space directions
NAll = (2 * miller_max + 1 )**3 ! from -miller_max to miller_max is (2*miller_max + 1), for 3 space directions
if(Skip000) then
NAll = NAll - 1 ! remove the (0, 0, 0) lattice vector
write(*,*) 'Skipping the (0,0,0) lattice vector...'
@ -184,17 +189,17 @@ implicit none
!
if(NUser.gt.0) then
NAll = NAll + NUser
write(*,*) "Creating ", NAll, " vectors from ", NMax, " indexes and ", NUser, " user-given vectors"
write(*,*) "Creating ", NAll, " vectors from ", miller_max, " indexes and ", NUser, " user-given vectors"
else
write(*,*) "Creating ", NAll, " vectors from ", NMax, " indexes"
write(*,*) "Creating ", NAll, " vectors from ", miller_max, " indexes"
end if
!
allocate ( VecAll(3,NAll), ModAll(NAll), MapAll(NAll) )
!
ivec = 0
do ii = -NMax, NMax
do jj= -NMax, NMax
do kk= -NMax, NMax
do ii = -miller_max, miller_max
do jj= -miller_max, miller_max
do kk= -miller_max, miller_max
if(Skip000.and.(ii.eq.0.and.jj.eq.0.and.kk.eq.0)) cycle
ivec = ivec + 1
VecAll(:,ivec) = dble(ii)*at(:,1) + dble(jj)*at(:,2) + dble(kk)*at(:,3)
@ -214,8 +219,8 @@ implicit none
end if
!
if(ivec.ne.NAll) then
write(*,*) "ERROR: wrong number of lattice vectors for a given NMax"
write(*,*) "NMax= ",NMax," ivec=",ivec," NAll=",NAll, " NUser=",NUser
write(*,*) "ERROR: wrong number of lattice vectors for a given miller_max"
write(*,*) "miller_max= ",miller_max," ivec=",ivec," NAll=",NAll, " NUser=",NUser
stop
endif
!
@ -483,12 +488,12 @@ implicit none
real(dp) :: rmod, rho
integer :: ic, iexp
!
rho = C(1)
if(NC.gt.1) then
rho = RoughC(1)
if(RoughN.gt.1) then
rmod = sqrt(dot_product(vec,vec))
do ic = 2, NC
do ic = 2, RoughN
iexp = 2*(ic - 1 )
rho = rho + C(ic) * rmod**iexp
rho = rho + RoughC(ic) * rmod**iexp
end do
end if
sqrt_rho = sqrt(rho)

243
PP/src/globalmod.f90
View File

@ -13,10 +13,7 @@
MODULE globalmod
USE kinds, ONLY : dp
implicit none
!INTEGER, PARAMETER :: DP = selected_real_kind(14,200)
!
! method
character(len=50) :: method
!
! band indexes
integer :: Nb
@ -42,30 +39,112 @@ implicit none
real(dp), allocatable :: Op_tmp(:,:,:) ! this is just a buffer to convert Op in cartesian units.
! quite redundant here, but useful to use s_axis_to_cart without modifications
!
logical :: trough = .false.
logical :: tuser = .false.
!
CONTAINS
!----------------------------------------------------------------------------
subroutine read_input ()
subroutine card_user_stars( input_line )
use fouriermod , only : NUser, VecUser
use parser, ONLY : read_line
implicit none
integer :: ivec
LOGICAL :: tend,terr
CHARACTER(len=256) :: input_line
!
IF ( tuser ) THEN
CALL errore( ' card_user_stars ', ' two occurrences', 2 )
ENDIF
!
! ... input Star vectors
!
CALL read_line( input_line, end_of_file = tend, error = terr )
IF (tend) GOTO 10
IF (terr) GOTO 20
READ(input_line, *, END=10, ERR=20) NUser
!
if(NUser.gt.0) then
allocate( VecUser(3,NUser) )
!
do ivec = 1, NUser
CALL read_line( input_line, end_of_file = tend, error = terr )
IF (tend) GOTO 10
IF (terr) GOTO 20
READ(input_line,*, END=10, ERR=20) VecUser(1:3,ivec)
end do
end if
!
tuser = .true.
!
RETURN
!
10 CALL errore ('card_user_stars', ' end of file while reading roughness function ', 1)
20 CALL errore ('card_user_stars', ' error while reading roughness function', 1)
!
end subroutine card_user_stars
!----------------------------------------------------------------------------
subroutine card_roughness( input_line )
use fouriermod , only : RoughN, RoughC
use parser, ONLY : read_line
implicit none
LOGICAL :: tend,terr
CHARACTER(len=256) :: input_line
!
IF ( trough ) THEN
CALL errore( ' card_roughness ', ' two occurrences', 2 )
ENDIF
!
! ... input coefficients for the roughness function
!
CALL read_line( input_line, end_of_file = tend, error = terr )
IF (tend) GOTO 10
IF (terr) GOTO 20
READ(input_line, *, END=10, ERR=20) RoughN
!
if(RoughN.gt.1) then
deallocate( RoughC )
allocate( RoughC(RoughN) )
end if
!
if(RoughN.gt.0) then
CALL read_line( input_line, end_of_file = tend, error = terr )
IF (tend) GOTO 10
IF (terr) GOTO 20
READ(input_line,*, END=10, ERR=20) RoughC(1:RoughN)
!
else
Call errore( ' card_roughness ', ' RoughN must be a positive integer ', 2 )
end if
!
trough = .true.
!
RETURN
!
10 CALL errore ('card_roughness', ' end of file while reading roughness function ', 1)
20 CALL errore ('card_roughness', ' error while reading roughness function', 1)
!
end subroutine card_roughness
!----------------------------------------------------------------------------
subroutine read_xml_input ()
!
! read the input file and make all allocations
! read the xml input file and make all allocations
!
use fouriermod, ONLY : NMax, allocate_fourier, NC, C, NUser, VecUser
use idwmod, ONLY : PMetric, ScaleSphere
use qes_read_module, only: qes_read
use qes_types_module, only: band_structure_type, atomic_structure_type, symmetries_type, basis_set_type
USE io_global, ONLY : stdout
use qes_read_module, ONLY : qes_read
use qes_types_module, ONLY : band_structure_type, atomic_structure_type, symmetries_type, basis_set_type
use fox_dom
implicit none
integer :: ik, iq, ib, ikl, ivec
integer :: ilines
integer :: iq
integer :: isym
character(len=50) :: string
!
type (node),pointer :: doc, pnode1, pnode2
type (band_structure_type) :: bandstr
type (node),pointer :: doc, pnode1, pnode2
type (band_structure_type) :: bandstr
type (atomic_structure_type) :: atstr
type (symmetries_type) :: symstr
type (basis_set_type) :: basisstr
!
! read data from the xml_file
!
doc => parsefile('pwscf.xml')
pnode1 => item(getElementsByTagname(doc, 'output'),0)
pnode2 => item(getElementsByTagname(pnode1, 'atomic_structure'),0)
@ -78,124 +157,14 @@ implicit none
call qes_read(pnode2, basisstr)
call destroy (doc)
nullify (doc, pnode1, pnode2)
read(*, *)
read(*, *) method
write(*,*) 'Interpolation method: ', method
if( TRIM(method).ne.'idw'.and.TRIM(method).ne.'idw-sphere'&
.and.TRIM(method).ne.'fourier'.and.TRIM(method).ne.'fourier-diff' ) then
write(*,*) 'ERROR: Wrong method ', method
stop
end if
!
! read specific parameters for the interpolation methods
!
if( TRIM(method).eq.'idw') THEN
!
! read parameters for IDW interpolation
!
PMetric = -1
ScaleSphere = -1.0d0
read(*,*)
read(*,*) PMetric
write(*,*) 'PMetric: ', PMetric
if(PMetric.lt.0) THEN
write(*,*) "ERROR: Wrong input for idw method"
stop
end if
!
elseif( TRIM(method).eq.'idw-sphere' ) THEN
!
! read parameters for IDW interpolation inside a sphere
!
PMetric = -1
ScaleSphere = -1.0d0
read(*,*)
read(*,*) PMetric, ScaleSphere
write(*,*) 'PMetric: ', PMetric, 'ScaleSphere ', ScaleSphere
if(PMetric.lt.0.or.ScaleSphere.lt.0.0d0) THEN
write(*,*) "ERROR: Wrong input for IDW-Sphere method"
stop
end if
!
elseif( TRIM(method).eq.'fourier'.or.TRIM(method).eq.'fourier-diff' ) THEN
!
! optionally add user-given star functions to the basis set
!
NUser = -1
read(*,*)
read(*,*) NUser
if(NUser.gt.0) then
write(*,*) NUser, ' user-given star functions found'
allocate( VecUser(3,NUser) )
VecUser = 0.0d0
do ivec = 1, NUser
read(*,*) VecUser(1:3,ivec)
write(*,'(3f12.6)') VecUser(1:3,ivec)
end do
elseif(NUser.eq.0) then
write(*,*) 'No user-given star functions provided'
else
write(*,*) 'ERROR: Wrong NUser'
write(*,*) ' Please provide non-negative NUser'
stop
end if
!
! read parameters for Fourier interpolation
!
NMax = -1
NC = -1
read(*, *)
read(*, *) string, NMax
if( NMax.le.0) then
write(*,*) 'Wrong NMax: ', NMax
write(*,*) 'NMax must be greater than 0 '
stop
end if
read(*, *) string, NC
if( NC.le.0) then
write(*,*) 'Wrong NC: ', NC
write(*,*) 'NC must be greater than 0 '
stop
end if
Call allocate_fourier( )
read(*, *) string, C(1:NC)
write(*,*) NC, ' coefficients read for rho expansion: ', C(:)
end if
!
! read the list of Nkl special points
!
read(*, *)
read(*, *) Nkl
!
write(*,*) Nkl, ' special points read'
!
! create the abscissa values for bands plotting
!
allocate( kl(3,Nkl), kln(Nkl ) )
!
do ikl = 1, Nkl
read(*, *) kl(:,ikl), kln(ikl)
write(*,'(3f12.6,I5)') kl(:,ikl), kln(ikl)
end do
!
Nlines = Nkl - 1
Nk = sum(kln(1:Nlines)) + 1
write(*,*) 'Creating ', Nlines, ' lines connecting ', Nkl, ' special points with ', Nk, ' k-points'
!
allocate( k(3,Nk), t(Nk) )
!
Call build_kpath()
!
deallocate( kl, kln )
!
! read the uniform grid of q-points from xml
!
Nq = bandstr%nks
Nb = bandstr%nbnd
!
write(*,*) Nq, ' points on the uniform grid, ', Nb, ' bands'
!write(*,'(A)') 'iq, q(iq, :), e(iq, :)'
write(stdout,'(2(I5,A))') Nq, ' points on the uniform grid, ', Nb, ' bands'
!write(stdout,'(A)') 'iq, q(iq, :), e(iq, :)'
!
allocate( q(3, Nq), eq(Nq, Nb), ek(Nk,Nb) )
!
@ -204,7 +173,7 @@ implicit none
end do
do iq = 1, Nq
eq(iq,:) = bandstr%ks_energies(iq)%eigenvalues%vector(:)*27.211386245988
!write(*,'(I5,11f12.6)') iq, q(iq, :), eq(iq, :)
!write(stdout, '(I5,11f12.6)') iq, q(iq, :), eq(iq, :)
end do
!
! read from xml crystalline group specifications
@ -213,19 +182,19 @@ implicit none
at(1:3,1) = atstr%cell%a1 / atstr%alat
at(1:3,2) = atstr%cell%a2 / atstr%alat
at(1:3,3) = atstr%cell%a3 / atstr%alat
write(*,*) ' Crystal lattice vectors found '
write(*,*) 'Ra: ' , at(:,1)
write(*,*) 'Rb: ' , at(:,2)
write(*,*) 'Rc: ' , at(:,3)
write(stdout,'(A)') ' Crystal lattice vectors found '
write(stdout,'(A,3f12.6)') 'Ra: ' , at(:,1)
write(stdout,'(A,3f12.6)') 'Rb: ' , at(:,2)
write(stdout,'(A,3f12.6)') 'Rc: ' , at(:,3)
bg(1:3,1) = basisstr%reciprocal_lattice%b1
bg(1:3,2) = basisstr%reciprocal_lattice%b2
bg(1:3,3) = basisstr%reciprocal_lattice%b3
write(*,*) ' Reciprocal lattice vectors found '
write(*,*) 'Ga: ' , bg(:,1)
write(*,*) 'Gb: ' , bg(:,2)
write(*,*) 'Gc: ' , bg(:,3)
write(stdout,'(A)') ' Reciprocal lattice vectors found '
write(stdout,'(A,3f12.6)') 'Ga: ' , bg(:,1)
write(stdout,'(A,3f12.6)') 'Gb: ' , bg(:,2)
write(stdout,'(A,3f12.6)') 'Gc: ' , bg(:,3)
Nsym = symstr%nsym
write(*,*) Nsym, ' symmetry operations found '
write(stdout,'(I5,A)') Nsym, ' symmetry operations found '
!
allocate( Op(3,3,Nsym), Op_tmp(3,3,Nsym) )
!
@ -237,7 +206,7 @@ implicit none
!
return
!
end subroutine read_input
end subroutine read_xml_input
!----------------------------------------------------------------------------
subroutine build_kpath ()
!

10
PP/src/idwmod.f90
View File

@ -22,9 +22,9 @@ save
!
integer, parameter :: dp = selected_real_kind(14,200)
!
integer :: PMetric ! metric for the (inverse) distance
integer :: p_metric ! metric for the (inverse) distance
!
real(dp) :: ScaleSphere ! scaling factor for the radius of the sphere for the modified method
real(dp) :: scale_sphere ! scaling factor for the radius of the sphere for the modified method
!
CONTAINS
!----------------------------------------------------------------------------
@ -71,7 +71,7 @@ implicit none
!
end do
end do
R = ScaleSphere * Rmin
R = scale_sphere * Rmin
write(*,*) 'Sphere radius: ', Rmin, ' Scaled sphere radius: ', R
end if
!
@ -98,14 +98,14 @@ implicit none
NCount(1) = NCount(1) + 1
if(iwhat.eq.1) then
! basic idw method
w = 1.0d0/(d**PMetric)
w = 1.0d0/(d**p_metric)
elseif(iwhat.eq.2) then
! search only inside the sphere R (idw-sphere)
!w = (max(0.0d0, (R-d))/(R*d))**2
w = 0.0d0
if(d.lt.R) then
NCount(2) = NCount(2) + 1
w = 1.0d0/(d**PMetric) !((R-d)/(R*d))**2
w = 1.0d0/(d**p_metric) !((R-d)/(R*d))**2
end if
end if
dsum = dsum + w