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

Clean up of the atomic code. Improved separation of pseudopotential 

generation and test.


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@3945 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
dalcorso 2007-05-12 11:01:56 +00:00
parent bfd169fce2
commit f352671639
22 changed files with 761 additions and 580 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
atomic/Makefile
View File

@ -8,17 +8,19 @@ LD1OBJS = \
add_exchange.o \
all_electron.o \
ascheq.o \
ascheqlocps.o \
ascheqps.o \
ascheqps_drv.o \
c6_dft.o \
c6_tfvw.o \
cfdsol.o \
chargeps.o \
compute_chi.o \
compute_chi_tm.o \
compute_det.o \
compute_phi.o \
compute_phipot.o \
compute_phi_tm.o \
compute_phius.o \
compute_potps.o \
compute_solution.o \
descreening.o \
dfx_new.o \
@ -65,6 +67,7 @@ scf.o \
seriebes.o \
series.o \
set_rho_core.o \
set_psi_in.o \
set_sl3.o \
sic_correction.o \
starting_potential.o \

89
atomic/ascheqps_drv.f90 Normal file
View File

@ -0,0 +1,89 @@
!
! Copyright (C) 2007 QUANTUM-ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
!--------------------------------------------------------------------------
subroutine ascheqps_drv(veff, ncom, thresh, flag_all)
!--------------------------------------------------------------------------
!
! This routine is a driver that calculates for the test
! configuration the solutions of the Kohn and Sham equation
! with a fixed pseudo-potential. The potentials are
! assumed to be screened. The effective potential veff is given
! in input.
! The output wavefunctions are written in phits and are normalized.
! If flag is .true. compute all wavefunctions, otherwise only
! the wavefunctions with positive occupation.
!
use ld1inc
implicit none
integer :: &
ncom ! number of components of the pseudopotential
real(DP) :: &
veff(ndm,ncom) ! work space for writing the potential
logical :: flag_all ! if true calculates all the wavefunctions
integer :: &
ns, & ! counter on pseudo functions
is, & ! counter on spin
nbf, & ! auxiliary nbeta
n, & ! index on r point
ind
real(DP) :: &
vaux(ndm,2) ! work space for writing the potential
real(DP) :: thresh ! threshold for selfconsistency
!
! compute the pseudowavefunctions in the test configuration
!
if (pseudotype.eq.1) then
nbf=0
else
nbf=nbeta
endif
do ns=1,nwfts
if (octs(ns).gt.0.0_dp.or.(octs(ns).gt.-1.0_dp .and. flag_all)) then
is=iswts(ns)
if (ncom==1.and.is==2) call &
errore('ascheqps_drv','uncompatible spin',1)
if (pseudotype ==1) then
if ( rel < 2 .or. llts(ns) == 0 .or. &
abs(jjts(ns)-llts(ns)+0.5_dp) < 0.001_dp) then
ind=1
else if ( rel == 2 .and. llts(ns) > 0 .and. &
abs(jjts(ns)-llts(ns)-0.5_dp) < 0.001_dp) then
ind=2
else
call errore('ascheqps_drv','something strange',1)
endif
do n=1,mesh
vaux(n,is)=veff(n,is)+vnl(n,llts(ns),ind)
enddo
else
do n=1,mesh
vaux(n,is)=veff(n,is)
enddo
endif
call ascheqps(nnts(ns),llts(ns),jjts(ns),enlts(ns), &
mesh,ndm,dx,r,r2,sqr,vaux(1,is),thresh,phits(1,ns), &
betas,ddd(1,1,is),qq,nbf,nwfsx,lls,jjs,ikk)
! write(6,*) ns, nnts(ns),llts(ns), jjts(ns), enlts(ns)
!
! normalize the wavefunctions
!
call normalize(phits(1,ns),llts(ns),jjts(ns))
endif
enddo
return
end subroutine ascheqps_drv

45
atomic/chargeps.f90
View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2004 PWSCF group
! Copyright (C) 2007 QUANTUM-ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
@ -7,22 +7,23 @@
!
!
!---------------------------------------------------------------
subroutine chargeps(nwff,lli,jji,oci,iswfi)
subroutine chargeps(rho_i,phi_i,nwf_i,ll_i,jj_i,oc_i,iswf_i)
!---------------------------------------------------------------
!
! calculate the (spherical) pseudo charge density and
! spin polarization
! calculate the (spherical) pseudo charge density
!
use ld1inc
integer :: &
nwff, & ! input: the number of wavefunctions
iswfi(nwff),& ! input: their spin
lli(nwff) ! input: their angular momentum
nwf_i, & ! input: the number of wavefunctions
iswf_i(nwfsx),& ! input: their spin
ll_i(nwfsx) ! input: their angular momentum
real(DP) :: &
jji(nwff), & ! input: their total angular momentum
oci(nwff) ! input: the occupation
jj_i(nwfsx), & ! input: their total angular momentum
oc_i(nwfsx), & ! input: the occupation
phi_i(ndm,nwfsx), & ! input: the functions to add
rho_i(ndm,2) ! output: the (nspin) components of the charge
integer :: &
is, & ! counter on spin
@ -35,15 +36,15 @@ subroutine chargeps(nwff,lli,jji,oci,iswfi)
gi(ndm) ! used to compute the integrals
rhos=0.0_dp
rho_i=0.0_dp
!
! compute the square modulus of the eigenfunctions
!
do ns=1,nwff
if (oci(ns).gt.0.0_dp) then
is=iswfi(ns)
do ns=1,nwf_i
if (oc_i(ns).gt.0.0_dp) then
is=iswf_i(ns)
do n=1,mesh
rhos(n,is)=rhos(n,is)+oci(ns)*phis(n,ns)**2
rho_i(n,is)=rho_i(n,is)+oc_i(ns)*phi_i(n,ns)**2
end do
endif
enddo
@ -51,16 +52,16 @@ subroutine chargeps(nwff,lli,jji,oci,iswfi)
! if US pseudopotential compute the augmentation part
!
if (pseudotype.eq.3) then
do ns=1,nwff
if (oci(ns).gt.0.0_dp) then
is=iswfi(ns)
do ns=1,nwf_i
if (oc_i(ns).gt.0.0_dp) then
is=iswf_i(ns)
do n1=1,nbeta
if (lli(ns).eq.lls(n1).and. &
abs(jji(ns)-jjs(n1)).lt.1.e-7_dp) then
nst=(lli(ns)+1)*2
if (ll_i(ns).eq.lls(n1).and. &
abs(jj_i(ns)-jjs(n1)).lt.1.e-7_dp) then
nst=(ll_i(ns)+1)*2
ikl=ikk(n1)
do n=1,ikl
gi(n)=betas(n,n1)*phis(n,ns)
gi(n)=betas(n,n1)*phi_i(n,ns)
enddo
work(n1)=int_0_inf_dr(gi,r,r2,dx,ikl,nst)
else
@ -73,7 +74,7 @@ subroutine chargeps(nwff,lli,jji,oci,iswfi)
do n1=1,nbeta
do n2=1,nbeta
do n=1,mesh
rhos(n,is)=rhos(n,is)+qvan(n,n1,n2)*oci(ns)* &
rho_i(n,is)=rho_i(n,is)+qvan(n,n1,n2)*oc_i(ns)* &
work(n1)*work(n2)
enddo
enddo

97
atomic/compute_chi.f90
View File

