README.md

@@ -3,16 +3,41 @@ Welcome to the source code of FLEUR
 
 Please note that the documentation of the
 code can be found at the [FLEUR Homepage]
-(http://www.flapw.de/pm/index.php?n=User-Documentation.Documentation).
+(http://www.flapw.de/).
 
 For further instructions on Installation/Usage,
 please check the [FLEUR Homepage]
-(http://www.flapw.de/pm/index.php?n=User-Documentation.Documentation).
+(http://www.flapw.de/).
 
-The rest of this document summarizes only the 
-structure of the FLEUR source code. Did
-we mention to check the Homepage?
 
+## Bugs in FLEUR
+
+You might experience bugs in FLEUR :-).
+
+If you find a bug you should:
+
+A)  Report this bug by generating an Issue. 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.
+
+or/and
+
+B) Provide a bugfix. If the bug is only present in the development branch/ is due
+to a new feature under development simply commit your fix to the development branch.
+If you are fixing a bug in a release-version, please:
+* check out the git release branch: ```git checkout --track origin/release```
+* create a bugfix branch: ```git checkout -b bugfix_SOME_NAME_HERE```
+* apply your changes, test them and commit them
+* push your bugfix branch to the server: ``` git push -u origin bugfix_SOME_NAME_HERE```
+* create a merge request on the gitlab to have you bugfix merged with the release branch
+* check out the develop branch: ```git checkout develop```
+* merge your fix into the develop branch: ```git merge bugfix_SOME_NAME_HERE```
+
+
+
+
+## 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

cdn/cdnval.F90

@@ -134,7 +134,7 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
    CALL denCoeffs%init(atoms,sphhar,jsp_start,jsp_end)
    ! The last entry in denCoeffsOffdiag%init is l_fmpl. It is meant as a switch to a plot of the full magnet.
    ! density without the atomic sphere approximation for the magnet. density. It is not completely implemented (lo's missing).
-   CALL denCoeffsOffdiag%init(atoms,noco,sphhar,noco%l_mtnocopot)
+   CALL denCoeffsOffdiag%init(atoms,noco,sphhar,noco%l_mtnocopot.OR.noco%l_mperp)
    CALL force%init1(input,atoms)
    CALL orb%init(atoms,noco,jsp_start,jsp_end)
 

cdn/m_perp.f90

@@ -36,7 +36,7 @@ CONTAINS
     !     .. Local Scalars ..
     INTEGER iri
     REAL b_xavh,scale,b_con_outx,b_con_outy,mx,my,mz,&
-         &     alphh,betah,mz_tmp,mx_mix,my_mix,mz_mix
+         &     alphh,betah,mz_tmp,mx_mix,my_mix,mz_mix,absmag
     REAL    rho11,rho22, alphdiff
     COMPLEX rho21
     !     ..
@@ -56,11 +56,12 @@ CONTAINS
     mx = 2*REAL(qa21(itype))
     my = 2*AIMAG(qa21(itype))
     mz = chmom(itype,1) - chmom(itype,2)
-    WRITE  (6,8025) mx,my
+    absmag=SQRT(mx*mx+my*my+mz*mz)
+    WRITE  (6,8025) mx,my,mz,absmag
     !---> determine the polar angles of the moment vector in the local frame
     CALL pol_angle(mx,my,mz,betah,alphh)
     WRITE  (6,8026) betah,alphh
-8025 FORMAT(2x,'--> local frame: ','mx=',f9.5,' my=',f9.5)
+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)
 

cdn/pwden.F90

@@ -409,7 +409,7 @@ CONTAINS
              length_zfft(1) = stars%kq1_fft
              length_zfft(2) = stars%kq2_fft
              length_zfft(3) = stars%kq3_fft
-             call fft_interface(3,length_zfft,zfft,forw)
+             call fft_interface(3,length_zfft,zfft,forw,iv1d(1:lapw%nv(jspin),jspin))
              psir = real(zfft)
              psii = aimag(zfft)
              !--------------------------------
@@ -443,7 +443,7 @@ CONTAINS
                    length_zfft(1) = stars%kq1_fft
                    length_zfft(2) = stars%kq2_fft
                    length_zfft(3) = stars%kq3_fft
-                   call fft_interface(3,length_zfft,zfft,forw)
+                   call fft_interface(3,length_zfft,zfft,forw,iv1d(1:lapw%nv(jspin),jspin))
                    kpsir = real(zfft)
                    kpsii = aimag(zfft)
                    !--------------------------------

cdn_mt/CMakeLists.txt

@@ -10,10 +10,13 @@ cdn_mt/abcof3.F90
 cdn_mt/abcrot2.f90
 cdn_mt/cdnmt.F90
 cdn_mt/cdncore.F90
+cdn_mt/denMultipoleExp.f90
 cdn_mt/magMoms.f90
+cdn_mt/magDiMom.f90
 cdn_mt/orbMagMoms.f90
 cdn_mt/orb_comp2.f90
 cdn_mt/radfun.f90
+cdn_mt/resMoms.f90
 cdn_mt/rhomt.f90
 cdn_mt/rhomt21.f90
 cdn_mt/rhomtlo.f90
@@ -24,5 +27,7 @@ cdn_mt/rhosphnlo.f90
 cdn_mt/setabc1locdn.f90
 cdn_mt/setabc1locdn1.f90
 cdn_mt/calcDenCoeffs.f90
+cdn_mt/magnMomFromDen.f90
+cdn_mt/alignSpinAxisMagn.f90
 )
 

cdn_mt/alignSpinAxisMagn.f90

