Commit c2f8cb88 authored by Daniel Wortmann's avatar Daniel Wortmann

Removed some old files

parent 91df57eb
MODULE m_hsmt_nonsph_GPU
#define CPP_BLOCKSIZE 4096
! USE m_juDFT
!$ USE omp_lib
!TODO:
! Check what can be done in l_noco=.true. case in terms of use of zgemm or aa_block
! Check what happens in case of CPP_INVERSION -> real matrix a
IMPLICIT NONE
CONTAINS
SUBROUTINE hsmt_nonsph_GPU(DIMENSION,atoms,sym,SUB_COMM, n_size,n_rank,input,isp,nintsp,&
hlpmsize,noco,l_socfirst, lapw, cell,tlmplm, fj,gj,gk,vk,oneD,l_real,aa_r,aa_c)
#include"cpp_double.h"
USE m_constants, ONLY : tpi_const
USE m_ylm
USE m_hsmt_spinor
USE m_hsmt_hlptomat
USE m_types
#if defined(_OPENACC)
! USE nvtx
USE cublas
#endif
IMPLICIT NONE
TYPE(t_dimension),INTENT(IN):: DIMENSION
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_input),INTENT(IN) :: input
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_lapw),INTENT(INOUT) :: lapw !lapw%nv_tot is updated
! ..
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: nintsp,isp
INTEGER, INTENT (IN) :: SUB_COMM,n_size,n_rank
INTEGER, INTENT (IN) :: hlpmsize
LOGICAL, INTENT (IN) :: l_socfirst
! ..
! .. Array Arguments ..
TYPE(t_tlmplm),INTENT(IN)::tlmplm
REAL, INTENT(IN) :: fj(:,0:,:,:),gj(:,0:,:,:)
REAL,INTENT(IN) :: gk(:,:,:),vk(:,:,:)
!-odim
!+odim
LOGICAL, INTENT(IN) :: l_real
REAL, ALLOCATABLE, INTENT (INOUT) :: aa_r(:)!(matsize)
COMPLEX,ALLOCATABLE, INTENT (INOUT) :: aa_c(:)
COMPLEX,PARAMETER :: one=CMPLX(1.0,0.0),zero=CMPLX(0.0,0.0)
! ..
! .. Local Scalars ..
INTEGER :: i,iii,ii,ij,im,in,k,ki,kj,l,ll1,lm,lmp,lp,jd,m
INTEGER :: mp,n,na,nn,np,kjmax,iintsp,jintsp
INTEGER :: kiStart, kiiSTart, kc
INTEGER :: nc ,kii,spin2,ab_dim,lnonsphd,bsize,bsize2,kb
REAL :: th,invsfct
COMPLEX :: term,chi11,chi21,chi22,chihlp
! ..
! .. Local Arrays ..
COMPLEX,ALLOCATABLE :: aa_block(:,:)
COMPLEX,ALLOCATABLE :: dtd(:,:),dtu(:,:),utd(:,:),utu(:,:)
REAL :: bmrot(3,3),gkrot(DIMENSION%nvd,3),vmult(3),v(3)
COMPLEX:: ylm( (atoms%lmaxd+1)**2 ),chi(2,2)
! ..
#if defined(_OPENACC)
COMPLEX,PINNED,ALLOCATABLE :: a(:,:,:),b(:,:,:)
#else
COMPLEX, ALLOCATABLE :: a(:,:,:),b(:,:,:)
#endif
COMPLEX, ALLOCATABLE :: ax(:,:),bx(:,:)
COMPLEX, ALLOCATABLE :: c_ph(:,:)
COMPLEX,ALLOCATABLE :: aahlp(:),aa_tmphlp(:)
INTEGER :: cublas_stream,istat
<<<<<<< HEAD
!call nvtxStartRange("hsmt_nonsph",1)
=======
#if defined(_OPENACC)
! call nvtxStartRange("hsmt_nonsph",1)
#endif
>>>>>>> hsmt_simple
lnonsphd=MAXVAL(atoms%lnonsph)*(MAXVAL(atoms%lnonsph)+2)
ALLOCATE(dtd(0:lnonsphd,0:lnonsphd),utd(0:lnonsphd,0:lnonsphd),dtu(0:lnonsphd,0:lnonsphd),utu(0:lnonsphd,0:lnonsphd))
!Decide how to distribute the work
IF ( noco%l_noco .AND. (.NOT. noco%l_ss) ) THEN
ALLOCATE ( aahlp(hlpmsize),aa_tmphlp(hlpmsize) )
ELSE
ALLOCATE ( aahlp(1),aa_tmphlp(1) )
END IF
ALLOCATE(aa_block(CPP_BLOCKSIZE,MAXVAL(lapw%nv)))
ab_dim=1
IF (noco%l_ss) ab_dim=2
ALLOCATE(a(DIMENSION%nvd,0:DIMENSION%lmd,ab_dim),b(DIMENSION%nvd,0:DIMENSION%lmd,ab_dim))
ALLOCATE(ax(DIMENSION%nvd,0:DIMENSION%lmd),bx(DIMENSION%nvd,0:DIMENSION%lmd))
ALLOCATE(c_ph(DIMENSION%nvd,ab_dim))
!$acc data copy(aa_r,aa_c) copyin(tlmplm, tlmplm%tdd, tlmplm%tdu, tlmplm%tud,tlmplm%tuu, tlmplm%ind, atoms,atoms%lnonsph,lapw,lapw%nv,noco ) create(utu,utd,dtu,dtd,ax,bx,a,b,aa_block,aahlp)
#if defined(_OPENACC)
! call nvtxStartRange("hsmt_nonsph_outer_data",2)
#endif
ntyploop: DO n=1,atoms%ntype
IF (noco%l_noco) THEN
IF (.NOT.noco%l_ss) aahlp=CMPLX(0.0,0.0)
IF (.NOT.noco%l_ss) aa_tmphlp=CMPLX(0.0,0.0)
CALL hsmt_spinor(isp,n, noco,input, chi, chi11, chi21, chi22)
ENDIF
DO nn = 1,atoms%neq(n)
a=0.0
b=0.0
na = SUM(atoms%neq(:n-1))+nn
IF (atoms%lnonsph(n)<0) CYCLE ntyploop
IF ((atoms%invsat(na)==0) .OR. (atoms%invsat(na)==1)) THEN
IF (atoms%invsat(na)==0) invsfct = 1
IF (atoms%invsat(na)==1) invsfct = 2
np = sym%invtab(atoms%ngopr(na))
IF (oneD%odi%d1) np = oneD%ods%ngopr(na)
!Using double buffering create_ab could be overlapped with following GPU work
#if defined(_OPENACC)
! call nvtxStartRange("create_ab",6)
#endif
!---> loop over interstitial spins
DO iintsp = 1,nintsp
IF (noco%l_constr.OR.l_socfirst) THEN
spin2=isp
ELSE
spin2=iintsp
ENDIF
!---> set up phase factors
DO k = 1,lapw%nv(iintsp)
th= DOT_PRODUCT((/lapw%k1(k,iintsp),lapw%k2(k,iintsp),lapw%k3(k,iintsp)/)+(iintsp-1.5)*noco%qss,atoms%taual(:,na))
c_ph(k,iintsp) = CMPLX(COS(tpi_const*th),-SIN(tpi_const*th))
END DO
IF (np==1) THEN
gkrot( 1:lapw%nv(iintsp),:) = gk( 1:lapw%nv(iintsp),:,iintsp)
ELSE
IF (oneD%odi%d1) THEN
bmrot=MATMUL(oneD%ods%mrot(:,:,np),cell%bmat)
ELSE
bmrot=MATMUL(1.*sym%mrot(:,:,np),cell%bmat)
END IF
DO k = 1,lapw%nv(iintsp)
!--> apply the rotation that brings this atom into the
!--> representative (this is the definition of ngopr(na)
!--> and transform to cartesian coordinates
v(:) = vk(k,:,iintsp)
gkrot(k,:) = MATMUL(TRANSPOSE(bmrot),v)
END DO
END IF
DO k = 1,lapw%nv(iintsp)
!--> generate spherical harmonics
vmult(:) = gkrot(k,:)
CALL ylm4(atoms%lnonsph(n),vmult,ylm)
!--> synthesize the complex conjugates of a and b
DO l = 0,atoms%lnonsph(n)
ll1 = l* (l+1)
DO m = -l,l
term = c_ph(k,iintsp)*ylm(ll1+m+1)
a(k,ll1+m,iintsp) = fj(k,l,n,spin2)*term
b(k,ll1+m,iintsp) = gj(k,l,n,spin2)*term
END DO
END DO
ENDDO !k-loop
!---> end loop over interstitial spin
ENDDO
!$acc update device(a,b)
#if defined(_OPENACC)
! call nvtxEndRange
#endif
!---> loops over the interstitial spin
DO iintsp = 1,nintsp
#if defined(_OPENACC)
! call nvtxStartRange("hsmt_nonsph_DO_jintsp",3)
#endif
DO jintsp = 1,iintsp
jd = 1 ; IF (noco%l_noco) jd = isp
!---> loop over l',m'
!$acc kernels
utu=0.0;utd=0.0;dtu=0.0;dtd=0.0
!!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(lp,mp,lmp,l,m,lm,in,utu,dtu,utd,dtd,im,k) &
!!$OMP SHARED(tlmplm,invsfct,lnonsph,nv,jintsp,jd,n)
!$acc loop collapse(2)
DO lmp=0,atoms%lnonsph(n)*(atoms%lnonsph(n)+2)
! lp=FLOOR(SQRT(1.0*lmp))
! mp=lmp-lp*(lp+1)
! IF (lp>atoms%lnonsph(n).OR.ABS(mp)>lp) STOP "BUG"
!---> loop over l,m
DO lm = 0,atoms%lnonsph(n)*(atoms%lnonsph(n)+2)
! DO l = 0,atoms%lnonsph(n)
! DO m = -l,l
! lm = l* (l+1) + m
in = tlmplm%ind(lmp,lm,n,jd)
IF (in/=-9999) THEN
IF (in>=0) THEN
utu(lm,lmp) =CONJG(tlmplm%tuu(in,n,jd))*invsfct
dtu(lm,lmp) =CONJG(tlmplm%tdu(in,n,jd))*invsfct
utd(lm,lmp) =CONJG(tlmplm%tud(in,n,jd))*invsfct
dtd(lm,lmp) =CONJG(tlmplm%tdd(in,n,jd))*invsfct
ELSE
im = -in
utu(lm,lmp) =tlmplm%tuu(im,n,jd)*invsfct
dtu(lm,lmp) =tlmplm%tud(im,n,jd)*invsfct
utd(lm,lmp) =tlmplm%tdu(im,n,jd)*invsfct
dtd(lm,lmp) =tlmplm%tdd(im,n,jd)*invsfct
END IF
!---> update ax, bx
END IF
! END DO
END DO
ENDDO
!$acc end loop
!$acc end kernels
!!$OMP END PARALLEL DO
lmp=atoms%lnonsph(n)*(atoms%lnonsph(n)+2)
!ax(:nv(jintsp),0:lmp)=(matmul(a(:nv(jintsp),0:lmp,jintsp),utu(0:lmp,0:lmp))+matmul(b(:nv(jintsp),0:lmp,jintsp),utd(0:lmp,0:lmp)))
!bx(:nv(jintsp),0:lmp)=(matmul(a(:nv(jintsp),0:lmp,jintsp),dtu(0:lmp,0:lmp))+matmul(b(:nv(jintsp),0:lmp,jintsp),dtd(0:lmp,0:lmp)))
#if defined(_OPENACC)
! call nvtxStartRange("hsmt_nonsph_zgemm_1",4)
#endif
!$acc host_data use_device(a,b,utu,utd,dtu,dtd,ax,bx )
CALL zgemm("N","N",lapw%nv(jintsp),lmp+1,lmp+1,one,a(1,0,jintsp),SIZE(a,1),utu(0,0),SIZE(utu,1),zero,ax,SIZE(ax,1))
CALL zgemm("N","N",lapw%nv(jintsp),lmp+1,lmp+1,one,b(1,0,jintsp),SIZE(a,1),utd(0,0),SIZE(utu,1),one,ax,SIZE(ax,1))
CALL zgemm("N","N",lapw%nv(jintsp),lmp+1,lmp+1,one,a(1,0,jintsp),SIZE(a,1),dtu(0,0),SIZE(utu,1),zero,bx,SIZE(ax,1))
CALL zgemm("N","N",lapw%nv(jintsp),lmp+1,lmp+1,one,b(1,0,jintsp),SIZE(a,1),dtd(0,0),SIZE(utu,1),one,bx,SIZE(ax,1))
!$acc end host_data
#if defined(_OPENACC)
! call nvtxEndRange
#endif
!
!---> update hamiltonian and overlap matrices
nc = 0
IF ( noco%l_noco .AND. (n_size>1) ) THEN
lapw%nv_tot = lapw%nv(1) + lapw%nv(2)
ELSE
lapw%nv_tot = lapw%nv(iintsp)
ENDIF
kii=n_rank
#if defined(_OPENACC)
! call nvtxStartRange("hsmt_nonsph_while_kii",5)
#endif
DO WHILE(kii<lapw%nv_tot)
!DO kii = n_rank, nv_tot-1, n_size
ki = MOD(kii,lapw%nv(iintsp)) + 1
bsize=MIN(SIZE(aa_block,1),(lapw%nv(iintsp)-ki)/n_size+1) !Either use maximal blocksize or number of rows left to calculate
IF (bsize<1) EXIT !nothing more to do here
bsize2=bsize*n_size
bsize2=min(bsize2,lapw%nv(iintsp)-ki+1)
!aa_block(:bsize,:ki+bsize2-1)=matmul(a(ki:ki+bsize2-1:n_size,0:lmp,iintsp),conjg(transpose(ax(:ki+bsize2-1,0:lmp))))+ &
! matmul(b(ki:ki+bsize2-1:n_size,0:lmp,iintsp),conjg(transpose(bx(:ki+bsize2-1,0:lmp))))
IF (n_size==1) THEN !Make this a special case to avoid copy-in of a array
!$acc host_data use_device( a,b,ax,bx,aa_block )
call zgemm("N","C",bsize,ki+bsize2-1,lmp+1,one,a(ki,0,iintsp),SIZE(a,1),ax(1,0),SIZE(ax,1),zero,aa_block,SIZE(aa_block,1))
call zgemm("N","C",bsize,ki+bsize2-1,lmp+1,one,b(ki,0,iintsp),SIZE(a,1),bx(1,0),SIZE(ax,1),one ,aa_block,SIZE(aa_block,1))
!$acc end host_data
ELSE
stop "Not implemented"
!$acc host_data use_device( a,b,ax,bx,aa_block )
!CALL zgemm("N","C",bsize,ki+bsize2-1,lmp+1,one,a(ki:ki+bsize2-1:n_size,0:lmp,iintsp), &
! SIZE(a(ki:ki+bsize2-1:n_size,0:lmp,iintsp),1),ax(1,0),SIZE(ax,1),zero,aa_block,SIZE(aa_block,1))
!CALL zgemm("N","C",bsize,ki+bsize2-1,lmp+1,one,b(ki:ki+bsize2-1:n_size,0:lmp,iintsp), &
! SIZE(a(ki:ki+bsize2-1:n_size,0:lmp,iintsp),1),bx(1,0),SIZE(ax,1),one,aa_block,SIZE(aa_block,1))
!$acc end host_data
ENDIF
IF ( noco%l_noco .AND. (.NOT. noco%l_ss) ) THEN
!$acc kernels
!$acc loop independent
DO kb=1,bsize
nc = 1+kii/n_size
ii = nc*(nc-1)/2*n_size-(nc-1)*(n_size-n_rank-1)
IF ( (n_size==1).OR.(kii+1<=lapw%nv(1)) ) THEN !
stop "not implemented"
!comments below must be reactivated
!aahlp(ii+1:ii+ki) = aahlp(ii+1:ii+ki)+MATMUL(CONJG(ax(:ki,:lmp)),a(ki,:,iintsp))+MATMUL(CONJG(bx(:ki,:lmp)),b(ki,:lmp,iintsp))
ELSE ! components for <2||2> block unused
!aa_tmphlp(:ki) = MATMUL(CONJG(ax(:ki,:lmp)),a(ki,:lmp,iintsp))+MATMUL(CONJG(bx(:ki,:DIMENSION%lmd)),b(ki,:lmp,iintsp))
!---> spin-down spin-down part
ij = ii + lapw%nv(1)
aa_c(ij+1:ij+ki)=aa_c(ij+1:ij+ki)+chi22*aa_tmphlp(:ki)
!---> spin-down spin-up part, lower triangle
ij = ii
aa_c(ij+1:ij+ki)=aa_c(ij+1:ij+ki)+chi21*aa_tmphlp(:ki)
ENDIF
!-||
ki=ki+n_size
kii=kii+n_size
ENDDO
!$acc end loop
!$acc end kernels
ELSEIF ( noco%l_noco .AND. noco%l_ss ) THEN
!$acc kernels
!$acc loop independent
DO kb=1,bsize
IF ( iintsp==1 .AND. jintsp==1 ) THEN
!---> spin-up spin-up part
kjmax = ki
chihlp = chi11
ii = (ki-1)*(ki)/2
ELSEIF ( iintsp==2 .AND. jintsp==2 ) THEN
!---> spin-down spin-down part
kjmax = ki
chihlp = chi22
ii = (lapw%nv(1)+atoms%nlotot+ki-1)*(lapw%nv(1)+atoms%nlotot+ki)/2+&
lapw%nv(1)+atoms%nlotot
ELSE
!---> spin-down spin-up part
kjmax = lapw%nv(1)
chihlp = chi21
ii = (lapw%nv(1)+atoms%nlotot+ki-1)*(lapw%nv(1)+atoms%nlotot+ki)/2
ENDIF
STOP "NOT implemented"
!next line must be reactivated
!aa_c(ii+1:ii+kjmax) = aa_c(ii+1:ii+kjmax) + chihlp*&
! (MATMUL(CONJG(ax(:kjmax,:lmp)),a(ki,:,iintsp))+MATMUL(CONJG(bx(:kjmax,:lmp)),b(ki,:lmp,iintsp)))
ki=ki+n_size
kii=kii+n_size
ENDDO
!$acc end loop
!$acc end kernels
ELSE
kiStart = ki
kiiSTart = kii
!$acc kernels
!$acc loop independent worker
DO kb=1,bsize
ki = kiStart + (kb-1)*n_size
kii = kiiStart + (kb-1)*n_size
nc = 1+kii/n_size
ii = nc*(nc-1)/2*n_size- (nc-1)*(n_size-n_rank-1)
if (l_real) THEN
! aa_r(ii+1:ii+ki) = aa_r(ii+1:ii+ki) + aa_block(kb,:ki)
!$acc loop independent vector
DO kc = 1, ki
aa_r(ii+kc) = aa_r(ii+kc) + aa_block(kb,kc)
END DO
ELSE
! aa_c(ii+1:ii+ki) = aa_c(ii+1:ii+ki) + aa_block(kb,:ki)
!$acc loop independent vector
DO kc = 1, ki
aa_c(ii+kc) = aa_c(ii+kc) + aa_block(kb,kc)
END DO
endif
!print*,ii,ki,kb
! IF (.not.apw(l)) THEN
!aa(ii+1:ii+ki) = aa(ii+1:ii+ki) + b(ki,lmp,iintsp)*bx(:ki)
! ENDIF
! ki=ki+n_size
! kii=kii+n_size
ENDDO
!$acc end loop
!$acc end kernels
ki = kiStart+bsize*n_size
kii = kiiSTart+bsize*n_size
ENDIF
!---> end loop over ki
END DO
#if defined(_OPENACC)
! call nvtxEndRange
#endif
!---> end loops over interstitial spin
ENDDO
#if defined(_OPENACC)
! call nvtxEndRange
#endif
ENDDO
ENDIF ! atoms%invsat(na) = 0 or 1
!---> end loop over equivalent atoms
END DO
IF ( noco%l_noco .AND. (.NOT. noco%l_ss) )THEN
!$acc update host(aahlp,aa_c)
CALL hsmt_hlptomat(atoms%nlotot,lapw%nv,sub_comm,chi11,chi21,chi22,aahlp,aa_c)
!$acc update device(aa_c)
ENDIF
!---> end loop over atom types (ntype)
ENDDO ntyploop
#if defined(_OPENACC)
! call nvtxEndRange
#endif
!$acc end data
#if defined(_OPENACC)
! call nvtxEndRange
#endif
RETURN
END SUBROUTINE hsmt_nonsph_GPU
END MODULE m_hsmt_nonsph_GPU
MODULE m_enpara
use m_juDFT
USE m_constants
! *************************************************************
! Module containing three subroutines
! r_enpara: read enpara file
! w_enpara: write enpara file
! mix_enpara: calculate new energy parameters
! *************************************************************
CONTAINS
SUBROUTINE w_enpara( &
& atoms,jspin,film,&
& enpara,&
& id)
!
! write enpara-file
!
USE m_types
IMPLICIT NONE
INTEGER, INTENT (IN) :: jspin,id
LOGICAL,INTENT(IN) :: film
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_enpara),INTENT(IN) :: enpara
INTEGER n,l,lo
LOGICAL l_opened
INQUIRE(unit=40,OPENED=l_opened)
if (.not.l_opened) return
WRITE (40,FMT=8035) jspin,enpara%enmix(jspin)
WRITE (40,FMT=8036)
8035 FORMAT (5x,'energy parameters for spin ',i1,' mix=',f10.6)
8036 FORMAT (t6,'atom',t15,'s',t24,'p',t33,'d',t42,'f')
DO n = 1,atoms%ntype
WRITE (6,FMT=8040) n, (enpara%el0(l,n,jspin),l=0,3),&
& (enpara%lchange(l,n,jspin),l=0,3),enpara%skiplo(n,jspin)
WRITE (id,FMT=8040) n, (enpara%el0(l,n,jspin),l=0,3),&
& (enpara%lchange(l,n,jspin),l=0,3),enpara%skiplo(n,jspin)
WRITE (40,FMT=8040) n, (enpara%el0(l,n,jspin),l=0,3),&
& (enpara%lchange(l,n,jspin),l=0,3),enpara%skiplo(n,jspin)
!---> energy parameters for the local orbitals
IF (atoms%nlo(n).GE.1) THEN
WRITE (6,FMT=8039) (enpara%ello0(lo,n,jspin),lo=1,atoms%nlo(n))
WRITE (6,FMT=8038) (enpara%llochg(lo,n,jspin),lo=1,atoms%nlo(n))
WRITE (id,FMT=8039) (enpara%ello0(lo,n,jspin),lo=1,atoms%nlo(n))
WRITE (id,FMT=8038) (enpara%llochg(lo,n,jspin),lo=1,atoms%nlo(n))
WRITE (40,FMT=8039) (enpara%ello0(lo,n,jspin),lo=1,atoms%nlo(n))
WRITE (40,FMT=8038) (enpara%llochg(lo,n,jspin),lo=1,atoms%nlo(n))
END IF
ENDDO
8038 FORMAT (' --> change ',60(l1,8x))
8039 FORMAT (' --> lo ',60f9.5)
8040 FORMAT (' -->',i3,1x,4f9.5,' change: ',4l1,' skiplo: ',i3)
IF (film) THEN
WRITE (40,FMT=8050) enpara%evac(1,jspin),enpara%lchg_v(1,jspin),enpara%evac(2,jspin)
WRITE (6,FMT=8050) enpara%evac(1,jspin),enpara%lchg_v(1,jspin),enpara%evac(2,jspin)
WRITE (id,FMT=8050) enpara%evac(1,jspin),enpara%lchg_v(1,jspin),enpara%evac(2,jspin)
8050 FORMAT (' vacuum parameter=',f9.5,' change: ',l1,&
& ' second vacuum=',f9.5)
ENDIF
RETURN
END SUBROUTINE w_enpara
!
!------------------------------------------------------------------
SUBROUTINE r_enpara(&
& atoms,input,jsp,&
& enpara)
!------------------------------------------------------------------
USE m_types
IMPLICIT NONE
INTEGER, INTENT (IN) :: jsp
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_input),INTENT(IN) :: input
TYPE(t_enpara),INTENT(INOUT):: enpara
INTEGER n,l,lo,skip_t,io_err
enpara%lchange(:,:,jsp) =.false.
enpara%el0(:,:,jsp) =0.0
enpara%ello0(:,:,jsp) =0.0
!--> first line contains mixing parameter!
enpara%enmix(jsp) = 0.0
READ (40,FMT ='(48x,f10.6)',iostat=io_err) enpara%enmix(jsp)
IF (io_err /= 0) THEN
!use defaults
enpara%lchange(:,:,jsp)=.false.
enpara%llochg(:,:,jsp)=.false.
! enpara%evac0(:,jsp) = eVac0Default_const
enpara%skiplo(:,jsp) = 0
enpara%enmix(jsp) = 0.0
enpara%lchg_v(:,jsp)=.false.
enpara%el0(:,:,jsp)=-99999.
CALL default_enpara(jsp,atoms,enpara)
WRITE (6,FMT=8001) jsp
WRITE (6,FMT = 8000)
DO n = 1,atoms%ntype
WRITE (6,FMT=8140) n,(enpara%el0(l,n,jsp),l=0,3),&
& (enpara%lchange(l,n,jsp),l=0,3),enpara%skiplo(n,jsp)
IF (atoms%nlo(n)>=1) THEN
WRITE (6,FMT = 8139) (enpara%ello0(lo,n,jsp),lo=1,atoms%nlo(n))
WRITE (6,FMT = 8138) (enpara%llochg(lo,n,jsp),lo = 1,atoms%nlo(n))
ENDIF
ENDDO
RETURN
ENDIF
READ (40,*) ! skip next line
IF (enpara%enmix(jsp).EQ.0.0) enpara%enmix(jsp) = 1.0
WRITE (6,FMT=8001) jsp
WRITE (6,FMT=8000)
skip_t = 0
DO n = 1,atoms%ntype
READ (40,FMT=8040,END=200) (enpara%el0(l,n,jsp),l=0,3),&
& (enpara%lchange(l,n,jsp),l=0,3),enpara%skiplo(n,jsp)
WRITE (6,FMT=8140) n,(enpara%el0(l,n,jsp),l=0,3),&
& (enpara%lchange(l,n,jsp),l=0,3),enpara%skiplo(n,jsp)
!
!---> energy parameters for the local orbitals
!
IF (atoms%nlo(n).GE.1) THEN
skip_t = skip_t + enpara%skiplo(n,jsp) * atoms%neq(n)
READ (40,FMT=8039,END=200) (enpara%ello0(lo,n,jsp),lo=1,atoms%nlo(n))
READ (40,FMT=8038,END=200) (enpara%llochg(lo,n,jsp),lo=1,atoms%nlo(n))
WRITE (6,FMT=8139) (enpara%ello0(lo,n,jsp),lo=1,atoms%nlo(n))
WRITE (6,FMT=8138) (enpara%llochg(lo,n,jsp),lo=1,atoms%nlo(n))
ELSEIF (enpara%skiplo(n,jsp).GT.0) THEN
WRITE (6,*) "for atom",n," no LO's were specified"
WRITE (6,*) 'but skiplo was set to',enpara%skiplo
CALL juDFT_error("No LO's but skiplo",calledby ="enpara",&
& hint="If no LO's are set skiplo must be 0 in enpara")
END IF
!
!---> set the energy parameters with l>3 to the value of l=3
!
DO l = 4,atoms%lmax(n)
enpara%el0(l,n,jsp) = enpara%el0(3,n,jsp)
ENDDO
ENDDO ! atoms%ntype
IF (input%film) THEN
enpara%lchg_v = .true.
READ (40,FMT=8050,END=200) enpara%evac0(1,jsp),enpara%lchg_v(1,jsp),enpara%evac0(2,jsp)
WRITE (6,FMT=8150) enpara%evac0(1,jsp),enpara%lchg_v(1,jsp),enpara%evac0(2,jsp)
ENDIF
IF (atoms%nlod.GE.1) THEN
WRITE (6,FMT=8090) jsp,skip_t
WRITE (6,FMT=8091)
END IF
! input formats
8038 FORMAT (14x,60(l1,8x))
8039 FORMAT (8x,60f9.5)
8040 FORMAT (8x,4f9.5,9x,4l1,9x,i3)
8050 FORMAT (19x,f9.5,9x,l1,15x,f9.5)
! output formats
8138 FORMAT (' --> change ',60(l1,8x))
8139 FORMAT (' --> lo ',60f9.5)
8140 FORMAT (' -->',i3,1x,4f9.5,' change: ',4l1,' skiplo: ',i3)
8150 FORMAT (' vacuum parameter=',f9.5,' change: ',l1,&
& ' second vacuum=',f9.5)
8001 FORMAT ('READING enpara for spin: ',i1)
8000 FORMAT (/,' energy parameters:',/,t10,'s',t20,&
& 'p',t30,'d',t37,'higher l - - -')
8090 FORMAT ('Spin: ',i1,' -- ',i3,'eigenvalues')
8091 FORMAT ('will be skipped for energyparameter computation')
RETURN
200 WRITE (6,*) 'the end of the file enpara has been reached while'
WRITE (6,*) 'reading the energy-parameters.'
WRITE (6,*) 'possible reason: energy parameters have not been'
WRITE (6,*) 'specified for all atom types.'
WRITE (6,FMT='(a,i4)')&
& 'the actual number of atom-types is: ntype=',atoms%ntype
CALL juDFT_error&
& ("unexpected end of file enpara reached while reading"&
& ,calledby ="enpara")
END SUBROUTINE r_enpara
!
!------------------------------------------------------------------
SUBROUTINE mix_enpara(&
& jsp,atoms,vacuum,&
& obsolete,input,enpara,&
& vr,vz,pvac,svac,&
& ener,sqal,enerlo,sqlo)
!------------------------------------------------------------------
USE m_types
IMPLICIT NONE
INTEGER,INTENT(IN) :: jsp
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_obsolete),INTENT(IN) :: obsolete
TYPE(t_input),INTENT(IN) :: input
TYPE(t_enpara),INTENT(INOUT) :: enpara
REAL, INTENT(IN) :: vr(atoms%jmtd,atoms%ntype)
REAL, INTENT(IN) :: ener(0:3,atoms%ntype),sqal(0:3,atoms%ntype)
REAL, INTENT(IN) :: enerlo(atoms%nlod,atoms%ntype),sqlo(atoms%nlod,atoms%ntype)
REAL, INTENT(IN) :: pvac(2),svac(2),vz(vacuum%nmzd,2)
INTEGER ityp,j,l,lo
REAl vbar,maxdist
INTEGER same(atoms%nlod)
LOGICAL l_int_enpara
l_int_enpara=all(enpara%el0==int(enpara%el0)) !test only enpara of lapw, no lo
maxdist=0.0
DO ityp = 1,atoms%ntype
! look for LO's energy parameters equal to the LAPW (and previous LO) ones
same = 0
DO lo = 1,atoms%nlo(ityp)
IF(enpara%el0(atoms%llo(lo,ityp),ityp,jsp).eq.enpara%ello0(lo,ityp,jsp)) same(lo)=-1
DO l = 1,lo-1
IF(atoms%llo(l,ityp).ne.atoms%llo(lo,ityp)) cycle
IF(enpara%ello0(l,ityp,jsp).eq.enpara%ello0(lo,ityp,jsp).and.same(lo).eq.0)&
& same(lo)=l
ENDDO
ENDDO
!
!---> change energy parameters
!
IF ( obsolete%lepr.EQ.1) THEN
j = atoms%jri(ityp) - (log(4.0)/atoms%dx(ityp)+1.51)
vbar = vr(j,ityp)/( atoms%rmt(ityp)*exp(atoms%dx(ityp)*(j-atoms%jri(ityp))) )
ELSE
vbar = 0.0
END IF
DO l = 0,3
IF ( enpara%lchange(l,ityp,jsp) ) THEN
write(6,*) 'Type:',ityp,' l:',l
write(6,FMT=777) enpara%el0(l,ityp,jsp),&
& (ener(l,ityp)/sqal(l,ityp) - vbar),&
& abs(enpara%el0(l,ityp,jsp)-(ener(l,ityp)/sqal(l,ityp) - vbar))
maxdist=max(maxdist,&
& abs(enpara%el0(l,ityp,jsp)-(ener(l,ityp)/sqal(l,ityp) - vbar)))
enpara%el0(l,ityp,jsp) =(1.0-enpara%enmix(jsp))*enpara%el0(l,ityp,jsp) + &
& enpara%enmix(jsp)*(ener(l,ityp)/sqal(l,ityp) - vbar)
ENDIF
ENDDO
DO l = 4, atoms%lmaxd
IF ( enpara%lchange(3,ityp,jsp) ) THEN
enpara%el0(l,ityp,jsp) = enpara%el0(3,ityp,jsp)
ENDIF
ENDDO
!
!---> determine and change local orbital energy parameters
!
DO lo = 1,atoms%nlo(ityp)
IF (atoms