@ -7,7 +7,7 @@
!
!
!--------------------------------------------------------------------------
subroutine compute_chi(lam,ik,ns,xc,lbes4)
subroutine compute_chi(lam,ik,ikk_in,phi_in,chi_out,xc,e,lbes4)
!--------------------------------------------------------------------------
!
! This routine computes the chi functions:
@ -19,41 +19,30 @@ subroutine compute_chi(lam,ik,ns,xc,lbes4)
implicit none
integer :: &
ik, & ! the point corresponding to rc
ikk_in,& ! the point after which the chi should be zero
ns, & ! the wavefunction
lam ! the angular momentum
logical :: &
lbes4
real(DP) :: &
xc(8)
e, & ! input: the energy
xc(8), & ! input: the coefficients of the bessel function
phi_in(ndm), & ! input: pseudo wavefunction
chi_out(ndm) ! output: the chi function
!
real(DP) :: &
j1(ndm),aux(ndm), &
j1(ndm), aux(ndm), gi(ndm),&
b(4),c(4), arow(ndm),brow(ndm),crow(ndm),drow(ndm), &
b0e, g0, g1, g2, &
ddx12, &
x4l6, &
x6l12, dpoly
real(DP), external :: pr, d2pr, dpr
ddx12, &
x4l6, &
x6l12, &
int_0_inf_dr, &
integral
integer :: &
n, nstart
!
! Troullier-Martins: use the analytic formula
!
if (tm) then
do n=1,ik
dpoly = dpr(xc,xc(7),r(n))
! dpr = first derivate of polynomial pr
! d2pr= second derivate of polynomial pr
chis(n,ns) = (enls(ns) + (2*lam+2)/r(n)*dpoly + &
d2pr(xc,xc(7),r(n)) + dpoly**2 - vpsloc(n))*phis(n,ns)
enddo
do n = ik+1,mesh
chis(n,ns) = (vpot(n,1) - vpsloc(n))*phis(n,ns)
enddo
return
end if
n, nstart, nst
!
! RRKJ: first expand in a taylor series the phis function
! Since we know that the phis functions are a sum of Bessel
@ -66,7 +55,7 @@ subroutine compute_chi(lam,ik,ns,xc,lbes4)
x6l12=6*lam+12
do n=1,6
j1(n)=phis(n,ns)/r(n)**(lam+1)
j1(n)=phi_in(n)/r(n)**(lam+1)
enddo
call seriesbes(j1,r,r2,6,c)
!
@ -127,14 +116,14 @@ subroutine compute_chi(lam,ik,ns,xc,lbes4)
!
! and compute the taylor expansion of the chis
!
b0e=(b(1)-enls(ns))
b0e=(b(1)-e)
g0=x4l6*c(3)-b0e*c(1)
g1=x6l12*c(4)-c(1)*b(2)
g2=-(b0e*c(3)+b(3)*c(1))
nstart=5
do n=1,nstart-1
chis(n,ns)= (g0+r(n)*(g1+g2*r(n)))*r(n)**(lam+3)/sqr(n)
chi_out(n)= (g0+r(n)*(g1+g2*r(n)))*r(n)**(lam+3)/sqr(n)
enddo
do n=1,mesh
aux(n)= (g0+r(n)*(g1+g2*r(n)))
@ -143,53 +132,52 @@ subroutine compute_chi(lam,ik,ns,xc,lbes4)
! set up the equation
!
do n=1,mesh
phis(n,ns)=phis(n,ns)/sqr(n)
gi(n)=phi_in(n)/sqr(n)
enddo
do n=1,mesh
j1(n)=r2(n)*(vpsloc(n)-enls(ns))+(lam+0.5_dp)**2
j1(n)=r2(n)*(vpsloc(n)-e)+(lam+0.5_dp)**2
j1(n)=1.0_dp-ddx12*j1(n)
enddo
do n=nstart,mesh-3
drow(n)= phis(n+1,ns)*j1(n+1) &
+ phis(n,ns)*(-12.0_dp+10.0_dp*j1(n))+ &
phis(n-1,ns)*j1(n-1)
drow(n)= gi(n+1)*j1(n+1) &
+ gi(n)*(-12.0_dp+10.0_dp*j1(n))+ &
gi(n-1)*j1(n-1)
brow(n)=10.0_dp*ddx12
crow(n)=ddx12
arow(n)=ddx12
enddo
drow(nstart)=drow(nstart)-ddx12*chis(nstart-1,ns)
chis(mesh-2,ns)=0.0_dp
chis(mesh-1,ns)=0.0_dp
chis(mesh,ns)=0.0_dp
drow(nstart)=drow(nstart)-ddx12*chi_out(nstart-1)
chi_out(mesh-2)=0.0_dp
chi_out(mesh-1)=0.0_dp
chi_out(mesh)=0.0_dp
!
! and solve it
!
call tridiag(arow(nstart),brow(nstart),crow(nstart), &
drow(nstart),chis(nstart,ns),mesh-3-nstart)
drow(nstart),chi_out(nstart),mesh-3-nstart)
!
! put the correct normalization and r dependence
!
do n=1,mesh
phis(n,ns)=phis(n,ns)*sqr(n)
chis(n,ns)=chis(n,ns)*sqr(n)/r2(n)
chi_out(n)=chi_out(n)*sqr(n)/r2(n)
! if(lam.eq.0)
! + write(*,'(5(e20.13,1x))')
! + r(n),chis(n,ns),chis(n,ns)/r(n)**(lam+1),
! + r(n),chi_out(n),chi_out(n)/r(n)**(lam+1),
! + aux(n),aux(n)*r(n)**(lam+1)
enddo
!
! smooth close to the origin with asymptotic expansion
!
do n=nstart,mesh
if (abs(chis(n,ns)/r(n)**(lam+1)-aux(n)) &
if (abs(chi_out(n)/r(n)**(lam+1)-aux(n)) &
.lt.1.e-3_dp*abs(aux(n)) ) goto 100
chis(n,ns)=aux(n)*r(n)**(lam+1)
chi_out(n)=aux(n)*r(n)**(lam+1)
enddo
100 if (n.eq.mesh+1.or.r(n).gt.0.05_dp)then
print*,lam,ns,n,mesh,r(n)
print*,lam,n,mesh,r(n)
call errore('compute_chi','n is too large',1)
endif
!
@ -197,9 +185,28 @@ subroutine compute_chi(lam,ik,ns,xc,lbes4)
!
do n=mesh,1,-1
if (r(n).lt.9.0_dp) goto 200
chis(n,ns)=0.0_dp
chi_out(n)=0.0_dp
enddo
200 continue
! check that the chi are zero beyond ikk
nst=0
gi=0.0_dp
do n=ikk_in+1,mesh
gi(n)=chi_out(n)**2
enddo
do n=min(ikk_in+20,mesh),mesh
chi_out(n)=0.0_dp
enddo
integral=int_0_inf_dr(gi,r,r2,dx,mesh,nst)
if (integral > 2.e-6_dp) then
write(6, '(5x,''ns='',i4,'' l='',i4, '' integral='',f15.9, &
& '' r(ikk) '',f15.9)') ns, lam, integral, r(ikk_in)
call infomsg ('gener_pseudo ','chi too large beyond r_c', -1)
do n=ikk_in,mesh
write(6,*) r(n),gi(n)
enddo
stop
endif
return
end subroutine compute_chi

52
atomic/compute_chi_tm.f90 Normal file
View File

@ -0,0 +1,52 @@
!
! Copyright (C) 2004 PWSCF group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
!--------------------------------------------------------------------------
subroutine compute_chi_tm(lam,ik,ikk_in,phi_in,chi_out,xc,e)
!--------------------------------------------------------------------------
!
! This routine computes the chi functions:
! |chi> = (\epsilon -T -V_{loc)) |psi>
!
use kinds, only : DP
use ld1inc
implicit none
integer :: &
ik, & ! the point corresponding to rc
ikk_in,& ! the point after which the chi should be zero
lam ! the angular momentum
real(DP) :: &
e, & ! input: the energy
xc(8), & ! input: the parameters of the fit
phi_in(ndm), & ! input: pseudo wavefunction
chi_out(ndm) ! output: the chi function
!
real(DP) :: &
dpoly
real(DP), external :: pr, d2pr, dpr
integer :: &
n
!
! Troullier-Martins: use the analytic formula
!
do n=1,ik
dpoly = dpr(xc,xc(7),r(n))
! dpr = first derivate of polynomial pr
! d2pr= second derivate of polynomial pr
chi_out(n) = (e + (2*lam+2)/r(n)*dpoly + &
d2pr(xc,xc(7),r(n)) + dpoly**2 - vpsloc(n))*phi_in(n)
enddo
do n = ik+1,mesh
chi_out(n) = (vpot(n,1) - vpsloc(n))*phi_in(n)
enddo
return
end subroutine compute_chi_tm

335
atomic/compute_phi.f90
View File

@ -6,13 +6,13 @@
! or http://www.gnu.org/copyleft/gpl.txt .
!
!--------------------------------------------------------------------------
subroutine compute_phi(lam,ik,nwf0,ns,xc,iflag,iok,occ)
subroutine compute_phi(lam,ik,chir,phi_out,xc,iflag,occ,e,els_in)
!--------------------------------------------------------------------------
!
! This routine computes the phi functions by pseudizing the
! all_electron chi functions. In input it receives, the point
! ik where the cut is done, the angular momentum lam,
! and the correspondence with the all electron wavefunction
! and the all electron wavefunction
!
!
!
@ -22,15 +22,20 @@ subroutine compute_phi(lam,ik,nwf0,ns,xc,iflag,iok,occ)
implicit none
integer :: &
lam, & ! the angular momentum
ik, & ! the point corresponding to rc
nwf0, & !
ns, & ! the function to pseudize
iflag,& ! if 1 print
iok ! output: number of nodes of the pseudized wavefunction
lam, & ! input: the angular momentum
ik, & ! input: the point corresponding to rc
iflag ! input: if 1 print the message
real(DP) :: &
xc(8), occ
xc(8), & ! output: the coefficients of the Bessel functions
occ, & ! input: the occupation of the level
e ! input: the energy of the level
real(DP) :: &
chir(ndm), & ! input: the all-electron function
phi_out(ndm) ! output: the phi function
character(len=2) :: els_in ! input: the label of the state
!
real(DP) :: &
fae, & ! the value of the all-electron function
@ -45,15 +50,16 @@ subroutine compute_phi(lam,ik,nwf0,ns,xc,iflag,iok,occ)
m, n, nst, nnode, nc, nc1, ij, imax, iq, i
real(DP) :: &
chir(ndm,nwfx), chi_dir(ndm,2), gi(ndm), j1(ndm,4), &
chi_dir(ndm,2), gi(ndm), j1(ndm,4), &
f1aep1, f1aem1, jnor, psnor, fact(4), &
cm(10), bm(4), ze2, cn(6), c2, &
delta, a, b, c, deter, gamma, &
lamda0, lamda3, lamda4, mu0, mu4, s0, s4, t0, t4
real(DP), external :: deriv_7pts, deriv2_7pts, int_0_inf_dr, pr
real(DP), external :: deriv_7pts, deriv2_7pts, int_0_inf_dr
integer :: &
nbes ! number of Bessel functions to be used
iok, & ! number of nodes
nbes ! number of Bessel functions to be used
!
! decide whether to use 4 Bessel functions or 3 (default)
@ -66,220 +72,173 @@ subroutine compute_phi(lam,ik,nwf0,ns,xc,iflag,iok,occ)
!
nst=(lam+1)*2
!
! compute the reference wavefunction
!
if (new(ns)) then
if (rel == 1) then
call lschps(3,zed,exp(dx),dx,mesh,mesh,mesh, &
1,lam,enls(ns),chir(1,ns),r,vpot)
elseif (rel == 2) then
call dir_outward(ndm,mesh,lam,jjs(ns),enls(ns),dx,chi_dir,r,rab,vpot)
chir(:,ns)=chi_dir(:,1)
else
ze2=-zed*2.0_dp
call intref(lam,enls(ns),mesh,dx,r,r2,sqr,vpot,ze2,chir(1,ns))
endif
!
! fix arbitrarily the norm at the cut-off radius equal to 0.5
!
jnor=chir(ik,ns)
do n=1,mesh
chir(n,ns)=chir(n,ns)*0.5_dp/jnor
enddo
else
do n=1,mesh
chir(n,ns)=psi(n,1,nwf0)
enddo
endif
!
! save the all-electron function for the PAW setup
!
if (lpaw) psipaw(1:mesh,ns) = chir(1:mesh,ns)
!
! compute the first and second derivative of all-electron function
!
fae=chir(ik,ns)
f1ae=deriv_7pts(chir(1,ns),ik,r(ik),dx)
f2ae=deriv2_7pts(chir(1,ns),ik,r(ik),dx)
fae=chir(ik)
f1ae=deriv_7pts(chir,ik,r(ik),dx)
f2ae=deriv2_7pts(chir,ik,r(ik),dx)
!
! compute the norm of the all-electron function
!
do n=1,ik+1
gi(n)=chir(n,ns)**2
gi(n)=chir(n)**2
enddo
faenor=int_0_inf_dr(gi,r,r2,dx,ik,nst)
!
!
if (tm) then
!
! TM: the pseudo-wavefunction is written as polynomial times exponential
!
call find_coefficients &
(ik, chir(1,ns), enls(ns), r, dx, faenor, vpot, cn, c2, lam, mesh)
do i=1,ik
phis(i,ns) =sign( r(i)**(lam+1)*exp(pr(cn,c2,r(i))), &
chir(ik,ns) )
end do
xc(1:6) = cn(:)
xc(7) = c2
else
! RRKJ: the pseudo-wavefunction is written as an expansion into (3 or 4)
! spherical bessel functions for r < r_c
! find q_i with the correct log derivatives
! RRKJ: the pseudo-wavefunction is written as an expansion into (3 or 4)
! spherical bessel functions for r < r_c
! find q_i with the correct log derivatives
call find_qi(f1ae/fae,xc(nbes+1),ik,lam,nbes,1,iok)
if (iok.ne.0) return
!
! compute the bessel functions
!
do nc=1,nbes
call sph_bes(ik+5,r,xc(nbes+nc),lam,j1(1,nc))
jnor=j1(ik,nc)*r(ik)
fact(nc)=chir(ik,ns)/jnor
do n=1,ik+5
j1(n,nc)=j1(n,nc)*r(n)*chir(ik,ns)/jnor
enddo
call find_qi(f1ae/fae,xc(nbes+1),ik,lam,nbes,1,iok)
if (iok.ne.0) &
call errore('compute phi','problem with the q_i coefficients',1)
!
! compute the bessel functions
!
do nc=1,nbes
call sph_bes(ik+5,r,xc(nbes+nc),lam,j1(1,nc))
jnor=j1(ik,nc)*r(ik)
fact(nc)=chir(ik)/jnor
do n=1,ik+5
j1(n,nc)=j1(n,nc)*r(n)*chir(ik)/jnor
enddo
enddo
!
! compute the bm functions (second derivative of the j1)
! and the integrals for cm
!
ij=0
do nc=1, nbes
bm(nc)=deriv2_7pts(j1(1,nc),ik,r(ik),dx)
do nc1=1,nc
ij=ij+1
do n=1,ik
gi(n)=j1(n,nc)*j1(n,nc1)
enddo
cm(ij)=int_0_inf_dr(gi,r,r2,dx,ik,nst)
ij=0
do nc=1, nbes
bm(nc)=deriv2_7pts(j1(1,nc),ik,r(ik),dx)
do nc1=1,nc
ij=ij+1
do n=1,ik
gi(n)=j1(n,nc)*j1(n,nc1)
enddo
cm(ij)=int_0_inf_dr(gi,r,r2,dx,ik,nst)
enddo
enddo
!
! solve the second order equation
!
if (nbes == 4) then
wmax=0.0_dp
do n=1,mesh
if (abs(chir(n,ns)) > wmax .and. r(n) < 4.0_dp) then
wmax=abs(chir(n,ns))
if(chir(n,ns).lt.0.0_dp)then
isign=-1
else
isign=+1
endif
if (nbes == 4) then
wmax=0.0_dp
do n=1,mesh
if (abs(chir(n)) > wmax .and. r(n) < 4.0_dp) then
wmax=abs(chir(n))
if(chir(n).lt.0.0_dp)then
isign=-1
else
isign=+1
endif
enddo
lamda0=(f2ae-bm(1))/(bm(2)-bm(1))
lamda3=(bm(3)-bm(1))/(bm(2)-bm(1))
lamda4=(bm(4)-bm(1))/(bm(2)-bm(1))
if (new(ns)) then
den=1.0_dp
else
den=abs(occ)
endif
mu0=(isign*sqrt(rho0*4.0_dp*pi/den) &
-fact(1)+fact(1)*lamda0-fact(2)*lamda0) &
/(fact(1)*lamda3-fact(1)-fact(2)*lamda3+fact(3))
mu4=(fact(1)*lamda4-fact(1)-fact(2)*lamda4+fact(4)) &
/(fact(1)*lamda3-fact(1)-fact(2)*lamda3+fact(3))
s0=lamda0-mu0*lamda3
s4=lamda4-mu4*lamda3
t0=1.0_dp-lamda0+mu0*lamda3-mu0
t4=1.0_dp+mu4*lamda3-mu4-lamda4
a=cm(1)*t4**2+cm(3)*s4**2+cm(6)*mu4**2+cm(10) &
+2.0_dp*cm(4)*t4*mu4+2.0_dp*cm(2)*t4*s4+2.0_dp*cm(5)*s4*mu4 &
-2.0_dp*cm(7)*t4-2.0_dp*cm(8)*s4-2.0_dp*cm(9)*mu4
b=-2.0_dp*cm(1)*t0*t4-2.0_dp*cm(3)*s0*s4-2.0_dp*cm(6)*mu0*mu4 &
-2.0_dp*cm(4)*t0*mu4-2.0_dp*cm(4)*t4*mu0-2.0_dp*cm(2)*t0*s4 &
-2.0_dp*cm(2)*t4*s0-2.0_dp*cm(5)*s0*mu4-2.0_dp*cm(5)*s4*mu0 &
+2.0_dp*cm(7)*t0+2.0_dp*cm(8)*s0+2.0_dp*cm(9)*mu0
c=cm(1)*t0**2+cm(3)*s0**2+cm(6)*mu0**2+2.0_dp*cm(4)*t0*mu0 &
+2.0_dp*cm(2)*t0*s0+2.0_dp*cm(5)*s0*mu0-faenor
deter=b**2-4.0_dp*a*c
if (deter < 0.0_dp) then
call infomsg ('compute phi','negative determinant', -1)
write(6,100) ns,f1ae/fae, f2ae, faenor
100 format(/5x,i5,' ld= ',f10.6,' f2ae',f10.6,' faenor',f10.6)
iok=1
return
endif
xc(4)=(-b+sqrt(deter))/(2.0_dp*a)
xc(3)=mu0-mu4*xc(4)
xc(2)=s0-s4*xc(4)
xc(1)=t0-t4*xc(4)
!
do n=1,ik
phis(n,ns)= xc(1)*j1(n,1) + xc(2)*j1(n,2) &
+ xc(3)*j1(n,3) + xc(4)*j1(n,4)
enddo
else
gamma=(bm(3)-bm(1))/(bm(2)-bm(1))
delta=(f2ae-bm(1))/(bm(2)-bm(1))
a=(gamma-1.0_dp)**2*cm(1)+gamma**2*cm(3) &
-2.0_dp*gamma*(gamma-1.0_dp)*cm(2) &
+2.0_dp*(gamma-1.0_dp)*cm(4)-2.0_dp*gamma*cm(5)+cm(6)
b=2.0_DP*(1.0_dp-delta)*(gamma-1.0_dp)*cm(1) &
-2.0_dp*gamma*delta*cm(3) &
-2.0_dp*gamma*(1.0_dp-delta)*cm(2) &
+2.0_dp*delta*(gamma-1.0_dp)*cm(2) &
+2.0_dp*(1.0_dp-delta)*cm(4)+2.0_dp*delta*cm(5)
c=-faenor+(1.0_dp-delta)**2*cm(1)+delta**2*cm(3) &
+2.0_dp*delta*(1.0_dp-delta)*cm(2)
deter=b**2-4.0_dp*a*c
if (deter < 0.0_dp) then
call infomsg ('compute phi','negative determinant', -1)
write(6,110) ns,f1ae/fae, f2ae, faenor
110 format (/5x,i5,' ld= ',f10.6,' f2ae',f10.6, ' faenor',f10.6)
iok=1
return
endif
xc(3)=(-b+sqrt(deter))/(2.0_dp*a)
xc(2)=-xc(3)*gamma+delta
xc(1)=1.0_dp-xc(2)-xc(3)
do n=1,ik
phis(n,ns)= xc(1)*j1(n,1) + xc(2)*j1(n,2) &
+ xc(3)*j1(n,3)
enddo
endif
do nc=1,nbes
xc(nc)=xc(nc)*fact(nc)
enddo
end if
lamda0=(f2ae-bm(1))/(bm(2)-bm(1))
lamda3=(bm(3)-bm(1))/(bm(2)-bm(1))
lamda4=(bm(4)-bm(1))/(bm(2)-bm(1))
den=abs(occ)
mu0=(isign*sqrt(rho0*4.0_dp*pi/den) &
-fact(1)+fact(1)*lamda0-fact(2)*lamda0) &
/(fact(1)*lamda3-fact(1)-fact(2)*lamda3+fact(3))
mu4=(fact(1)*lamda4-fact(1)-fact(2)*lamda4+fact(4)) &
/(fact(1)*lamda3-fact(1)-fact(2)*lamda3+fact(3))
s0=lamda0-mu0*lamda3
s4=lamda4-mu4*lamda3
t0=1.0_dp-lamda0+mu0*lamda3-mu0
t4=1.0_dp+mu4*lamda3-mu4-lamda4
a=cm(1)*t4**2+cm(3)*s4**2+cm(6)*mu4**2+cm(10) &
+2.0_dp*cm(4)*t4*mu4+2.0_dp*cm(2)*t4*s4+2.0_dp*cm(5)*s4*mu4 &
-2.0_dp*cm(7)*t4-2.0_dp*cm(8)*s4-2.0_dp*cm(9)*mu4
b=-2.0_dp*cm(1)*t0*t4-2.0_dp*cm(3)*s0*s4-2.0_dp*cm(6)*mu0*mu4 &
-2.0_dp*cm(4)*t0*mu4-2.0_dp*cm(4)*t4*mu0-2.0_dp*cm(2)*t0*s4 &
-2.0_dp*cm(2)*t4*s0-2.0_dp*cm(5)*s0*mu4-2.0_dp*cm(5)*s4*mu0 &
+2.0_dp*cm(7)*t0+2.0_dp*cm(8)*s0+2.0_dp*cm(9)*mu0
c=cm(1)*t0**2+cm(3)*s0**2+cm(6)*mu0**2+2.0_dp*cm(4)*t0*mu0 &
+2.0_dp*cm(2)*t0*s0+2.0_dp*cm(5)*s0*mu0-faenor
deter=b**2-4.0_dp*a*c
if (deter < 0.0_dp) then
call infomsg ('compute phi','negative determinant', -1)
write(6,100) f1ae/fae, f2ae, faenor
100 format(/5x,' ld= ',f10.6,' f2ae',f10.6,' faenor',f10.6)
iok=1
return
endif
xc(4)=(-b+sqrt(deter))/(2.0_dp*a)
xc(3)=mu0-mu4*xc(4)
xc(2)=s0-s4*xc(4)
xc(1)=t0-t4*xc(4)
!
do n=1,ik
phi_out(n)= xc(1)*j1(n,1) + xc(2)*j1(n,2) &
+ xc(3)*j1(n,3) + xc(4)*j1(n,4)
enddo
else
gamma=(bm(3)-bm(1))/(bm(2)-bm(1))
delta=(f2ae-bm(1))/(bm(2)-bm(1))
a=(gamma-1.0_dp)**2*cm(1)+gamma**2*cm(3) &
-2.0_dp*gamma*(gamma-1.0_dp)*cm(2) &
+2.0_dp*(gamma-1.0_dp)*cm(4)-2.0_dp*gamma*cm(5)+cm(6)
b=2.0_DP*(1.0_dp-delta)*(gamma-1.0_dp)*cm(1) &
-2.0_dp*gamma*delta*cm(3) &
-2.0_dp*gamma*(1.0_dp-delta)*cm(2) &
+2.0_dp*delta*(gamma-1.0_dp)*cm(2) &
+2.0_dp*(1.0_dp-delta)*cm(4)+2.0_dp*delta*cm(5)
c=-faenor+(1.0_dp-delta)**2*cm(1)+delta**2*cm(3) &
+2.0_dp*delta*(1.0_dp-delta)*cm(2)
deter=b**2-4.0_dp*a*c
if (deter < 0.0_dp) then
call infomsg ('compute phi','negative determinant', -1)
write(6,110) f1ae/fae, f2ae, faenor
110 format (/5x,' ld= ',f10.6,' f2ae',f10.6, ' faenor',f10.6)
iok=1
return
endif
xc(3)=(-b+sqrt(deter))/(2.0_dp*a)
xc(2)=-xc(3)*gamma+delta
xc(1)=1.0_dp-xc(2)-xc(3)
do n=1,ik
phi_out(n)= xc(1)*j1(n,1) + xc(2)*j1(n,2) &
+ xc(3)*j1(n,3)
enddo
endif
do nc=1,nbes
xc(nc)=xc(nc)*fact(nc)
enddo
!
! for r > r(ik) the pseudo and all-electron psi(r) coincide
!
do n=ik+1,mesh
phis(n,ns)= chir(n,ns)
phi_out(n)= chir(n)
enddo
!
! check for the norm of the pseudo wavefunction
!
do n=1,ik
gi(n)=phis(n,ns)**2
gi(n)=phi_out(n)**2
enddo
psnor=int_0_inf_dr(gi,r,r2,dx,ik,nst)
if (iflag == 1) then
if (tm) then
write(6,120) els(ns), rcut(ns)
write(6,120) els_in, r(ik)
120 format (/ /5x, ' Wfc ',a3,' rcut=',f6.3, &
' Using Troullier-Martins method ')
else
write(6,130) els(ns),rcut(ns),2.0_dp*xc(6)**2
write(6,130) els_in,r(ik),2.0_dp*xc(6)**2
130 format (/ /5x, ' Wfc ',a3,' rcut=',f6.3, &
' Estimated cut-off energy= ', f11.2,' Ry')
if (nbes == 4) write(6,140) rho0
@ -293,9 +252,9 @@ subroutine compute_phi(lam,ik,nwf0,ns,xc,iflag,iok,occ)
!
nnode=0
do n=1,ik+1
if ( phis(n,ns) .ne. sign(phis(n,ns),phis(n+1,ns)) ) then
if (iflag==1) write(6,150) lam,ns,r(n)
150 format (5x,'l=',i4,' ns=',i4,' Node at ',f10.8)
if ( phi_out(n) .ne. sign(phi_out(n),phi_out(n+1)) ) then
if (iflag==1) write(6,150) lam,r(n)
150 format (5x,'l=',i4,' Node at ',f10.8)
nnode=nnode+1
endif
enddo

93
atomic/compute_phi_tm.f90 Normal file
View File

@ -0,0 +1,93 @@
!
! Copyright (C) 2004 PWSCF group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!--------------------------------------------------------------------------
subroutine compute_phi_tm(lam,ik,chir,phi_out,iflag,xc,e,els_in)
!--------------------------------------------------------------------------
!
! This routine computes the phi functions by pseudizing the
! all_electron chi functions. In input it receives, the point
! ik where the cut is done, the angular momentum lam,
! and the all electron wavefunction
!
!
!
use kinds, only : DP
use constants, only: pi
use ld1inc
implicit none
integer :: &
lam, & ! input: the angular momentum
ik, & ! input: the point corresponding to rc
iflag ! input: if 1 print the message
real(DP) :: &
e ! input: the energy of the level
real(DP) :: &
xc(8), & ! output: parameters of the fit
chir(ndm), & ! input: the all-electron function
phi_out(ndm) ! output: the phi function
character(len=2) :: els_in ! input: the label of the state
!
real(DP) :: &
faenor ! the norm of the function
integer :: &
n, nst
real(DP) :: &
gi(ndm), &
psnor, &
cn(6), c2
real(DP), external :: deriv_7pts, deriv2_7pts, int_0_inf_dr, pr
!
!
! compute the norm of the all-electron function
!
nst=(lam+1)*2
do n=1,ik+1
gi(n)=chir(n)**2
enddo
faenor=int_0_inf_dr(gi,r,r2,dx,ik,nst)
!
!
! TM: the pseudo-wavefunction is written as polynomial times exponential
!
call find_coefficients &
(ik, chir, e, r, dx, faenor, vpot, cn, c2, lam, mesh)
do n=1,ik
phi_out(n)=sign(r(n)**(lam+1)*exp(pr(cn,c2,r(n))),chir(ik))
end do
xc(1:6) = cn(:)
xc(7) = c2
!
! for r > r(ik) the pseudo and all-electron psi(r) coincide
!
do n=ik+1,mesh
phi_out(n)= chir(n)
enddo
!
! check for the norm of the pseudo wavefunction
!
! do n=1,ik
! gi(n)=phi_out(n)**2
! enddo
! psnor=int_0_inf_dr(gi,r,r2,dx,ik,nst)
! write(6,'(5x," AE norm = ",f6.3," PS norm = ",f6.3)') faenor, psnor
if (iflag == 1) then
write(6,120) els_in, r(ik)
120 format (/ /5x, ' Wfc ',a3,' rcut=',f6.3, &
' Using Troullier-Martins method ')
end if
return
end subroutine compute_phi_tm

45
atomic/compute_phius.f90
View File

@ -1,3 +1,4 @@
!
! Copyright (C) 2004 PWSCF group
! This file is distributed under the terms of the
@ -6,7 +7,7 @@
! or http://www.gnu.org/copyleft/gpl.txt .
!
!--------------------------------------------------------------------------
subroutine compute_phius(lam,ik,ns,xc,iflag)
subroutine compute_phius(lam,ik,psi_in,phi_out,xc,iflag,els_in)
!--------------------------------------------------------------------------
!
! This routine computes the phi functions by pseudizing the
@ -19,6 +20,12 @@
use ld1inc
implicit none
real(DP) :: &
psi_in(ndm), & ! input: the norm conserving wavefunction
phi_out(ndm) ! output: the us wavefunction
character(len=2) :: els_in
real(DP) :: &
fae, & ! the value of the all-electron function
f1ae, & ! its first derivative
@ -28,20 +35,19 @@ use ld1inc
integer :: &
ik, & ! the point corresponding to rc
ns, & ! the function to pseudize
iflag,& ! if 1 print
iok, & ! if 0 there are no problem
iflag, & ! if 1 print
iok, & ! if 0 there are no problem
lam ! the angular momentum
real(DP) :: &
f1aep1,f1aem1,jnor, & ! auxilairy quantities
bm(2), & ! the derivative of the bessel
ff, & ! contain deltah corrections
fact(2), & ! factor of normalization
j1(ndm,8) ! the bessel functions
real(DP) :: &
deriv_7pts, deriv2_7pts
deriv_7pts, deriv2_7pts, p1aep1, p1aem1
integer :: &
@ -51,18 +57,9 @@ use ld1inc
!
! compute first and second derivative
!
ff=1-dx**2/48.0_dp
! fae=(psipsus(ik+1,ns)+psipsus(ik,ns))*0.5_dp
! f1aep1=psipsus(ik+1,ns)*ff/sqr(ik+1)
! f1ae=psipsus(ik,ns)*(-12.0_dp+10.0_dp*ff)/sqr(ik)
! f1aem1=psipsus(ik-1,ns)*ff/sqr(ik-1)
! f2ae=(f1aep1+f1ae+f1aem1)/dx**2
! f1ae=(psipsus(ik+1,ns)-psipsus(ik,ns))/(r(ik+1)-r(ik))
fae=psipsus(ik,ns)
f1ae=deriv_7pts(psipsus(1,ns),ik,r(ik),dx)
f2ae=deriv2_7pts(psipsus(1,ns),ik,r(ik),dx)
fae=psi_in(ik)
f1ae=deriv_7pts(psi_in,ik,r(ik),dx)
f2ae=deriv2_7pts(psi_in,ik,r(ik),dx)
!
! find the q_i of the bessel functions
!
@ -74,7 +71,7 @@ use ld1inc
!
do nc=1,2
call sph_bes(ik+5,r,xc(3+nc),lam,j1(1,nc))
fact(nc)=psipsus(ik,ns)/(j1(ik,nc)*r(ik))
fact(nc)=psi_in(ik)/(j1(ik,nc)*r(ik))
do n=1,ik+5
j1(n,nc)=j1(n,nc)*r(n)*fact(nc)
enddo
@ -85,17 +82,13 @@ use ld1inc
!
do nc=1,2
! f1aep1=j1(ik+1,nc)*ff/sqr(ik+1)
! f1ae=j1(ik,nc)*(-12.0_dp+10.0_dp*ff)/sqr(ik)
! f1aem1=j1(ik-1,nc)*ff/sqr(ik-1)
! bm(nc)=(f1aep1+f1ae+f1aem1)/dx**2
bm(nc)=deriv2_7pts(j1(1,nc),ik,r(ik),dx)
bm(nc)=deriv2_7pts(j1(1,nc),ik,r(ik),dx)
enddo
xc(2)=(f2ae-bm(1))/(bm(2)-bm(1))
xc(1)=1.0_dp-xc(2)
if (iflag.eq.1) then
write(6,110) els(ns),rcutus(ns),2.0_dp*xc(5)**2
write(6,110) els_in,r(ik),2.0_dp*xc(5)**2
110 format (5x, ' Wfc-us ',a3,' rcutus=',f6.3, &
' Estimated cut-off energy= ', f8.2,' Ry')
endif
@ -103,11 +96,11 @@ use ld1inc
! define the phis function
!
do n=1,ik
phis(n,ns)=xc(1)*j1(n,1)+xc(2)*j1(n,2)
phi_out(n)=xc(1)*j1(n,1)+xc(2)*j1(n,2)
enddo
do n=ik+1,mesh
phis(n,ns)=psipsus(n,ns)
phi_out(n)=psi_in(n)
enddo
do nc=1,2

91
atomic/compute_potps.f90 Normal file
View File

@ -0,0 +1,91 @@
!
! Copyright (C) 2007 QUANTUM-ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!--------------------------------------------------------------------------
SUBROUTINE compute_potps(ik,v_in,v_out,xc)
!--------------------------------------------------------------------------
!
! This routine computes the pseudized pseudopotential by pseudizing the
! all_electron potential. In input it receives, the point
! ik where the cut is done.
!
use ld1inc
IMPLICIT NONE
REAL(DP) :: &
v_in(ndm), & ! input: the potential to pseudize
v_out(ndm),& ! output: the pseudized potential
xc(8) ! output: the coefficients of the fit
INTEGER :: &
ik ! input: the point corresponding to rc
REAL(DP) :: &
fae, & ! the value of the all-electron function
f1ae, & ! its first derivative
f2ae ! the second derivative
REAL(DP) :: &
f1aep1,f1aem1,jnor, & ! auxilairy quantities
bm(2), & ! the derivative of the bessel
fact(2), & ! factor of normalization
j1(ndm,8) ! the bessel functions
REAL(DP) :: &
deriv_7pts, deriv2_7pts, p1aep1, p1aem1
INTEGER :: &
iok, & ! if 0 there are no problem
n, & ! counter on mesh points
nc ! counter on bessel
!
! compute first and second derivative
!
fae=v_in(ik)
f1ae=deriv_7pts(v_in,ik,r(ik),dx)
f2ae=deriv2_7pts(v_in,ik,r(ik),dx)
!
! find the q_i of the bessel functions
!
CALL find_qi(f1ae/fae,xc(4),ik,0,2,0,iok)
IF (iok.NE.0) &
CALL errore('compute_potps','problems with find_qi',1)
!
! compute the functions
!
DO nc=1,2
call sph_bes(ik+1,r,xc(3+nc),0,j1(1,nc))
fact(nc)=v_in(ik)/j1(ik,nc)
DO n=1,ik+1
j1(n,nc)=j1(n,nc)*fact(nc)
ENDDO
ENDDO
!
! compute the second derivative and impose continuity of zero,
! first and second derivative
!
DO nc=1,2
p1aep1=(j1(ik+1,nc)-j1(ik,nc))/(r(ik+1)-r(ik))
p1aem1=(j1(ik,nc)-j1(ik-1,nc))/(r(ik)-r(ik-1))
bm(nc)=(p1aep1-p1aem1)*2.0_dp/(r(ik+1)-r(ik-1))
ENDDO
xc(2)=(f2ae-bm(1))/(bm(2)-bm(1))
xc(1)=1.0_dp-xc(2)
!
! define the v_out function
!
DO n=1,ik
v_out(n)=xc(1)*j1(n,1)+xc(2)*j1(n,2)
ENDDO
DO n=ik+1,mesh
v_out(n)=v_in(n)
ENDDO
RETURN
END SUBROUTINE compute_potps

74
atomic/descreening.f90
View File

@ -18,7 +18,6 @@ subroutine descreening
use ld1inc
implicit none
integer :: &
ns, & ! counter on pseudo functions
ns1, & ! counter on pseudo functions
@ -45,61 +44,23 @@ subroutine descreening
! a test configuration equal to the one used for pseudopotential
! generation is strongly suggested
!
nc=1
nwfts=nwftsc(nc)
! nc=1
! nwfts=nwftsc(nc)
do n=1,nwfts
nnts(n)=nntsc(n,nc)
llts(n)=lltsc(n,nc)
elts(n)=eltsc(n,nc)
! rcutts(n)=rcut(n)
! rcutusts(n)=rcutus(n)
jjts(n) = jjtsc(n,nc)
iswts(n)=iswtsc(n,nc)
octs(n)=octsc(n,nc)
nstoae(n)=nstoaec(n,nc)
enlts(n)=enl(nstoae(n))
new(n)=.false.
enddo
do ns=1,nwfs
do n=1,mesh
phist(n,ns)=phis(n,ns)
enddo
! nnts(n)=nntsc(n,nc)
! llts(n)=lltsc(n,nc)
! elts(n)=eltsc(n,nc)
! jjts(n) = jjtsc(n,nc)
! iswts(n)=iswtsc(n,nc)
! octs(n)=octsc(n,nc)
! nstoaets(n)=nstoaec(n,nc)
enlts(n)=enl(nstoaets(n))
! new(n)=.false.
enddo
!
! compute the pseudowavefunctions in the test configuration
!
if (pseudotype.eq.1) then
nbf=0
else
nbf=nbeta
endif
do ns=1,nwfts
if (octs(ns).gt.0.0_dp) then
is=iswts(ns)
if (pseudotype ==1) then
if ( rel < 2 .or. llts(ns) == 0 .or. &
abs(jjts(ns)-llts(ns)+0.5_dp) < 0.001_dp) then
ind=1
else if ( rel == 2 .and. llts(ns) > 0 .and. &
abs(jjts(ns)-llts(ns)-0.5_dp) < 0.001_dp) then
ind=2
endif
do n=1,mesh
vaux(n,1)=vpsloc(n)+vnl(n,llts(ns),ind)
enddo
else
do n=1,mesh
vaux(n,1)=vpsloc(n)
enddo
endif
call ascheqps(nnts(ns),llts(ns),jjts(ns),enlts(ns), &
mesh,ndm,dx,r,r2,sqr,vaux,thresh,phis(1,ns), &
betas,bmat,qq,nbf,nwfsx,lls,jjs,ikk)
! write(6,*) ns, nnts(ns),llts(ns), jjts(ns), enlts(ns)
endif
enddo
call ascheqps_drv(vpsloc, 1, thresh, .false.)
!
! descreening the D coefficients
!
@ -127,8 +88,7 @@ subroutine descreening
! descreening the local pseudopotential
!
iwork=1
call normalize
call chargeps(nwfts,llts,jjts,octs,iwork)
call chargeps(rhos,phits,nwfts,llts,jjts,octs,iwork)
call new_potential(ndm,mesh,r,r2,sqr,dx,0.0_dp,vxt,lsd,nlcc,latt,enne, &
rhoc,rhos,vh,vaux)
@ -148,14 +108,6 @@ subroutine descreening
enddo
close(20)
endif
!
! copy the phis used to construct the pseudopotential
!
do ns=1,nwfs
do n=1,mesh
phis(n,ns)=phist(n,ns)
enddo
enddo
return
end subroutine descreening

8
atomic/elsdps.f90
View File

@ -28,7 +28,7 @@ implicit none
exc_t, & ! the exchange correlation energy
vxcp(2), &
rho_tot, & !
gi(ndm), & !
gi(ndm), & !
work(nwfsx) !
real(DP),allocatable :: &
@ -139,7 +139,7 @@ implicit none
lam=llts(ns)
if (octs(ns).gt.0.0_DP) then
do n=1, mesh
f1(n,1) = f1(n,1) + phis(n,ns)**2 * octs(ns)
f1(n,1) = f1(n,1) + phits(n,ns)**2 * octs(ns)
enddo
endif
do n=1,mesh
@ -157,7 +157,7 @@ implicit none
nst=(llts(ns)+1)*2
ikl=ikk(n1)
do n=1,ikl
gi(n)=betas(n,n1)*phis(n,ns)
gi(n)=betas(n,n1)*phits(n,ns)
enddo
work(n1)=int_0_inf_dr(gi,r,r2,dx,ikl,nst)
else
@ -176,7 +176,7 @@ implicit none
ekin = int_0_inf_dr(f2,r,r2,dx,mesh,2) - epseu
do ns=1,nwfts
if (octs(ns).gt.0.0_DP) then
ekin=ekin+octs(ns)*enls(ns)
ekin=ekin+octs(ns)*enlts(ns)
endif
end do

55
atomic/gener_pseudo.f90
View File

@ -45,7 +45,9 @@ subroutine gener_pseudo
real(DP) :: &
xc(8), & ! parameters of bessel functions
psi_in(ndm), & ! the all_electron wavefunction
gi(ndm,2), & ! auxiliary to compute the integrals
occ, &
sum, db, work(nwfsx) ! work space
real(DP), allocatable :: &
@ -131,47 +133,40 @@ subroutine gener_pseudo
else
ikk(ns)=max(ik+10,ikloc+5)
endif
if (new(ns)) then
call set_psi_in(ik,lam,jjs(ns),enls(ns),psi_in)
occ=1.d0
else
psi_in(:)=psi(:,1,nwf0)
occ=ocs(ns)
endif
!
! save the all-electron function for the PAW setup
!
if (lpaw) psipaw(1:mesh,ns) = psi_in(1:mesh)
!
! compute the phi functions
!
nnode=0
call compute_phi(lam,ik,nwf0,ns,xc,1,nnode,ocs(ns))
if (nnode.ne.0) call errore('gener_pseudo','too many nodes',1)
if (tm) then
call compute_phi_tm(lam,ik,psi_in,phis(1,ns),1,xc,enls(ns),els(ns))
else
call compute_phi(lam,ik,psi_in,phis(1,ns),xc,1,occ,enls(ns),els(ns))
endif
!
! US only on the components where ikus <> ik
!
do n=1,mesh
psipsus(n,ns)=phis(n,ns)
enddo
psipsus(:,ns)=phis(:,ns)
if (ikus.ne.ik) then
call compute_phius(lam,ikus,ns,xc,1)
call compute_phius(lam,ikus,psipsus(1,ns),phis(1,ns),xc,1,els(ns))
lbes4=.true.
else
lbes4=.false.
endif
call compute_chi(lam,ik,ns,xc,lbes4)
!
! check that the chi are zero beyond ikk
nst=0
do n=1,mesh
gi(n,1)=0.0_dp
enddo
! do n=ikk(ns)+1,min(ikk(ns)+20,mesh)
do n=ikk(ns)+1,mesh
gi(n,1)=chis(n,ns)**2
enddo
do n=min(ikk(ns)+20,mesh),mesh
chis(n,ns)=0.0_dp
enddo
sum=int_0_inf_dr(gi,r,r2,dx,mesh,nst)
if (sum > 2.e-6_dp) then
write(6, '(5x,''ns='',i4,'' l='',i4, '' sum='',f15.9, &
& '' r(ikk) '',f15.9)') ns, lam, sum, r(ikk(ns))
call infomsg ('gener_pseudo ','chi too large beyond r_c', -1)
do n=ikk(ns),mesh
write(6,*) r(n),gi(n,1)
enddo
stop
if (tm) then
call compute_chi_tm(lam,ik,ikk(ns),phis(1,ns),chis(1,ns),xc,enls(ns))
else
call compute_chi(lam,ik,ikk(ns),phis(1,ns),chis(1,ns),xc,enls(ns),lbes4)
endif
enddo

3
atomic/ld1.f90
View File

@ -56,10 +56,7 @@ program ld1
!
call all_electron(.false.)
call gener_pseudo ( )
! save energies used to generate PP (run_test overwrites them)
enlts(1:nwfs) = enls(1:nwfs)
call run_test ( )
enls(1:nwfs) = enlts(1:nwfs)
call ld1_writeout ( )
!
else

15
atomic/ld1_setup.f90
View File

@ -64,6 +64,21 @@ subroutine ld1_setup
& //eltsc(n,nc)//' not found',nc)
enddo
enddo
!
! set the test configuration for descreening
!
if (iswitch == 3) then
nwfts=nwftsc(1)
do n=1,nwfts
nnts(n)=nntsc(n,1)
llts(n)=lltsc(n,1)
elts(n)=eltsc(n,1)
jjts(n) = jjtsc(n,1)
iswts(n)=iswtsc(n,1)
octs(n)=octsc(n,1)
nstoaets(n)=nstoaec(n,1)
enddo
endif
endif
!
new(:)=.false.

1
atomic/ld1inc.f90
View File

@ -136,6 +136,7 @@ module ld1inc
nnts(nwfsx), & ! the main quantum number of pseudopotential
llts(nwfsx), & ! the angular momentum of pseudopotential
iswts(nwfsx), & ! spin of the wfc. if(.not.lsd) all 1 (default)
nstoaets(nwfsx), & ! for each test wavefunction the all-electron
nwfts ! the number of pseudo wavefunctions
real(DP) :: &

81
atomic/normalize.f90
View File

@ -7,69 +7,68 @@
!
!
!---------------------------------------------------------------
subroutine normalize
subroutine normalize(phi,l,j)
!---------------------------------------------------------------
!
! normalize the US wavefunction so that <phis|S|phis>=1
!
!
!
use ld1inc
implicit none
integer :: &
real(DP) :: &
phi(ndm), & ! function to normalize
j ! total angular momentum
integer :: &
l ! orbital angular momentum
integer :: &
n,n1,n2, & ! counters on beta and mesh function
ns,nst,ikl ! counter on wavefunctions
nst,ikl ! counter on wavefunctions
real(DP) :: &
work(nwfsx), & ! auxiliary variable for becp
work1, & ! the norm
int_0_inf_dr,& ! integration function
gi(ndm) ! used to compute the integrals
gi(ndm) ! used to compute the integrals
if (pseudotype.ne.3) return
!
! if US pseudopotential compute the augmentation part
!
do ns=1,nwfts
if (octs(ns).gt.0.0_dp) then
nst=(llts(ns)+1)*2
do n1=1,nbeta
if (llts(ns).eq.lls(n1).and.abs(jjts(ns)-jjs(n1)).lt.1.e-7_dp) then
ikl=ikk(n1)
do n=1,ikl
gi(n)=betas(n,n1)*phis(n,ns)
enddo
work(n1)=int_0_inf_dr(gi,r,r2,dx,ikl,nst)
else
work(n1)=0.0_dp
endif
enddo
do n=1,mesh
gi(n)=phis(n,ns)*phis(n,ns)
enddo
work1=int_0_inf_dr(gi,r,r2,dx,mesh,nst)
!
! and adding to the charge density
!
do n1=1,nbeta
do n2=1,nbeta
work1=work1+qq(n1,n2)*work(n1)*work(n2)
enddo
enddo
if (abs(work1) < 1e-10_dp) then
call infomsg('normalize','zero norm: not a true US PP ?',-1)
work1=1.0_dp
else if (work1 < -1e-10_dp) then
call errore('normalize','negative norm?',1)
end if
work1=sqrt(work1)
do n=1,mesh
phis(n,ns)=phis(n,ns)/work1
nst=(l+1)*2
do n1=1,nbeta
if (l.eq.lls(n1).and.abs(j-jjs(n1)).lt.1.e-7_dp) then
ikl=ikk(n1)
do n=1,ikl
gi(n)=betas(n,n1)*phi(n)
enddo
work(n1)=int_0_inf_dr(gi,r,r2,dx,ikl,nst)
else
work(n1)=0.0_dp
endif
enddo
do n=1,mesh
gi(n)=phi(n)*phi(n)
enddo
work1=int_0_inf_dr(gi,r,r2,dx,mesh,nst)
!
! and adding to the charge density
!
do n1=1,nbeta
do n2=1,nbeta
work1=work1+qq(n1,n2)*work(n1)*work(n2)
enddo
enddo
if (abs(work1) < 1e-10_dp) then
call infomsg('normalize','zero norm: not a true US PP ?',-1)
work1=1.0_dp
else if (work1 < -1e-10_dp) then
call errore('normalize','negative norm?',1)
end if
work1=sqrt(work1)
do n=1,mesh
phi(n)=phi(n)/work1
enddo
return
end subroutine normalize

72
atomic/pseudovloc.f90
View File

@ -36,6 +36,7 @@ subroutine pseudovloc
xc(8), & ! the coefficients of the fit
bm(2), & ! the derivative of the bessel
vaux(ndm,2), & ! keeps the potential
psi_in(ndm), & ! auxiliary
j1(ndm,4) ! the bessel functions
real(DP) :: &
@ -44,6 +45,7 @@ subroutine pseudovloc
integer :: &
n, & ! counter on mesh points
ns, & ! auxiliary
indi,rep, & ! auxiliary
indns(0:1), & ! auxiliary
nc ! counter on bessel
@ -63,55 +65,14 @@ subroutine pseudovloc
if (mod(ik,2) == 0) ik=ik+1
if (ik <= 1 .or. ik > mesh) &
call errore('pseudovloc','wrong matching point',1)
!
! compute first and second derivative
!
fae=vpot(ik,1)
f1ae=deriv_7pts(vpot,ik,r(ik),dx)
f2ae=deriv2_7pts(vpot,ik,r(ik),dx)
!
! find the q_i of the bessel functions
!
call find_qi(f1ae/fae,xc(3),ik,0,2,0,iok)
if (iok /= 0) &
call errore('pseudovloc','problems with find_qi',1)
!
! compute the functions
!
do nc=1,2
call sph_bes(ik+1,r,xc(2+nc),0,j1(1,nc))
jnor=j1(ik,nc)
do n=1,ik+1
j1(n,nc)=j1(n,nc)*vpot(ik,1)/jnor
enddo
enddo
!
! compute the second derivative and impose continuity of zero,
! first and second derivative
!
!
! smooth the potential before ik.
!
do nc=1,2
p1aep1=(j1(ik+1,nc)-j1(ik,nc))/(r(ik+1)-r(ik))
p1aem1=(j1(ik,nc)-j1(ik-1,nc))/(r(ik)-r(ik-1))
bm(nc)=(p1aep1-p1aem1)*2.0_dp/(r(ik+1)-r(ik-1))
enddo
xc(2)=(f2ae-bm(1))/(bm(2)-bm(1))
xc(1)=1.0_dp-xc(2)
write(6, 110) rcloc,xc(4)**2
call compute_potps(ik,vpot,vpsloc,xc)
write(6, 110) r(ik),xc(5)**2
110 format (/5x, ' Local pseudo, rcloc=',f6.3, &
' Estimated cut-off energy= ', f8.2,' Ry')
!
! define the local pseudopotential
!
do n=1,ik
vpsloc(n)=xc(1)*j1(n,1)+xc(2)*j1(n,2)
enddo
do n=ik+1,mesh
vpsloc(n)=vpot(n,1)
enddo
else
!
! if a given angular momentum gives the local component this is done
@ -134,8 +95,7 @@ subroutine pseudovloc
vaux=0.0_dp
do indi=0,rep
nwf0=nstoae(nsloc+indi)
if (enls(nsloc+indi) == 0.0_dp) &
enls(nsloc+indi)=enl(nstoae(nsloc+indi))
if (enls(nsloc+indi) == 0.0_dp) enls(nsloc+indi)=enl(nwf0)
!
! compute the ik closer to r_cut
!
@ -154,7 +114,7 @@ subroutine pseudovloc
if (rel==2) then
write(6,"(/,5x,' Generating local pot.: lloc=',i1, &
&', j=',f5.2,', matching radius rcloc = ',f8.4)") &
lloc, lloc-0.5+indi, rcloc
lloc, lloc-0.5d0+indi, rcloc
else
write(6,"(/,5x,' Generating local pot.: lloc=',i1, &
&', spin=',i1,', matching radius rcloc = ',f8.4)") &
@ -164,7 +124,17 @@ subroutine pseudovloc
!
! compute the phi functions
!
call compute_phipot(lloc,ik,nwf0,indns(indi),xc)
ns=indns(indi)
if (new(ns)) then
call set_psi_in(ik,lloc,jjs(ns),enls(ns),psi_in)
else
psi_in(:)=psi(:,1,nwf0)
endif
!
! compute the phi and chi functions
!
call compute_phi_tm(lloc,ik,psi_in,phis(1,ns),0,xc,enls(ns),els(ns))
call compute_chi_tm(lloc,ik,ik+10,phis(1,ns),chis(1,ns),xc,enls(ns))
!
! set the local potential equal to the all-electron one at large r
!
@ -172,7 +142,7 @@ subroutine pseudovloc
if (r(n) > rcloc) then
vaux(n,indi+1)=vpot(n,1)
else
vaux(n,indi+1)=chis(n,indns(indi))/phis(n,indns(indi))
vaux(n,indi+1)=chis(n,ns)/phis(n,ns)
endif
enddo
enddo

89
atomic/run_pseudo.f90
View File

@ -50,13 +50,8 @@ subroutine run_pseudo
! initial estimate of the eigenvalues
!
do ns=1,nwfts
enls(ns)=enl(nstoae(ns))
enlts(ns)=enl(nstoaets(ns))
enddo
if (pseudotype.eq.1) then
nbf=0
else
nbf=nbeta
endif
!
! compute an initial estimate of the potential
!
@ -89,40 +84,11 @@ subroutine run_pseudo
! iterate to self-consistency
!
do iter=1,itmax
do ns=1,nwfts
if (octs(ns).gt.0.0_dp) then
is=iswts(ns)
if (pseudotype == 1) then
if ( rel < 2 .or. llts(ns) == 0 .or. &
abs(jjts(ns)-llts(ns)+0.5_dp) < 0.001_dp) then
ind=1
else if ( rel==2 .and. llts(ns)>0 .and. &
abs(jjts(ns)-llts(ns)-0.5_dp) < 0.001_dp) then
ind=2
else
call errore('run-pseudo','spin-orbit?!?',3)
endif
do n=1,mesh
vaux(n)=vpstot(n,is)+vnl (n,llts(ns),ind)
enddo
else
do n=1,mesh
vaux(n)=vpstot(n,is)
enddo
endif
!
call ascheqps(nnts(ns), llts(ns), jjts(ns), enls(ns), &
mesh, ndm, dx, r, r2, sqr, vaux(1), thresh, phis(1,ns), &
betas, ddd(1,1,is), qq, nbf, nwfsx, lls, jjs, ikk)
! write(6,*) 'run_pseu',ns, nwfts, nnts(ns),llts(ns), jjts(ns),enls(ns)
endif
enddo
call normalize ( )
call ascheqps_drv(vpstot, nspin, thresh, .false.)
if (.not.lpaw) then
!
call chargeps(nwfts,llts,jjts,octs,iswts)
call chargeps(rhos,phits,nwfts,llts,jjts,octs,iswts)
call new_potential(ndm,mesh,r,r2,sqr,dx,0.0_dp,vxt,lsd, &
nlcc,latt,enne,rhoc,rhos,vh,vnew)
@ -148,7 +114,7 @@ subroutine run_pseudo
else
!
call new_paw_hamiltonian (vnew, dddnew, etots, &
pawsetup, nwfts, llts, nspin, iswts, octs, phis, enls)
pawsetup, nwfts, llts, nspin, iswts, octs, phits, enlts)
do is=1,nspin
vnew(1:mesh,is)=vnew(1:mesh,is)-pawsetup%psloc(1:mesh)
enddo
@ -161,7 +127,7 @@ subroutine run_pseudo
!
endif
! write(6,*) 'iteration number',iter, eps0
! write(6,*) 'iteration number',iter, eps0
if (conv) goto 900
enddo
call infomsg('run_pseudo','convergence not achieved', -1)
@ -171,54 +137,13 @@ subroutine run_pseudo
!
900 continue
do ns=1,nwfts
is=iswts(ns)
if (octs(ns).gt.-1.0_dp) then
if (pseudotype == 1) then
if ( rel < 2 .or. llts(ns) == 0 .or. &
abs(jjts(ns)-llts(ns)+0.5_dp) < 0.001_dp) then
ind=1
else if ( rel == 2 .and. llts(ns) > 0 .and. &
abs(jjts(ns)-llts(ns)-0.5_dp) < 0.001_dp) then
ind=2
endif
do n=1,mesh
vaux(n)=vpstot(n,is)+vnl (n,llts(ns),ind)
enddo
else
do n=1,mesh
vaux(n)=vpstot(n,is)
enddo
endif
!
call ascheqps(nnts(ns), llts(ns), jjts(ns), enls(ns), &
mesh, ndm, dx, r, r2, sqr, vaux(1), thresh, phis(1,ns), &
betas, ddd(1,1,is), qq, nbf, nwfsx, lls, jjs, ikk)
call ascheqps_drv(vpstot, nspin, thresh, .true.)
! write(6,*) ns, nnts(ns),llts(ns), enls(ns)
endif
enddo
call normalize ( )
if (.not.lpaw) then
call elsdps ( )
else
call new_paw_hamiltonian (vnew, dddnew, etots, &
pawsetup, nwfts, llts, nspin, iswts, octs, phis, enls)
endif
!
! if iswitch=3 we write on the pseudopotential file the calculated
! selfconsistent wavefunctions
!
if (iswitch.eq.3) then
phits=0.0_dp
do nst=1,nwfts
if (octs(nst).ge.0.0_dp) then
phits(:,nst)=phis(:,nst)
endif
enddo
else
phits=phis
pawsetup, nwfts, llts, nspin, iswts, octs, phits, enlts)
endif
if (file_recon.ne.' ') call write_paw_recon ( )

18
atomic/run_test.f90
View File

@ -67,8 +67,8 @@ subroutine run_test
iswts(n)=iswtsc(n,nc)
octs(n)=octsc(n,nc)
if (rel==2) jjts(n)=jjtsc(n,nc)
nstoae(n)=nstoaec(n,nc)
oc(nstoae(n))=octs(n)
nstoaets(n)=nstoaec(n,nc)
oc(nstoaets(n))=octs(n)
enddo
call all_electron(.true.)
if (nc.eq.1) etot0=etot
@ -87,27 +87,27 @@ subroutine run_test
dum=0.0_dp
im=2
do ir=1,mesh-1
dum=abs(psi(ir+1,1,nstoae(n)))
if(dum.gt.abs(psi(ir,1,nstoae(n)))) im=ir+1
dum=abs(psi(ir+1,1,nstoaets(n)))
if(dum.gt.abs(psi(ir,1,nstoaets(n)))) im=ir+1
enddo
if (pseudotype.lt.3) then
rcutts(n)=r(im)*1.1_dp
rcutusts(n)=r(im)*1.1_dp
else
if (ll(nstoae(n)).eq.0) then
if (ll(nstoaets(n)).eq.0) then
rcutts(n)=r(im)*1.6_dp
rcutusts(n)=r(im)*1.7_dp
elseif (ll(nstoae(n)).eq.1) then
elseif (ll(nstoaets(n)).eq.1) then
rcutts(n)=r(im)*1.6_dp
rcutusts(n)=r(im)*1.7_dp
if (el(nstoae(n)).eq.'2P') then
if (el(nstoaets(n)).eq.'2P') then
rcutts(n)=r(im)*1.7_dp
rcutusts(n)=r(im)*1.8_dp
endif
elseif (ll(nstoae(n)).eq.2) then
elseif (ll(nstoaets(n)).eq.2) then
rcutts(n)=r(im)*2.0_dp
rcutusts(n)=r(im)*2.2_dp
if (el(nstoae(n)).eq.'3D') then
if (el(nstoaets(n)).eq.'3D') then
rcutts(n)=r(im)*2.5_dp
rcutusts(n)=r(im)*3.0_dp
endif