+!--------------------------------------------------------------------------------
+! Copyright (c) 2018 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
+! This file is part of FLEUR and avhttps://gcc.gnu.org/onlinedocs/gfortran/SQRT.htmlailable as free software under the conditions
+! of the MIT license as expressed in the LICENSE file in more detail.
+!------------------------------------------------------------------------------
+!  This routine allows to rotate the cdn in a way that the direction of magnetization aligns with the direction of the spin quantization axis. 
+!  This routine also allows to reverse the rotation by using the angles stored in atoms (phi_mt_avg,theta_mt_avg) which are generated by the 
+!  routine magnMomFromDen. 
+!
+! Robin Hilgers, Nov '19
+MODULE m_alignSpinAxisMagn
+
+USE m_magnMomFromDen
+USE m_types
+USE m_flipcdn
+USE m_constants
+USE m_polangle
+
+CONTAINS
+SUBROUTINE rotateMagnetToSpinAxis(vacuum,sphhar,stars&
+,sym,oneD,cell,noco,input,atoms,den)
+   TYPE(t_input), INTENT(INOUT)  :: input
+   TYPE(t_atoms), INTENT(INOUT)  :: atoms
+   TYPE(t_noco), INTENT(INOUT)   :: noco
+   TYPE(t_stars),INTENT(IN)      :: stars
+   TYPE(t_vacuum),INTENT(IN)     :: vacuum
+   TYPE(t_sphhar),INTENT(IN)     :: sphhar
+   TYPE(t_sym),INTENT(IN)        :: sym
+   TYPE(t_oneD),INTENT(IN)       :: oneD
+   TYPE(t_cell),INTENT(IN)       :: cell
+   TYPE(t_potden), INTENT(INOUT) :: den 
+
+   REAL                          :: moments(atoms%ntype,3)
+   REAL                          :: phiTemp(atoms%ntype),thetaTemp(atoms%ntype),pi
+   INTEGER                       :: i
+   pi=pimach()   
+!!TEMP
+!   REAL :: x,y,z
+   
+   phiTemp=noco%alph
+   thetaTemp=noco%beta
+   CALL magnMomFromDen(input,atoms,noco,den,moments)
+  ! DO i=1, atoms%ntype
+     ! IF (abs(atoms%theta_mt_avg(i)).LE. 0.001) THEN 
+     !    atoms%phi_mt_avg(i)=0.0
+     !    atoms%theta_mt_avg(i)=0.0
+     ! END IF
+  ! END DO
+   !write(*,*) "mx1"
+   !write(*,*) moments(1,1)
+   !write(*,*) "mz1"
+   !write(*,*) moments(1,3)
+   !write(*,*) "mx2"
+   !write(*,*) moments(2,1)
+   !write(*,*) "mz2"
+   !write(*,*) moments(2,3)
+   CALL flipcdn(atoms,input,vacuum,sphhar,stars,sym,noco,oneD,cell,-atoms%phi_mt_avg,-atoms%theta_mt_avg,den)
+  !write (*,*)"mx                my                     mz"
+  !CALL sphericaltocart(SQRT(moments(1,1)**2+moments(1,2)**2+moments(1,3)**2),thetaTemp(1),phiTemp(1),x,y,z)
+   !write(*,*) x,y,z
+   !CALL sphericaltocart(SQRT(moments(2,1)**2+moments(2,2)**2+moments(2,3)**2),thetaTemp(2),phiTemp(2),x,y,z)
+   !write(*,*) x,y,z
+   !write(*,*) "atoms%phi_mt_avg"
+   !!write(*,*) atoms%phi_mt_avg
+   !write(*,*) "atoms%theta_mt_avg"
+   !write(*,*) atoms%theta_mt_avg
+   noco%alph=mod(atoms%phi_mt_avg+phiTemp,2*pimach())
+   noco%beta=mod(atoms%theta_mt_avg+thetaTemp,2*pimach())
+   atoms%phi_mt_avg=noco%alph
+   atoms%theta_mt_avg=noco%beta
+ 
+   DO i=1, atoms%ntype
+      IF(noco%alph(i)<0) noco%alph(i)=noco%alph(i)+2*pi
+      IF(noco%beta(i)<0) THEN
+         noco%beta(i)=-noco%beta(i) 
+         noco%alph=noco%alph+pi 
+END IF
+      IF(noco%beta(i)>pi) THEN 
+         noco%beta(i)=pi-mod(noco%beta(i),pi)
+         noco%alph(i)=noco%alph(i)+pi
+      END IF
+      noco%alph=mod(noco%alph,2*pi)
+   End Do
+   write(*,*) "Noco Phi"
+   write(*,*) noco%alph
+   write(*,*) "Noco Theta"
+   write(*,*) noco%beta
+END SUBROUTINE rotateMagnetToSpinAxis
+
+
+SUBROUTINE rotateMagnetFromSpinAxis(noco,vacuum,sphhar,stars&
+,sym,oneD,cell,input,atoms,den,inDen)
+   TYPE(t_input), INTENT(INOUT)  :: input
+   TYPE(t_atoms), INTENT(INOUT)  :: atoms
+   TYPE(t_noco), INTENT(INOUT)	 :: noco
+   TYPE(t_stars),INTENT(IN)	 :: stars
+   TYPE(t_vacuum),INTENT(IN)     :: vacuum
+   TYPE(t_sphhar),INTENT(IN)     :: sphhar
+   TYPE(t_sym),INTENT(IN)        :: sym
+   TYPE(t_oneD),INTENT(IN)	 :: oneD
+   TYPE(t_cell),INTENT(IN)	 :: cell
+   TYPE(t_potden), INTENT(INOUT) :: den, inDen
+
+
+   CALL flipcdn(atoms,input,vacuum,sphhar,stars,sym,noco,oneD,cell,atoms%phi_mt_avg,atoms%theta_mt_avg,den)
+   CALL flipcdn(atoms,input,vacuum,sphhar,stars,sym,noco,oneD,cell,atoms%phi_mt_avg,atoms%theta_mt_avg,inDen)
+  
+   atoms%flipSpinPhi=0
+   atoms%flipSpinTheta=0
+   noco%alph=0
+   noco%beta=0
+
+END SUBROUTINE rotateMagnetFromSpinAxis
+
+
+END MODULE m_alignSpinAxisMagn
+

cdn_mt/cdnmt.F90

@@ -38,7 +38,7 @@ CONTAINS
     TYPE (t_denCoeffsOffdiag), INTENT(IN) :: denCoeffsOffdiag
     !     ..
     !     .. Local Scalars ..
-    INTEGER itype,na,nd,l,lp,llp ,lh,j,ispin,noded,nodeu
+    INTEGER itype,na,nd,l,lp,llp ,lh,j,ispin,noded,nodeu,llpb
     INTEGER ilo,ilop,i
     REAL s,wronk,sumlm,qmtt
     COMPLEX cs
@@ -63,11 +63,11 @@ CONTAINS
        ENDIF
     ENDIF
 
-    !$OMP PARALLEL DEFAULT(none) &
-    !$OMP SHARED(usdus,rho,moments,qmtl) &
-    !$OMP SHARED(atoms,jsp_start,jsp_end,enpara,vr,denCoeffs,sphhar)&
-    !$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,cs)
+!    !$OMP PARALLEL DEFAULT(none) &
+!    !$OMP SHARED(usdus,rho,moments,qmtl) &
+!    !$OMP SHARED(atoms,jsp_start,jsp_end,enpara,vr,denCoeffs,sphhar)&
+!    !$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
        ALLOCATE ( f(atoms%jmtd,2,0:atoms%lmaxd,jspd),g(atoms%jmtd,2,0:atoms%lmaxd,jspd) )
     ELSE
@@ -77,7 +77,7 @@ CONTAINS
 
     qmtl = 0
     
-    !$OMP DO
+!    !$OMP DO
     DO itype = 1,atoms%ntype
        na = 1
        DO i = 1, itype - 1
@@ -88,11 +88,13 @@ CONTAINS
           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)
+             llp = (l* (l+1))/2 + l
              DO j = 1,atoms%jri(itype)
                 s = denCoeffs%uu(l,itype,ispin)*( f(j,1,l,ispin)*f(j,1,l,ispin)+f(j,2,l,ispin)*f(j,2,l,ispin) )&
                      +   denCoeffs%dd(l,itype,ispin)*( g(j,1,l,ispin)*g(j,1,l,ispin)+g(j,2,l,ispin)*g(j,2,l,ispin) )&
                      + 2*denCoeffs%du(l,itype,ispin)*( f(j,1,l,ispin)*g(j,1,l,ispin)+f(j,2,l,ispin)*g(j,2,l,ispin) )
                 rho(j,0,itype,ispin) = rho(j,0,itype,ispin)+ s/(atoms%neq(itype)*sfp_const)
+                moments%rhoLRes(j,0,llp,itype,ispin) = moments%rhoLRes(j,0,llp,itype,ispin)+ s/(atoms%neq(itype)*sfp_const)
              ENDDO
           ENDDO
 
@@ -110,7 +112,7 @@ CONTAINS
                denCoeffs%cclo(1,1,itype,ispin),denCoeffs%acnmt(0,1,1,itype,ispin),&
                denCoeffs%bcnmt(0,1,1,itype,ispin),denCoeffs%ccnmt(1,1,1,itype,ispin),&
                f(1,1,0,ispin),g(1,1,0,ispin),&
-               rho(:,0:,itype,ispin),qmtllo)
+               rho(:,0:,itype,ispin),qmtllo,moments%rhoLRes(:,0:,:,itype,ispin))
 
 
           !--->       l-decomposed density for each atom type
@@ -146,6 +148,7 @@ CONTAINS
                            + denCoeffs%dunmt(llp,lh,itype,ispin)*(g(j,1,l,ispin)*f(j,1,lp,ispin)&
                            + g(j,2,l,ispin)*f(j,2,lp,ispin) )
                       rho(j,lh,itype,ispin) = rho(j,lh,itype,ispin)+ s/atoms%neq(itype)
+                      moments%rhoLRes(j,lh,llp,itype,ispin) = moments%rhoLRes(j,lh,llp,itype,ispin)+ s/atoms%neq(itype)
                    ENDDO
                 ENDDO
              ENDDO
@@ -183,6 +186,7 @@ CONTAINS
              !--->        calculate off-diagonal part of the density matrix
              !--->        spherical component
              DO l = 0,atoms%lmax(itype)
+                llp = (l* (l+1))/2 + l
                 DO j = 1,atoms%jri(itype)
                    cs = denCoeffsOffdiag%uu21(l,itype)*( f(j,1,l,2)*f(j,1,l,1) +f(j,2,l,2)*f(j,2,l,1) )&
                         + denCoeffsOffdiag%ud21(l,itype)*( f(j,1,l,2)*g(j,1,l,1) +f(j,2,l,2)*g(j,2,l,1) )&
@@ -191,7 +195,9 @@ CONTAINS
                    !rho21(j,0,itype) = rho21(j,0,itype)+ conjg(cs)/(atoms%neq(itype)*sfp_const)
                    rho21=CONJG(cs)/(atoms%neq(itype)*sfp_const)
                    rho(j,0,itype,3)=rho(j,0,itype,3)+REAL(rho21)
-                   rho(j,0,itype,4)=rho(j,0,itype,4)+imag(rho21)
+                   rho(j,0,itype,4)=rho(j,0,itype,4)+aimag(rho21)
+                   moments%rhoLRes(j,0,llp,itype,3) = moments%rhoLRes(j,0,llp,itype,3)+ REAL(conjg(cs)/(atoms%neq(itype)*sfp_const))
+                   moments%rhoLRes(j,0,llp,itype,4) = moments%rhoLRes(j,0,llp,itype,4)+ AIMAG(conjg(cs)/(atoms%neq(itype)*sfp_const))
                 ENDDO
              ENDDO
 
@@ -201,6 +207,7 @@ CONTAINS
                 DO l = 0,atoms%lmax(itype)
                    DO lp = 0,atoms%lmax(itype)
                       llp = lp*(atoms%lmax(itype)+1)+l+1
+                      llpb = (MAX(l,lp)* (MAX(l,lp)+1))/2 + MIN(l,lp)
                       DO j = 1,atoms%jri(itype)
                          cs = denCoeffsOffdiag%uunmt21(llp,lh,itype)*(f(j,1,lp,2)*f(j,1,l,1)&
                               + f(j,2,lp,2)*f(j,2,l,1) )+ denCoeffsOffdiag%udnmt21(llp,lh,itype)*(f(j,1,lp,2)*g(j,1,l,1)&
@@ -210,7 +217,9 @@ CONTAINS
                          !rho21(j,lh,itype)= rho21(j,lh,itype)+ CONJG(cs)/atoms%neq(itype)
                          rho21=CONJG(cs)/atoms%neq(itype)
                          rho(j,lh,itype,3)=rho(j,lh,itype,3)+REAL(rho21)
-                         rho(j,lh,itype,4)=rho(j,lh,itype,4)+imag(rho21)
+                         rho(j,lh,itype,4)=rho(j,lh,itype,4)+aimag(rho21)
+                         moments%rhoLRes(j,lh,llpb,itype,3)= moments%rhoLRes(j,lh,llpb,itype,3) + REAL(conjg(cs)/atoms%neq(itype))
+                         moments%rhoLRes(j,lh,llpb,itype,4)= moments%rhoLRes(j,lh,llpb,itype,4) + AIMAG(conjg(cs)/atoms%neq(itype))
                       ENDDO
                    ENDDO
                 ENDDO
@@ -220,9 +229,9 @@ CONTAINS
        ENDIF ! noco%l_mperp
 
     ENDDO ! end of loop over atom types
-    !$OMP END DO
+!    !$OMP END DO
     DEALLOCATE ( f,g)
-    !$OMP END PARALLEL
+!    !$OMP END PARALLEL
 
     WRITE (6,FMT=8000)
 8000 FORMAT (/,5x,'l-like charge',/,t6,'atom',t15,'s',t24,'p',&

cdn_mt/denMultipoleExp.f90

+MODULE m_denMultipoleExp
+
+IMPLICIT NONE
+
+CONTAINS
+
+SUBROUTINE denMultipoleExp(input, mpi, atoms, sphhar, stars, sym, cell, oneD, den)
+
+   USE m_types
+   USE m_constants
+   USE m_mpmom
+
+   TYPE(t_input),  INTENT(IN) :: input
+   TYPE(t_mpi),    INTENT(IN) :: mpi
+   TYPE(t_atoms),  INTENT(IN) :: atoms
+   TYPE(t_sphhar), INTENT(IN) :: sphhar
+   TYPE(t_stars),  INTENT(IN) :: stars
+   TYPE(t_sym),    INTENT(IN) :: sym
+   TYPE(t_cell),   INTENT(IN) :: cell
+   TYPE(t_oneD),   INTENT(IN) :: oneD
+   TYPE(t_potden), INTENT(IN) :: den
+
+   type(t_potden)             :: workDen
+   COMPLEX                    :: qlm(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
+
+   IF(input%jspins == 2) THEN
+      IF(mpi%irank.EQ.0) THEN
+         WRITE(6,*) 'Multipole expansion for spin-up density:'
+         WRITE(6,*) '======================================='
+      END IF
+      qlm = CMPLX(0.0,0.0)
+      workDen = den
+      CALL mpmom(input,mpi,atoms,sphhar,stars,sym,cell,oneD,workDen%pw(1:,1),workDen%mt(:,0:,1:,1),POTDEN_TYPE_DEN,qlm,.FALSE.)
+      IF(mpi%irank.EQ.0) THEN
+         WRITE(6,*) '======================================='
+      END IF
+
+      IF(mpi%irank.EQ.0) THEN
+         WRITE(6,*) 'Multipole expansion for spin-down density:'
+         WRITE(6,*) '======================================='
+      END IF
+      qlm = CMPLX(0.0,0.0)
+      CALL mpmom(input,mpi,atoms,sphhar,stars,sym,cell,oneD,workDen%pw(1:,2),workDen%mt(:,0:,1:,2),POTDEN_TYPE_DEN,qlm,.FALSE.)
+      IF(mpi%irank.EQ.0) THEN
+         WRITE(6,*) '======================================='
+      END IF
+   END IF
+
+   IF(mpi%irank.EQ.0) THEN
+      WRITE(6,*) 'Multipole expansion for charge density:'
+      WRITE(6,*) '======================================='
+   END IF
+   qlm = CMPLX(0.0,0.0)
+   workDen = den
+   IF(input%jspins == 2) CALL workDen%SpinsToChargeAndMagnetisation()
+   CALL mpmom(input,mpi,atoms,sphhar,stars,sym,cell,oneD,workDen%pw(1:,1),workDen%mt(:,0:,1:,1),POTDEN_TYPE_DEN,qlm,.FALSE.)
+   IF(mpi%irank.EQ.0) THEN
+      WRITE(6,*) '======================================='
+   END IF
+
+   IF(input%jspins == 2) THEN
+      IF(mpi%irank.EQ.0) THEN
+         WRITE(6,*) 'Multipole expansion for magnetization density:'
+         WRITE(6,*) '======================================='
+      END IF
+      qlm = CMPLX(0.0,0.0)
+      CALL mpmom(input,mpi,atoms,sphhar,stars,sym,cell,oneD,workDen%pw(1:,2),workDen%mt(:,0:,1:,2),POTDEN_TYPE_DEN,qlm,.FALSE.)
+      IF(mpi%irank.EQ.0) THEN
+         WRITE(6,*) '======================================='
+      END IF
+   END IF
+
+END SUBROUTINE denMultipoleExp
+
+END MODULE m_denMultipoleExp

cdn_mt/magnMomFromDen.f90

+!--------------------------------------------------------------------------------
+! Copyright (c) 2018 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
+! This file is part of FLEUR and avhttps://gcc.gnu.org/onlinedocs/gfortran/SQRT.htmlailable as free software under the conditions
+! of the MIT license as expressed in the LICENSE file in more detail.
+!------------------------------------------------------------------------------
+! This routine calculates the magnetic moments and the corresponding directions 
+! (angles) according to the Atoms in the system. 
+! 
+!
+! Robin Hilgers, Nov '19
+
+MODULE m_magnMomFromDen
+CONTAINS
+SUBROUTINE magnMomFromDen(input,atoms,noco,den,moments)
+   USE m_constants
+   USE m_types
+   USE m_intgr
+   USE m_juDFT
+   USE m_polangle
+   IMPLICIT NONE
+
+
+   TYPE(t_input), INTENT(IN)     ::  input
+   TYPE(t_atoms), INTENT(INOUT)  ::  atoms
+   TYPE(t_noco), INTENT(IN)      ::  noco
+   TYPE(t_potden),INTENT(IN)     ::  den
+   REAL, INTENT(OUT)             ::  moments(atoms%ntype,3)
+
+   INTEGER                       ::  jsp,i,j
+   REAL                          ::  mx,my,mz
+
+   REAL, ALLOCATABLE             ::  dummyResults(:,:)
+
+  ALLOCATE(dummyResults(SIZE(den%mt,3),SIZE(den%mt,4)))
+
+  IF(noco%l_mtNocoPot) THEN
+     jsp=4
+  ELSE 
+     jsp=input%jspins
+  END IF
+!!Loop over Spins and Atoms
+   DO i=1, atoms%ntype 
+      DO j=1, jsp
+!!Integration over r
+           CALL intgr3(den%mt(:,0,i,j), atoms%rmsh(:,i),atoms%dx(i),atoms%jri(i),dummyResults(i,j))
+!!Considering Lattice harmonics integral (Only L=0 component does not vanish and has a factor of sqrt(4*Pi))
+          dummyResults(i,j)=dummyResults(i,j)*sfp_const
+      END DO
+   END DO 
+!!Assign results
+   DO i=1 , atoms%ntype
+   IF (noco%l_mtNocoPot) THEN
+      moments(i,1:2)=2*dummyResults(i,3:4)
+   END IF
+      moments(i,3)=dummyResults(i,1)-dummyResults(i,2)
+   END DO
+DEALLOCATE(dummyResults)
+
+!!Calculation of Angles
+   DO i=1 , atoms%ntype
+      mx=moments(i,1)
+      my=moments(i,2)
+      mz=moments(i,3)
+      CALL pol_angle(mx,my,mz,atoms%theta_mt_avg(i),atoms%phi_mt_avg(i))
+      IF(mx<0) atoms%theta_mt_avg(i)=-atoms%theta_mt_avg(i)
+   ENDDO
+
+END SUBROUTINE magnMomFromDen
+END MODULE m_magnMomFromDen

cdn_mt/resMoms.f90

+! Copyright (c) 2018 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
+! This file is part of FLEUR and available as free software under the conditions
+! of the MIT license as expressed in the LICENSE file in more detail.
+!--------------------------------------------------------------------------------
+
+MODULE m_resMoms
+
+CONTAINS
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! This subroutine calculates and writes out intraatomic electric and magnetic dipole 
+! moments resolved with respect to their orbital (angular momentum) origins.
+!
+!                                           GM'2018
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+SUBROUTINE resMoms(input,atoms,sphhar,noco,den,rhoLRes)
+
+   USE m_constants
+   USE m_types
+   USE m_juDFT
+   USE m_magDiMom
+
+   IMPLICIT NONE
+
+   TYPE(t_input),         INTENT(IN)    :: input
+   TYPE(t_atoms),         INTENT(IN)    :: atoms
+   TYPE(t_sphhar),        INTENT(IN)    :: sphhar
+   TYPE(t_noco),          INTENT(IN)    :: noco
+   TYPE(t_potden),        INTENT(IN)    :: den
+   REAL,                  INTENT(IN)    :: rhoLRes(:,0:,0:,:,:)
+
+   REAL,    ALLOCATABLE :: rhoTemp(:,:,:,:)
+
+   REAL    :: t_op(3,atoms%ntype), elecDip(3,atoms%ntype)
+   REAL    :: res_T_op(3,atoms%ntype,0:(atoms%lmaxd*(atoms%lmaxd+1))/2+atoms%lmaxd)
+   REAL    :: resElecDip(3,atoms%ntype,0:(atoms%lmaxd*(atoms%lmaxd+1))/2+atoms%lmaxd)
+
+   INTEGER :: iType, l, lp, llp
+
+   IF(input%jspins.EQ.1) RETURN
+   IF(.NOT.noco%l_noco) RETURN
+
+   t_op = 0.0
+   res_T_op = 0.0
+   elecDip = 0.0
+   resElecDip = 0.0
+   ALLOCATE(rhoTemp(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,4))
+
+   rhoTemp = 0.0
+
+   rhoTemp(:,:,:,1) = den%mt(:,:,:,1)
+   rhoTemp(:,:,:,2) = den%mt(:,:,:,2)
+
+   IF (noco%l_mperp) THEN
+      rhoTemp(:,:,:,3) = den%mt(:,:,:,3)
+      rhoTemp(:,:,:,4) = den%mt(:,:,:,4)
+!      WRITE(5000,'(f15.8)') den%mt(:,:,:,3)
+!      WRITE(5000,'(f15.8)') den%mt(:,:,:,4)
+   END IF
+
+   CALL magDiMom(input,atoms,sphhar,noco,noco%l_mperp,rhoTemp,t_op,elecDip)
+
+   DO l = 0, atoms%lmaxd
+      DO lp = 0, l
+         llp = (l* (l+1))/2 + lp
+         rhoTemp = 0.0
+         rhoTemp(:,:,:,1) = rhoLRes(:,:,llp,:,1)
+         rhoTemp(:,:,:,2) = rhoLRes(:,:,llp,:,2)
+         rhoTemp(:,:,:,3) = rhoLRes(:,:,llp,:,3)
+         rhoTemp(:,:,:,4) = rhoLRes(:,:,llp,:,4)
+         CALL magDiMom(input,atoms,sphhar,noco,noco%l_mperp,rhoTemp,res_T_op(:,:,llp),resElecDip(:,:,llp))
+      END DO
+   END DO
+
+   DO iType = 1, atoms%ntype
+      WRITE(6,*) 'Intraatomic electric and magnetic dipole moments for atom type ', iType,':'
+      WRITE(6,'(a)')        '             lowL  largeL      p_x            p_y            p_z            t_x            t_y            t_z'
+      WRITE(6,'(a,6f15.8)') 'Overall:              ', elecDip(:,iType), t_op(:,iType)
+      DO l = 0, atoms%lmax(iType)
+         DO lp = 0, l
+            llp = (l* (l+1))/2 + lp
+            IF(ALL(ABS(res_T_op(:,iType,llp)).LT.1.0e-8).AND.&
+               ALL(ABS(resElecDip(:,iType,llp)).LT.1.0e-8)) CYCLE
+            WRITE(6,'(a,2i6,6f15.8)') '          ', lp, l, resElecDip(:,iType,llp),res_T_op(:,iType,llp)
+         END DO
+      END DO
+   END DO
+
+END SUBROUTINE resMoms
+
+END MODULE m_resMoms

cdn_mt/rhosphnlo.f90

@@ -13,7 +13,7 @@ MODULE m_rhosphnlo
   !***********************************************************************
 CONTAINS
   SUBROUTINE rhosphnlo(itype,atoms,sphhar, uloulopn,dulon,uulon,&
-       ello,vr, aclo,bclo,cclo,acnmt,bcnmt,ccnmt,f,g, rho,qmtllo)
+       ello,vr, aclo,bclo,cclo,acnmt,bcnmt,ccnmt,f,g, rho,qmtllo,rhoLRes)
 
     USE m_constants, ONLY : c_light,sfp_const
     USE m_radsra
@@ -36,11 +36,12 @@ CONTAINS
     REAL,    INTENT (IN) :: f(atoms%jmtd,2,0:atoms%lmaxd),g(atoms%jmtd,2,0:atoms%lmaxd)
     REAL,    INTENT (INOUT) :: qmtllo(0:atoms%lmaxd)
     REAL,    INTENT (INOUT) :: rho(atoms%jmtd,0:sphhar%nlhd)
+    REAL,    INTENT (INOUT) :: rhoLRes(atoms%jmtd,0:sphhar%nlhd,0:(atoms%lmaxd*(atoms%lmaxd+1))/2+atoms%lmaxd)
     !     ..
     !     .. Local Scalars ..
     REAL dsdum,usdum ,c_1,c_2
-    INTEGER j,l,lh,lo,lop,lp,nodedum
-    REAL dus,ddn,c
+    INTEGER j,l,lh,lo,lop,lp,nodedum,llp
+    REAL dus,ddn,c,temp
     !     ..
     !     .. Local Arrays ..
     REAL,    ALLOCATABLE :: flo(:,:,:),glo(:,:)
@@ -90,16 +91,21 @@ CONTAINS
 
     DO lo = 1,atoms%nlo(itype)
        l = atoms%llo(lo,itype)
+       llp = (l* (l+1))/2 + l
        DO j = 1,atoms%jri(itype)
-          rho(j,0) = rho(j,0) + c_2 *&
-               (aclo(lo) * ( f(j,1,l)*flo(j,1,lo) +f(j,2,l)*flo(j,2,lo) ) +&
-               bclo(lo) * ( g(j,1,l)*flo(j,1,lo) +g(j,2,l)*flo(j,2,lo) ) )
+          temp = c_2 *&
+                 (aclo(lo) * ( f(j,1,l)*flo(j,1,lo) +f(j,2,l)*flo(j,2,lo) ) +&
+                 bclo(lo) * ( g(j,1,l)*flo(j,1,lo) +g(j,2,l)*flo(j,2,lo) ) )
+          rho(j,0) = rho(j,0) + temp
+          rhoLRes(j,0,llp) = rhoLRes(j,0,llp) + temp
        END DO
        DO lop = 1,atoms%nlo(itype)
           IF (atoms%llo(lop,itype).EQ.l) THEN
              DO j = 1,atoms%jri(itype)
-                rho(j,0) = rho(j,0) + c_2 * cclo(lop,lo) *&
+                temp = c_2 * cclo(lop,lo) *&
                      ( flo(j,1,lop)*flo(j,1,lo) +flo(j,2,lop)*flo(j,2,lo) )
+                rho(j,0) = rho(j,0) + temp
+                rhoLRes(j,0,llp) = rhoLRes(j,0,llp) + temp
              END DO
           END IF
        END DO
@@ -111,18 +117,27 @@ CONTAINS
     DO lh = 1,sphhar%nlh(atoms%ntypsy(atoms%nat))
        DO lp = 0,atoms%lmax(itype)
           DO lo = 1,atoms%nlo(itype)
+             l = atoms%llo(lo,itype)
+             llp = (MAX(l,lp)* (MAX(l,lp)+1))/2 + MIN(l,lp)
              DO j = 1,atoms%jri(itype)
-                rho(j,lh) = rho(j,lh) + c_1 * (&
+                temp = c_1 * (&
                      acnmt(lp,lo,lh) * (f(j,1,lp)*flo(j,1,lo) +f(j,2,lp)*flo(j,2,lo) ) +&
                      bcnmt(lp,lo,lh) * (g(j,1,lp)*flo(j,1,lo) +g(j,2,lp)*flo(j,2,lo) ) )
+                rho(j,lh) = rho(j,lh) + temp
+                rhoLRes(j,lh,llp) = rhoLRes(j,lh,llp) + temp
              END DO
           END DO
        END DO
        DO lo = 1,atoms%nlo(itype)
+          l = atoms%llo(lo,itype)
           DO lop = 1,atoms%nlo(itype)
+             lp = atoms%llo(lop,itype)
+             llp = (MAX(l,lp)* (MAX(l,lp)+1))/2 + MIN(l,lp)
              DO j = 1,atoms%jri(itype)
-                rho(j,lh) = rho(j,lh) + c_1 * ccnmt(lop,lo,lh) *&
+                temp = c_1 * ccnmt(lop,lo,lh) *&
                      ( flo(j,1,lop)*flo(j,1,lo) +flo(j,2,lop)*flo(j,2,lo) )
+                rho(j,lh) = rho(j,lh) + temp
+                rhoLRes(j,lh,llp) = rhoLRes(j,lh,llp) + temp
              END DO
           END DO
        END DO

cmake/CompilerConfig.txt

@@ -8,6 +8,7 @@ include("cmake/tests/test_XML.cmake")
 include("cmake/tests/test_LAPACK.cmake")
 include("cmake/tests/test_MPI.cmake")
 include("cmake/tests/test_FFTMKL.cmake")
+include("cmake/tests/test_SpFFT.cmake")
 include("cmake/tests/test_HDF5.cmake")
 include("cmake/tests/test_Wannier90.cmake")
 #include("cmake/tests/test_Wannier4.cmake")

cmake/Files_and_Targets.txt

@@ -21,6 +21,7 @@ include(cdn/CMakeLists.txt)
 include(diagonalization/CMakeLists.txt)
 include(eigen_soc/CMakeLists.txt)
 include(math/CMakeLists.txt)
+include(fft/CMakeLists.txt)
 include(orbdep/CMakeLists.txt)
 include(cdn_mt/CMakeLists.txt)
 include(dos/CMakeLists.txt)

cmake/ReportConfig.txt

@@ -20,6 +20,7 @@ message("${Green}  LAPACK Library found    : ${CReset}    ${FLEUR_USE_LAPACK}")
 message("${Red}These Libraries are optional:${CReset}")
 message("${Green}  ELPA (one node) Library found : ${CReset}    ${FLEUR_USE_ELPA_ONENODE}")
 message("${Green}  FFT from MKL found            : ${CReset}    ${FLEUR_USE_FFTMKL}")
+message("${Green}  SpFFT found                   : ${CReset}    ${FLEUR_USE_SPFFT}")
 message("${Green}  LibXC Library found           : ${CReset}    ${FLEUR_USE_LIBXC}")
 message("${Green}  HDF5 Library found            : ${CReset}    ${FLEUR_USE_HDF5}")
 message("${Green}  Wannier90 1.2 Library found   : ${CReset}    ${FLEUR_USE_WANN}")

cmake/process_arguments.sh

@@ -90,7 +90,7 @@ the test and might lead to the configuration to fail.
 
 Command line option to compile external libraries:
   -external # : download and compile external libraries before building FLEUR
-                currently 'hdf5' and 'libxc' are possible options. The switch 
+                currently 'hdf5','libxc' and 'chase' are possible options. The switch 
                 can be specified multiple times
 
 Options to specify Fortran/Linker flags:

cmake/tests/test_SpFFT.cmake

+#First check if we can compile with SpFFT
+try_compile(FLEUR_USE_SPFFT ${CMAKE_BINARY_DIR} ${CMAKE_SOURCE_DIR}/cmake/tests/test_SpFFT.f90
+LINK_LIBRARIES ${FLEUR_LIBRARIES}
+            )
+
+message("SpFFT found:${FLEUR_USE_SPFFT}")
+if (FLEUR_USE_SPFFT)
+   set(FLEUR_MPI_DEFINITIONS ${FLEUR_MPI_DEFINITIONS} "CPP_SPFFT")
+   set(FLEUR_DEFINITIONS ${FLEUR_DEFINITIONS} "CPP_SPFFT")
+endif()
+

cmake/tests/test_SpFFT.f90

+! The contents of this test program are adapted from one of the examples
+! of the SpFFT library. This was done with permission from the SpFFT
+! developer Simon Frasch. Thank you for that! - G.M.
+
+! Note the following copyright and license aspects of the SpFFT library:
+
+! Copyright (c) 2019 ETH Zurich, Simon Frasch
+! 
+! Redistribution and use in source and binary forms, with or without
+! modification, are permitted provided that the following conditions are met:
+! 
+! 1. Redistributions of source code must retain the above copyright notice,
+!    this list of conditions and the following disclaimer.
+! 2. Redistributions in binary form must reproduce the above copyright
+!    notice, this list of conditions and the following disclaimer in the
+!    documentation and/or other materials provided with the distribution.
+! 3. Neither the name of the copyright holder nor the names of its contributors
+!    may be used to endorse or promote products derived from this software
+!    without specific prior written permission.
+! 
+! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+! POSSIBILITY OF SUCH DAMAGE.
+
+program test
+    use iso_c_binding
+    use spfft
+    implicit none
+    integer :: i, j, k, counter
+    integer, parameter :: dimX = 2
+    integer, parameter :: dimY = 2
+    integer, parameter :: dimZ = 2
+    integer, parameter :: maxNumLocalZColumns = dimX * dimY
+    integer, parameter :: processingUnit = 1
+    integer, parameter :: maxNumThreads = -1
+    type(c_ptr) :: grid = c_null_ptr
+    type(c_ptr) :: transform = c_null_ptr
+    integer :: errorCode = 0
+    integer, dimension(dimX * dimY * dimZ * 3):: indices = 0
+    complex(C_DOUBLE_COMPLEX), dimension(dimX * dimY * dimZ):: frequencyElements
+    complex(C_DOUBLE_COMPLEX), pointer :: spaceDomain(:,:,:)
+    type(c_ptr) :: realValuesPtr
+
+
+    counter = 0
+    do k = 1, dimZ
+        do j = 1, dimY
+            do i = 1, dimX
+             frequencyElements(counter + 1) = cmplx(counter, -counter)
+             indices(counter * 3 + 1) = i - 1
+             indices(counter * 3 + 2) = j - 1
+             indices(counter * 3 + 3) = k - 1
+             counter = counter + 1
+            end do
+        end do
+    end do
+
+    ! print input
+    print *, "Input:"
+    do i = 1, size(frequencyElements)
+         print *, frequencyElements(i)
+    end do
+
+
+    ! create grid and transform
+    errorCode = spfft_grid_create(grid, dimX, dimY, dimZ, maxNumLocalZColumns, processingUnit, maxNumThreads);
+    if (errorCode /= SPFFT_SUCCESS) error stop
+    errorCode = spfft_transform_create(transform, grid, processingUnit, 0, &
+       dimX, dimY, dimZ, dimZ, size(frequencyElements), 0, indices)
+    if (errorCode /= SPFFT_SUCCESS) error stop
+
+    ! grid can be safely destroyed after creating all required transforms
+    errorCode = spfft_grid_destroy(grid)
+    if (errorCode /= SPFFT_SUCCESS) error stop
+
+    ! set space domain array to use memory allocted by the library
+    errorCode = spfft_transform_get_space_domain(transform, processingUnit, realValuesPtr)
+    if (errorCode /= SPFFT_SUCCESS) error stop
+
+    ! transform backward
+    errorCode = spfft_transform_backward(transform, frequencyElements, processingUnit)
+    if (errorCode /= SPFFT_SUCCESS) error stop
+
+
+    call c_f_pointer(realValuesPtr, spaceDomain, [dimX,dimY,dimZ])
+
+    print *, ""
+    print *, "After backward transform:"
+    do k = 1, size(spaceDomain, 3)
+        do j = 1, size(spaceDomain, 2)
+            do i = 1, size(spaceDomain, 1)
+             print *, spaceDomain(i, j, k)
+            end do
+        end do
+    end do
+
+    ! transform forward (will invalidate space domain data)
+    errorCode = spfft_transform_forward(transform, processingUnit, frequencyElements, 0)
+    if (errorCode /= SPFFT_SUCCESS) error stop
+
+    print *, ""
+    print *, "After forward transform (without scaling):"
+    do i = 1, size(frequencyElements)
+             print *, frequencyElements(i)
+    end do
+
+    ! destroying the final transform will free the associated memory
+    errorCode = spfft_transform_destroy(transform)
+    if (errorCode /= SPFFT_SUCCESS) error stop
+
+end program test

cmake/tests/test_Wannier90.cmake

@@ -22,7 +22,9 @@ foreach(ADD_String "-lwannier;-lmkl_intel_lp64;-lmkl_sequential;-lmkl_core" )
             LINK_LIBRARIES ${TEST_LIBRARIES} OUTPUT_VARIABLE TELL_ME
             )
 
-     message("TELL_ME=${TELL_ME}")
+     if(DEFINED ENV{VERBOSE})
+        message("TELL_ME=${TELL_ME}")
+     endif()
 
 
      if (FLEUR_USE_WANN)

core/setcor.f90

@@ -30,6 +30,9 @@ CONTAINS
     CHARACTER(len=13) :: fname
     !     ..
 
+    kappa(:) = -100
+    nprnc(:) = -100
+
     l_clf = .FALSE.  
     WRITE(fname,"('corelevels.',i2.2)") NINT(atoms%zatom(itype))
     INQUIRE (file=fname, exist=l_clf)
@@ -206,7 +209,7 @@ CONTAINS
        ENDIF
     ENDIF
 
-    ! modify default electron configuration according to explicitely provided setting in inp.xml
+    ! modify default electron configuration according to explicitly provided setting in inp.xml
     IF(input%l_inpXML) THEN
        nst = max(nst,atoms%numStatesProvided(itype))
        IF (atoms%numStatesProvided(itype).NE.0) THEN

diagonalization/chase_diag.F90

@@ -89,9 +89,7 @@ CONTAINS
 
 #ifdef CPP_CHASE
     SUBROUTINE init_chase(mpi,DIMENSION,input,atoms,kpts,noco,l_real)
-    USE m_types_mpimat
-    USE m_types_setup
-    USE m_types_mpi
+    USE m_types
     USE m_judft
     USE m_eig66_io
 

diagonalization/eigen_diag.F90

@@ -20,7 +20,7 @@ CONTAINS
     USE m_elpa_onenode
     USE m_scalapack
     USE m_elemental
-!    USE m_chase_diag
+    USE m_chase_diag
     USE m_types_mpimat
     USE m_types_gpumat
 !    USE m_matrix_copy
@@ -71,7 +71,7 @@ CONTAINS
        CALL lapack_diag(hmat,smat,ne,eig,ev)
     CASE (diag_chase)
        IF (.NOT.(PRESENT(ikpt).AND.PRESENT(jsp).AND.PRESENT(iter))) CALL judft_error("Optional arguments must be present for chase in eigen_diag")
-!       CALL chase_diag(hmat,smat,ikpt,jsp,iter,ne,eig,ev)
+       CALL chase_diag(hmat,smat,ikpt,jsp,iter,ne,eig,ev)
     CASE (diag_debugout)
        CALL diag_writeout(smat,hmat)
     CASE default

eigen/hsmt_ab.F90

@@ -198,7 +198,7 @@ CONTAINS
     REAL,ALLOCATABLE   :: gkrot(:,:)
     LOGICAL :: l_apw
    
-    ALLOCATE(c_ph(maxval(lapw%nv),MERGE(2,1,noco%l_ss)))
+    ALLOCATE(c_ph(maxval(lapw%nv),MERGE(2,1,noco%l_ss.or.noco%l_mtNocoPot)))
     ALLOCATE(gkrot(3,maxval(lapw%nv)))
 
     lmax=MERGE(atoms%lnonsph(n),atoms%lmax(n),l_nonsph)

eigen/hsmt_fjgj.F90

@@ -169,7 +169,7 @@ CONTAINS
                 fj(k,l,ispin,intspin) = 1.0*con1 * ff / usdus%us(l,n,ispin)
                 gj(k,l,ispin,intspin) = 0.0
              ELSE
-                IF (noco%l_constr.or.l_socfirst) THEN
+                IF (noco%l_constr.or.l_socfirst.OR.noco%l_mtNocoPot) THEN
                    DO jspin = 1, input%jspins
                       fj(k,l,jspin,intspin) = ws(jspin) * ( usdus%uds(l,n,jspin)*gg - usdus%duds(l,n,jspin)*ff )
                       gj(k,l,jspin,intspin) = ws(jspin) * ( usdus%dus(l,n,jspin)*ff - usdus%us(l,n,jspin)*gg )

eigen/hsmt_mtNocoPot_offdiag.F90

@@ -41,7 +41,7 @@ CONTAINS
     CALL hmat_tmp%clear()
     !The spin1,2 matrix is calculated(real part of potential)
     CALL hsmt_nonsph(n,mpi,sym,atoms,3,1,1,chi_one,noco,cell,lapw,td,&
-         fj(:,0:,1,:),gj(:,0:,2,:),hmat_tmp)
+         fj(:,0:,1,:),gj(:,0:,1,:),hmat_tmp)
 
     CALL hsmt_spinor(3,n,noco,chi) !spinor for off-diagonal part
     CALL hsmt_distspins(chi,hmat_tmp,hmat)
@@ -56,7 +56,7 @@ CONTAINS
     !The spin1,2 matrix is calculated(imag part of potential)
     chi_one=CMPLX(0.,1.)
     CALL hsmt_nonsph(n,mpi,sym,atoms,4,1,1,chi_one,noco,cell,lapw,td,&
-         fj(:,0:,1,:),gj(:,0:,2,:),hmat_tmp)
+         fj(:,0:,1,:),gj(:,0:,1,:),hmat_tmp)
 
     CALL hsmt_spinor(3,n,noco,chi) 
     CALL hsmt_distspins(chi,hmat_tmp,hmat)

