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

688 lines
21 KiB
Fortran
Raw Permalink Blame History

!
! Copyright (C) 2020-2023 Quantum ESPRESSO Foundation
! This file is distributed under the terms of the
! GNU General Public License. See the file 'License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
MODULE cpmd_module
!
! Read and convert a pseudopotential written in the CPMD format, TYPE:
! NORMCONSERVING [ NUMERIC, LOGARITHMIC, CAR ], single radial grid
! to unified pseudopotential format (v.2)
!
USE upf_kinds, ONLY: dp
USE upf_io, ONLY: stdout
!
IMPLICIT NONE
PRIVATE
PUBLIC :: read_cpmd
!
CHARACTER (len=80) title
!
! All variables read from CPMD file format
!
INTEGER :: ixc, pstype = 0
REAL(dp) :: alphaxc
REAL(dp) :: z, zv
! Grid variables
INTEGER :: mesh
REAL(dp) :: amesh, rmax, xmin
REAL(dp), ALLOCATABLE :: r(:)
! PP variables
INTEGER, PARAMETER :: lmaxx=3
INTEGER ::lmax, nwfc=0
! Car PP variables
REAL(dp) :: alphaloc, alpha(0:lmaxx), a(0:lmaxx), b(0:lmaxx)
! Goedecker PP variables
INTEGER, PARAMETER :: ncmax=4, nlmax=3
INTEGER :: nc, nl(0:lmaxx)
REAL(dp) :: rc, rl(0:lmaxx), c(ncmax), h(0:lmaxx, nlmax*(nlmax+1)/2 )
! Numeric PP variables
REAL(dp), ALLOCATABLE :: vnl(:,:)
REAL(dp), ALLOCATABLE :: chi(:,:)
! Core correction variables
LOGICAL :: nlcc=.false.
REAL(dp), ALLOCATABLE :: rho_atc(:)
! Variables used for reading and for checks
INTEGER :: maxinfo, info_lines
PARAMETER (maxinfo = 100)
CHARACTER (len=80), ALLOCATABLE :: info_sect(:)
!------------------------------
!
CONTAINS
!
! ----------------------------------------------------------
SUBROUTINE read_cpmd(iunps, upf, ierr)
! ----------------------------------------------------------
!
USE pseudo_types, ONLY : pseudo_upf
USE upf_utils, ONLY : matches
!
IMPLICIT NONE
INTEGER, INTENT(in) :: iunps
INTEGER, INTENT(out) :: ierr
TYPE (pseudo_upf) :: upf
!
INTEGER :: found = 0, closed = 0, unknown = 0
INTEGER :: i, l, dum, ios
CHARACTER (len=256) line
REAL (dp) :: amesh_, vnl0(0:3)
LOGICAL :: grid_read = .false., wfc_read=.false.
!
ierr = 1
!
info_lines = 0
10 READ (iunps,'(A)',end=20,err=20) line
IF (matches ("&ATOM", trim(line)) ) THEN
found = found + 1
! Z
READ (iunps,'(a)',end=200,err=200) line
l = len_trim(line)
i = INDEX(line,'=')
READ (line(i+1:l),*) z
! Zv
READ (iunps,'(a)',end=200,err=200) line
l = len_trim(line)
i = INDEX(line,'=')
READ (line(i+1:l),*) zv
! XC
READ (iunps,'(a)',end=200,err=200) line
l = len_trim(line)
i = INDEX(line,'=')
READ (line(i+1:l),*) ixc, alphaxc
! TYPE
READ (iunps,'(a)',end=200,err=200) line
IF ( matches("NORMCONSERVING",line) ) THEN
IF ( matches("NUMERIC",line) .or. matches("LOGARITHMIC",line) ) THEN
pstype = 1
ELSEIF ( matches("CAR",line) ) THEN
pstype = 2
ELSEIF ( matches("GOEDECKER",line) ) THEN
pstype = 3
ENDIF
ENDIF
IF (pstype == 0 ) then
write(stdout,'(5x,"read_cpmd: unknown type ",a)') trim(line)
return
END IF
ELSEIF (matches ("&INFO", trim(line)) ) THEN
found = found + 1
! read (iunps,'(a)') title
! store info section for later perusal
ALLOCATE (info_sect(maxinfo))
DO i=1,maxinfo
READ (iunps,'(a)',end=20,err=20) title
IF (matches ("&END", trim(title)) ) THEN
closed = closed + 1
GOTO 10
ELSE
info_sect(i) = trim(title)
info_lines = i
ENDIF
ENDDO
ELSEIF (matches ("&POTENTIAL", trim(line)) ) THEN
found = found + 1
READ (iunps,'(a)') line
IF ( pstype == 1 ) THEN
!
! NORMCONSERVING NUMERIC
!
READ (line,*,iostat=ios) mesh, amesh_
IF ( ios /= 0) THEN
READ (line,*,iostat=ios) mesh
amesh_ = -1.0d0
ENDIF
IF ( .not. grid_read ) ALLOCATE (r(mesh))
!
! determine the number of angular momenta
!
READ (iunps, '(a)') line
ios = 1
lmax= 4
DO WHILE (ios /= 0)
lmax = lmax - 1
READ(line,*,iostat=ios) r(1),(vnl0(l),l=0,lmax)
ENDDO
ALLOCATE (vnl(mesh,0:lmax))
vnl(1,0:lmax) = vnl0(0:lmax)
DO i=2,mesh
READ(iunps, *) r(i),(vnl(i,l),l=0,lmax)
ENDDO
IF ( .not.grid_read ) THEN
CALL check_radial_grid ( amesh_, mesh, r, amesh, ios )
if (ios /= 0 ) return
grid_read = .true.
ENDIF
ELSEIF ( pstype == 2 ) THEN
!
! NORMCONSERVING CAR
!
READ(iunps, *) alphaloc
! convert r_c's written in file to alpha's: alpha = 1/r_c^2
alphaloc = 1.d0/alphaloc**2
DO lmax=-1,2
READ(iunps, '(A)') line
IF (matches ("&END", trim(line)) ) THEN
closed = closed + 1
exit
ENDIF
READ(line, *) alpha(lmax+1), a(lmax+1), b(lmax+1)
alpha(lmax+1) = 1.d0/alpha(lmax+1)**2
ENDDO
ELSEIF ( pstype == 3 ) THEN
!
! NORMCONSERVING GOEDECKER
!
c(:) = 0.d0
rl(:) = 0.d0
nl(:) = 0
h(:,:) = 0.d0
READ(iunps, *) lmax
lmax = lmax - 1
IF ( lmax > lmaxx ) THEN
write(stdout,'(5x,"read_cpmd: incorrect lmax read")')
return
END IF
READ(iunps, *) rc
READ(iunps, '(A)') line
READ(line, *) nc
IF ( nc > ncmax ) THEN
! I am not sure if it is possible to use nc in the same line
! where it is read. Just in case, better to read twice
write(stdout,'(5x,"read_cpmd: incorrect nc read")')
return
END IF
READ(line, *) dum, (c(i), i=1,nc)
DO l=0,lmax+1
READ(iunps, '(A)') line
IF ( matches ("&END", trim(line)) ) THEN
closed = closed + 1
exit
ENDIF
READ(line, *) rl(l), nl(l)
IF ( nl(l) > nlmax ) THEN
write(stdout,'(5x,"read_cpmd: incorrect nl read")')
return
END IF
IF ( nl(l) > 0 ) &
READ(line, *) rl(l), dum, ( h(l,i), i=1,nl(l)*(nl(l)+1)/2)
ENDDO
ENDIF
ELSEIF (matches ("&WAVEFUNCTION", trim(line)) ) THEN
wfc_read=.true.
found = found + 1
READ (iunps,'(A)') line
READ (line,*,iostat=ios) mesh, amesh_
IF ( ios /= 0) THEN
READ (line,*,iostat=ios) mesh
amesh_ = -1.0d0
ENDIF
IF ( .not. grid_read ) ALLOCATE(r(mesh))
! find number of atomic wavefunctions
READ (iunps,'(A)') line
DO nwfc = lmax+1,1,-1
READ(line,*,iostat=ios) r(1),(vnl0(l),l=0,nwfc-1)
IF ( ios == 0 ) exit
ENDDO
IF ( ios /= 0 ) THEN
write(stdout,'(5x,"read_cpmd: at least one atomic wvfct should be present")')
return
END IF
ALLOCATE(chi(mesh,nwfc))
chi(1,1:nwfc) = vnl0(0:nwfc-1)
DO i=2,mesh
READ(iunps, *) r(i),(chi(i,l),l=1,nwfc)
ENDDO
IF ( .not.grid_read ) THEN
CALL check_radial_grid ( amesh_, mesh, r, amesh, ios )
IF ( ios /= 0 ) return
grid_read = .true.
ENDIF
ELSEIF (matches ("&NLCC", trim(line)) ) THEN
found = found + 1
READ (iunps, '(a)') line
READ(iunps, *) mesh
nlcc = ( mesh > 0 )
IF (nlcc) THEN
IF ( .not. matches ("NUMERIC", trim(line)) ) THEN
write(stdout,'(5x,"read_cpmd: only NUMERIC core correction supported")')
return
END IF
ALLOCATE (rho_atc(mesh))
READ(iunps, * ) (r(i), rho_atc(i), i=1,mesh)
ENDIF
ELSEIF (matches ("&ATDENS", trim(line)) ) THEN
! skip over &ATDENS section, add others here, if there are more.
DO WHILE(.not. matches("&END", trim(line)))
READ (iunps,'(a)') line
ENDDO
ELSEIF (matches ("&END", trim(line)) ) THEN
closed = closed + 1
ELSE
WRITE(stdout, '(5x,"read_cpmd: line ",a," ignored")') trim(line)
unknown = unknown + 1
ENDIF
GOTO 10
20 CONTINUE
IF ( pstype /= 3 ) THEN
IF (nlcc .and. found /= 5 .or. .not.nlcc .and. found /= 4) THEN
write(stdout,'(5x,"read_cpmd: some &FIELD card missing")')
return
end if
ELSE
IF (found /= 3) THEN
write(stdout,'(5x,"read_cpmd: some &FIELD card missing")')
return
end if
ENDIF
IF (closed /= found) THEN
write(stdout,'(5x,"read_cpmd: some &END card missing")')
return
end if
IF (unknown /= 0 ) WRITE(stdout,'("WARNING: ",i3," cards not read")') unknown
!
IF ( .not. grid_read ) THEN
xmin = -7.0d0
amesh=0.0125d0
rmax =100.0d0
WRITE(stdout,'("A radial grid must be provided. We use the following one:")')
WRITE(stdout,'("r_i = e^{xmin+(i-1)*dx}/Z, i=1,mesh, with parameters:")')
WRITE(stdout,'("Z=",f6.2,", xmin=",f6.2," dx=",f8.4," rmax=",f6.1)') &
z, xmin, amesh, rmax
mesh = 1 + (log(z*rmax)-xmin)/amesh
mesh = (mesh/2)*2+1 ! mesh is odd (for historical reasons?)
ALLOCATE (r(mesh))
DO i=1, mesh
r(i) = exp (xmin+(i-1)*amesh)/z
ENDDO
WRITE(stdout,'(I4," grid points, rmax=",f8.4)') mesh, r(mesh)
grid_read = .true.
ENDIF
rmax = r(mesh)
xmin = log(z*r(1))
!
IF ( .not. wfc_read ) Write(stdout,'("Notice: atomic wfcs not found")')
!
IF ( pstype == 2 ) THEN
ALLOCATE (vnl(mesh,0:lmax))
DO l=0, lmax
DO i=1, mesh
vnl(i,l) = ( a(l) + b(l)*r(i)**2 ) * exp (-alpha(l)*r(i)**2) - &
zv * erf (sqrt(alphaloc)*r(i))/r(i)
ENDDO
ENDDO
ENDIF
ierr = 0
CALL convert_cpmd(upf)
RETURN
200 WRITE(stdout, '("read_cpmd: error in reading file")')
END SUBROUTINE read_cpmd
! ----------------------------------------------------------
SUBROUTINE convert_cpmd(upf)
! ----------------------------------------------------------
!
USE pseudo_types, ONLY : pseudo_upf
USE upf_const, ONLY : e2
USE upf_utils, ONLY : spdf_to_l, l_to_spdf
!
IMPLICIT NONE
!
TYPE(pseudo_upf) :: upf
!
REAL(dp), ALLOCATABLE :: aux(:)
REAL(dp) :: x, x2, vll, rcloc, fac
CHARACTER (len=2):: label
CHARACTER (len=2), EXTERNAL :: atom_name
INTEGER :: lloc, my_lmax, l, i, j, ij, ir, iv, jv
!
! NOTE: many CPMD pseudopotentials created with the 'Hamann' code
! from Juerg Hutter's homepage have additional (bogus) entries for
! pseudo-potential and wavefunction. In the 'report' they have
! the same rc and energy eigenvalue than the previous angular momentum.
! we need to be able to ignore that part or the resulting UPF file
! will be useless. so we first print the info section and ask
! for the LMAX to really use. AK 2005/03/30.
!
DO i=1,info_lines
WRITE(stdout,'(A)') info_sect(i)
ENDDO
IF ( pstype == 3 ) THEN
! not actually used, except by write_upf to write a meaningful message
lloc = -3
rcloc=0.0
ELSE
WRITE(*,'("max L to use ( <= ",I1," ) > ")', advance="NO") lmax
READ (5,*) my_lmax
IF ((my_lmax <= lmax) .and. (my_lmax >= 0)) lmax = my_lmax
WRITE(*,'("local L ( <= ",I1," ), Rc for local pot (au) > ")', advance="NO") lmax
READ (5,*) lloc, rcloc
ENDIF
!
IF ( pstype == 3 ) THEN
upf%generated= "Generated in analytical, separable form"
upf%author = "Goedecker/Hartwigsen/Hutter/Teter"
upf%date = "Phys.Rev.B58, 3641 (1998); B54, 1703 (1996)"
ELSE
upf%generated= "Generated using unknown code"
upf%author = "unknown"
upf%date = "unknown"
ENDIF
upf%nv = "2.0.1"
upf%comment = "Info: automatically converted from CPMD format"
upf%psd = atom_name ( nint(z) )
! reasonable assumption
IF (z > 18) THEN
upf%rel = 'no'
ELSE
upf%rel = 'scalar'
ENDIF
IF ( pstype == 3 ) THEN
upf%typ = 'NC'
ELSE
upf%typ = 'SL'
ENDIF
upf%tvanp = .false.
upf%tpawp = .false.
upf%tcoulombp=.false.
upf%nlcc = nlcc
!
IF (ixc==900) THEN
WRITE(stdout,'("Pade approx. not implemented! assuming Perdew-Zunger LDA")')
upf%dft='SLA-PZ-NOGX-NOGC'
ELSEIF (ixc==1100) THEN
upf%dft='SLA-PZ-NOGX-NOGC'
ELSEIF (ixc==1111) THEN
upf%dft='SLA-PZ-B86-P88'
ELSEIF (ixc==1134 .or. ixc==1434) THEN
upf%dft='SLA-PW-PBX-PBC'
ELSEIF (ixc==1134) THEN
upf%dft='revPBE'
ELSEIF (ixc==1197) THEN
upf%dft='PBESOL'
ELSEIF (ixc==1312) THEN
upf%dft='BLYP'
ELSEIF (ixc==362) THEN
upf%dft='OLYP'
ELSEIF (ixc==1372) THEN
upf%dft='XLYP'
ELSEIF (ixc==55) THEN
upf%dft='HCTH'
ELSE
WRITE(*,'("Unknown DFT ixc=",i4,". Please provide a DFT name > ")', advance="NO") ixc
READ *, upf%dft
ENDIF
WRITE(stdout,'("Assuming DFT: ",A," . Please check this is what you want")') &
trim(upf%dft)
!
upf%zp = zv
upf%etotps =0.0d0
upf%ecutrho=0.0d0
upf%ecutwfc=0.0d0
IF ( lmax == lloc) THEN
upf%lmax = lmax-1
ELSE
upf%lmax = lmax
ENDIF
upf%lloc = lloc
upf%lmax_rho = 0
upf%nwfc = nwfc
!
ALLOCATE( upf%els(upf%nwfc) )
ALLOCATE( upf%oc(upf%nwfc) )
ALLOCATE( upf%epseu(upf%nwfc) )
ALLOCATE( upf%lchi(upf%nwfc) )
ALLOCATE( upf%nchi(upf%nwfc) )
ALLOCATE( upf%rcut_chi (upf%nwfc) )
ALLOCATE( upf%rcutus_chi(upf%nwfc) )
DO i=1, upf%nwfc
10 WRITE(*,'("Wavefunction # ",i1,": label (e.g. 4s), occupancy > ")', advance="NO") i
READ (5,*) label, upf%oc(i)
READ (label(1:1),*, err=10) l
upf%els(i) = label
upf%nchi(i) = l
l =spdf_to_l(label(2:2))
upf%lchi(i) = l
upf%rcut_chi(i) = 0.0d0
upf%rcutus_chi(i)= 0.0d0
upf%epseu(i) = 0.0d0
ENDDO
upf%mesh = mesh
upf%dx = amesh
upf%rmax = rmax
upf%xmin = xmin
upf%zmesh= z
ALLOCATE(upf%rab(upf%mesh))
ALLOCATE(upf%r(upf%mesh))
upf%r(:) = r(1:upf%mesh)
upf%rab(:)=upf%r(:)*amesh
ALLOCATE (upf%rho_atc(upf%mesh))
IF (upf%nlcc) upf%rho_atc(:) = rho_atc(1:upf%mesh)
upf%rcloc = rcloc
ALLOCATE (upf%vloc(upf%mesh))
!
! the factor e2=2 converts from Hartree to Rydberg
!
IF ( upf%typ == "SL" ) THEN
upf%vloc(:) = vnl(1:upf%mesh,upf%lloc)*e2
ALLOCATE(upf%vnl(upf%mesh,0:upf%lmax,1))
upf%vnl(:,:,1) = vnl(1:upf%mesh,0:upf%lmax)*e2
upf%nbeta= lmax
ELSE
DO i=1,upf%mesh
x = upf%r(i)/rc
x2=x**2
upf%vloc(i) = e2 * ( -upf%zp*erf(x/sqrt(2.d0))/upf%r(i) + &
exp ( -0.5d0*x2 ) * (c(1) + x2*( c(2) + x2*( c(3) + x2*c(4) ) ) ) )
ENDDO
upf%nbeta=0
DO l=0,upf%lmax
upf%nbeta = upf%nbeta + nl(l)
ENDDO
ENDIF
IF (upf%nbeta > 0) THEN
ALLOCATE(upf%els_beta(upf%nbeta) )
ALLOCATE(upf%lll(upf%nbeta))
ALLOCATE(upf%kbeta(upf%nbeta))
IF ( pstype == 3 ) THEN
upf%kbeta(:) = upf%mesh
ELSE
iv=0 ! counter on beta functions
DO i=1,upf%nwfc
l=upf%lchi(i)
IF (l/=lloc) THEN
iv=iv+1
upf%kbeta(iv)=upf%mesh
DO ir = upf%mesh,1,-1
IF ( abs ( vnl(ir,l) - vnl(ir,lloc) ) > 1.0E-6 ) THEN
! include points up to the last with nonzero value
upf%kbeta(iv)=ir+1
exit
ENDIF
ENDDO
ENDIF
ENDDO
! the number of points used in the evaluation of integrals
! should be even (for simpson integration)
DO i=1,upf%nbeta
IF ( mod (upf%kbeta(i),2) == 0 ) upf%kbeta(i)=upf%kbeta(i)+1
upf%kbeta(i)=MIN(upf%mesh,upf%kbeta(i))
ENDDO
upf%kkbeta = maxval(upf%kbeta(:))
ENDIF
ALLOCATE(upf%beta(upf%mesh,upf%nbeta))
ALLOCATE(upf%dion(upf%nbeta,upf%nbeta))
upf%beta(:,:) =0.d0
upf%dion(:,:) =0.d0
ALLOCATE(upf%rcut (upf%nbeta))
ALLOCATE(upf%rcutus(upf%nbeta))
IF ( pstype == 3 ) THEN
iv=0 ! counter on beta functions
DO l=0,upf%lmax
ij = 0
DO i=1, nl(l)
iv = iv+1
upf%lll(iv)=l
WRITE (upf%els_beta(iv), '(I1,A1)' ) i, l_to_spdf(l)
DO j=i, nl(l)
jv = iv+j-i
ij=ij+1
upf%dion(iv,jv) = h(l,ij)/e2
IF ( j > i ) upf%dion(jv,iv) = upf%dion(iv,jv)
ENDDO
fac= sqrt(2d0*rl(l)) / ( rl(l)**(l+2*i) * sqrt(mygamma(l+2*i)) )
DO ir=1,upf%mesh
x2 = (upf%r(ir)/rl(l))**2
upf%beta(ir,iv) = upf%r(ir)**(l+2*(i-1)) * &
exp ( -0.5d0*x2 ) * fac * e2
! ...remember: the beta functions in the UPF format
! ...have to be multiplied by a factor r !!!
upf%beta(ir,iv) = upf%beta(ir,iv)*upf%r(ir)
!
ENDDO
! look for index kbeta such that v(i)=0 if i>kbeta
DO ir=upf%mesh,1,-1
IF ( abs(upf%beta(ir,iv)) > 1.D-12 ) exit
ENDDO
IF ( ir < 2 ) THEN
WRITE(stdout,'("cpmd2upf: zero beta function?!?")')
return
ELSEIF ( mod(ir,2) /= 0 ) THEN
! even index
upf%kbeta(iv) = ir
ELSEIF ( ir < upf%mesh .and. mod(ir,2) == 0 ) THEN
! odd index
upf%kbeta(iv) = ir+1
ELSE
upf%kbeta(iv) = upf%mesh
ENDIF
! not really the same thing as rc in PP generation
upf%rcut (iv) = upf%r(upf%kbeta(iv))
upf%rcutus(iv) = 0.0
ENDDO
ENDDO
upf%kkbeta = maxval(upf%kbeta(:))
ELSE
ALLOCATE(aux(upf%kkbeta))
iv=0 ! counter on beta functions
DO i=1,upf%nwfc
l=upf%lchi(i)
IF (l/=lloc) THEN
iv=iv+1
upf%lll(iv)=l
upf%els_beta(iv)=upf%els(i)
DO ir=1,upf%kbeta(iv)
! the factor e2 converts from Hartree to Rydberg
upf%beta(ir,iv) = e2 * chi(ir,l+1) * &
( vnl(ir,l) - vnl(ir,lloc) )
aux(ir) = chi(ir,l+1) * upf%beta(ir,iv)
ENDDO
upf%rcut (iv) = upf%r(upf%kbeta(iv))
upf%rcutus(iv) = 0.0
CALL simpson(upf%kbeta(iv),aux,upf%rab,vll)
upf%dion(iv,iv) = 1.0d0/vll
ENDIF
ENDDO
DEALLOCATE(aux)
ENDIF
ELSE
! prevents funny errors when writing file
ALLOCATE(upf%dion(upf%nbeta,upf%nbeta))
ENDIF
ALLOCATE (upf%chi(upf%mesh,upf%nwfc))
upf%chi(:,:) = chi(1:upf%mesh,1:upf%nwfc)
ALLOCATE (upf%rho_at(upf%mesh))
upf%rho_at(:) = 0.d0
DO i=1,upf%nwfc
upf%rho_at(:) = upf%rho_at(:) + upf%oc(i) * upf%chi(:,i) ** 2
ENDDO
! ----------------------------------------------------------
WRITE (6,'(a)') 'Pseudopotential successfully converted'
! ----------------------------------------------------------
RETURN
END SUBROUTINE convert_cpmd
!
! ------------------------------------------------------------------
SUBROUTINE check_radial_grid ( amesh_, mesh, r, amesh, ierr )
! ------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER, INTENT (in) :: mesh
INTEGER, INTENT (out) :: ierr
REAL(dp), INTENT (in) :: amesh_, r(mesh)
REAL(dp), INTENT (out):: amesh
INTEGER :: i
!
ierr = 1
! get amesh if not available directly, check its value otherwise
WRITE(stdout, "('Radial grid r(i) has ',i4,' points')") mesh
WRITE(stdout, "('Assuming log radial grid: r(i)=exp[(i-1)*amesh]*r(1), with:')")
IF (amesh_ < 0.0d0) THEN
amesh = log (r(mesh)/r(1))/(mesh-1)
WRITE(stdout,"('amesh = log (r(mesh)/r(1))/(mesh-1) = ',f10.6)") amesh
ELSE
! not clear whether the value of amesh read from file
! matches the above definition, or if it is exp(amesh) ...
amesh = log (r(mesh)/r(1))/(mesh-1)
IF ( abs ( amesh - amesh_ ) > 1.0d-5 ) THEN
IF ( abs ( amesh - exp(amesh_) ) < 1.0d-5 ) THEN
amesh = log(amesh_)
WRITE(stdout, "('amesh = log (value read from file) = ',f10.6)") amesh
ELSE
WRITE(stdout,'("check_radial_grid: unknown real-space grid")')
return
ENDIF
ELSE
amesh = amesh_
WRITE(stdout,"('amesh = value read from file = ',f10.6)") amesh
ENDIF
ENDIF
! check if the grid is what we expect
DO i=2,mesh
IF ( abs(r(i) - exp((i-1)*amesh)*r(1)) > 1.0d-5) THEN
WRITE(stdout,"('check_radial_grid: grid point ',i4,', found ', &
& f10.6,', expected ',f10.6)") i, r(i), exp((i-1)*amesh)*r(1)
return
ENDIF
ENDDO
ierr = 0
!
END SUBROUTINE check_radial_grid
! ------------------------------------------------------------------
FUNCTION mygamma ( n )
!------------------------------------------------------------------
!
! mygamma(n) = \Gamma(n-1/2) = sqrt(pi)*(2n-3)!!/2**(n-1)
!
USE upf_const, ONLY : pi
IMPLICIT NONE
REAL(dp) :: mygamma
INTEGER, INTENT(in) :: n
INTEGER :: i
!
IF ( n < 2 ) write(stdout,'("mygamma: unexpected input argument")')
mygamma = sqrt(pi) / 2.d0**(n-1)
! next factor is (2n-3)!!
do i = 2*n-3, 1, -2
mygamma = i*mygamma
end do
!
END FUNCTION mygamma
END MODULE cpmd_module
Reference in New Issue View Git Blame Copy Permalink