Initial release of FLEUR v0.27

Based on v0.26 hybrid branch

CMakeLists.txt

+cmake_minimum_required(VERSION 2.8)
+project(FLEUR LANGUAGES NONE)
+
+include("cmake/Architectures.txt")
+include("cmake/ReportConfig.txt")
+
+#Here the targets and the files are defined
+include("cmake/Files_and_Targets.txt")
+
+#install(TARGETS fleur inpgen DESTINATION bin)
+install(PROGRAMS ${CMAKE_BINARY_DIR}/fleur 
+        CONFIGURATIONS Debug
+        DESTINATION bin
+        RENAME fleur_dbg)

README

-Currently, there are no releases or stable commits. Please checkout the development branch by running:
+Welcome to the source code of FLEUR
 
-git checkout -t origin/develop
+Outline:
+1. Short Overview of the code
+2. Installation
+
+----------------------------------------------
+1. Overview of the 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
+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 (see Installation)
+
+----------------------------------------------
+2. Installation
+----------------------------------------------
+
+For the compilation of FLEUR you will need:
+
+- Fortran compiler (should support Fortran2003)
+- Compatible C compiler for the XML-IO
+- cmake for controlling the make-process
+- Libraries:
+  * LAPACK/BLAS
+  * MPI for parallelization
+  * SCALAPACK or ELPA or ELEMENTAL if you want to use parallel diagonalization
+    	      (without you can only exploit k-point parallelism)
+  * HDF5 for parallel IO (if you have lots of memory, you might not need to do IO :-)
+
+You should probably create a directory structure like this:
+
+somewhere/src -this is where the source code lives, i.e. the location of this README
+somewhere/build -here you can build the code
+
+change to the build directory and type:
+
+cmake ../src
+make 
+
+This might generate the FLEUR executables. If you get errors, you have to adjust cmake
+so that your build system is correctly set-up. To do so, please adjust cmake.config. You might
+find inspiration in the cmake.???.config files provided for several systems.
+
+If your build environment is recognized correctly you can obtain the following executables:
+inpgen      - the input generator of FLEUR
+fleur       - general version of FLEUR
+fleur_INVS  - version of FLEUR suitable for systems with inversion symmetry leading to real matrices
+fleur_SOC   - version of FLEUR with spin-orbit coupling
+
+All codes might also have a _MPI attached to indicate parallel versions.
\ No newline at end of file

TODO

+List of BUGS/ things not tested or implemented:
+-omp bug in vmtxc oder xcall
+-MPI version not compiling
+-DOS IO missing
+-Wannier part
+-Hybrid functionals
+-vdW
+-one-dimensional code
+

cdn/CMakeLists.txt

+set(fleur_F77 ${fleur_F77}
+cdn/diflgr.f
+cdn/od_chirot.f
+cdn/rcerf.f
+cdn/rot_den_mat.f
+)
+set(fleur_F90 ${fleur_F90}
+cdn/cdnovlp.F90
+cdn/cdntot.f90
+cdn/cdnval.F90
+cdn/eparas.f90
+cdn/int_21.f90
+cdn/int_21lo.f90
+cdn/m_perp.f90
+cdn/n_mat.f90
+cdn/od_abvac.f90
+cdn/prp_qfft_map.f90
+cdn/pwden.F90
+cdn/pwint.f90
+cdn/pwint_sl.f90
+cdn/q_int_sl.F90
+cdn/q_mt_sl.f90
+cdn/qal_21.f90
+cdn/qpw_to_nmt.f90
+cdn/slab_dim.f90
+cdn/slabgeom.f90
+cdn/vacden.F90
+)

cdn/cdntot.f90

