Commit efbd6bd3 authored by Daniel Wortmann's avatar Daniel Wortmann
Browse files

Started to implement banddos switch from inp.xml

parent f63498cb
......@@ -48,6 +48,9 @@ MODULE m_types_banddos
LOGICAL :: l_slab=.false.
LOGICAL,ALLOCATABLE :: dos_atom(:) ! for each atom (not type) switch on DOS
INTEGER,ALLOCATABLE :: dos_typelist(:) !list of types for which DOS is calculated
INTEGER,ALLOCATABLE :: dos_atomlist(:) !list of atoms for which DOS is calculated
INTEGER,ALLOCATABLE :: map_atomtype(:) !map an atomtype to corresponding entry in DOS
!INTEGER :: ndir =0
!INTEGER :: orbCompAtom=0
!INTEGER :: jDOSAtom=0
......@@ -99,6 +102,10 @@ CONTAINS
CALL mpi_bc(this%locy(2),rank,mpi_comm)
CALL mpi_bc(this%starcoeff,rank,mpi_comm)
CALL mpi_bc(this%izlay,rank,mpi_comm)
CALL mpi_bc(this%dos_atom,rank,mpi_comm)
CALL mpi_bc(this%dos_atomlist,rank,mpi_comm)
CALL mpi_bc(this%dos_typelist,rank,mpi_comm)
CALL mpi_bc(this%map_atomtype,rank,mpi_comm)
END SUBROUTINE mpi_bc_banddos
SUBROUTINE read_xml_banddos(this,xml)
......@@ -107,8 +114,10 @@ CONTAINS
TYPE(t_xml),INTENT(INOUT)::xml
CHARACTER(len=300) :: xPathA, xPathB,str
INTEGER::numberNodes,iType,i,na,n
LOGICAL::l_orbcomp,l_jDOS,all_atoms
INTEGER::numberNodes,iType,i,na,n,n_dos_atom,n_dos_type
LOGICAL::l_orbcomp,l_jDOS,all_atoms,dos_atom_found
integer,allocatable:: dos_atomlist(:),dos_typelist(:),neq(:)
this%band = evaluateFirstBoolOnly(xml%GetAttributeValue('/fleurInput/output/@band'))
this%dos = evaluateFirstBoolOnly(xml%GetAttributeValue('/fleurInput/output/@dos'))
!this%l_slab = evaluateFirstBoolOnly(xml%GetAttributeValue('/fleurInput/output/@slab'))
......@@ -131,14 +140,18 @@ CONTAINS
END IF
allocate(this%dos_atom(xml%get_nat()))
this%dos_atom=all_atoms
allocate(neq(xml%get_ntype()))
this%dos_atom=(all_atoms.and.(this%dos.or.this%band))
na = 0
IF (xml%versionNumber > 31) then
DO iType = 1, xml%GetNumberOfNodes('/fleurInput/atomGroups/atomGroup')
WRITE(xPathA,*) '/fleurInput/atomGroups/atomGroup[',iType,']'
DO i = 1, xml%GetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/relPos')
neq(itype)= xml%GetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/relPos')
DO i = 1, neq(itype)
na = na + 1
WRITE(xPathB,*) TRIM(ADJUSTL(xPathA))//'/relPos[',i,']'
if (xml%GetNumberOfNodes(TRIM(ADJUSTL(xPathB))//'/@banddos')==1) &
this%dos_atom(na) = evaluateFirstBoolOnly(xml%GetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@banddos'))
ENDDO
ENDDO
......@@ -176,6 +189,39 @@ CONTAINS
allocate(this%izlay(0,2))
END IF
!Create a list of all atoms and all types for which the DOS is calculated
ALLOCATE(dos_atomlist(xml%get_nat()))
allocate(dos_typelist(xml%get_ntype()))
na=0
n_dos_atom=0
n_dos_type=0
DO itype=1,xml%get_ntype()
dos_atom_found=.false.
DO n=1,neq(itype)
na=na+1
if (this%dos_atom(na)) then
dos_atom_found=.true.
n_dos_atom=n_dos_atom+1
dos_atomlist(n_dos_atom)=na
endif
enddo
if (dos_atom_found)then
n_dos_type=n_dos_type+1
dos_typelist(n_dos_type)=iType
ENDIF
ENDDO
this%dos_atomlist=dos_atomlist(:n_dos_atom)
this%dos_typelist=dos_typelist(:n_dos_type)
!create map
allocate(this%map_atomtype(size(neq)))
this%map_atomtype=0
DO itype=1,size(neq)
DO n=1,size(dos_typelist)
if (dos_typelist(n)==itype) this%map_atomtype(itype)=n
ENDDO
ENDDO
END SUBROUTINE read_xml_banddos
......
......@@ -464,9 +464,8 @@
<xsd:extension base="String3DVecType">
<xsd:attribute default=" " name="label" type="xsd:string" use="optional"/>
<xsd:attribute default="F" name="wannier" type="FleurBool" use="optional"/>
<xsd:attribute default="F" name="orbcomp" type="FleurBool" use="optional"/>
<xsd:attribute default="F" name="jDOS" type="FleurBool" use="optional"/>
<xsd:attribute default="F" name="banddos" type="FleurBool" use="optional"/>
<xsd:attribute name="banddos" type="FleurBool" use="optional"/>
</xsd:extension>
</xsd:simpleContent>
</xsd:complexType>
......
MODULE m_mcdinit
CONTAINS
SUBROUTINE mcd_init(atoms,input,vr,g,f,mcd,itype,jspin)
SUBROUTINE mcd_init(atoms,banddos,input,vr,g,f,mcd,itype,jspin)
!-----------------------------------------------------------------------
!
......@@ -22,6 +22,7 @@ CONTAINS
TYPE(t_input),INTENT(IN) :: input
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_banddos),INTENT(IN) :: banddos
TYPE(t_mcd),INTENT(INOUT) :: mcd
INTEGER, PARAMETER :: l_max = 3
......@@ -36,7 +37,7 @@ CONTAINS
! Locals ...
INTEGER kap,mue,iri,l,ispin,i,icore,korb,nst,n_core,ierr
INTEGER kap,mue,iri,l,ispin,i,icore,korb,nst,n_core,ierr,n_dos
REAL c,t2,e,fj,fl,fn ,d,ms,rn
INTEGER kappa(maxval(atoms%econf%num_states)),nprnc(maxval(atoms%econf%num_states)),l_core(maxval(atoms%econf%num_states))
REAL vrd(atoms%msh),occ(maxval(atoms%econf%num_states),2),a(atoms%msh),b(atoms%msh),j_core(maxval(atoms%econf%num_states)),e_mcd1(maxval(atoms%econf%num_states))
......@@ -46,13 +47,15 @@ CONTAINS
!-----------------------------------------------------------------------
if (.not.any(banddos%dos_atom(sum(atoms%neq(:itype-1))+1:sum(atoms%neq(:itype))))) return
c = c_light(1.0)
ALLOCATE ( gc(atoms%jri(itype),atoms%econf(itype)%num_core_states,input%jspins) )
ALLOCATE ( fc(atoms%jri(itype),atoms%econf(itype)%num_core_states,input%jspins) )
! core setup
mcd%ncore(itype) = 0
n_dos=banddos%map_atomtype(itype)
mcd%ncore(n_dos) = 0
CALL atoms%econf(itype)%get_core(nst,nprnc,kappa,occ)
DO ispin = jspin, jspin
......@@ -105,7 +108,7 @@ CONTAINS
ALLOCATE (dgv(atoms%jri(itype),0:l_max,input%jspins,2) )
ALLOCATE ( fv(atoms%jri(itype),0:l_max,input%jspins,2) )
DO i = 1, 2
DO iri = 3*(itype-1)+1 , 3*(itype-1)+3
DO iri = 3*(n_dos-1)+1 , 3*(n_dos-1)+3
DO l = 1, (l_max+1)**2
DO icore = 1, maxval(atoms%econf%num_states)
mcd%m_mcd(icore,l,iri,i) = CMPLX(0.0,0.0)
......@@ -139,15 +142,15 @@ CONTAINS
DO i = 1, 2
! write(*,*) j_core(icore),l_core(icore),l_max,ms
CALL nabla(itype,icore,atoms%jri(itype),atoms%dx(itype),maxval(atoms%econf%num_states),atoms%ntype,&
CALL nabla(n_dos,icore,atoms%jri(itype),atoms%dx(itype),maxval(atoms%econf%num_states),atoms%ntype,&
j_core(icore),l_core(icore),l_max,ms,atoms%rmsh(:,itype),gc(:,icore,ispin),&
gv(:,0:,ispin,i),dgv(:,0:,ispin,i), mcd%m_mcd(:,:,:,i) )
ENDDO
DO i = 1, 2*icore*l_core(icore)
mcd%ncore(itype) = mcd%ncore(itype) + 1
IF (mcd%ncore(itype)>maxval(atoms%econf%num_states)) CALL juDFT_error("maxval(atoms%econf%num_states) too small" ,calledby ="mcd_init")
mcd%e_mcd(itype,ispin,mcd%ncore(itype)) = e_mcd1(icore)
mcd%ncore(n_dos) = mcd%ncore(n_dos) + 1
IF (mcd%ncore(n_dos)>maxval(atoms%econf%num_states)) CALL juDFT_error("maxval(atoms%econf%num_states) too small" ,calledby ="mcd_init")
mcd%e_mcd(n_dos,ispin,mcd%ncore(n_dos)) = e_mcd1(icore)
ENDDO
ENDDO
ENDDO
......
MODULE m_nabla
use m_juDFT
CONTAINS
SUBROUTINE nabla(
......@@ -22,13 +22,13 @@
! psi(r) ... core wavefunction
! phi(r,l) ... valence wavefunction
! dphi(r,l) ... radial derivative of valence wavefunction
!
!
!----------------------------------------------------------------
USE m_constants
USE m_clebsch
USE m_intgr, ONLY : intgr3
IMPLICIT NONE
IMPLICIT NONE
INTEGER, INTENT(IN) :: ispecies, number_of_j1, grid_size
INTEGER, INTENT(IN) :: l1, lmax, nstd, ntypd
......@@ -58,80 +58,80 @@
ALLOCATE ( f(grid_size),STAT=alloc_error )
IF (alloc_error /= 0) CALL juDFT_error("Couldn't allocate f",
+ calledby ="nabla")
lmn1 = 2 * (number_of_j1 - 1) * l1
mu = -j1
DO WHILE (mu <= j1)
lmn1 = lmn1 + 1
m1 = INT(mu - ms)
m1 = INT(mu - ms)
lmn2 = 0
DO l2 = 0, lmax
DO m2 = -l2, l2
lmn2 = lmn2 + 1
IF(l1 == l2 + 1)THEN
IF(l1 == l2 + 1)THEN
total_result = 0.00
result = 0.00
result1 = 0.00
!
!
! (l+1)/srqt[(2l+1)(2l+3)] < phi_core | (d phi_valence / dr ) >
!
f(:) = psi(:) * dphi(:,l2) ! assumed to be already multiplied with * ri(:) * ri(:)
CALL intgr3(f,ri,delta_x,grid_size,result)
result = result * (l2 + 1.00) /
+ sqrt((2.00*l2 +1.00)*(2.*l2+3.00))
!
! - l(l+1)/srqt[(2l+1)(2l+3)] < phi_core | (1/r) phi_valence >
!
f(:) = psi(:) * phi(:,l2) ! assumed to be already multiplied with 1G* ri(:)
CALL intgr3(f,ri,delta_x,grid_size,result1)
result1 = - result1 * l2 * (l2 + 1.0) /
+ sqrt((2.0 * l2 + 1.00) * (2. * l2 + 3.00))
result1 = - result1 * l2 * (l2 + 1.0) /
+ sqrt((2.0 * l2 + 1.00) * (2. * l2 + 3.00))
!
! Sum up and decorate with Clebsch-Gordon coefficients
!
result = result + result1
total_result = result*clebsch(real(l1),spin,mu-ms,ms,j1,mu)
index = (ispecies - 1) * 3 + 1
psi_phi(lmn1,lmn2,index)= cgc(l2,1,l1,m2,1,m1) * ! left polarization
+ total_result / cgc(l2,1,l1,0,0,0)
index = index + 1
psi_phi(lmn1,lmn2,index)= cgc(l2,1,l1,m2,-1,m1) * ! right polarization
+ total_result / cgc(l2,1,l1,0,0,0)
index = index + 1
psi_phi(lmn1,lmn2,index)= cgc(l2,1,l1,m2,0,m1) * ! z-polarization
+ total_result / cgc(l2,1,l1,0,0,0)
ELSEIF(l1== l2-1)THEN
result = cnst_zero
ELSEIF(l1== l2-1)THEN
result = cnst_zero
result1 = cnst_zero
!
! l/srqt[(2l-1)(2l+1)] < phi_core | (d phi_valence / dr ) >
!
f(:) = psi(:)* dphi(:,l2) * ri(:) * ri(:)
CALL intgr3(f,ri,delta_x,grid_size,result)
result = result * l2 / sqrt((2.0*l2 - 1.0)*(2.0*l2 + 1.0))
result = result * l2 / sqrt((2.0*l2 - 1.0)*(2.0*l2 + 1.0))
!
! l(l+1)/srqt[(2l-1)(2l+1)] < phi_core | (1/r) phi_valence >
!
f(:) = psi(:)* phi(:,l2) * ri(:)
CALL intgr3(f,ri,delta_x,grid_size,result1)
result1 = result1 * l2 * (l2 + 1.0) /
+ sqrt((2.00 * l2 - 1.00) * (2.00 * l2 + 1.00))
+ sqrt((2.00 * l2 - 1.00) * (2.00 * l2 + 1.00))
!
! Sum up and decorate with Clebsch-Gordon coefficients
!
result = result + result1
total_result= result*clebsch(real(l1),spin,mu-ms,ms,j1,mu)
index = (ispecies - 1) * 3 + 1
psi_phi(lmn1,lmn2,index) = cgc(l2,1,l1,m2,1,m1) * ! left polarization
+ total_result / cgc(l2,1,l1,0,0,0)
......@@ -141,7 +141,7 @@
index = index + 1
psi_phi(lmn1,lmn2,index) = cgc(l2,1,l1,m2,0,m1) * ! z-polarization
+ total_result / cgc(l2,1,l1,0,0,0)
ENDIF
ENDIF
ENDDO
ENDDO
mu = mu + 1.00
......@@ -153,37 +153,37 @@
FUNCTION cgc(l1,l2,l3,m1,m2,m3)
IMPLICIT NONE
IMPLICIT NONE
INTEGER :: l1, l2, l3, m1, m2, m3
REAL :: two_l1p1, two_l1p2, l1pm3, l1pm3p1, l1mm3p1, l1mm3, cgc
REAL :: two_l1p1, two_l1p2, l1pm3, l1pm3p1, l1mm3p1, l1mm3, cgc
IF (m3 /= m1 + m2) THEN
cgc = 0.0
RETURN
END IF
END IF
! gb m3 = m1 + m2
two_l1p1 = 2 * l1 + 1
two_l1p2 = 2 * l1 + 2
l1pm3 = l1 + m3
l1pm3p1 = l1 + m3 + 1
l1mm3p1 = l1 - m3 + 1
l1mm3 = l1 - m3
cgc = 0.0
l1mm3 = l1 - m3
cgc = 0.0
IF (l3 == l1 + 1) THEN
IF (m2 == 1) then
cgc = sqrt( (l1pm3 * l1pm3p1) / (two_l1p1 * two_l1p2))
cgc = sqrt( (l1pm3 * l1pm3p1) / (two_l1p1 * two_l1p2))
ELSEIF (m2 == 0) THEN
cgc = sqrt( (l1mm3p1 * l1pm3p1) / (two_l1p1 * (l1 + 1)))
ELSEIF (m2 == -1) THEN
cgc = sqrt( (l1mm3 * l1mm3p1) / (two_l1p1 * two_l1p2))
cgc = sqrt( (l1mm3 * l1mm3p1) / (two_l1p1 * two_l1p2))
END IF
ELSE IF(l3 == l1 -1) THEN
IF (m2 == 1) then
cgc = sqrt( (l1mm3 * l1mm3p1) / (2.0 * l1 * two_l1p1))
cgc = sqrt( (l1mm3 * l1mm3p1) / (2.0 * l1 * two_l1p1))
ELSEIF (m2 == 0) THEN
cgc = -sqrt( (l1mm3 * l1pm3) / (l1 * (two_l1p1)))
cgc = -sqrt( (l1mm3 * l1pm3) / (l1 * (two_l1p1)))
ELSEIF (m2 == -1) THEN
cgc = sqrt( (l1pm3p1 * l1pm3) / (2.0 * l1 * two_l1p1))
cgc = sqrt( (l1pm3p1 * l1pm3) / (2.0 * l1 * two_l1p1))
END IF
END IF
END FUNCTION cgc
......
......@@ -95,26 +95,26 @@ SUBROUTINE dos_init(thisDOS,input,atoms,kpts,banddos,eig)
INTEGER :: ntype,l,i,ind
character :: spdfg(0:4)=["s","p","d","f","g"]
thisDOS%name_of_dos="Local"
thisDOS%neq=atoms%neq
thisDOS%neq=atoms%neq(banddos%dos_typelist)
thisDOS%eig=eig
ALLOCATE(thisDOS%jsym(input%neig,kpts%nkpt,input%jspins))
ALLOCATE(thisDOS%ksym(input%neig,kpts%nkpt,input%jspins))
ALLOCATE(thisDOS%qis(input%neig,kpts%nkpt,input%jspins))
ALLOCATE(thisDOS%qal(0:3,atoms%ntype,input%neig,kpts%nkpt,input%jspins))
ALLOCATE(thisDOS%qal(0:3,size(banddos%dos_typelist),input%neig,kpts%nkpt,input%jspins))
thisDOS%jsym = 0
thisDOS%ksym = 0
thisDOS%qis = 0.0
thisDOS%qal = 0.0
allocate(thisDOS%weight_names(2+4*atoms%ntype))
allocate(thisDOS%weight_names(2+4*size(banddos%dos_typelist)))
thisDOS%weight_names(1)="Total"
thisDOS%weight_names(2)="INT"
ind=2
DO ntype=1,atoms%ntype
DO ntype=1,size(banddos%dos_typelist)
DO l=0,3
ind=ind+1
write(thisDOS%weight_names(ind),"(a,i0,a)") "MT:",ntype,spdfg(l)
write(thisDOS%weight_names(ind),"(a,i0,a)") "MT:",banddos%dos_typelist(ntype),spdfg(l)
ENDDO
ENDDO
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment