Commit d55a9801 authored by Stefan Rost's avatar Stefan Rost

Merge branch 'develop' of iffgit.fz-juelich.de:fleur/fleur into develop

parents f62f823c 755deeac
......@@ -175,4 +175,4 @@ gfortran-coverage:
url: https://fleur.iffgit.fz-juelich.de/fleur/coverage_html
only:
- web
\ No newline at end of file
......@@ -10,3 +10,6 @@
[submodule "external/wannier90"]
path = external/wannier90
url = https://github.com/wannier-developers/wannier90.git
[submodule "edsolver-library"]
path = external/edsolver-library
url = gitlab@iffgit.fz-juelich.de:janssen/edsolver-library.git
<div align="center">
<img src="https://www.flapw.de/site/img/fleur.gif" width="220">
Welcome to the source code of FLEUR
===================================
=====================
[Report bug](https://iffgit.fz-juelich.de/fleur/fleur/issues/new?template=Bug.md)
.
[Request feature](https://iffgit.fz-juelich.de/fleur/fleur/issues/new?template=FeatureRequest.md&labels=feature)
[Homepage and Documentation](https://www.flapw.de)
</div>
## Table of contents
- [Using the FLEUR git repository](#fleur-git-repository)
- [Dealing with Bugs and problems](#bugs-and-problems)
- [Installation of FLEUR](#installation-of-FLEUR)
- [Contributing](#contributing)
Please note that the documentation of the
code can be found at the [FLEUR Homepage]
(http://www.flapw.de/).
## FLEUR git repository
For further instructions on Installation/Usage,
please check the [FLEUR Homepage]
(http://www.flapw.de/).
The primary git-repository of FLEUR can be found on the [iffgit-Server at FZ-Jülich](https://iffgit.fz-juelich.de/fleur/fleur/).
You can clone the repository by using
```
git clone https://iffgit.fz-juelich.de/fleur/fleur.git
```
If you are a FLEUR developer you should use
```
git clone gitlab@iffgit.fz-juelich.de:fleur/fleur.git
```
to be able to push changes back to the server. If you are not a developer yet but want to contribute, please contact [Gregor](g.michalicek@fz-juelich.de) or [Daniel](d.wortmann@fz-juelich.de).
## Bugs in FLEUR
Please note, that the default branch you will see after cloning the repository is the 'develop' branch. In general you might find
the following branches on the server.
* develop: this is the default branch with the most up-to-date version of FLEUR. Small changes and developments should be committed
directly into this branch. When doing so you should try to keep the code operational. It should still compile and the test should run.
* release: this branch collects the official releases. You cannot commit to this branch and bugfixes should be handled as [described below](#bugs-and-problems).
* stable: this branch contains snapshots of the development branch considered "stable".
In addition several other branches can/will be present corresponding to features currently under development. If you start your own larger development
it can be advisable to create your own branch. In this case you should try to follow changes in 'develop' by frequently merging 'develop' into your branch
and you should create a merge request with 'develop' as soon as you are finished or reached some usefull state in your development.
## Bugs and Problems
You might experience bugs in FLEUR :-).
If you find a bug you should:
A) Report this bug by generating an Issue. Please describe in
A) [Report this bug by generating an Issue](https://iffgit.fz-juelich.de/fleur/fleur/issues/new?template=Bug.md). Please describe in
detail the relevant input and what happens. You should consider using
the bug-template for your issue as this will help you providing us with
the relevant information.
......@@ -35,102 +71,15 @@ If you are fixing a bug in a release-version, please:
* merge your fix into the develop branch: ```git merge bugfix_SOME_NAME_HERE```
## Installation of FLEUR
To install and use FLEUR, please check the [Documentation](https://www.flapw.de).
## Structure of the source code
The source of FLEUR is organized in several
subdirectories. Some of them collect code
specific for particular features, others code
relevant for crutial steps in the code or simply
code that is usually executed together.
Here a short description of the directories:
* main: contains the main-program and several core subroutines
* init: stuff for the initialization (called from fleur_init in main)
* vgen: potential generation (called from vgen in main)
* eigen: setup of the eigenproblem
* diagonalization: various methods to diagonalize the Hamiltonian
* cdn: general code for the generation of charge
* cdn_mt: charge generation in MT-spheres
* force: code related to the evaluation of forces
* mix: charge/potential mixing routines
* ldau: routines needed in case of LDA+U calculations
* inpgen: code for the input generator (seperate executable inpgen)
* fermi: determination of the fermi-level
* eigen_secvar: second variational solution of the Hamiltonian
* eigen_soc: Spin-orbit related code
* core: Core states
* dos: Code for Density of states, bandstructures
* orbdep: Code for quantities depending on orbitals
* optional: code that is used in special cases like inital charge generation
* wannier: wannier related code
* xc-pot: various exchange-correlation potential routines
* mpi: code for parallel execution
* io: subroutines doing IO
* juDFT: timing, error handling, etc
* eels: code for electron-energy loss spectroscopy
* hybrid: code for hybrid functionals
* math: code providing math functionality
* include: c-type include files
* global: code used everywhere (here you find types.F90 with the data-types)
* cmake: definitions used by cmake
If you modify FLEUR please do so in the develop branch by running
'git checkout -t origin/develop'
after cloning the git repository. For larger changes you might want to
create your own branch.
## Compiling Fleur on Jureca
To compile Fleur on Jureca the following modules need to be loaded:
```bash
module load intel-para CMake HDF5 libxml2/.2.9.7 ELPA/2017.11.001-hybrid
```
which is then configured with
```bash
./configure.sh JURECA_INTEL
```
or
```bash
./configure.sh -external libxc JURECA_INTEL
```
if libXC should be supported aswell.
## Compilling Fleur on Jureca Booster
To compile on the booster you need to first switch to the KNL environment
## Contributing
```bash
module purge
ml Architecture/KNL
module load intel-para CMake HDF5 libxml2/.2.9.7 ELPA/2017.11.001-hybrid
```
Then you can compile as discrebed above.
## Running Fleur on Jureca Booster
To run on the Boosters you need to switch the architecture in each job script aswell
```bash
#!/bin/bash -x
#SBATCH --nodes=1
##SBATCH --ntasks=10
##SBATCH --ntasks-per-node=10
#SBATCH --cpus-per-task=68
#SBATCH --output=mpi-%j.out
#SBATCH --error=mpi-%j.err
#SBATCH --time=1:00:00
#SBATCH --partition=develbooster
#SBATCH --gres=mem96
##SBATCH --mail-user=your.name@fz-juelich.de
##SBATCH --mail-type=END
export OMP_NUM_THREADS=${SLURM_CPUS_PER_TASK}
module purge
ml Architecture/KNL
module load intel-para CMake HDF5 libxml2/.2.9.7 ELPA/2017.11.001-hybrid
/work/ias-1/s.rost/fleur_booster/fleur/build/fleur
```
FLEUR is an open source code under the [MIT-license](https://iffgit.fz-juelich.de/fleur/fleur/blob/develop/LICENSE).
## Developing Fleur
Your are very welcome to contribute to its development. If you need help or access to the git repository,
please contact [Gregor](g.michalicek@fz-juelich.de) or [Daniel](d.wortmann@fz-juelich.de).
We agreed to use a unified indentation-width of 3.
Hint: [vim](http://vim.wikia.com/wiki/Converting_tabs_to_spaces) [emacs](https://www.gnu.org/software/emacs/manual/html_node/efaq/Changing-the-length-of-a-Tab.html)
\ No newline at end of file
Please also use the [Wiki](https://iffgit.fz-juelich.de/fleur/fleur/wikis/home) for sharing information relevant for developers.
\ No newline at end of file
......@@ -13,6 +13,7 @@ cdn/int_21.f90
cdn/int_21lo.f90
cdn/m_perp.f90
cdn/n_mat.f90
cdn/n_mat21.f90
cdn/od_abvac.f90
cdn/prp_qfft_map.f90
cdn/pwden.F90
......
......@@ -12,7 +12,7 @@ CONTAINS
SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,stars,&
vacuum,dimension,sphhar,sym,vTot,oneD,cdnvalJob,den,regCharges,dos,results,&
moments,coreSpecInput,mcd,slab,orbcomp)
moments,hub1,coreSpecInput,mcd,slab,orbcomp,greensfCoeffs,angle)
!************************************************************************************
! This is the FLEUR valence density generator
......@@ -34,10 +34,13 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
USE m_qal21 ! off-diagonal part of partial charges
USE m_abcof
USE m_nmat ! calculate density matrix for LDA + U
USE m_nmat21
USE m_vacden
USE m_pwden
USE m_forcea8
USE m_checkdopall
USE m_tetrahedronInit
USE m_gfcalc
USE m_cdnmt ! calculate the density and orbital moments etc.
USE m_orbmom ! coeffd for orbital moments
USE m_qmtsl ! These subroutines divide the input%film into vacuum%layers
......@@ -75,10 +78,13 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
TYPE(t_regionCharges), INTENT(INOUT) :: regCharges
TYPE(t_dos), INTENT(INOUT) :: dos
TYPE(t_moments), INTENT(INOUT) :: moments
TYPE(t_hub1ham), OPTIONAL, INTENT(INOUT) :: hub1
TYPE(t_coreSpecInput), OPTIONAL, INTENT(IN) :: coreSpecInput
TYPE(t_mcd), OPTIONAL, INTENT(INOUT) :: mcd
TYPE(t_slab), OPTIONAL, INTENT(INOUT) :: slab
TYPE(t_orbcomp), OPTIONAL, INTENT(INOUT) :: orbcomp
TYPE(t_greensfCoeffs), OPTIONAL, INTENT(INOUT) :: greensfCoeffs
REAL, OPTIONAL, INTENT(IN) :: angle(:)
! Scalar Arguments
INTEGER, INTENT(IN) :: eig_id, jspin
......@@ -88,7 +94,7 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
#endif
! Local Scalars
INTEGER :: ikpt,ikpt_i,jsp_start,jsp_end,ispin,jsp
INTEGER :: ikpt,ikpt_i,jsp_start,jsp_end,ispin,jsp,ne
INTEGER :: iErr,nbands,noccbd,iType
INTEGER :: skip_t,skip_tt,nbasfcn
LOGICAL :: l_orbcomprot, l_real, l_dosNdir, l_corespec
......@@ -97,6 +103,8 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
REAL,ALLOCATABLE :: we(:),eig(:)
INTEGER,ALLOCATABLE :: ev_list(:)
REAL, ALLOCATABLE :: f(:,:,:,:),g(:,:,:,:),flo(:,:,:,:) ! radial functions
REAL, ALLOCATABLE :: dosWeights(:,:),resWeights(:,:),eMesh(:)
INTEGER, ALLOCATABLE :: dosBound(:,:)
TYPE (t_lapw) :: lapw
TYPE (t_orb) :: orb
......@@ -110,10 +118,11 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
CALL timestart("cdnval")
l_real = sym%invs.AND.(.NOT.noco%l_soc).AND.(.NOT.noco%l_noco)
l_dosNdir = banddos%dos.AND.(banddos%ndir.EQ.-3)
IF (noco%l_mperp) THEN
IF (noco%l_mperp.OR.input%l_gfmperp) THEN
! when the off-diag. part of the desinsity matrix, i.e. m_x and
! m_y, is calculated inside the muffin-tins (l_mperp = T), cdnval
! is called only once. therefore, several spin loops have been
......@@ -160,9 +169,9 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
DO iType = 1, atoms%ntype
DO ispin = jsp_start, jsp_end
CALL genMTBasis(atoms,enpara,vTot,mpi,iType,ispin,usdus,f(:,:,0:,ispin),g(:,:,0:,ispin),flo(:,:,:,ispin))
CALL genMTBasis(atoms,enpara,vTot,mpi,iType,ispin,usdus,f(:,:,0:,ispin),g(:,:,0:,ispin),flo(:,:,:,ispin),input%l_dftspinpol)
END DO
IF (noco%l_mperp) CALL denCoeffsOffdiag%addRadFunScalarProducts(atoms,f,g,flo,iType)
IF (noco%l_mperp.OR.input%l_gfmperp) CALL denCoeffsOffdiag%addRadFunScalarProducts(atoms,f,g,flo,iType)
IF (banddos%l_mcd) CALL mcd_init(atoms,input,dimension,vTot%mt(:,0,:,:),g,f,mcd,iType,jspin)
IF (l_coreSpec) CALL corespec_rme(atoms,input,iType,dimension%nstd,input%jspins,jspin,results%ef,&
dimension%msh,vTot%mt(:,0,:,:),f,g)
......@@ -174,12 +183,15 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
jsp = MERGE(1,jspin,noco%l_noco)
!Get the energy mesh for the tetrahedron method
IF(atoms%n_gf>0) CALL greensfCoeffs%eMesh(ne,eMesh)
DO ikpt_i = 1,size(cdnvalJob%k_list)
ikpt=cdnvalJob%k_list(ikpt_i)
CALL lapw%init(input,noco, kpts,atoms,sym,ikpt,cell,.false., mpi)
skip_t = skip_tt
ev_list=cdnvaljob%compact_ev_list(ikpt_i,banddos%dos)
ev_list=cdnvaljob%compact_ev_list(ikpt_i,banddos%dos.OR.input%l_gf)
noccbd = SIZE(ev_list)
we = cdnvalJob%weights(ev_list,ikpt)
eig = results%eig(ev_list,ikpt,jsp)
......@@ -218,12 +230,26 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
! valence density in the atomic spheres
CALL eigVecCoeffs%init(input,DIMENSION,atoms,noco,jspin,noccbd)
IF (atoms%n_gf.GT.0.AND.(input%tria.OR.input%gfTet)) THEN
CALL timestart("TetrahedronWeights")
ALLOCATE(dosWeights(greensfCoeffs%ne,noccbd))
ALLOCATE(resWeights(greensfCoeffs%ne,noccbd))
resWeights = 0.0
ALLOCATE(dosBound(noccbd,2))
CALL tetrahedronInit(kpts,ikpt,results%eig(ev_list,:,jsp),noccbd,eMesh,ne,input%film,dosWeights,bounds=dosBound,dos=.TRUE.)
CALL timestop("TetrahedronWeights")
ENDIF
DO ispin = jsp_start, jsp_end
IF (input%l_f) CALL force%init2(noccbd,input,atoms)
CALL abcof(input,atoms,sym,cell,lapw,noccbd,usdus,noco,ispin,oneD,&
eigVecCoeffs%acof(:,0:,:,ispin),eigVecCoeffs%bcof(:,0:,:,ispin),&
eigVecCoeffs%ccof(-atoms%llod:,:,:,:,ispin),zMat,eig,force)
IF (atoms%n_u.GT.0) CALL n_mat(atoms,sym,noccbd,usdus,ispin,we,eigVecCoeffs,den%mmpMat(:,:,:,jspin))
IF (atoms%n_u.GT.0) CALL n_mat(atoms,sym,noccbd,usdus,ispin,we,eigVecCoeffs,den%mmpMat(:,:,:,ispin))
IF (atoms%n_u.GT.0.AND.noco%l_mperp.AND.(ispin==jsp_end)) CALL n_mat21(atoms,sym,noccbd,denCoeffsOffdiag,&
we,eigVecCoeffs,angle,den%mmpMat(:,:,:,3))
IF (atoms%n_gf.GT.0) CALL bzIntegrationGF(atoms,sym,input,angle,ispin,noccbd,dosWeights,resWeights,dosBound,kpts%wtkpt(ikpt),&
eig,denCoeffsOffdiag,usdus,eigVecCoeffs,greensfCoeffs,ispin==jsp_end)
! perform Brillouin zone integration and summation over the
! bands in order to determine the energy parameters for each atom and angular momentum
......@@ -232,6 +258,7 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
IF (noco%l_mperp.AND.(ispin==jsp_end)) CALL qal_21(dimension,atoms,input,noccbd,ev_list,noco,eigVecCoeffs,denCoeffsOffdiag,ikpt,dos)
! layer charge of each valence state in this k-point of the SBZ from the mt-sphere region of the film
IF (l_dosNdir) THEN
IF (PRESENT(slab)) CALL q_mt_sl(ispin,atoms,noccbd,ev_list,ikpt,noccbd,skip_t,noccbd,eigVecCoeffs,usdus,slab)
......@@ -249,6 +276,7 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
IF(l_coreSpec) CALL corespec_dos(atoms,usdus,ispin,dimension%lmd,kpts%nkpt,ikpt,dimension%neigd,&
noccbd,results%ef,banddos%sig_dos,eig,we,eigVecCoeffs)
END DO ! end loop over ispin
IF (atoms%n_gf.GT.0.AND.(input%tria.OR.input%gfTet)) DEALLOCATE(dosWeights,resWeights,dosBound)
IF (noco%l_mperp) CALL denCoeffsOffdiag%calcCoefficients(atoms,sphhar,sym,eigVecCoeffs,we,noccbd)
IF ((banddos%dos.OR.banddos%vacdos.OR.input%cdinf).AND.(banddos%ndir.GT.0)) THEN
......@@ -260,12 +288,13 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
#ifdef CPP_MPI
DO ispin = jsp_start,jsp_end
CALL mpi_col_den(mpi,sphhar,atoms,oneD,stars,vacuum,input,noco,ispin,regCharges,dos,&
results,denCoeffs,orb,denCoeffsOffdiag,den,den%mmpMat(:,:,:,jspin),mcd,slab,orbcomp)
results,denCoeffs,orb,denCoeffsOffdiag,den,mcd,slab,orbcomp,greensfCoeffs)
END DO
#endif
CALL cdnmt(mpi,input%jspins,atoms,sphhar,noco,jsp_start,jsp_end,&
enpara,vTot%mt(:,0,:,:),denCoeffs,usdus,orb,denCoeffsOffdiag,moments,den%mt)
enpara,vTot%mt(:,0,:,:),denCoeffs,usdus,orb,denCoeffsOffdiag,moments,den%mt,hub1,input%l_dftspinpol)
IF (mpi%irank==0) THEN
IF (l_coreSpec) CALL corespec_ddscs(jspin,input%jspins)
DO ispin = jsp_start,jsp_end
......
......@@ -62,8 +62,8 @@ CONTAINS
CALL pol_angle(mx,my,mz,betah,alphh)
WRITE (6,8026) betah,alphh
8025 FORMAT(2x,'--> local frame: ','mx=',f9.5,' my=',f9.5,' mz=',f9.5,' |m|=',f9.5)
8026 FORMAT(2x,'-->',10x,' delta beta=',f9.5,&
& ' delta alpha=',f9.5)
8026 FORMAT(2x,'-->',10x,' local beta=',f9.5,&
& ' local alpha=',f9.5)
IF (noco%l_relax(itype)) THEN
!---> rotate the (total (integrated) density matrix to obtain
......
MODULE m_nmat21
! ************************************************************
! This subroutine calculates the density matrix n^{s}_{m,m'}
! for a given atom 'n' and l-quantum number 'l'. The l's for
! all atoms are stored in lda_u(), if lda_u()<0, no +U is used.
! For details see Eq.(12) of Shick et al. PRB 60, 10765 (1999)
! Part of the LDA+U package G.B., Oct. 2000
! ************************************************************
CONTAINS
SUBROUTINE n_mat21(atoms,sym, ne,denCoeffsOffdiag,we,eigVecCoeffs,angle,n_mmp)
!
USE m_types
USE m_constants
IMPLICIT NONE
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_eigVecCoeffs),INTENT(IN) :: eigVecCoeffs
TYPE(t_denCoeffsOffDiag), INTENT(IN) :: denCoeffsOffDiag
REAL, INTENT(IN) :: angle(sym%nop)
! ..
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: ne
! ..
! .. Array Arguments ..
REAL, INTENT (IN) :: we(:)!(dimension%neigd)
COMPLEX, INTENT (INOUT) :: n_mmp(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_u)
! ..
! .. Local Scalars ..
COMPLEX c_0,phase
INTEGER i,j,k,l,m,mp,n,it,is,isi,natom,n_ldau,lp
INTEGER ilo,ilop,ll1,nn,lmp,lm,i_u,natomTemp
REAL fac
! ..
! .. Local Arrays ..
COMPLEX n_tmp(-3:3,-3:3),nr_tmp(-3:3,-3:3),d_tmp(-3:3,-3:3)
COMPLEX n1_tmp(-3:3,-3:3)
! ..
! calculate n_mat:
!
natom = 0
i_u = 1
DO n = 1,atoms%ntype
DO WHILE (i_u.LE.atoms%n_u)
IF (atoms%lda_u(i_u)%atomType.GT.n) EXIT
natomTemp = natom
n_tmp(:,:) = cmplx(0.0,0.0)
l = atoms%lda_u(i_u)%l
ll1 = (l+1)*l
DO nn = 1, atoms%neq(n)
natomTemp = natomTemp + 1
!
! prepare n_mat in local frame (in noco-calculations this depends
! also on alpha(n) and beta(n) )
!
DO m = -l,l
lm = ll1+m
DO mp = -l,l
lmp = ll1+mp
c_0 = cmplx(0.0,0.0)
DO i = 1,ne
c_0 = c_0 + we(i) * ( &
+ conjg(eigVecCoeffs%acof(i,lmp,natomTemp,2))*eigVecCoeffs%acof(i,lm,natomTemp,1) * denCoeffsOffdiag%uu21n(l,n) &
+ conjg(eigVecCoeffs%acof(i,lmp,natomTemp,2))*eigVecCoeffs%bcof(i,lm,natomTemp,1) * denCoeffsOffdiag%ud21n(l,n) &
+ conjg(eigVecCoeffs%bcof(i,lmp,natomTemp,2))*eigVecCoeffs%acof(i,lm,natomTemp,1) * denCoeffsOffdiag%du21n(l,n) &
+ conjg(eigVecCoeffs%bcof(i,lmp,natomTemp,2))*eigVecCoeffs%bcof(i,lm,natomTemp,1) * denCoeffsOffdiag%dd21n(l,n))
ENDDO
n_tmp(m,mp) = c_0
ENDDO
ENDDO
!
! add local orbital contribution (if there is one) (untested so far)
!
DO ilo = 1, atoms%nlo(n)
IF (atoms%llo(ilo,n).EQ.l) THEN
DO m = -l,l
lm = ll1+m
DO mp = -l,l
lmp = ll1+mp
c_0 = cmplx(0.0,0.0)
DO i = 1,ne
c_0 = c_0 + we(i) * ( &
conjg(eigVecCoeffs%acof(i,lmp,natomTemp,2))*eigVecCoeffs%ccof(m,i,ilo,natomTemp,1) * denCoeffsOffdiag%uulo21n(l,n) + &
conjg(eigVecCoeffs%ccof(mp,i,ilo,natomTemp,2))*eigVecCoeffs%acof(i,lm,natomTemp,1) * denCoeffsOffdiag%ulou21n(l,n) + &
conjg(eigVecCoeffs%bcof(i,lmp,natomTemp,2))*eigVecCoeffs%ccof(m,i,ilo,natomTemp,1) * denCoeffsOffdiag%dulo21n(l,n) + &
conjg(eigVecCoeffs%ccof(mp,i,ilo,natomTemp,2))*eigVecCoeffs%bcof(i,lm,natomTemp,1) * denCoeffsOffdiag%ulod21n(l,n))
ENDDO
DO ilop = 1, atoms%nlo(n)
IF (atoms%llo(ilop,n).EQ.l) THEN
DO i = 1,ne
c_0 = c_0 + we(i) * denCoeffsOffdiag%uloulop21n(ilo,ilop,n) *conjg(eigVecCoeffs%ccof(mp,i,ilop,natomTemp,2)) *eigVecCoeffs%ccof(m ,i,ilo ,natomTemp,1)
ENDDO
ENDIF
ENDDO
n_tmp(m,mp) = n_tmp(m,mp) + c_0
ENDDO
ENDDO
ENDIF
ENDDO
!
! n_mmp should be rotated by D_mm' ; compare force_a21
!
DO it = 1, sym%invarind(natomTemp)
fac = 1.0 / ( sym%invarind(natomTemp) * atoms%neq(n) )
is = sym%invarop(natomTemp,it)
isi = sym%invtab(is)
d_tmp(:,:) = cmplx(0.0,0.0)
DO m = -l,l
DO mp = -l,l
d_tmp(m,mp) = sym%d_wgn(m,mp,l,isi)
ENDDO
ENDDO
nr_tmp = matmul( transpose( conjg(d_tmp) ) , n_tmp)
n1_tmp = matmul( nr_tmp, d_tmp )
phase = exp((0.0,1.0)*angle(isi))
DO m = -l,l
DO mp = -l,l
n_mmp(m,mp,i_u) = n_mmp(m,mp,i_u) + conjg(n1_tmp(m,mp)) *fac*phase
ENDDO
ENDDO
ENDDO
ENDDO ! sum over equivalent
i_u = i_u +1
ENDDO
natom = natom + atoms%neq(n)
ENDDO ! loop over atom types
! do m=-l,l
! write(*,'(14f12.6)') (n_mmp21(m,mp,1),mp=-l,l)
! enddo
RETURN
END SUBROUTINE n_mat21
END MODULE m_nmat21
......@@ -10,11 +10,13 @@
! Robin Hilgers, Nov '19
MODULE m_alignSpinAxisMagn
USE m_magnMomFromDen
USE m_types
USE m_flipcdn
USE m_constants
USE m_polangle
IMPLICIT NONE
CONTAINS
SUBROUTINE rotateMagnetToSpinAxis(vacuum,sphhar,stars&
......
......@@ -11,7 +11,7 @@ MODULE m_cdnmt
!***********************************************************************
CONTAINS
SUBROUTINE cdnmt(mpi,jspd,atoms,sphhar,noco,jsp_start,jsp_end,enpara,&
vr,denCoeffs,usdus,orb,denCoeffsOffdiag,moments,rho)
vr,denCoeffs,usdus,orb,denCoeffsOffdiag,moments,rho,hub1,l_dftspinpol)
use m_constants,only: sfp_const
USE m_rhosphnlo
USE m_radfun
......@@ -26,6 +26,8 @@ CONTAINS
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_enpara), INTENT(IN) :: enpara
TYPE(t_moments), INTENT(INOUT) :: moments
TYPE(t_hub1ham), OPTIONAL, INTENT(INOUT) :: hub1
LOGICAL, INTENT(IN) :: l_dftspinpol
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: jsp_start,jsp_end,jspd
......@@ -39,12 +41,13 @@ CONTAINS
! ..
! .. Local Scalars ..
INTEGER itype,na,nd,l,lp,llp ,lh,j,ispin,noded,nodeu,llpb
INTEGER ilo,ilop,i
INTEGER ilo,ilop,i,i_hia,i_exc
REAL s,wronk,sumlm,qmtt
COMPLEX cs
LOGICAL l_hia
! ..
! .. Local Arrays ..
REAL qmtl(0:atoms%lmaxd,jspd,atoms%ntype),qmtllo(0:atoms%lmaxd)
REAL qmtl(0:atoms%lmaxd,jspd,atoms%ntype),qmtllo(0:atoms%lmaxd),vrTmp(atoms%jmtd)
CHARACTER(LEN=20) :: attributes(6)
! ..
......@@ -64,8 +67,8 @@ CONTAINS
ENDIF
! !$OMP PARALLEL DEFAULT(none) &
! !$OMP SHARED(usdus,rho,moments,qmtl) &
! !$OMP SHARED(atoms,jsp_start,jsp_end,enpara,vr,denCoeffs,sphhar)&
! !$OMP SHARED(usdus,rho,moments,qmtl,hub1) &
! !$OMP SHARED(atoms,jsp_start,jsp_end,enpara,vr,denCoeffs,sphhar,l_dftspinpol)&
! !$OMP SHARED(orb,noco,denCoeffsOffdiag,jspd)&
! !$OMP PRIVATE(itype,na,ispin,l,rho21,f,g,nodeu,noded,wronk,i,j,s,qmtllo,qmtt,nd,lh,lp,llp,llpb,cs)
IF (noco%l_mperp) THEN
......@@ -86,8 +89,25 @@ CONTAINS
!---> spherical component
DO ispin = jsp_start,jsp_end
DO l = 0,atoms%lmax(itype)
CALL radfun(l,itype,ispin,enpara%el0(l,itype,ispin),vr(1,itype,ispin),atoms,&
f(1,1,l,ispin),g(1,1,l,ispin),usdus, nodeu,noded,wronk)
!Check if the orbital is treated with Hubbard 1
l_hia=.FALSE.
DO i = atoms%n_u+1, atoms%n_u+atoms%n_hia
IF(atoms%lda_u(i)%atomType.EQ.itype.AND.atoms%lda_u(i)%l.EQ.l) THEN
l_hia=.TRUE.
ENDIF
ENDDO
!In the case of a spin-polarized calculation with Hubbard 1 we want to treat
!the correlated orbitals with a non-spin-polarized basis
IF(l_hia.AND.jspd.EQ.2.AND..NOT.l_dftspinpol) THEN
vrTmp = (vr(:,itype,1) + vr(:,itype,2))/2.0
ELSE
vrTmp = vr(:,itype,ispin)
ENDIF
CALL radfun(l,itype,ispin,enpara%el0(l,itype,ispin),vrTmp,atoms,&
f(1,1,l,ispin),g(1,1,l,ispin),usdus, nodeu,noded,wronk)