Commit c74e6d87 authored by Daniel Wortmann's avatar Daniel Wortmann

Merge branch 'develop' into release

parents af9fa5b2 8429d14c

Too many changes to show.

To preserve performance only 260 of 260+ files are displayed.

init/compileinfo.h
*~
#*
build
*.o
*.mod
*.x
*.swp
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)
The MIT License (MIT)
Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"),
to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute,
sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
IN THE SOFTWARE.
Currently, there are no releases or stable commits. Please checkout the development branch by running:
Welcome to the source code of FLEUR
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)
If you modify FLEUR please do so in the develop branch by running
git checkout -t origin/develop
after cloning the git repository.
----------------------------------------------
2. Installation
----------------------------------------------
For the compilation of FLEUR you will need:
- Fortran compiler (should support Fortran2008)
- 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. It most probable will only work on systems we know as
only for those there will be specific configurations in the cmake directory.
If you get errors, you have to configure your build system yourself.
There are two options to do that.
- Either set the environment variable FC to either ifort,pgfortran or gfortran.
Then the configuration in the cmake directory will be used from cmake.[ifort|gfortran|pgfortran].config will be used that can again work for your system or can be adjusted.
- Or you run cmake with the option -DFleur_custom_toolchain and you adjust cmake.config to your needs.
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
List of BUGS/ things not tested or implemented:
-omp bug in vmtxc oder xcall
-Wannier part
-Hybrid functionals
-one-dimensional code
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
)
!--------------------------------------------------------------------------------
! Copyright (c) 2016 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_cdnovlp
USE m_juDFT
IMPLICIT NONE
PRIVATE
PUBLIC :: cdnovlp
CONTAINS
SUBROUTINE cdnovlp(mpi,&
& sphhar,stars,atoms,sym,&
& DIMENSION,vacuum,cell,&
& input,oneD,l_st,&
& jspin,rh,&
& qpw,rhtxy,rho,rht)
!*****************************************************************
! calculates the overlapping core tail density and adds
! its contribution to the corresponging components of
! valence density.
!
! OLD VERSION:
! The previous version to calculate the overlap of the
! core tails was done in real space:
! A three dimensional real space grid was set up. On each
! of these grid points the charge density of all atoms inside
! the unit cell and neigboring unit cells was calculated.
! This calculation required a lagrange fit since the core
! charge is on a radial mesh. The same was done in the vacuum
! region, except on each z-plane we worked on a two dimension
! grid. After the charge density was generated on a equidistant
! grid the charge density was FFTed into G-space. The set up
! of the charge density on the real space grid is rather time
! consuming. 3-loops are required for the 3D grid
! 1-loop over all atoms in the unit cells
! Larange interpolation
! 3D and 2D FFTs
! In order to save time the non-spherical contributions inside
! the sphere had been ignored. It turns out that the later
! approximation is pure in the context of force calculations.
!
! PRESENT VERSION:
! The present version is written from scratch. It is based on the
! idea that the FFT of an overlap of spherically symmetric
! charges can be expressed by the product of
!
! sum_natype{ F(G,ntype) * sum_atom(atype) {S(\vec{G},atom)}}
!
! of form factor F and structure factor S. The Form factor
! depends only G while the structure factor depends on \vec{G}
! and can build up in G-space. F of a gaussian chargedensity can
! be calculated analytically.
! The core-tails to the vacuum region are described by an
! exponentially decaying function into the vacuum:
! rho(r_||,z,ivac)= sum_n{ rho(n,ivac) * exp(-kappa*(z-z_v))
! * exp(iG(n)r_||) }
! And the plane waves are expanded into lattice harmonics
! up to a l_cutoff. Tests of the accuracy inside the sphere
! have shown that a reduction of the l_cutoff inside the
! in order to save time leads to sizable errors and should
! be omitted.
! rho_L(r) = 4 \pi i^l \sum_{g =|= 0} \rho_int(g) r_i^{2} \times
! j_{l} (gr_i) \exp{ig\xi_i} Y^*_{lm} (g)
! Tests have shown that the present version is about 10 times
! faster than the previous one also all nonspherical terms are
! included up to l=8 and the previous one included only l=0.
! For l=16 it is still faster by factor 2.
! coded Stefan Bl"ugel, IFF Nov. 1997
! tested RObert Abt , IFF Dez. 1997
!*****************************************************************
!
USE m_constants
USE m_qpwtonmt
USE m_cylbes
USE m_dcylbs
USE m_od_cylbes
USE m_diflgr
USE m_types
#ifdef CPP_MPI
USE m_mpi_bc_st
#endif
!
! .. Parameters ..
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_dimension),INTENT(IN)::DIMENSION
TYPE(t_vacuum),INTENT(in):: vacuum
TYPE(t_input),INTENT(in)::input
INTEGER method1, method2
PARAMETER (method1 = 2, method2 = 1)
! ..
! .. Scalar Arguments ..
INTEGER,INTENT (IN) :: jspin
LOGICAL,INTENT (IN) :: l_st
! ..
! .. Array Arguments ..
COMPLEX,INTENT (INOUT) :: qpw(stars%n3d,DIMENSION%jspd)
COMPLEX,INTENT (INOUT) :: rhtxy(vacuum%nmzxyd,oneD%odi%n2d-1,2,DIMENSION%jspd)
REAL, INTENT (INOUT) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntypd,DIMENSION%jspd)
REAL, INTENT (INOUT) :: rht(vacuum%nmzd,2,DIMENSION%jspd)
REAL, INTENT (INOUT) :: rh(DIMENSION%msh,atoms%ntypd)
! ..
! .. Local Scalars ..
COMPLEX czero,carg,VALUE,slope,ci
REAL dif,dxx,g,gz,dtildh,&
& rkappa,sign,signz,tol_14,z,zero,zvac,&
& g2,phi,gamma,qq
INTEGER ig3,imz,ivac,j,j1,k,kz,k1,k2,l_cutoff,m0,&
& n,nz,nrz,nzvac,&
& irec2,irec3,irec1,m,gzi
! ..
! .. Local Arrays ..
COMPLEX, ALLOCATABLE :: qpwc(:)
REAL acoff(atoms%ntypd),alpha(atoms%ntypd),rho_out(2)
REAL rat(DIMENSION%msh,atoms%ntypd)
INTEGER mshc(atoms%ntypd)
REAL fJ(-oneD%odi%M:oneD%odi%M),dfJ(-oneD%odi%M:oneD%odi%M)
! ..
DATA czero /(0.0,0.0)/, zero /0.0/, tol_14 /1.0e-10/!-14
#ifdef CPP_MPI
EXTERNAL MPI_BCAST
INTEGER ierr
#include "cpp_double.h"
INCLUDE "mpif.h"
#endif
!
!----> Abbreviation
!
! l_cutoff : cuts off the l-expansion of the nonspherical charge
! density due to the core-tails of neigboring atoms
! mshc : maximal radial meshpoint for which the radial coretail
! density is larger than tol_14
! method1 : two different ways to calculate the derivative of the
! charge density at the sphere boundary.
! (1) use subroutine diflgr based on lagrange interpol.
! (2) use two point formular in real space,
! see notes of SB.
! Tests have shown that (2) is more accurate.
! method2 : two different integration routines to calculate form
! factor of coretails outside the sphere.
! (1) use subroutine intgrz to integrate the tails from
! outside to inside.
! (2) use subroutine intgr3 to integrate the tails from
! muffin-tin radius to outside and include correction
! for start up.
! Tests have shown that (1) is more accurate.
!
!
ci = CMPLX(0.0,1.0)
ALLOCATE (qpwc(stars%n3d))
!
!----> prepare local array to store pw-expansion of pseudo core charge
!
DO k = 1 , stars%ng3
qpwc(k) = czero
ENDDO
!
!----> (1) set up radial mesh beyond muffin-tin radius
! (2) cut_off core tails from noise
!
#ifdef CPP_MPI
CALL MPI_BCAST(rh,DIMENSION%msh*atoms%ntypd,CPP_MPI_REAL,0,mpi%mpi_comm,ierr)
#endif
nloop: DO n = 1 , atoms%ntype
IF ((atoms%ncst(n).GT.0).OR.l_st) THEN
DO j = 1 , atoms%jri(n)
rat(j,n) = atoms%rmsh(j,n)
ENDDO
dxx = EXP(atoms%dx(n))
DO j = atoms%jri(n) + 1 , DIMENSION%msh
rat(j,n) = rat(j-1,n)*dxx
ENDDO
DO j = atoms%jri(n) - 1 , DIMENSION%msh
rh(j,n) = rh(j,n)/ (fpi_const*rat(j,n)*rat(j,n))
ENDDO
DO j = DIMENSION%msh , atoms%jri(n) , -1
IF ( rh(j,n) .GT. tol_14 ) THEN
mshc(n) = j
CYCLE nloop
END IF
ENDDO
mshc(n) = atoms%jri(n)
ENDIF
ENDDO nloop
!
!-----> the core density inside the spheres is replaced by a
! gaussian-like pseudo density : n(r) = acoff*exp(-alpha*r*r)
! acoff and alpha determined to obtain a continous and
! differentiable density at the sphere boundary.
! IF mshc = jri either core tail too small or no core (i.e. H)
!
DO n = 1,atoms%ntype
IF ((mshc(n).GT.atoms%jri(n)).AND.((atoms%ncst(n).GT.0).OR.l_st)) THEN
j1 = atoms%jri(n) - 1
IF ( method1 .EQ. 1) THEN
dif = diflgr(rat(j1,n),rh(j1,n))
WRITE (6,FMT=8000) n,rh(atoms%jri(n),n),dif
alpha(n) = -0.5 * dif / ( rh(atoms%jri(n),n)*atoms%rmt(n) )
ELSEIF ( method1 .EQ. 2) THEN
alpha(n) = LOG( rh(j1,n) / rh(atoms%jri(n),n) )
alpha(n) = alpha(n)&
& / ( atoms%rmt(n)*atoms%rmt(n)*( 1.0-EXP( -2.0*atoms%dx(n) ) ) )
ELSE
WRITE (6,'('' error in choice of method1 in cdnovlp '')')
CALL juDFT_error("error in choice of method1 in cdnovlp"&
& ,calledby ="cdnovlp")
ENDIF
acoff(n) = rh(atoms%jri(n),n) * EXP( alpha(n)*atoms%rmt(n)*atoms%rmt(n) )
WRITE (6,FMT=8010) alpha(n),acoff(n)
DO j = 1,atoms%jri(n) - 1
rh(j,n) = acoff(n) * EXP( -alpha(n)*rat(j,n)**2 )
ENDDO
ELSE
alpha(n) = 0.0
ENDIF
ENDDO
!
IF (mpi%irank ==0) THEN
8000 FORMAT (/,10x,'core density and its first derivative ',&
& 'at sph. bound. for atom type',&
& i2,' is',3x,2e15.7)
8010 FORMAT (/,10x,'alpha=',f10.5,5x,'acoff=',f10.5)
END IF
!
!=====> calculate the fourier transform of the core-pseudocharge
CALL ft_of_CorePseudocharge(mpi,DIMENSION,atoms,mshc,alpha,tol_14,rh, &
acoff,stars,method2,rat,cell,oneD,sym,qpwc)
DO k = 1 , stars%ng3
qpw(k,jspin) = qpw(k,jspin) + qpwc(k)
ENDDO
IF (mpi%irank ==0) THEN
!
!=====> calculate core-tails to the vacuum region
! Coretails expanded in exponentially decaying functions.
! Describe vacuum by:
! rho(r_||,z,ivac)= sum_n{ rho(n,ivac) * exp(-kappa*(z-z_v))
! * exp(iG(n)r_||) }
IF (input%film .AND. .NOT.oneD%odi%d1) THEN
!+gu
dtildh = 0.5 * tpi_const / cell%bmat(3,3)
IF (vacuum%nvac.EQ.1) THEN
rho_out(1) = qpwc(1)*cell%z1
DO k = 2,stars%ng3
IF ((stars%kv3(1,k).EQ.0).AND.(stars%kv3(2,k).EQ.0)) THEN
nz = 1
IF (sym%invs.OR.sym%zrfs) nz = 2
g = stars%kv3(3,k) * cell%bmat(3,3)
rho_out(1) = rho_out(1) + nz*qpwc(k)*SIN(g*cell%z1)/g
ENDIF
ENDDO
rho_out(1) = qpwc(1) * dtildh - rho_out(1)
ELSE
DO ivac = 1, vacuum%nvac
carg = CMPLX(0.0,0.0)
DO k = 2,stars%ng3
IF ((stars%kv3(1,k).EQ.0).AND.(stars%kv3(2,k).EQ.0)) THEN
g = stars%kv3(3,k) * cell%bmat(3,3) * (3. - 2.*ivac)
carg = carg -qpwc(k)*(EXP(ci*g*dtildh)-EXP(ci*g*cell%z1))/g
ENDIF
ENDDO
rho_out(ivac) = qpwc(1) * ( dtildh-cell%z1 ) - AIMAG(carg)
ENDDO
ENDIF
!-gu
! nzvac = min(50,nmz)
m0 = -stars%mx3
IF (sym%zrfs) m0 = 0
!
!---> loop over 2D stars
!
DO 280 k = 1,stars%ng2
k1 = stars%kv2(1,k)
k2 = stars%kv2(2,k)
DO 270 ivac = 1,vacuum%nvac
VALUE = czero
slope = czero
sign = 3. - 2.*ivac
!
! ---> sum over gz-stars
DO 250 kz = m0,stars%mx3
ig3 = stars%ig(k1,k2,kz)
! ----> use only stars within the g_max sphere (oct.97 shz)
IF (ig3.NE.0) THEN
nz = 1
IF (sym%zrfs) nz = stars%nstr(ig3)/stars%nstr2(k)
gz = kz*cell%bmat(3,3)
DO 240 nrz = 1,nz
signz = 3. - 2.*nrz
carg = ci*sign*signz*gz
VALUE = VALUE + qpwc(ig3)* EXP(carg*cell%z1)
slope = slope + carg*qpwc(ig3)* EXP(carg*cell%z1)
240 ENDDO
END IF
250 ENDDO
! roa work-around
IF ( ABS(REAL(VALUE)).GT.zero ) THEN
! roa work-around
! gb works also around
rkappa = - REAL( slope/VALUE )
IF (k.EQ.1) rkappa = VALUE/rho_out(ivac)
! rkappa = - sign * real( slope/value )
IF (rkappa.LE.zero) rkappa=MIN(rkappa,-tol_14)
IF (rkappa.GT.zero) rkappa=MAX(rkappa,tol_14)
! gb works also around
zvac = - LOG( tol_14/cabs(VALUE) ) / rkappa
nzvac = INT( zvac/vacuum%delz ) + 1
! IF ( rkappa.GT.zero .AND. real(value).GT.zero ) THEN
IF ( rkappa.GT.zero ) THEN
z = 0.
IF ( k.EQ.1 ) THEN
DO imz = 1 , MIN( nzvac,vacuum%nmz )
rht(imz,ivac,jspin) = &
& rht(imz,ivac,jspin) + VALUE*EXP(-rkappa*z)
z = z + vacuum%delz
220 ENDDO
ELSE
DO imz = 1 , MIN( nzvac,vacuum%nmzxy )
rhtxy(imz,k-1,ivac,jspin) = &
& rhtxy(imz,k-1,ivac,jspin) + VALUE*EXP(-rkappa*z)
z = z + vacuum%delz
230 ENDDO
END IF
END IF
! roa work-around
END IF
! roa work-around
270 ENDDO
280 ENDDO
ELSEIF (oneD%odi%d1) THEN
!-odim
!---> rho_out is the charge lost between the interstitial and the
!---> rectangular unit cell
rho_out(1) = cell%vol*qpwc(1)
DO k = 2,stars%ng3
IF (stars%kv3(3,k).EQ.0) THEN
irec3 = stars%ig(stars%kv3(1,k),stars%kv3(2,k),stars%kv3(3,k))
IF (irec3.NE.0) THEN
irec2 = stars%ig2(irec3)
IF (irec2.NE.1) THEN
g2 = stars%sk2(irec2)
CALL cylbes(oneD%odi%M,g2*cell%z1,fJ)
rho_out(1) = rho_out(1) +&
& qpwc(k)*2.*cell%vol*fJ(1)/(cell%z1*g2)
END IF
END IF
END IF
END DO
rho_out(1) = cell%bmat(3,3)*(qpwc(1)*cell%omtil - rho_out(1))/(tpi_const*tpi_const)
! then we are constructing our radial components of the vacuum
! charge density so that the so that they are continuous in the
! value and slope on the boundary
! value = sum{gpar}[qpw(gz,gpar)exp{-i*m*phi(gpar))J_m(gpar*z1)]
! slope = .......the same...................*gpar*dJ_m(gpar*z1)]
! for determining the rht we need only continuity, because the
! second condition is the charge neutrality of the system
! Y.Mokrousov
DO irec1 = 1,oneD%odi%nq2
VALUE = czero
slope = czero
gzi = oneD%odi%kv(1,irec1)
m = oneD%odi%kv(2,irec1)
DO irec2 = 1,stars%ng2
irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),gzi)
IF (irec3.NE.0) THEN
g2 = stars%sk2(irec2)
phi = stars%phi2(irec2)
CALL cylbes(oneD%odi%M,g2*cell%z1,fJ)
CALL dcylbs(oneD%odi%M,g2*cell%z1,fJ,dfJ)
VALUE = VALUE + (ci**m)*qpwc(irec3)*&
& EXP(CMPLX(0.,-m*phi))*fJ(m)
slope = slope + (ci**m)*g2*qpwc(irec3)*&
& EXP(CMPLX(0.,-m*phi))*dfJ(m)