eigen/hsmt_nonsph.F90

@@ -201,6 +201,7 @@ CONTAINS
                 CALL ZHERK("U","N",lapw%nv(iintsp),ab_size,Rchi,CONJG(ab1),SIZE(ab1,1),1.0,hmat%data_c,SIZE(hmat%data_c,1))
              ELSE !This is the case of a local off-diagonal contribution.
                   !It is not Hermitian, so we need to USE zgemm CALL
+                CALL hsmt_ab(sym,atoms,noco,isp,2,n,na,cell,lapw,fj,gj,ab,ab_size,.TRUE.)
                 CALL zgemm("N","T",lapw%nv(iintsp),lapw%nv(jintsp),ab_size,chi,CONJG(ab),SIZE(ab,1),&
                      ab1,SIZE(ab1,1),CMPLX(1.0,0.0),hmat%data_c,SIZE(hmat%data_c,1))
              ENDIF

eigen/mt_setup.F90

@@ -41,7 +41,7 @@ CONTAINS
 
     DO jsp=1,MERGE(4,input%jspins,noco%l_mtNocoPot)
        !CALL tlmplm_cholesky(sphhar,atoms,DIMENSION,enpara, jsp,1,mpi,vTot%mt(:,0,1,jsp),input,vTot%mmpMat, td,ud)
-       CALL tlmplm_cholesky(sphhar,atoms,noco,enpara,jsp,jsp,mpi,vTot,input,td,ud)
+       CALL tlmplm_cholesky(sphhar,atoms,noco,enpara,jsp,mpi,vTot,input,td,ud)
        IF (input%l_f) CALL write_tlmplm(td,vTot%mmpMat,atoms%n_u>0,jsp,jsp,input%jspins)
     END DO
     CALL timestop("tlmplm")

eigen/tlmplm.F90

 MODULE m_tlmplm
-
+  use m_judft
   IMPLICIT NONE
   !*********************************************************************
   !     sets up the local Hamiltonian, i.e. the Hamiltonian in the
@@ -7,7 +7,7 @@ MODULE m_tlmplm
   !*********************************************************************
 CONTAINS
   SUBROUTINE tlmplm(n,sphhar,atoms,enpara,&
-       jspin,jsp,mpi,v,input,td,ud)
+       jspin,mpi,v,input,td,ud)
     USE m_constants
     USE m_intgr, ONLY : intgr3
     USE m_genMTBasis
@@ -25,7 +25,7 @@ CONTAINS
     TYPE(t_tlmplm),INTENT(INOUT) :: td
     TYPE(t_usdus),INTENT(INOUT)  :: ud
 
-    INTEGER, INTENT (IN) :: n,jspin,jsp !atom index,physical spin&spin index for data
+    INTEGER, INTENT (IN) :: n,jspin !atom index,physical spin&spin index for data
 
     REAL, ALLOCATABLE   :: dvd(:,:),dvu(:,:),uvd(:,:),uvu(:,:),f(:,:,:,:),g(:,:,:,:),x(:),flo(:,:,:)
     INTEGER,ALLOCATABLE :: indt(:)
@@ -35,7 +35,7 @@ CONTAINS
     COMPLEX  :: cil
     REAL     :: temp
     INTEGER i,l,l2,lamda,lh,lm,lmin,lmin0,lmp,lmpl,lmplm,lmx,lmxx,lp,info,in
-    INTEGER lp1,lpl ,mem,mems,mp,mu,nh,na,m,nsym,s,i_u,jspin1,jspin2
+    INTEGER lp1,lpl ,mem,mems,mp,mu,nh,na,m,nsym,s,i_u,jspin1,jspin2,jsp
 
     ALLOCATE( dvd(0:atoms%lmaxd*(atoms%lmaxd+3)/2,0:sphhar%nlhd ))
     ALLOCATE( dvu(0:atoms%lmaxd*(atoms%lmaxd+3)/2,0:sphhar%nlhd ))
