Commit 9fe1447a authored by Daniel Wortmann's avatar Daniel Wortmann

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

parents 8c9bdbc4 9cdd9786
......@@ -116,7 +116,7 @@
REAL, INTENT (INOUT) :: rh(DIMENSION%msh,atoms%ntype)
! ..
! .. Local Scalars ..
COMPLEX czero,carg,VALUE,slope,ci
COMPLEX czero,carg,VALUE,slope,ci,c_ph
REAL dif,dxx,g,gz,dtildh,&
& rkappa,sign,signz,tol_14,z,zero,zvac,&
& g2,phi,gamma,qq
......@@ -300,6 +300,7 @@
! ---> sum over gz-stars
DO 250 kz = m0,stars%mx3
ig3 = stars%ig(k1,k2,kz)
c_ph = stars%rgphs(k1,k2,kz) ! phase factor for invs=T & zrfs=F
! ----> use only stars within the g_max sphere (oct.97 shz)
IF (ig3.NE.0) THEN
nz = 1
......@@ -308,8 +309,8 @@
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)
VALUE = VALUE + c_ph*qpwc(ig3)* EXP(carg*cell%z1)
slope = slope + c_ph*carg*qpwc(ig3)* EXP(carg*cell%z1)
240 ENDDO
END IF
250 ENDDO
......
......@@ -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,orbcomp,coreSpecInput,mcd,slab)
moments,coreSpecInput,mcd,slab,orbcomp)
!************************************************************************************
! This is the FLEUR valence density generator
......@@ -75,10 +75,10 @@ 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_orbcomp), INTENT(INOUT) :: orbcomp
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
! Scalar Arguments
INTEGER, INTENT(IN) :: eig_id, jspin
......@@ -249,7 +249,7 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
INQUIRE (file='orbcomprot',exist=l_orbcomprot)
IF (l_orbcomprot) CALL abcrot2(atoms,noccbd,eigVecCoeffs,ispin) ! rotate ab-coeffs
CALL orb_comp(ispin,ikpt,noccbd,atoms,noccbd,usdus,eigVecCoeffs,orbcomp)
IF (PRESENT(orbcomp)) CALL orb_comp(ispin,ikpt,noccbd,atoms,noccbd,usdus,eigVecCoeffs,orbcomp)
END IF
CALL calcDenCoeffs(atoms,sphhar,sym,we,noccbd,eigVecCoeffs,ispin,denCoeffs)
......@@ -270,8 +270,8 @@ 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,orbcomp,&
results,denCoeffs,orb,denCoeffsOffdiag,den,den%mmpMat(:,:,:,jspin),mcd,slab)
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)
END DO
#endif
......
......@@ -60,16 +60,13 @@ CONTAINS
!
IF ((ikpt.LE.mpi%isize).AND..NOT.l_evp) THEN
IF (l_mcd) THEN
mcd%mcd(:,:,:,ikpt,jsp) = 0.0
ENDIF
regCharges%ener(:,:,jsp) = 0.0
regCharges%sqal(:,:,jsp) = 0.0
regCharges%enerlo(:,:,jsp) = 0.0
regCharges%sqlo(:,:,jsp) = 0.0
dos%qal(:,:,:,ikpt,jsp) = 0.0
END IF
!
!---> l-decomposed density for each occupied state
!
! DO 140 i = (skip_t+1),ne ! this I need for all states
......
......@@ -18,6 +18,7 @@ include("cmake/tests/test_LibXC.cmake")
if (FLEUR_USE_MPI)
include("cmake/tests/test_SCALAPACK.cmake")
include("cmake/tests/test_ELPA.cmake")
include("cmake/tests/test_ChASE.cmake")
endif()
include("cmake/compileenv.txt")
......@@ -29,9 +29,11 @@ message("${Green} MPI Library found : ${CReset} ${FLEUR_USE_MPI}")
if (FLEUR_USE_MPI)
message("${Green} SCALAPACK Library found : ${CReset} ${FLEUR_USE_SCALAPACK}")
message("${Green} ELPA Library found : ${CReset} ${FLEUR_USE_ELPA}")
message("${Green} ChASE Library found : ${CReset} ${FLEUR_USE_CHASE}")
else()
message("${Green} SCALAPACK Library found : ${CReset} ---")
message("${Green} ELPA Library found : ${CReset} ---")
message("${Green} ChASE Library found : ${CReset} ---")
endif()
message("${Green} Compile GPU version : ${CReset} ${FLEUR_USE_GPU}")
message("\n")
......
#First check if we can compile with ChASE
try_compile(FLEUR_USE_CHASE ${CMAKE_BINARY_DIR} ${CMAKE_SOURCE_DIR}/cmake/tests/test_ChASE.f90 LINK_LIBRARIES ${FLEUR_LIBRARIES})
#if (NOT FLEUR_USE_CHASE)
# find_package(chase)
# find_package(chase-fleur)
#
# set(TEST_LIBRARIES ${FLEUR_LIBRARIES} ${CHASE_LIBRARIES} ${CHASE-FLEUR_LIBRARIES})
# try_compile(FLEUR_USE_CHASE ${CMAKE_BINARY_DIR} ${CMAKE_SOURCE_DIR}/cmake/tests/test_ChASE.f90 LINK_LIBRARIES ${TEST_LIBRARIES})
#
# if (FLEUR_USE_CHASE)
# set(FLEUR_MPI_LIBRARIES ${FLEUR_MPI_LIBRARIES} ${CHASE_LIBRARIES} ${CHASE-FLEUR_LIBRARIES})
# endif()
#endif()
message("ChASE Library found:${FLEUR_USE_CHASE}")
if (FLEUR_USE_CHASE)
# set(FLEUR_DEFINITIONS ${FLEUR_DEFINITIONS} "CPP_CHASE")
set(FLEUR_MPI_DEFINITIONS ${FLEUR_MPI_DEFINITIONS} "CPP_CHASE")
endif()
!--------------------------------------------------------------------------------
! Copyright (c) 2017 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.
!--------------------------------------------------------------------------------
! Author: Miriam Hinzen
PROGRAM test
use iso_c_binding
interface
subroutine chase_r( h, n, v, ritzv, nev, nex, deg, tol, mode, opt ) bind( c, name = 'dchase_' )
use iso_c_binding
real(c_double) :: h(n,*), v(n,*)
integer(c_int) :: n, deg, nev, nex
real(c_double) :: ritzv(*), tol
character(len=1,kind=c_char) :: mode, opt
end subroutine chase_r
end interface
real(c_double), ALLOCATABLE :: h(:,:), v(:,:), ritzv(:)
integer(c_int) :: n, deg, nev, nex
real(c_double) :: tol
character(len=1,kind=c_char) :: mode, opt
CALL chase_r( h, n, v, ritzv, nev, nex, deg, tol, mode, opt )
END
......@@ -6,4 +6,5 @@ diagonalization/lapack_diag.F90
diagonalization/magma.F90
diagonalization/elpa.F90
diagonalization/scalapack.F90
diagonalization/chase_diag.F90
diagonalization/elemental.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.
!
! @authors: Miriam Hinzen, Gregor Michalicek
!--------------------------------------------------------------------------------
MODULE m_chase_diag
#ifdef CPP_CHASE
IMPLICIT NONE
interface
subroutine chase_c( h, n, v, ritzv, nev, nex, deg, tol, mode, opt ) bind( c, name = 'zchase_' )
use, intrinsic :: iso_c_binding
complex(c_double_complex) :: h(n,*), v(n,*)
integer(c_int) :: n, deg, nev, nex
real(c_double) :: ritzv(*), tol
character(len=1,kind=c_char) :: mode, opt
end subroutine chase_c
end interface
interface
subroutine chase_r( h, n, v, ritzv, nev, nex, deg, tol, mode, opt ) bind( c, name = 'dchase_' )
use, intrinsic :: iso_c_binding
real(c_double_complex) :: h(n,*), v(n,*)
integer(c_int) :: n, deg, nev, nex
real(c_double) :: ritzv(*), tol
character(len=1,kind=c_char) :: mode, opt
end subroutine chase_r
end interface
CONTAINS
SUBROUTINE chase_diag(hmat,smat,ne,eig,zmat)
USE m_types
USE m_judft
USE iso_c_binding
!Simple driver to solve Generalized Eigenvalue Problem using the ChASE library
IMPLICIT NONE
TYPE(t_mat), INTENT(INOUT) :: hmat,smat
INTEGER, INTENT(INOUT) :: ne
CLASS(t_mat), ALLOCATABLE, INTENT(OUT) :: zmat
REAL, INTENT(OUT) :: eig(:)
INTEGER :: i, j, nev, nex
INTEGER :: info
REAL(c_double), ALLOCATABLE :: eigenvalues(:)
REAL(c_double), ALLOCATABLE :: eigenvectors_r(:,:)
COMPLEX(c_double_complex), ALLOCATABLE :: eigenvectors_c(:,:)
ALLOCATE(t_mat::zmat)
CALL zmat%alloc(hmat%l_real,hmat%matsize1,ne)
IF (hmat%l_real) THEN
! --> start with Cholesky factorization of b ( so that b = l * l^t)
! --> b is overwritten by l
CALL dpotrf('U',smat%matsize1,smat%data_r,SIZE(smat%data_r,1),info)
IF (info.NE.0) THEN
WRITE (*,*) 'Error in dpotrf: info =',info
CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
ENDIF
! --> now reduce a * z = eig * b * z to the standard form a' * z' = eig * z'
! --> where a' = (l)^-1 * a * (l^t)^-1 and z' = l^t * z
CALL dsygst(1,'U',smat%matsize1,hmat%data_r,SIZE(hmat%data_r,1),smat%data_r,SIZE(smat%data_r,1),info)
IF (info.NE.0) THEN
WRITE (6,*) 'Error in dsygst: info =',info
CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
ENDIF
! --> solve a' * z' = eig * z' for eigenvalues eig between lb und ub
nev = min(ne,hmat%matsize1)
nex = min(max(nev/4, 45), hmat%matsize1-nev) !dimensioning for workspace
ALLOCATE(eigenvectors_r(smat%matsize1,nev+nex))
ALLOCATE(eigenvalues(nev+nex))
eigenvectors_r = 0.0
eigenvalues = 0.0
do j = 1, hmat%matsize1
do i = 1, j
hmat%data_r(j,i) = hmat%data_r(i,j)
end do
end do
! if(first_entry_franza) then
call chase_r(hmat%data_r, hmat%matsize1, eigenvectors_r, eigenvalues, nev, nex, 25, 1e-10, 'R', 'S' )
! else
! call chase_r(hmat%data_r, hmat%matsize1, eigenvectors_r, eigenvalues, nev, nex, 25, 1e-10, 'A', 'S' )
! end if
ne = nev
! --> recover the generalized eigenvectors z by solving z' = l^t * z
CALL dtrtrs('U','N','N',hmat%matsize1,nev,smat%data_r,eigenvectors_r,zmat%matsize1,info)
IF (info.NE.0) THEN
WRITE (6,*) 'Error in dtrtrs: info =',info
CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
ENDIF
DO i = 1, ne
DO j = 1, hmat%matsize1
zmat%data_r(j,i) = eigenvectors_r(j,i)
END DO
eig(i) = eigenvalues(i)
END DO
!TODO: Store eigenvectors array to reuse it in next iteration
DEALLOCATE(eigenvalues)
DEALLOCATE(eigenvectors_r)
ELSE
! --> start with Cholesky factorization of b ( so that b = l * l^t)
! --> b is overwritten by l
CALL zpotrf('U',smat%matsize1,smat%data_c,SIZE(smat%data_c,1),info)
IF (info.NE.0) THEN
WRITE (*,*) 'Error in zpotrf: info =',info
CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
ENDIF
! --> now reduce a * z = eig * b * z to the standard form a' * z' = eig * z'
! --> where a' = (l)^-1 * a * (l^t)^-1 and z' = l^t * z
CALL zhegst(1,'U',smat%matsize1,hmat%data_c,SIZE(hmat%data_c,1),smat%data_c,SIZE(smat%data_c,1),info)
IF (info.NE.0) THEN
WRITE (6,*) 'Error in zhegst: info =',info
CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
ENDIF
! --> solve a' * z' = eig * z' for eigenvalues eig between lb und ub
nev = min(ne,hmat%matsize1)
nex = min(max(nev/4, 45), hmat%matsize1-nev) !dimensioning for workspace
ALLOCATE(eigenvectors_c(smat%matsize1,nev+nex))
ALLOCATE(eigenvalues(nev+nex))
eigenvectors_c = CMPLX(0.0,0.0)
eigenvalues = 0.0
do j = 1, hmat%matsize1
do i = 1, j
hmat%data_c(j,i) = conjg(hmat%data_c(i,j))
end do
end do
! if(first_entry_franza) then
call chase_c(hmat%data_c, hmat%matsize1, eigenvectors_c, eigenvalues, nev, nex, 25, 1e-10, 'R', 'S' )
! else
! call chase_c(hmat%data_c, hmat%matsize1, eigenvectors_c, eigenvalues, nev, nex, 25, 1e-10, 'A', 'S' )
! end if
ne = nev
! --> recover the generalized eigenvectors z by solving z' = l^t * z
CALL ztrtrs('U','N','N',hmat%matsize1,nev,smat%data_c,eigenvectors_c,zmat%matsize1,info)
IF (info.NE.0) THEN
WRITE (6,*) 'Error in ztrtrs: info =',info
CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
ENDIF
DO i = 1, ne
DO j = 1, hmat%matsize1
zmat%data_c(j,i) = eigenvectors_c(j,i)
END DO
eig(i) = eigenvalues(i)
END DO
!TODO: Store eigenvectors array to reuse it in next iteration
DEALLOCATE(eigenvalues)
DEALLOCATE(eigenvectors_c)
ENDIF
IF (info.NE.0) CALL judft_error("Diagonalization via ChASE failed", calledby = 'chase_diag')
END SUBROUTINE chase_diag
#endif
END MODULE m_chase_diag
......@@ -30,6 +30,7 @@ MODULE m_eigen_diag
INTEGER,PARAMETER:: diag_magma=-6
#endif
INTEGER,PARAMETER:: diag_lapack=4
INTEGER,PARAMETER:: diag_chase=7
INTEGER,PARAMETER:: diag_debugout=99
PUBLIC eigen_diag,parallel_solver_available
CONTAINS
......@@ -44,9 +45,10 @@ CONTAINS
USE m_elpa
USE m_scalapack
USE m_elemental
use m_types_mpimat
USE m_chase_diag
USE m_types_mpimat
IMPLICIT NONE
#ifdef CPP_MPI
#ifdef CPP_MPI
include 'mpif.h'
#endif
CLASS(t_mat),INTENT(INOUT) :: smat,hmat
......@@ -76,6 +78,12 @@ CONTAINS
!CALL magma_diag(hmat,smat,ne,eig,ev)
CASE (diag_lapack)
CALL lapack_diag(hmat,smat,ne,eig,ev)
CASE (diag_chase)
#ifdef CPP_CHASE
CALL chase_diag(hmat,smat,ne,eig,ev)
#else
CALL juDFT_error('ChASE eigensolver selected but not available', calledby = 'eigen_diag')
#endif
CASE (diag_debugout)
CALL priv_debug_out(smat,hmat)
END SELECT
......@@ -147,6 +155,7 @@ CONTAINS
IF (juDFT_was_argument("-diag:elemental")) diag_solver=diag_elemental
IF (juDFT_was_argument("-diag:lapack")) diag_solver=diag_lapack
IF (juDFT_was_argument("-diag:magma")) diag_solver=diag_magma
IF (juDFT_was_argument("-diag:chase")) diag_solver=diag_chase
IF (juDFT_was_argument("-diag:debugout")) diag_solver=diag_debugout
!Check if solver is possible
......
......@@ -37,7 +37,7 @@ CONTAINS
USE m_eig66_io, ONLY : open_eig, write_eig, close_eig,read_eig
USE m_xmlOutput
#ifdef CPP_MPI
USE m_mpi_bc_pot
USE m_mpi_bc_potden
#endif
IMPLICIT NONE
......@@ -110,8 +110,7 @@ CONTAINS
IF (mpi%n_size > 1) l_zref = .FALSE.
#ifdef CPP_MPI
CALL mpi_bc_pot(mpi,stars,sphhar,atoms,input,vacuum,&
v%iter,v%mt,v%pw,v%vacz,v%vacxy)
CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,v)
#endif
!IF (mpi%irank.EQ.0) CALL openXMLElementFormPoly('iteration',(/'numberForCurrentRun','overallNumber '/),(/it,v%iter/),&
! RESHAPE((/19,13,5,5/),(/2,2/)))
......
......@@ -68,7 +68,7 @@ CONTAINS
CALL timestart("Interstitial part")
!Generate interstitial part of Hamiltonian
CALL hs_int(input,noco,stars,lapw,mpi,cell,isp,v%pw,smat,hmat)
CALL hs_int(input,noco,stars,lapw,mpi,cell,isp,v%pw_w,smat,hmat)
CALL timestop("Interstitial part")
CALL timestart("MT part")
!MT-part of Hamiltonian. In case of noco, we need an loop over the local spin of the atoms
......
......@@ -8,10 +8,10 @@ MODULE m_mt_setup
CONTAINS
SUBROUTINE mt_setup(atoms,sym,sphhar,input,noco,enpara,inden,vTot,mpi,results,DIMENSION,td,ud)
USE m_types
USE m_usetup
USE m_tlmplm_cholesky
USE m_tlmplm_store
USE m_types
USE m_spnorb
IMPLICIT NONE
TYPE(t_results),INTENT(INOUT):: results
......
......@@ -21,7 +21,7 @@ MODULE m_eigenso
!
CONTAINS
SUBROUTINE eigenso(eig_id,mpi,DIMENSION,stars,vacuum,atoms,sphhar,&
obsolete,sym,cell,noco,input,kpts,oneD,vTot,enpara)
obsolete,sym,cell,noco,input,kpts,oneD,vTot,enpara,results)
USE m_eig66_io, ONLY : read_eig,write_eig
USE m_spnorb
......@@ -33,21 +33,22 @@ CONTAINS
#endif
IMPLICIT NONE
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_dimension),INTENT(IN) :: DIMENSION
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_obsolete),INTENT(IN) :: obsolete
TYPE(t_input),INTENT(IN) :: input
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_kpts),INTENT(IN) :: kpts
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_potden),INTENT(IN) :: vTot
TYPE(t_enpara),INTENT(IN) :: enpara
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_dimension),INTENT(IN) :: DIMENSION
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_obsolete),INTENT(IN) :: obsolete
TYPE(t_input),INTENT(IN) :: input
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_kpts),INTENT(IN) :: kpts
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_potden),INTENT(IN) :: vTot
TYPE(t_enpara),INTENT(IN) :: enpara
TYPE(t_results),INTENT(INOUT) :: results
! ..
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: eig_id
......@@ -155,6 +156,19 @@ CONTAINS
ENDIF ! (input%eonly) ELSE
deallocate(zso)
ENDDO ! DO nk
! Sorry for the following strange workaround to fill the results%neig and results%eig arrays.
! At some point someone should have a closer look at how the eigenvalues are
! distributed and fill the arrays without using the eigenvalue-IO.
DO jspin = 1, wannierspin
DO nk = 1,kpts%nkpt
CALL read_eig(eig_id,nk,jspin,results%neig(nk,jspin),results%eig(:,nk,jspin))
#ifdef CPP_MPI
CALL MPI_BARRIER(mpi%MPI_COMM,ierr)
#endif
END DO
END DO
RETURN
END SUBROUTINE eigenso
END MODULE m_eigenso
......@@ -64,10 +64,13 @@ CONTAINS
! some praparations to match array sizes
!
nv1(1) = lapw%nv(1) ; nv1(DIMENSION%jspd) = lapw%nv(1)
ALLOCATE ( g1(DIMENSION%nvd,DIMENSION%jspd),g2(DIMENSION%nvd,DIMENSION%jspd),g3(DIMENSION%nvd,DIMENSION%jspd) )
g1(:,1) = lapw%k1(:,1) ; g1(:,DIMENSION%jspd) = lapw%k1(:,1)
g2(:,1) = lapw%k2(:,1) ; g2(:,DIMENSION%jspd) = lapw%k2(:,1)
g3(:,1) = lapw%k3(:,1) ; g3(:,DIMENSION%jspd) = lapw%k3(:,1)
ALLOCATE (g1(DIMENSION%nvd,DIMENSION%jspd))
ALLOCATE (g2(DIMENSION%nvd,DIMENSION%jspd))
ALLOCATE (g3(DIMENSION%nvd,DIMENSION%jspd))
g1 = 0 ; g2 = 0 ; g3 = 0
g1(:SIZE(lapw%k1,1),1) = lapw%k1(:SIZE(lapw%k1,1),1) ; g1(:SIZE(lapw%k1,1),DIMENSION%jspd) = lapw%k1(:SIZE(lapw%k1,1),1)
g2(:SIZE(lapw%k1,1),1) = lapw%k2(:SIZE(lapw%k1,1),1) ; g2(:SIZE(lapw%k1,1),DIMENSION%jspd) = lapw%k2(:SIZE(lapw%k1,1),1)
g3(:SIZE(lapw%k1,1),1) = lapw%k3(:SIZE(lapw%k1,1),1) ; g3(:SIZE(lapw%k1,1),DIMENSION%jspd) = lapw%k3(:SIZE(lapw%k1,1),1)
chelp(:,:,:,:,input%jspins) = CMPLX(0.0,0.0)
......
......@@ -34,10 +34,10 @@ CONTAINS
! istepnow . steps to be done in this run
!
! *********************************************************************
USE m_types
USE m_rwinp
USE m_bfgs
USE m_bfgs0
USE m_types
USE m_constants
USE m_rinpXML
USE m_winpXML
......
......@@ -60,7 +60,7 @@ SUBROUTINE checkDOPAll(input,dimension,sphhar,stars,atoms,sym,vacuum,oneD,&
!---> m.t. boundaries
nat = 1
DO n = 1, atoms%ntype
CALL sphpts(xp,dimension%nspd,atoms%rmt(n),atoms%pos(1,atoms%nat))
CALL sphpts(xp,dimension%nspd,atoms%rmt(n),atoms%pos(1,nat))
CALL checkdop(xp,dimension%nspd,n,nat,0,-1,ispin,&
dimension,atoms,sphhar,stars,sym,vacuum,cell,oneD,potden)
nat = nat + atoms%neq(n)
......
......@@ -23,10 +23,9 @@ c
INTEGER neq,n,na,nlo,line,i,n2spg,nqpt
CHARACTER*4 namgrp,namex
CHARACTER*3 latnam
LOGICAL invs,zrfs,l_J
LOGICAL invs,zrfs
l_kpts = .false.
l_J=.false.
l_qpts=.true.
OPEN (5,file='inp',form='formatted',status='old')
......@@ -44,9 +43,9 @@ c
line = line + 1
READ (5,*)
line = line + 1
READ (5,8010) latnam,namgrp,invs,zrfs,l_J
READ (5,8010) latnam,namgrp,invs,zrfs
line = line + 1
8010 FORMAT (a3,1x,a4,6x,l1,6x,l1,31x,l1)
8010 FORMAT (a3,1x,a4,6x,l1,6x,l1)
READ (5,*)
line = line + 1
......@@ -122,15 +121,10 @@ c
line = line + 1
READ (5,'(16x,i2)',END=99,ERR=99) layerd
line = line + 1
c
c In the case of J-constants calculation check if
c the qpts-file is present:
c
IF (l_J) INQUIRE (file='qpts',exist=l_qpts)
c
c if (.not.l_kpts) then the kpts-file is missing and we have
c to find out how many k-points to generate...
c
INQUIRE (file='QGpsi',exist=l_kpts)
IF (.not.l_kpts) INQUIRE (file='kpts',exist=l_kpts)
IF (.not.l_kpts) THEN
......@@ -160,33 +154,9 @@ c
STOP
ENDIF
96 CONTINUE
c
c If (.not.l_qpts) reading in the number of q-points to be generated
c
IF (.not.l_qpts) THEN
IF(l_kpts)THEN
DO line = 1,7
READ (5,*)
ENDDO
ENDIF
WRITE (6,*) 'No qpts-file exists, trying to generate it'
READ (5,'(5x,i5,3(4x,i2))',END=102,ERR=102)
+ nqpt,nmopq(1),nmopq(2),nmopq(3)
GOTO 101
102 BACKSPACE (5)
READ (5,'(5x,i5,3(4x,i2))',END=103,ERR=103) nqpt
nmopq(1) = 0 ; nmopq(2) = 0 ; nmopq(3) = 0
GOTO 101
103 WRITE (6,*) 'Please give a q-mesh information at the end of'
WRITE (6,*) 'the inp-file'
CLOSE (6)
STOP
ENDIF
101 CONTINUE
!
! determine the number of symmetry operations
!
IF (namgrp.EQ.'any ') THEN
nop = 48
CLOSE (5)
......
......@@ -37,7 +37,7 @@ CONTAINS
REAL gmi,gla,eps
REAL gfx,gfy,pon,pon2
INTEGER j,k,k1,k2,k3,m0,mxx1,mxx2,n
INTEGER ned1,nint,kdone,i
INTEGER ned1,nint,kdone,i,i_sym,n_sym
LOGICAL NEW,l_cdn1,l_xcExtended, l_error
INTEGER kfx,kfy,kfz,kidx,nfftx,nffty,nfftz,kfft
INTEGER nfftxy,norm,n1,kidx2,k2i
......@@ -70,6 +70,12 @@ CONTAINS
GOTO 270
END IF
IF (input%film.AND.sym%invs.AND.(.not.sym%zrfs).AND.(.not.sym%symor)) THEN
n_sym = 2 ! needs reordering of 2d-stars
ELSE
n_sym = 1 ! as before...
ENDIF
mxx1 = 0
mxx2 = 0
stars%ng2 = 0
......@@ -80,45 +86,53 @@ CONTAINS
kv(1) = k1
k2_loop:DO k2 = stars%mx2,-stars%mx2,-1
kv(2) = k2
DO j = 1,2
g(j) = kv(1)*cell%bmat(1,j) + kv(2)*cell%bmat(2,j)
ENDDO
s = SQRT(g(1)**2+g(2)**2)
!---> determine the angle of the G_{||} needed in odim calculations
!+odim YM cute little 'angle' is written by me to determine the phi2
phi = angle(g(1),g(2))
!-odim
!
!---> check if generated vector belongs to a star already
!---> stars should be within the g_max-sphere ! (Oct.97) sbluegel
!
IF (s.LT.stars%gmax) THEN
CALL spgrot(&
& sym%nop,sym%symor,sym%mrot,sym%tau,sym%invtab,&
& kv,&
& kr)
DO n = 1,sym%nop2
IF (mxx1.LT.kr(1,n)) mxx1 = kr(1,n)
IF (mxx2.LT.kr(2,n)) mxx2 = kr(2,n)