+      MODULE m_cdntot
+!     ********************************************************
+!     calculate the total charge density in the interstial.,
+!     vacuum, and mt regions      c.l.fu
+!     ********************************************************
+      CONTAINS
+      SUBROUTINE cdntot(&
+     &                  stars,atoms,sym,&
+     &                  vacuum,input,cell,oneD,&
+     &                  qpw,rho,rht,&
+     &                  qtot,qistot)
+
+      USE m_intgr, ONLY : intgr3
+      USE m_constants
+      USE m_qsf
+      USE m_pwint
+      USE m_types
+      use m_juDFT
+      IMPLICIT NONE
+!     ..
+!     .. Scalar Arguments ..
+      TYPE(t_stars),INTENT(IN) :: stars
+      TYPE(t_atoms),INTENT(IN) :: atoms
+      TYPE(t_sym),INTENT(IN)   :: sym
+      TYPE(t_vacuum),INTENT(IN):: vacuum
+      TYPE(t_input),INTENT(IN) :: input
+      TYPE(t_oneD),INTENT(IN)  :: oneD
+      TYPE(t_cell),INTENT(IN)  :: cell
+      REAL,    INTENT (OUT):: qtot,qistot
+!     ..
+!     .. Array Arguments ..
+      COMPLEX, INTENT (IN) :: qpw(stars%n3d,input%jspins)
+      REAL,    INTENT (IN) :: rho(:,0:,:,:) !(atoms%jmtd,0:sphhar%nlhd,atoms%ntypd,input%jspins)
+      REAL,    INTENT (IN) :: rht(vacuum%nmzd,2,input%jspins)
+!-odim
+!+odim
+!     ..
+!     .. Local Scalars ..
+    ! COMPLEX x
+      COMPLEX x(stars%ng3)
+      REAL q,qis,w
+      INTEGER i,ivac,j,jspin,n,nz
+!     ..
+!     .. Local Arrays ..
+      REAL qmt(atoms%ntypd),qvac(2),q2(vacuum%nmz),rht1(vacuum%nmzd,2,input%jspins)
+!     ..
+!     .. Intrinsic Functions ..
+      INTRINSIC real
+!     ..
+!
+      CALL timestart("cdntot")
+      qtot = 0.e0
+      qistot = 0.e0
+      DO 40 jspin = 1,input%jspins
+         q = 0.e0
+!     -----mt charge
+         CALL timestart("MT")
+         DO 10 n = 1,atoms%ntype
+            CALL intgr3(rho(:,0,n,jspin),atoms%rmsh(:,n),atoms%dx(n),atoms%jri(n),w)
+            qmt(n) = w*sfp_const
+            q = q + atoms%neq(n)*qmt(n)
+   10    CONTINUE
+         CALL timestop("MT")
+!     -----vacuum region
+         IF (input%film) THEN
+            DO 20 ivac = 1,vacuum%nvac
+               DO nz = 1,vacuum%nmz
+                  IF (oneD%odi%d1) THEN
+                     rht1(nz,ivac,jspin) = (cell%z1+(nz-1)*vacuum%delz)*&
+     &                    rht(nz,ivac,jspin)
+                  ELSE
+                     rht1(nz,ivac,jspin) =  rht(nz,ivac,jspin)
+                  END IF
+               END DO
+               CALL qsf(vacuum%delz,rht1(1,ivac,jspin),q2,vacuum%nmz,0)
+               qvac(ivac) = q2(1)*cell%area
+               IF (.NOT.oneD%odi%d1) THEN
+                  q = q + qvac(ivac)*2./real(vacuum%nvac)
+               ELSE
+                  q = q + cell%area*q2(1)
+               END IF
+   20       CONTINUE
+         END IF
+!     -----is region
+         qis = 0.
+!         DO 30 j = 1,nq3
+!            CALL pwint(
+!     >                 k1d,k2d,k3d,n3d,ntypd,natd,nop,invtab,odi,
+!     >                 ntype,neq,volmts,taual,z1,vol,volint,
+!     >                 symor,tau,mrot,rmt,sk3,bmat,ig2,ig,
+!     >                 kv3(1,j),
+!     <                 x)
+!            qis = qis + qpw(j,jspin)*x*nstr(j)
+!   30    CONTINUE
+         CALL pwint_all(&
+     &                 stars,atoms,sym,oneD,&
+     &                 cell,&
+     &                 x)
+         DO j = 1,stars%ng3
+             qis = qis + qpw(j,jspin)*x(j)*stars%nstr(j)
+         ENDDO
+         qistot = qistot + qis
+         q = q + qis
+         WRITE (6,FMT=8000) jspin,q,qis, (qmt(n),n=1,atoms%ntype)
+         IF (input%film) WRITE (6,FMT=8010) (i,qvac(i),i=1,vacuum%nvac)
+         WRITE (16,FMT=8000) jspin,q,qis, (qmt(n),n=1,atoms%ntype)
+         IF (input%film) WRITE (16,FMT=8010) (i,qvac(i),i=1,vacuum%nvac)
+         qtot = qtot + q
+   40 CONTINUE
+      WRITE (6,FMT=8020) qtot
+      WRITE (16,FMT=8020) qtot
+ 8000 FORMAT (/,10x,'total charge for spin',i3,'=',f12.6,/,10x,&
+     &       'interst. charge =   ',f12.6,/,&
+     &       (10x,'mt charge=          ',4f12.6,/))
+ 8010 FORMAT (10x,'vacuum ',i2,'  charge=  ',f12.6)
+ 8020 FORMAT (/,10x,'total charge  =',f12.6)
+
+      CALL timestop("cdntot")
+      END SUBROUTINE cdntot
+      END MODULE m_cdntot