34
atomic/set_psi_in.f90 Normal file
View File

@ -0,0 +1,34 @@
SUBROUTINE set_psi_in(ik,l,j,e,psi_out)
!
! This subroutine calculates the all electron wavefunction psi at the
! input energy e
!
USE ld1inc
IMPLICIT NONE
INTEGER :: l, ik ! input: angular momentum and index of the cut-off radius
REAL(DP) :: e, j ! input: energy and total angular momentum
REAL(DP) :: psi_out(ndm) ! output: the function psi.
REAL(DP) :: psi_dir(ndm,2) ! auxiliary function.
REAL(DP) :: ze2, jnor
integer :: n
IF (rel == 1) THEN
CALL lschps(3,zed,exp(dx),dx,mesh,mesh,mesh,1,l,e,psi_out,r,vpot)
ELSEIF (rel == 2) THEN
CALL dir_outward(ndm,mesh,l,j,e,dx,psi_dir,r,rab,vpot)
psi_out(:)=psi_dir(:,1)
ELSE
ze2=-zed*2.0_dp
CALL intref(l,e,mesh,dx,r,r2,sqr,vpot,ze2,psi_out)
ENDIF
!
! fix arbitrarily the norm at the cut-off radius equal to 0.5
!
jnor=psi_out(ik)
DO n=1,mesh
psi_out(n)=psi_out(n)*0.5_dp/jnor
ENDDO
RETURN
END SUBROUTINE set_psi_in

