in the release notes of version 4.0.15 there is mentioned a "serious bug", which was fixed. But I have to stick to version 3.1.12 (because of the interface to DGrid 4.6) and DFT+U. So my question is: How "serious" is this bug, i.e. how reliable are my results, since I use DFT+U for atoms with symmetric equivalents?
Further information: I want to calculate the ground-state energy of LuPO4 (xenotime) in the primitive cell. DFT+U is applied for Lu.
Thanks a lot for a reply!
Best regards,
Roger
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
Dear Roger,
the problem was located already in your version it dated back to 2.3.16. I found it by comparing the magnetic energies of different configuation which exhibited a way to hight dependence on the U value. You can easily fix the thing by editing symdmat.f90
see my code here
! Copyright (C) 2007 F. Bultmark, F. Cricchio, L. Nordstrom and J. K. Dewhurst.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.
subroutine symdmat(lmax,ld,dmat)
use modmain
implicit none
! arguments
integer, intent(in) :: lmax
integer, intent(in) :: ld
complex(8), intent(inout) :: dmat(ld,nspinor,ld,nspinor,natmtot)
! local variables
integer is,ia,ja,ias,jas
integer isym,lspl,lspn,lmmax
real(8) t1
! automatic arrays
logical done(natmmax)
! allocatable arrays
complex(8), allocatable :: dm(:,:,:,:,:)
lmmax=(lmax+1)*2
! allocate local arrays
allocate(dm(ld,nspinor,ld,nspinor,natmmax))
t1=1.d0/dble(nsymcrys)
do is=1,nspecies
! make copy of the density matrices
do ia=1,natoms(is)
ias=idxas(ia,is)
dm(1:lmmax,:,1:lmmax,:,ia)=dmat(1:lmmax,:,1:lmmax,:,ias)
end do
done(:)=.false.
do ia=1,natoms(is)
if (done(ia)) cycle
ias=idxas(ia,is)
dmat(:,:,:,:,ias)=0.d0
do isym=1,nsymcrys
lspl=lsplsymc(isym)
lspn=lspnsymc(isym)
! equivalent atom index (symmetry rotates atom ja into atom ia)
ja=ieqatom(ia,is,isym)
call rotdmat(symlatc(:,:,lspl),symlatc(:,:,lspn),lmax,nspinor,ld, &
dm(:,:,:,:,ja),dmat(:,:,:,:,ias))
! end loop over crystal symmetries
end do
! normalise
dmat(:,:,:,:,ias)=t1dmat(:,:,:,:,ias)
done(ia)=.true.
! rotate into equivalent atoms
do isym=1,nsymcrys
ja=ieqatom(ia,is,isym)
if (.not.done(ja)) then
jas=idxas(ja,is)
! inverse symmetry (which rotates atom ia into atom ja)
lspl=isymlat(lsplsymc(isym))
lspn=isymlat(lspnsymc(isym))
dmat(:,:,:,:,jas)=0.d0
call rotdmat(symlatc(:,:,lspl),symlatc(:,:,lspn),lmax,nspinor,ld, &
dmat(:,:,:,:,ias),dmat(:,:,:,:,jas))
done(ja)=.true.
end if
end do
! end loop over atoms and species
end do
end do
deallocate(dm)
return
end subroutine
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
As Michael pointed out, there was a bug in the rotation of the density matrix used for finding the DFT+U potential (in fact, Michael was the one who discovered it).
You should be able to simply cut-and-paste either Michael's code or the same routine from one of the later released versions into 3.1.12, and it should work OK.
Regards,
Kay.
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
Dear Elk community,
in the release notes of version 4.0.15 there is mentioned a "serious bug", which was fixed. But I have to stick to version 3.1.12 (because of the interface to DGrid 4.6) and DFT+U. So my question is: How "serious" is this bug, i.e. how reliable are my results, since I use DFT+U for atoms with symmetric equivalents?
Further information: I want to calculate the ground-state energy of LuPO4 (xenotime) in the primitive cell. DFT+U is applied for Lu.
Thanks a lot for a reply!
Best regards,
Roger
Dear Roger,
the problem was located already in your version it dated back to 2.3.16. I found it by comparing the magnetic energies of different configuation which exhibited a way to hight dependence on the U value. You can easily fix the thing by editing symdmat.f90
see my code here
! Copyright (C) 2007 F. Bultmark, F. Cricchio, L. Nordstrom and J. K. Dewhurst.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.
subroutine symdmat(lmax,ld,dmat)
use modmain
implicit none
! arguments
integer, intent(in) :: lmax
integer, intent(in) :: ld
complex(8), intent(inout) :: dmat(ld,nspinor,ld,nspinor,natmtot)
! local variables
integer is,ia,ja,ias,jas
integer isym,lspl,lspn,lmmax
real(8) t1
! automatic arrays
logical done(natmmax)
! allocatable arrays
complex(8), allocatable :: dm(:,:,:,:,:)
lmmax=(lmax+1)*2
! allocate local arrays
allocate(dm(ld,nspinor,ld,nspinor,natmmax))
t1=1.d0/dble(nsymcrys)
do is=1,nspecies
! make copy of the density matrices
do ia=1,natoms(is)
ias=idxas(ia,is)
dm(1:lmmax,:,1:lmmax,:,ia)=dmat(1:lmmax,:,1:lmmax,:,ias)
end do
done(:)=.false.
do ia=1,natoms(is)
if (done(ia)) cycle
ias=idxas(ia,is)
dmat(:,:,:,:,ias)=0.d0
do isym=1,nsymcrys
lspl=lsplsymc(isym)
lspn=lspnsymc(isym)
! equivalent atom index (symmetry rotates atom ja into atom ia)
ja=ieqatom(ia,is,isym)
call rotdmat(symlatc(:,:,lspl),symlatc(:,:,lspn),lmax,nspinor,ld, &
dm(:,:,:,:,ja),dmat(:,:,:,:,ias))
! end loop over crystal symmetries
end do
! normalise
dmat(:,:,:,:,ias)=t1dmat(:,:,:,:,ias)
done(ia)=.true.
! rotate into equivalent atoms
do isym=1,nsymcrys
ja=ieqatom(ia,is,isym)
if (.not.done(ja)) then
jas=idxas(ja,is)
! inverse symmetry (which rotates atom ia into atom ja)
lspl=isymlat(lsplsymc(isym))
lspn=isymlat(lspnsymc(isym))
dmat(:,:,:,:,jas)=0.d0
call rotdmat(symlatc(:,:,lspl),symlatc(:,:,lspn),lmax,nspinor,ld, &
dmat(:,:,:,:,ias),dmat(:,:,:,:,jas))
done(ja)=.true.
end if
end do
! end loop over atoms and species
end do
end do
deallocate(dm)
return
end subroutine
Dear Roger,
As Michael pointed out, there was a bug in the rotation of the density matrix used for finding the DFT+U potential (in fact, Michael was the one who discovered it).
You should be able to simply cut-and-paste either Michael's code or the same routine from one of the later released versions into 3.1.12, and it should work OK.
Regards,
Kay.