cdn/diflgr.f

+      MODULE m_diflgr
+      CONTAINS
+      REAL FUNCTION diflgr(x,f)
+c     **********************************************************
+c     differentiate the function f,
+c     given at the points x0,x1,x2,
+c     at the point x1 by lagrange interpolation
+c     e.wimmer
+c     ***********************************************************
+C     .. Array Arguments ..
+
+      IMPLICIT NONE
+      REAL f(0:2),x(0:2)
+C     ..
+      diflgr = (x(1)-x(2))*f(0)/ ((x(0)-x(1))* (x(0)-x(2))) +
+     +         (2.0e0*x(1)-x(2)-x(0))*f(1)/ ((x(1)-x(0))* (x(1)-x(2))) +
+     +         (x(1)-x(0))*f(2)/ ((x(2)-x(0))* (x(2)-x(1)))
+c     print *,'x=',x
+c     print *,'f=',f
+c     print *,'diflgr=',diflgr
+      RETURN
+      END FUNCTION
+      END

cdn/eparas.f90

+MODULE m_eparas
+  !***********************************************************************
+  ! Calculates qlo, enerlo and sqlo, which are needed to determine the 
+  ! new energy parameters.
+  ! Philipp Kurz 99/04
+  !***********************************************************************
+  ! also the 'normal' energy parameters are now included...
+  !
+  ! if (l_mcd) then mcd contains mcd spectrum: first index = polarization
+  ! second = core level ; third = band index                  gb.2001
+  ! corrected to work also for multiple LO's of same l at the same atom
+  !                                                           gb.2005
+  !*************** ABBREVIATIONS *****************************************
+  ! qlo     : charge density of one local orbital at the current k-point
+  ! sqlo    : qlo integrated over the Brillouin zone
+  ! enerlo  : qlo*energy integrated over the Brillouin zone
+  !***********************************************************************
+  !
+CONTAINS
+  SUBROUTINE eparas(jsp,atoms,noccbd, mpi,ikpt,ne,we,eig,ccof, skip_t,l_evp,acof,bcof,&
+       usdus, ncore,l_mcd,m_mcd, enerlo,sqlo,ener,sqal,qal,mcd)
+    USE m_types
+    IMPLICIT NONE
+    TYPE(t_usdus),INTENT(IN)   :: usdus
+    TYPE(t_mpi),INTENT(IN)     :: mpi
+    TYPE(t_atoms),INTENT(IN)   :: atoms
+    !     ..
+    !     .. Scalar Arguments ..
+    INTEGER, INTENT (IN) :: noccbd,jsp     
+    INTEGER, INTENT (IN) :: ne,ikpt  ,skip_t
+    LOGICAL, INTENT (IN) :: l_mcd,l_evp
+    !     ..
+    !     .. Array Arguments ..
+    INTEGER, INTENT (IN)  :: ncore(atoms%ntypd)
+    REAL,    INTENT (IN)  :: eig(:)!(dimension%neigd),
+    REAL,    INTENT (IN)  :: we(noccbd) 
+    COMPLEX, INTENT (IN)  :: ccof(-atoms%llod:atoms%llod,noccbd,atoms%nlod,atoms%natd)
+    COMPLEX, INTENT (IN)  :: acof(:,0:,:)!(noccbd,0:dimension%lmd,atoms%natd)
+    COMPLEX, INTENT (IN)  :: bcof(:,0:,:)!(noccbd,0:dimension%lmd,atoms%natd)
+    COMPLEX, INTENT (IN)  :: m_mcd(:,:,:,:)!(dimension%nstd,(3+1)**2,3*ntypd ,2)
+    REAL,    INTENT (OUT) :: enerlo(atoms%nlod,atoms%ntypd),sqlo(atoms%nlod,atoms%ntypd)
+    REAL,    INTENT (OUT) :: ener(0:3,atoms%ntypd),sqal(0:3,atoms%ntypd)
+    REAL,    INTENT (OUT) :: qal(0:,:,:)!(0:3,atoms%ntypd,dimension%neigd)
+    REAL,    INTENT (OUT) :: mcd(:,:,:)!(3*atoms%ntypd,dimension%nstd,dimension%neigd)
+
+    !     ..
+    !     .. Local Scalars ..
+    INTEGER i,l,lo,lop ,natom,nn,ntyp,m
+    INTEGER nt1,nt2,lm,n,ll1,ipol,icore,index
+    REAL fac
+    COMPLEX suma,sumb,sumab,sumba
+    !     ..
+    !     .. Local Arrays ..
+    REAL qlo(noccbd,atoms%nlod,atoms%nlod,atoms%ntypd)
+    REAL qaclo(noccbd,atoms%nlod,atoms%ntypd),qbclo(noccbd,atoms%nlod,atoms%ntypd)
+    !     ..
+    !
+    !---> initialize ener, sqal, enerlo and sqlo on first call
+    !
+
+    IF ((ikpt.LE.mpi%isize).AND..NOT.l_evp) THEN
+       IF (l_mcd) THEN
+          mcd(:,:,:) = 0.0
+       ENDIF
+       ener(:,:) = 0.0
+       sqal(:,:) = 0.0
+       qal(:,:,:) = 0.0
+       enerlo(:,:) = 0.0
+       sqlo(:,:) = 0.0
+    END IF
+    !
+    !--->    l-decomposed density for each occupied state
+    !
+    !         DO 140 i = (skip_t+1),ne    ! this I need for all states
+    DO i = 1,ne              ! skip in next loop
+       nt1 = 1
+       DO n = 1,atoms%ntype
+          fac = 1./atoms%neq(n)
+          nt2 = nt1 + atoms%neq(n) - 1
+          DO l = 0,3
+             suma = CMPLX(0.,0.)
+             sumb = CMPLX(0.,0.)
+             ll1 = l* (l+1)
+             DO m = -l,l
+                lm = ll1 + m
+                IF ( .NOT.l_mcd ) THEN
+                   DO natom = nt1,nt2
+                      suma = suma + acof(i,lm,natom)*CONJG(acof(i,lm,natom))
+                      sumb = sumb + bcof(i,lm,natom)*CONJG(bcof(i,lm,natom))
+                   ENDDO
+                ELSE
+                   suma = CMPLX(0.,0.) ; sumab = CMPLX(0.,0.) 
+                   sumb = CMPLX(0.,0.) ; sumba = CMPLX(0.,0.)
+                   DO natom = nt1,nt2
+                      suma = suma + acof(i,lm,natom)*CONJG(acof(i,lm,natom))
+                      sumb = sumb + bcof(i,lm,natom)*CONJG(bcof(i,lm,natom))
+                      sumab= sumab + acof(i,lm,natom) *CONJG(bcof(i,lm,natom))
+                      sumba= sumba + bcof(i,lm,natom) *CONJG(acof(i,lm,natom))
+                   ENDDO
+                   DO icore = 1, ncore(n)
+                      DO ipol = 1, 3
+                         index = 3*(n-1) + ipol
+                         mcd(index,icore,i)=mcd(index,icore,i) + fac*(&
+                              suma * CONJG(m_mcd(icore,lm+1,index,1))*m_mcd(icore,lm+1,index,1)  +&
+                              sumb * CONJG(m_mcd(icore,lm+1,index,2))*m_mcd(icore,lm+1,index,2)  +&
+                              sumab* CONJG(m_mcd(icore,lm+1,index,2))*m_mcd(icore,lm+1,index,1)  +&
+                              sumba* CONJG(m_mcd(icore,lm+1,index,1))*m_mcd(icore,lm+1,index,2)  ) 
+                      ENDDO
+                   ENDDO
+                ENDIF     ! end MCD
+             ENDDO
+             qal(l,n,i) = (suma+sumb*usdus%ddn(l,n,jsp))/atoms%neq(n)
+          ENDDO
+          nt1 = nt1 + atoms%neq(n)
+       ENDDO
+    ENDDO
+    !
+    !--->    perform Brillouin zone integration and summation over the
+    !--->    bands in order to determine the energy parameters for each
+    !--->    atom and angular momentum
+    !
+    DO l = 0,3
+       DO n = 1,atoms%ntype
+          DO i = (skip_t+1),noccbd
+             ener(l,n) = ener(l,n) + qal(l,n,i)*we(i)*eig(i)
+             sqal(l,n) = sqal(l,n) + qal(l,n,i)*we(i)
+          ENDDO
+       ENDDO
+    ENDDO
+
+    !---> initialize qlo
+
+    qlo(:,:,:,:) = 0.0
+    qaclo(:,:,:) = 0.0
+    qbclo(:,:,:) = 0.0
+
+    !---> density for each local orbital and occupied state
+
+    natom = 0
+    DO ntyp = 1,atoms%ntype
+       DO nn = 1,atoms%neq(ntyp)
+          natom = natom + 1
+          DO lo = 1,atoms%nlo(ntyp)
+             l = atoms%llo(lo,ntyp)
+             ll1 = l* (l+1)
+             DO m = -l,l
+                lm = ll1 + m
+                DO i = 1,ne
+                   qbclo(i,lo,ntyp) = qbclo(i,lo,ntyp) +REAL(&
+                        bcof(i,lm,natom)*CONJG(ccof(m,i,lo,natom))+ccof(m,i,lo,natom)*CONJG(bcof(i,lm,natom)) )
+                   qaclo(i,lo,ntyp) = qaclo(i,lo,ntyp) + REAL(&
+                        acof(i,lm,natom)*CONJG(ccof(m,i,lo,natom))+ccof(m,i,lo,natom)*CONJG(acof(i,lm,natom)) )
+                ENDDO
+             ENDDO
+             DO lop = 1,atoms%nlo(ntyp)
+                IF (atoms%llo(lop,ntyp).EQ.l) THEN
+                   DO m = -l,l
+                      DO i = 1,ne
+                         qlo(i,lop,lo,ntyp) = qlo(i,lop,lo,ntyp) +  REAL(&
+                              CONJG(ccof(m,i,lop,natom))*ccof(m,i,lo,natom))
+                      ENDDO
+                   ENDDO
+                ENDIF
+             ENDDO
+          ENDDO
+       ENDDO
+    ENDDO
+
+    !---> perform brillouin zone integration and sum over bands
+
+    DO ntyp = 1,atoms%ntype
+       DO lo = 1,atoms%nlo(ntyp)
+          l = atoms%llo(lo,ntyp)
+          ! llo > 3 used for unoccupied states only
+          IF( l .GT. 3 ) CYCLE
+          DO i = 1,ne
+             qal(l,ntyp,i)= qal(l,ntyp,i)  + ( 1.0/atoms%neq(ntyp) )* (&
+                  qaclo(i,lo,ntyp)*usdus%uulon(lo,ntyp,jsp)+qbclo(i,lo,ntyp)*usdus%dulon(lo,ntyp,jsp)     )
+          END DO
+          DO lop = 1,atoms%nlo(ntyp)
+             IF (atoms%llo(lop,ntyp).EQ.l) THEN
+                DO i = 1,ne
+                   enerlo(lo,ntyp) = enerlo(lo,ntyp) +qlo(i,lop,lo,ntyp)*we(i)*eig(i)
+                   sqlo(lo,ntyp) = sqlo(lo,ntyp) + qlo(i,lop,lo,ntyp)*we(i)
+                   qal(l,ntyp,i)= qal(l,ntyp,i)  + ( 1.0/atoms%neq(ntyp) ) *&
+                        qlo(i,lop,lo,ntyp)*usdus%uloulopn(lop,lo,ntyp,jsp)
+                ENDDO
+             ENDIF
+          ENDDO
+       END DO
+    END DO
+
+  END SUBROUTINE eparas
+END MODULE m_eparas