21
atomic/start_potps.f90
View File

@ -13,8 +13,6 @@ subroutine start_potps
! This routine computes an initial estimate of the screening
! potential
!
!
!
use ld1inc
implicit none
@ -38,7 +36,7 @@ subroutine start_potps
do ns=1,nwfts
if (octs(ns).gt.0.0_dp) then
lam=llts(ns)
nwf0=nstoae(ns)
nwf0=nstoaets(ns)
!
! compute the ik closer to r_cut
!
@ -55,7 +53,13 @@ subroutine start_potps
!
! compute the phi functions
!
call compute_phi(lam,ik,nwf0,ns,xc,0,nnode,octs(ns))
if (tm) then
call compute_phi_tm(lam,ik,psi(1,1,nwf0),phis(1,ns),0,xc, &
enls(ns),els(ns))
else
call compute_phi(lam,ik,psi(1,1,nwf0),phis(1,ns),xc,0,octs(ns), &
enl(nwf0),' ')
endif
if (pseudotype.eq.3) then
!
! US only on the components where ikus <> ik
@ -63,20 +67,21 @@ subroutine start_potps
do n=1,mesh
psipsus(n,ns)=phis(n,ns)
enddo
if (ikus.ne.ik) call compute_phius(lam,ikus,ns,xc,0)
if (ikus.ne.ik) call compute_phius(lam,ikus,psipsus(1,ns), &
phis(1,ns),xc,0,' ')
endif
call normalize(phis(1,ns),llts(ns),jjts(ns))
endif
enddo
call normalize
call chargeps(nwfts,llts,jjts,octs,iswts)
call chargeps(rhos,phis,nwfts,llts,jjts,octs,iswts)
call new_potential(ndm,mesh,r,r2,sqr,dx,0.0_dp,vxt,lsd,nlcc, &
latt,enne,rhoc,rhos,vh,vnew)
do is=1,nspin
do n=1,mesh
vpstot(n,is)=vpsloc(n)+vnew(n,is)
! if (is.eq.1.and.nspin.eq.1.and.n.lt.420.and.n.gt.410) &
! if (is.eq.1) &
! write(6,'(3f25.16)') r(n), rhos(n,1),vpstot(n,1)
! if (is.eq.2.and.nspin.eq.2) &
! write(6,'(3f25.16)') 2.0_dp*rhos(n,1),vpstot(n,1),vpstot(n,2)

16
atomic/write_resultsps.f90
View File

@ -37,12 +37,12 @@ subroutine write_resultsps
' e PS (Ry) De AE-PS (Ry) ')
write(6,1100) &
(nnts(n),llts(n),elts(n),iswts(n),octs(n), &
enl(nstoae(n)),enls(n), &
enl(nstoae(n))-enls(n), n=1,nwfts)
enl(nstoaets(n)),enlts(n), &
enl(nstoaets(n))-enlts(n), n=1,nwfts)
write(13,1100) &
(nnts(n),llts(n),elts(n),iswts(n),octs(n), &
enl(nstoae(n)),enls(n), &
enl(nstoae(n))-enls(n), n=1,nwfts)
enl(nstoaets(n)),enlts(n), &
enl(nstoaets(n))-enlts(n), n=1,nwfts)
1100 format(4x,2i2,5x,a2,i4,'(',f5.2,')',f15.5,f15.5,f15.5)
else
write(6,1001)
@ -50,12 +50,12 @@ subroutine write_resultsps
' e PS (Ry) De AE-PS (Ry) ')
write(6,1101) &
(nnts(n),llts(n),jjts(n),elts(n),iswts(n),octs(n), &
enl(nstoae(n)),enls(n), &
enl(nstoae(n))-enls(n), n=1,nwfts)
enl(nstoaets(n)),enlts(n), &
enl(nstoaets(n))-enlts(n), n=1,nwfts)
write(13,1101) &
(nnts(n),llts(n),jjts(n),elts(n),iswts(n),octs(n), &
enl(nstoae(n)),enls(n), &
enl(nstoae(n))-enls(n), n=1,nwfts)
enl(nstoaets(n)),enlts(n), &
enl(nstoaets(n))-enlts(n), n=1,nwfts)
1101 format(4x,2i2,f4.1,1x,a2,i4,'(',f5.2,')',f15.5,f15.5,f15.5)
endif
write(6,1200) eps0,iter