Menu

DFT+U in version 3.1.12

Elk Users
2017-03-21
2017-03-21
  • Roger Kloditz

    Roger Kloditz - 2017-03-21

    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

     
  • mfechner

    mfechner - 2017-03-24

    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)=t1
    dmat(:,:,:,:,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

     
  • J. K. Dewhurst

    J. K. Dewhurst - 2017-03-28

    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.

     

Log in to post a comment.