cdn/int_21.f90

+MODULE m_int21
+  !-----------------------------------------------------------
+  !
+  ! Integrates f(:,:,l,s) and g(:,:,l,s) (s=1,2) for given l
+  ! and itype with different spins (s).
+  ! Output is ..n21(l,itype), where .. is a (u,d) combination 
+  ! dependet on the (f,g) combination used. 
+  !
+  !-----------------------------------------------------------
+CONTAINS
+  SUBROUTINE int_21(f,g,atoms,ityp,l, uun21,udn21,dun21,ddn21)
+    
+    USE m_intgr, ONLY : intgr3
+    USE m_types
+    IMPLICIT NONE
+    TYPE(t_atoms),INTENT(IN)   :: atoms
+    !     ..
+    !     .. Scalar Arguments ..
+    INTEGER, INTENT (IN) :: l,ityp
+    REAL,    INTENT (OUT):: uun21,udn21,dun21,ddn21
+    !     ... Array Arguments
+    REAL,    INTENT (IN) :: f(:,:,0:,:)!(atoms%jmtd,2,0:atoms%lmaxd,dimension%jspd)
+    REAL,    INTENT (IN) :: g(:,:,0:,:)!(atoms%jmtd,2,0:atoms%lmaxd,dimension%jspd)
+    !     ...local arrays
+    REAL        uu_tmp(atoms%jri(ityp))
+
+
+    uu_tmp(:atoms%jri(ityp)) = f(:atoms%jri(ityp),1,l,2)*f(:atoms%jri(ityp),1,l,1)&
+         + f(:atoms%jri(ityp),2,l,2)*f(:atoms%jri(ityp),2,l,1)
+    CALL intgr3(uu_tmp,atoms%rmsh(:,ityp),atoms%dx(ityp),atoms%jri(ityp),uun21)
+    
+    uu_tmp(:atoms%jri(ityp)) = f(:atoms%jri(ityp),1,l,2)*g(:atoms%jri(ityp),1,l,1)&
+         + f(:atoms%jri(ityp),2,l,2)*g(:atoms%jri(ityp),2,l,1)
+    CALL intgr3(uu_tmp,atoms%rmsh(:,ityp),atoms%dx(ityp),atoms%jri(ityp),udn21)
+    
+    uu_tmp(:atoms%jri(ityp)) = g(:atoms%jri(ityp),1,l,2)*f(:atoms%jri(ityp),1,l,1)&
+         + g(:atoms%jri(ityp),2,l,2)*f(:atoms%jri(ityp),2,l,1)
+    CALL intgr3(uu_tmp,atoms%rmsh(:,ityp),atoms%dx(ityp),atoms%jri(ityp),dun21)
+    
+    uu_tmp(:atoms%jri(ityp)) = g(:atoms%jri(ityp),1,l,2)*g(:atoms%jri(ityp),1,l,1)&
+         + g(:atoms%jri(ityp),2,l,2)*g(:atoms%jri(ityp),2,l,1)
+    CALL intgr3(uu_tmp,atoms%rmsh(:,ityp),atoms%dx(ityp),atoms%jri(ityp),ddn21)
+    
+  END SUBROUTINE int_21
+END MODULE m_int21

cdn/int_21lo.f90

+MODULE m_int21lo
+  !-----------------------------------------------------------
+  !
+  ! Integrates combinations of flo(:,:,llo,s) with itself, f(:,:,l,s) 
+  ! and g(:,:,l,s) (s=1,2) for given llo and itype with different 
+  ! spins (s).
+  ! Output is ..n21(l,itype), where .. is a combination of (u,d) and
+  ! ulo dependent on the (f,g) combination used. Also ..n12 and 
+  ! uloulopn21 are calculated.
+  !
+  !-----------------------------------------------------------
+CONTAINS
+  SUBROUTINE int_21lo(f,g,atoms,n, flo,ilo, uulon21,dulon21,uulon12,dulon12,uloulopn21)
+
+    USE m_intgr, ONLY : intgr3
+    USE m_types
+    USE m_types
+    IMPLICIT NONE
+    TYPE(t_atoms),INTENT(IN)   :: atoms
+    !     ..
+    !     .. Scalar Arguments ..
+    INTEGER, INTENT (IN) :: ilo,n
+    !     ... Array Arguments
+    REAL,    INTENT (IN) :: f(:,:,0:,:)!(atoms%jmtd,2,0:atoms%lmaxd,dimension%jspd)
+    REAL,    INTENT (IN) :: g(:,:,0:,:)!(atoms%jmtd,2,0:atoms%lmaxd,dimension%jspd)
+    REAL,    INTENT (IN) :: flo(:,:,:,:)!(atoms%jmtd,2,atoms%nlod,dimension%jspd)
+    REAL,    INTENT (OUT):: uulon21,uulon12
+    REAL,    INTENT (OUT):: dulon21,dulon12
+    REAL,    INTENT (OUT):: uloulopn21(atoms%nlod,atoms%nlod)
+
+    !     ...local scalars
+    INTEGER iri,l,lp,ilop
+    !     ...local arrays
+    REAL    uu_tmp(atoms%jri(n))
+
+    !
+    ! --> norm of product of u and ulo:
+    !
+    l = atoms%llo(ilo,n)
+    DO iri = 1, atoms%jri(n)
+       uu_tmp(iri) = f(iri,1,l,2)*flo(iri,1,ilo,1)+ f(iri,2,l,2)*flo(iri,2,ilo,1)
+    ENDDO
+    CALL intgr3(uu_tmp,atoms%rmsh(:,n),atoms%dx(n),atoms%jri(n),uulon21)
+    DO iri = 1, atoms%jri(n)
+       uu_tmp(iri) = f(iri,1,l,1)*flo(iri,1,ilo,2)+ f(iri,2,l,1)*flo(iri,2,ilo,2)
+    ENDDO
+    CALL intgr3(uu_tmp,atoms%rmsh(:,n),atoms%dx(n),atoms%jri(n),uulon12)
+    !
+    ! --> norm of product of du and ulo:
+    !
+    DO iri = 1, atoms%jri(n)
+       uu_tmp(iri) = g(iri,1,l,2)*flo(iri,1,ilo,1) + g(iri,2,l,2)*flo(iri,2,ilo,1)