@@ -48,7 +48,7 @@ CONTAINS
     ALLOCATE( vr0(SIZE(v%mt,1),0:SIZE(v%mt,2)-1))
 
 
-    
+    jsp=jspin
     vr0=v%mt(:,:,n,jsp)
     IF (jsp<3) vr0(:,0)=0.0
 
@@ -184,8 +184,9 @@ CONTAINS
     !--->   set up the t-matrices for the local orbitals,
     !--->   if there are any
     IF (atoms%nlo(n).GE.1) THEN
-       CALL tlo(atoms,sphhar,jspin,jsp,n,enpara,1,input,v%mt(1,0,n,jsp),&
-            na,flo,f(:,:,:,jspin),g(:,:,:,jspin),ud, ud%uuilon(:,:,jspin),ud%duilon(:,:,jspin),ud%ulouilopn(:,:,:,jspin), td)
+       IF (jspin>3) call judft_error("l_mtnocoPot=T and LOs not supported.")
+          CALL tlo(atoms,sphhar,jspin,jsp,n,enpara,1,input,v%mt(1,0,n,jsp),&
+               na,flo,f(:,:,:,jspin),g(:,:,:,jspin),ud, ud%uuilon(:,:,jspin),ud%duilon(:,:,jspin),ud%ulouilopn(:,:,:,jspin), td)
 
     ENDIF
   END SUBROUTINE tlmplm

eigen/tlmplm_cholesky.F90

@@ -9,7 +9,7 @@ MODULE m_tlmplm_cholesky
   !*********************************************************************
 CONTAINS
   SUBROUTINE tlmplm_cholesky(sphhar,atoms,noco,enpara,&
-       jspin,jsp,mpi,v,input,td,ud)
+       jspin,mpi,v,input,td,ud)
     USE m_tlmplm
     USE m_types
     IMPLICIT NONE
@@ -21,7 +21,7 @@ CONTAINS
     TYPE(t_enpara),INTENT(IN)   :: enpara
     !     ..
     !     .. Scalar Arguments ..
-    INTEGER, INTENT (IN) :: jspin,jsp !physical spin&spin index for data
+    INTEGER, INTENT (IN) :: jspin !physical spin&spin index for data
     !     ..
     TYPE(t_potden),INTENT(IN)    :: v
     TYPE(t_tlmplm),INTENT(INOUT) :: td
@@ -31,7 +31,7 @@ CONTAINS
     !     .. Local Scalars ..
     REAL temp
     INTEGER i,l,lm,lmin,lmin0,lmp,lmplm,lp,info,in
-    INTEGER lpl ,mp,n,m,s,i_u
+    INTEGER lpl ,mp,n,m,s,i_u, jsp
     LOGICAL OK
     !     ..
     !     .. Local Arrays ..
@@ -42,6 +42,7 @@ CONTAINS
 
 
     !     ..e_shift
+    jsp=jspin
     td%e_shift(:,jsp)=0.0
     IF (jsp<3) td%e_shift(:,jsp)=e_shift_min
 
@@ -59,7 +60,7 @@ CONTAINS
     !$OMP PRIVATE(OK,s,in,info)&
     !$OMP SHARED(atoms,jspin,jsp,sphhar,enpara,td,ud,v,mpi,input)
     DO  n = 1,atoms%ntype
-       CALL tlmplm(n,sphhar,atoms,enpara,jspin,jsp,mpi,v,input,td,ud)
+       CALL tlmplm(n,sphhar,atoms,enpara,jspin,mpi,v,input,td,ud)
        OK=.FALSE.
        cholesky_loop:DO WHILE(.NOT.OK)
           td%h_loc(:,:,n,jsp)=0.0