Commit e0eefd71 authored by Daniel Wortmann's avatar Daniel Wortmann

Merge branch 'develop' of fleur-git:fleur into develop

parents ff41dc65 95185571
...@@ -7,7 +7,7 @@ ...@@ -7,7 +7,7 @@
MODULE m_chani MODULE m_chani
CONTAINS CONTAINS
SUBROUTINE chani(M,N,Neigd, Myid,Np,Sub_comm,mpi_comm, & SUBROUTINE chani(M,N,Neigd, Myid,Np,Sub_comm,mpi_comm, &
Eig,Num,A_r,B_r,Z_r,A_c,B_c,Z_c) Eig,Num,hamOvlp,zMat)
! !
!---------------------------------------------------- !----------------------------------------------------
!- Parallel eigensystem solver - driver routine; gb99 !- Parallel eigensystem solver - driver routine; gb99
...@@ -33,6 +33,9 @@ CONTAINS ...@@ -33,6 +33,9 @@ CONTAINS
!#include"cpp_arch.h" !#include"cpp_arch.h"
#include"cpp_double.h" #include"cpp_double.h"
USE m_juDFT USE m_juDFT
USE m_types
USE m_subredist1
USE m_subredist2
IMPLICIT NONE IMPLICIT NONE
INCLUDE 'mpif.h' INCLUDE 'mpif.h'
...@@ -40,11 +43,9 @@ CONTAINS ...@@ -40,11 +43,9 @@ CONTAINS
INTEGER, INTENT (IN) :: SUB_COMM,np,myid,mpi_comm INTEGER, INTENT (IN) :: SUB_COMM,np,myid,mpi_comm
INTEGER, INTENT (INOUT) :: num INTEGER, INTENT (INOUT) :: num
REAL, INTENT (OUT) :: eig(neigd) REAL, INTENT (OUT) :: eig(neigd)
REAL,OPTIONAL, INTENT (INOUT) :: a_r(m,n),b_r(m,n) TYPE(t_hamOvlp),INTENT(INOUT) :: hamOvlp
REAL,OPTIONAL, INTENT (OUT) :: z_r(m,neigd) TYPE(t_zMat),INTENT(INOUT) :: zMat
COMPLEX, OPTIONAL,INTENT (INOUT) :: a_c(m,n),b_c(m,n)
COMPLEX, OPTIONAL,INTENT (OUT) :: z_c(m,neigd)
!... Local variables !... Local variables
! !
INTEGER nc,ic,ir,n_sym,jsym,num_j,icm,n_bound INTEGER nc,ic,ir,n_sym,jsym,num_j,icm,n_bound
...@@ -74,13 +75,11 @@ CONTAINS ...@@ -74,13 +75,11 @@ CONTAINS
REAL, ALLOCATABLE :: asca_r(:,:), bsca_r(:,:),work2_r(:) REAL, ALLOCATABLE :: asca_r(:,:), bsca_r(:,:),work2_r(:)
COMPLEX, ALLOCATABLE :: asca_c(:,:), bsca_c(:,:),work2_c(:) COMPLEX, ALLOCATABLE :: asca_c(:,:), bsca_c(:,:),work2_c(:)
EXTERNAL subredist1, subredist2, iceil, numroc EXTERNAL iceil, numroc
EXTERNAL CPP_LAPACK_slamch, descinit EXTERNAL CPP_LAPACK_slamch, descinit
EXTERNAL blacs_pinfo, blacs_gridinit EXTERNAL blacs_pinfo, blacs_gridinit
EXTERNAL MPI_COMM_DUP EXTERNAL MPI_COMM_DUP
LOGICAL :: l_real
l_real=present(a_r)
! !
! determine actual number of columns of input matrices A and B ! determine actual number of columns of input matrices A and B
! nc is number of columns the local processor will get, must be <=n ! nc is number of columns the local processor will get, must be <=n
...@@ -132,17 +131,17 @@ CONTAINS ...@@ -132,17 +131,17 @@ CONTAINS
mycolssca=(m-1)/(nb*npcol)*nb+ & mycolssca=(m-1)/(nb*npcol)*nb+ &
MIN(MAX(m-(m-1)/(nb*npcol)*nb*npcol-nb*mycol,0),nb) MIN(MAX(m-(m-1)/(nb*npcol)*nb*npcol-nb*mycol,0),nb)
! Number of columns the local process gets in ScaLAPACK distribution ! Number of columns the local process gets in ScaLAPACK distribution
if (l_real) THEN if (hamOvlp%l_real) THEN
ALLOCATE ( asca_r(myrowssca,mycolssca), stat=err ) ALLOCATE ( asca_r(myrowssca,mycolssca), stat=err )
else else
ALLOCATE ( asca_c(myrowssca,mycolssca), stat=err ) ALLOCATE ( asca_c(myrowssca,mycolssca), stat=err )
end if endif
IF (err.NE.0) THEN IF (err.NE.0) THEN
WRITE (*,*) 'In chani an error occured during the allocation of' WRITE (*,*) 'In chani an error occured during the allocation of'
WRITE (*,*) 'asca: ',err,' size: ',myrowssca*mycolssca WRITE (*,*) 'asca: ',err,' size: ',myrowssca*mycolssca
CALL juDFT_error("allocte in chani",calledby ="chani") CALL juDFT_error("allocte in chani",calledby ="chani")
ENDIF ENDIF
if (l_real) THEN if (hamOvlp%l_real) THEN
asca_r=0.0 asca_r=0.0
ALLOCATE ( bsca_r(myrowssca,mycolssca), stat=err ) ALLOCATE ( bsca_r(myrowssca,mycolssca), stat=err )
else else
...@@ -154,14 +153,14 @@ CONTAINS ...@@ -154,14 +153,14 @@ CONTAINS
WRITE (*,*) 'bsca :',err WRITE (*,*) 'bsca :',err
CALL juDFT_error("allocate in chani",calledby ="chani") CALL juDFT_error("allocate in chani",calledby ="chani")
ENDIF ENDIF
if (l_real) THEN if (hamOvlp%l_real) THEN
bsca_r=0.0 bsca_r=0.0
CALL subredist1(m,myrowssca,SUB_COMM, nprow, npcol, myid, ierr, nb ,a_r,asca_r) CALL subredist1(m,myrowssca,SUB_COMM, nprow, npcol, myid, ierr, nb ,achi_r=hamovlp%a_r,asca_r=asca_r)
CALL subredist1(m,myrowssca,SUB_COMM, nprow, npcol, myid, ierr, nb ,b_r,bsca_r) CALL subredist1(m,myrowssca,SUB_COMM, nprow, npcol, myid, ierr, nb ,achi_r=hamovlp%b_r,asca_r=bsca_r)
else else
bsca_c=0.0 bsca_c=0.0
CALL subredist1(m,myrowssca,SUB_COMM, nprow, npcol, myid, ierr, nb ,a_c,asca_c) CALL subredist1(m,myrowssca,SUB_COMM, nprow, npcol, myid, ierr, nb ,achi_c=hamovlp%a_c,asca_c=asca_c)
CALL subredist1(m,myrowssca,SUB_COMM, nprow, npcol, myid, ierr, nb ,b_c,bsca_c) CALL subredist1(m,myrowssca,SUB_COMM, nprow, npcol, myid, ierr, nb ,achi_c=hamovlp%b_c,asca_c=bsca_c)
end if end if
CALL BLACS_PINFO(iamblacs,npblacs) ! iamblacs = local process rank (e.g. myid) CALL BLACS_PINFO(iamblacs,npblacs) ! iamblacs = local process rank (e.g. myid)
! npblacs = number of available processes ! npblacs = number of available processes
...@@ -243,7 +242,7 @@ CONTAINS ...@@ -243,7 +242,7 @@ CONTAINS
CALL juDFT_error('Failed to allocated "eig2"', calledby ='chani') CALL juDFT_error('Failed to allocated "eig2"', calledby ='chani')
ENDIF ENDIF
! write(*,*) 'c :',myrowssca,mycolssca,desceigv ! write(*,*) 'c :',myrowssca,mycolssca,desceigv
if (l_real) THEN if (hamovlp%l_real) THEN
ALLOCATE ( eigvec_r(myrowssca,mycolssca), stat=err ) ! Eigenvectors for ScaLAPACK ALLOCATE ( eigvec_r(myrowssca,mycolssca), stat=err ) ! Eigenvectors for ScaLAPACK
else else
ALLOCATE ( eigvec_c(myrowssca,mycolssca), stat=err ) ! Eigenvectors for ScaLAPACK ALLOCATE ( eigvec_c(myrowssca,mycolssca), stat=err ) ! Eigenvectors for ScaLAPACK
...@@ -257,7 +256,7 @@ CONTAINS ...@@ -257,7 +256,7 @@ CONTAINS
nn=MAX(MAX(m,nb),2) nn=MAX(MAX(m,nb),2)
np0=numroc(nn,nb,0,0,nprow) np0=numroc(nn,nb,0,0,nprow)
mq0=numroc(MAX(MAX(neigd,nb),2),nb,0,0,npcol) mq0=numroc(MAX(MAX(neigd,nb),2),nb,0,0,npcol)
if (l_real) THEN if (hamovlp%l_real) THEN
lwork2=5*m+MAX(5*nn,np0*mq0+2*nb*nb)+ iceil(neigd,nprow*npcol)*nn lwork2=5*m+MAX(5*nn,np0*mq0+2*nb*nb)+ iceil(neigd,nprow*npcol)*nn
ALLOCATE ( work2_r(lwork2+10*m), stat=err ) ! Allocate more in case of clusters ALLOCATE ( work2_r(lwork2+10*m), stat=err ) ! Allocate more in case of clusters
else else
...@@ -293,7 +292,7 @@ CONTAINS ...@@ -293,7 +292,7 @@ CONTAINS
! !
! Compute size of workspace ! Compute size of workspace
! !
if (l_real) THEN if (hamovlp%l_real) THEN
uplo='U' uplo='U'
CALL CPP_LAPACK_pdsygvx(1,'V','I','U',m,asca_r,1,1,desca,bsca_r,1,1, desca,& CALL CPP_LAPACK_pdsygvx(1,'V','I','U',m,asca_r,1,1,desca,bsca_r,1,1, desca,&
0.0,1.0,1,num,abstol,num1,num2,eig2,orfac,eigvec_r,1,1,& 0.0,1.0,1,num,abstol,num1,num2,eig2,orfac,eigvec_r,1,1,&
...@@ -352,7 +351,7 @@ endif ...@@ -352,7 +351,7 @@ endif
! !
! Now solve generalized eigenvalue problem ! Now solve generalized eigenvalue problem
! !
if (l_real) THEN if (hamovlp%l_real) THEN
CALL CPP_LAPACK_pdsygvx(1,'V','I','U',m,asca_r,1,1,desca,bsca_r,1,1, desca,& CALL CPP_LAPACK_pdsygvx(1,'V','I','U',m,asca_r,1,1,desca,bsca_r,1,1, desca,&
1.0,1.0,1,num,abstol,num1,num2,eig2,orfac,eigvec_r,1,1,& 1.0,1.0,1,num,abstol,num1,num2,eig2,orfac,eigvec_r,1,1,&
desceigv,work2_r,lwork2,iwork,liwork,ifail,iclustr,& desceigv,work2_r,lwork2,iwork,liwork,ifail,iclustr,&
...@@ -425,10 +424,10 @@ endif ...@@ -425,10 +424,10 @@ endif
! Redistribute eigvec from ScaLAPACK distribution to each process ! Redistribute eigvec from ScaLAPACK distribution to each process
! having all eigenvectors corresponding to his eigenvalues as above ! having all eigenvectors corresponding to his eigenvalues as above
! !
if (l_real) THEN if (hamovlp%l_real) THEN
CALL subredist2(m,num2,myrowssca,SUB_COMM,nprow,npcol, myid,ierr,nb,z_r,eigvec_r) CALL subredist2(m,num2,myrowssca,SUB_COMM,nprow,npcol, myid,ierr,nb,zmat%z_r,eigvec_r)
ELSE ELSE
CALL subredist2(m,num2,myrowssca,SUB_COMM,nprow,npcol, myid,ierr,nb,z_c,eigvec_c) CALL subredist2(m,num2,myrowssca,SUB_COMM,nprow,npcol, myid,ierr,nb,achi_c=zmat%z_c,asca_c=eigvec_c)
end if end if
! !
!DEALLOCATE ( eigvec) !DEALLOCATE ( eigvec)
......
...@@ -162,7 +162,11 @@ CONTAINS ...@@ -162,7 +162,11 @@ CONTAINS
SELECT CASE (priv_select_solver(parallel)) SELECT CASE (priv_select_solver(parallel))
#ifdef CPP_ELPA #ifdef CPP_ELPA
CASE (diag_elpa) CASE (diag_elpa)
CALL elpa(lapw%nmat,n,SUB_COMM,eig,ne_found,hamOvlp%a_r,hamOvlp%b_r,zMat%z_r,hamOvlp%a_c,hamOvlp%b_c,zMat%z_c) IF (hamovlp%l_real) THEN
CALL elpa(lapw%nmat,n,SUB_COMM,hamOvlp%a_r,hamOvlp%b_r,zMat%z_r,eig,ne_found)
ELSE
CALL elpa(lapw%nmat,n,SUB_COMM,hamOvlp%a_c,hamOvlp%b_c,zMat%z_c,eig,ne_found)
ENDIF
#endif #endif
#ifdef CPP_ELEMENTAL #ifdef CPP_ELEMENTAL
CASE (diag_elemental) CASE (diag_elemental)
...@@ -175,7 +179,7 @@ CONTAINS ...@@ -175,7 +179,7 @@ CONTAINS
#endif #endif
#ifdef CPP_SCALAPACK #ifdef CPP_SCALAPACK
CASE (diag_scalapack) CASE (diag_scalapack)
CALL chani(lapw%nmat,dimension%nbasfcn/n_size,ndim, n_rank,n_size,SUB_COMM,mpi%mpi_comm,eig,ne_found,hamOvlp%a_r,hamOvlp%b_r,zMat%z_r,hamOvlp%a_c,hamOvlp%b_c,zMat%z_c) CALL chani(lapw%nmat,dimension%nbasfcn/n_size,ndim, n_rank,n_size,SUB_COMM,mpi%mpi_comm,eig,ne_found,hamOvlp,zMat)
#endif #endif
#ifdef CPP_MAGMA #ifdef CPP_MAGMA
CASE (diag_magma) CASE (diag_magma)
......
...@@ -26,6 +26,7 @@ CONTAINS ...@@ -26,6 +26,7 @@ CONTAINS
#define CPP_ONE 1.0 #define CPP_ONE 1.0
#define CPP_ZERO 0.0 #define CPP_ZERO 0.0
#define CPP_mult mult_at_b_real #define CPP_mult mult_at_b_real
#define CPP_REAL
SUBROUTINE elpa_r(m,n, SUB_COMM, a,b, z,eig,num) SUBROUTINE elpa_r(m,n, SUB_COMM, a,b, z,eig,num)
! !
!---------------------------------------------------- !----------------------------------------------------
...@@ -61,7 +62,7 @@ CONTAINS ...@@ -61,7 +62,7 @@ CONTAINS
#define CPP_ONE cmplx(1.,0.) #define CPP_ONE cmplx(1.,0.)
#define CPP_ZERO cmplx(0.,0.) #define CPP_ZERO cmplx(0.,0.)
#define CPP_mult mult_ah_b_complex #define CPP_mult mult_ah_b_complex
#undef CPP_REAL
SUBROUTINE elpa_c(m,n, SUB_COMM, a,b, z,eig,num) SUBROUTINE elpa_c(m,n, SUB_COMM, a,b, z,eig,num)
! !
!---------------------------------------------------- !----------------------------------------------------
......
...@@ -5,6 +5,8 @@ ...@@ -5,6 +5,8 @@
!-------------------------------------------------------------------------------- !--------------------------------------------------------------------------------
USE m_juDFT USE m_juDFT
USE elpa1 USE elpa1
USE m_subredist1
USE m_subredist2
#ifdef CPP_ELPA2 #ifdef CPP_ELPA2
USE elpa2 USE elpa2
#endif #endif
...@@ -19,9 +21,6 @@ ...@@ -19,9 +21,6 @@
CPP_REALCOMPLEX, INTENT (OUT) :: z(:,:) CPP_REALCOMPLEX, INTENT (OUT) :: z(:,:)
! !
#include "elpa_cpp.F90"
END SUBROUTINE elpa_r
!... Local variables !... Local variables
! !
INTEGER :: sc_desc(9) !SCALAPACK DESCRIPTOR INTEGER :: sc_desc(9) !SCALAPACK DESCRIPTOR
...@@ -56,9 +55,13 @@ ...@@ -56,9 +55,13 @@
! block-size nb ! block-size nb
! !
!print *,"Before subredist" !print *,"Before subredist"
CALL subredist1(a,m,asca,myrowssca,SUB_COMM, nprow, npcol, myid, ierr, nb ) #ifdef CPP_REAL
CALL subredist1(b,m,bsca,myrowssca,SUB_COMM, nprow, npcol, myid, ierr, nb ) CALL subredist1(m,myrowssca,SUB_COMM, nprow, npcol, myid, ierr, nb, achi_r=a, asca_r=asca )
CALL subredist1(m,myrowssca,SUB_COMM, nprow, npcol, myid, ierr, nb, achi_r=b, asca_r=bsca)
#else
CALL subredist1(m,myrowssca,SUB_COMM, nprow, npcol, myid, ierr, nb, achi_c=a, asca_c=asca )
CALL subredist1(m,myrowssca,SUB_COMM, nprow, npcol, myid, ierr, nb, achi_c=b, asca_c=bsca)
#endif
! for saving memory one might deallocate(a,b) ! for saving memory one might deallocate(a,b)
!print *,"Before Allocate" !print *,"Before Allocate"
ALLOCATE ( eig2(m), stat=err ) ! The eigenvalue array for ScaLAPACK ALLOCATE ( eig2(m), stat=err ) ! The eigenvalue array for ScaLAPACK
...@@ -181,7 +184,11 @@ ...@@ -181,7 +184,11 @@
! Redistribute eigvec from ScaLAPACK distribution to each process ! Redistribute eigvec from ScaLAPACK distribution to each process
! having all eigenvectors corresponding to his eigenvalues as above ! having all eigenvectors corresponding to his eigenvalues as above
! !
CALL subredist2(z,m,num2,asca,myrowssca,SUB_COMM,nprow,npcol, myid,ierr,nb) #ifdef CPP_REAL
CALL subredist2(m,num2,myrowssca,SUB_COMM,nprow,npcol, myid,ierr,nb,z,asca)
#else
CALL subredist2(m,num2,myrowssca,SUB_COMM,nprow,npcol, myid,ierr,nb,achi_c=z,asca_c=asca)
#endif
! !
DEALLOCATE ( asca ) DEALLOCATE ( asca )
......
...@@ -3,7 +3,9 @@ ...@@ -3,7 +3,9 @@
! This file is part of FLEUR and available as free software under the conditions ! 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. ! of the MIT license as expressed in the LICENSE file in more detail.
!-------------------------------------------------------------------------------- !--------------------------------------------------------------------------------
MODULE m_subredist1
CONTAINS
!DEC$ FREEFORM !DEC$ FREEFORM
SUBROUTINE subredist1(n,lda,SUB_COMM,nprow,npcol,iam,ierr,nb,achi_r,asca_r,achi_c,asca_c) SUBROUTINE subredist1(n,lda,SUB_COMM,nprow,npcol,iam,ierr,nb,achi_r,asca_r,achi_c,asca_c)
USE m_juDFT USE m_juDFT
...@@ -29,9 +31,9 @@ SUBROUTINE subredist1(n,lda,SUB_COMM,nprow,npcol,iam,ierr,nb,achi_r,asca_r,achi_ ...@@ -29,9 +31,9 @@ SUBROUTINE subredist1(n,lda,SUB_COMM,nprow,npcol,iam,ierr,nb,achi_r,asca_r,achi_
! !
! Parameters, I/O channel numbers ! Parameters, I/O channel numbers
REAL,OPTIONAL :: achi_r(:) ! Input, matrix in chani-distribution REAL,OPTIONAL :: achi_r(:) ! Input, matrix in chani-distribution
REAL,OPTIONAL :: asca_r(lda,*) ! Output, matrix in ScaLAPACK distribution REAL,OPTIONAL :: asca_r(:,:)!(lda,*) ! Output, matrix in ScaLAPACK distribution
COMPLEX,OPTIONAL :: achi_c(*) ! Input, matrix in chani-distribution COMPLEX,OPTIONAL :: achi_c(:) ! Input, matrix in chani-distribution
COMPLEX,OPTIONAL :: asca_c(lda,*) ! Output, matrix in ScaLAPACK distribution COMPLEX,OPTIONAL :: asca_c(:,:)!(lda,*) ! Output, matrix in ScaLAPACK distribution
! Matrix might be real or complex ! Matrix might be real or complex
INTEGER :: n,lda ! Global matrix size, local leading dimension of asca INTEGER :: n,lda ! Global matrix size, local leading dimension of asca
...@@ -666,3 +668,4 @@ SUBROUTINE subredist1(n,lda,SUB_COMM,nprow,npcol,iam,ierr,nb,achi_r,asca_r,achi_ ...@@ -666,3 +668,4 @@ SUBROUTINE subredist1(n,lda,SUB_COMM,nprow,npcol,iam,ierr,nb,achi_r,asca_r,achi_
DEALLOCATE(icommcol) DEALLOCATE(icommcol)
RETURN RETURN
END SUBROUTINE subredist1 END SUBROUTINE subredist1
END
MODULE m_subredist2
CONTAINS
!-------------------------------------------------------------------------------- !--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany ! 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 ! This file is part of FLEUR and available as free software under the conditions
...@@ -516,3 +518,4 @@ SUBROUTINE subredist2(n,m,lda,SUB_COMM,nprow,npcol,& ...@@ -516,3 +518,4 @@ SUBROUTINE subredist2(n,m,lda,SUB_COMM,nprow,npcol,&
DEALLOCATE(icommcol) DEALLOCATE(icommcol)
RETURN RETURN
END SUBROUTINE subredist2 END SUBROUTINE subredist2
END
...@@ -42,8 +42,8 @@ CONTAINS ...@@ -42,8 +42,8 @@ CONTAINS
TYPE(t_tlmplm),INTENT(INOUT) :: tlmplm TYPE(t_tlmplm),INTENT(INOUT) :: tlmplm
LOGICAL,INTENT(IN) :: l_real LOGICAL,INTENT(IN) :: l_real
REAL, OPTIONAL,INTENT (INOUT) :: aa_r(:)!(matsize) REAL, OPTIONAL,ALLOCATABLE,INTENT (INOUT) :: aa_r(:)!(matsize)
COMPLEX, OPTIONAL,INTENT (INOUT) :: aa_c(:) COMPLEX, OPTIONAL,ALLOCATABLE,INTENT (INOUT) :: aa_c(:)
! .. ! ..
! .. Local Scalars .. ! .. Local Scalars ..
COMPLEX axx,bxx,cxx,dtd,dtu,dtulo,ulotd,ulotu,ulotulo,utd,utu, utulo,chihlp COMPLEX axx,bxx,cxx,dtd,dtu,dtulo,ulotd,ulotu,ulotulo,utd,utu, utulo,chihlp
......
...@@ -60,8 +60,8 @@ ...@@ -60,8 +60,8 @@
INTEGER, PARAMETER :: neig12=12 INTEGER, PARAMETER :: neig12=12
!===> Local Variables !===> Local Variables
LOGICAL lerr,err_setup,invsym,inversionOp LOGICAL lerr,err_setup,invsym
INTEGER i,j,k,n,m,na,nt,mdet,mtr,nop0,fh INTEGER i,j,k,n,m,na,nt,mdet,mtr,nop0,fh,inversionOp
REAL t,volume,eps7,eps12 REAL t,volume,eps7,eps12
INTEGER optype(nop48) INTEGER optype(nop48)
......
...@@ -39,7 +39,7 @@ CONTAINS ...@@ -39,7 +39,7 @@ CONTAINS
TYPE(t_data_MPI),POINTER :: d TYPE(t_data_MPI),POINTER :: d
CALL priv_find_data(id,d) CALL priv_find_data(id,d)
CALL eig66_data_storedefault(d,jspins,nkpts,nmat,neig,lmax,nlotot,nlo,ntype,l_real,l_soc,l_dos,l_mcd,l_orb) CALL eig66_data_storedefault(d,jspins,nkpts,nmat,neig,lmax,nlotot,nlo,ntype,l_real.and..not.l_soc,l_soc,l_dos,l_mcd,l_orb)
IF (PRESENT(n_size_opt)) d%n_size=n_size_opt IF (PRESENT(n_size_opt)) d%n_size=n_size_opt
IF (ALLOCATED(d%pe_ev)) THEN IF (ALLOCATED(d%pe_ev)) THEN
...@@ -48,7 +48,7 @@ CONTAINS ...@@ -48,7 +48,7 @@ CONTAINS
d%eig_data=1E99 d%eig_data=1E99
d%int_data=9999999 d%int_data=9999999
d%real_data=1E99 d%real_data=1E99
if (l_real.and..not.l_soc) THEN if (d%l_real.and..not.l_soc) THEN
d%zr_data=0.0 d%zr_data=0.0
else else
d%zc_data=0.0 d%zc_data=0.0
...@@ -193,10 +193,10 @@ CONTAINS ...@@ -193,10 +193,10 @@ CONTAINS
IF (d%irank==0) THEN IF (d%irank==0) THEN
tmp_id=eig66_data_newid(DA_mode) tmp_id=eig66_data_newid(DA_mode)
IF (d%l_dos) CPP_error("Could not read DOS data") IF (d%l_dos) CPP_error("Could not read DOS data")
CALL open_eig_DA(tmp_id,nmat,neig,nkpts,jspins,lmax,nlo,ntype,nlotot,.FALSE.,.FALSE.,l_real,l_soc,.FALSE.,.FALSE.,filename) CALL open_eig_DA(tmp_id,nmat,neig,nkpts,jspins,lmax,nlo,ntype,nlotot,.FALSE.,.FALSE.,d%l_real,l_soc,.FALSE.,.FALSE.,filename)
DO jspin=1,jspins DO jspin=1,jspins
DO nk=1,nkpts DO nk=1,nkpts
if (l_real) THEN if (d%l_real) THEN
CALL read_eig_DA(tmp_id,nk,jspin,nv,i,k1,k2,k3,bk3,wk,ii,eig,el,ello,evac,kveclo,z=z_r) CALL read_eig_DA(tmp_id,nk,jspin,nv,i,k1,k2,k3,bk3,wk,ii,eig,el,ello,evac,kveclo,z=z_r)
CALL write_eig(id,nk,jspin,ii,ii,nv,nmat,k1,k2,k3,bk3,wk,eig,el,ello,evac,nlotot,kveclo,z=z_r) CALL write_eig(id,nk,jspin,ii,ii,nv,nmat,k1,k2,k3,bk3,wk,eig,el,ello,evac,nlotot,kveclo,z=z_r)
ELSE ELSE
......
...@@ -89,20 +89,20 @@ ...@@ -89,20 +89,20 @@
! transform. At the moment it is designed to use fftw3. ! transform. At the moment it is designed to use fftw3.
SUBROUTINE inplfft( fftin, idir ) SUBROUTINE inplfft( fftin, idir )
USE m_juDFT
USE param, ONLY: jnlout USE param, ONLY: jnlout
USE nonlocal_data,ONLY: nx,ny,nz,n_grid USE nonlocal_data,ONLY: nx,ny,nz,n_grid
implicit none implicit none
include 'fftw3.f' ! some definitions for fftw3
complex, intent(inout) :: fftin(n_grid) complex, intent(inout) :: fftin(n_grid)
integer :: idir integer :: idir
integer*8 :: plan integer*8 :: plan
#ifdef CPP_FFTW
include 'fftw3.f' ! some definitions for fftw3
if ( idir == -1 ) then if ( idir == -1 ) then
...@@ -125,6 +125,9 @@ ...@@ -125,6 +125,9 @@
STOP 'Error in FFT' STOP 'Error in FFT'
endif endif
#else
CALL judft_error("VdW functionals only available if compiled with CPP_FFTW")
#endif
END SUBROUTINE inplfft END SUBROUTINE inplfft
! !
END MODULE driver_fft END MODULE driver_fft
...@@ -203,30 +206,30 @@ ...@@ -203,30 +206,30 @@
! we also need LDA correlation per particle so I repeat here the formulas from calc_q0 ! we also need LDA correlation per particle so I repeat here the formulas from calc_q0
sqrt_r_s = dsqrt(r_s) sqrt_r_s = sqrt(r_s)
LDA_dummy_1 = (8.0*pi/3.0)*A*(1.0+ alpha_1*r_s) LDA_dummy_1 = (8.0*pi/3.0)*A*(1.0+ alpha_1*r_s)
LDA_dummy_2 = 2.0*A*(beta_1*sqrt_r_s + beta_2*r_s + & LDA_dummy_2 = 2.0*A*(beta_1*sqrt_r_s + beta_2*r_s + &
beta_3*sqrt_r_s*r_s + beta_4*r_s*r_s) beta_3*sqrt_r_s*r_s + beta_4*r_s*r_s)
eLDA_c = (-3.0/(4.0*pi))*LDA_dummy_1*dlog(1.0+1.0/LDA_dummy_2) eLDA_c = (-3.0/(4.0*pi))*LDA_dummy_1*log(1.0+1.0/LDA_dummy_2)
! now start calculation of PBE correlation. first check wether density is very small (approx ! now start calculation of PBE correlation. first check wether density is very small (approx
! zero) or negative which is also not correct ! zero) or negative which is also not correct
! !
k_F = (3.0*pi*pi*n(i1))**(1.0/3.0) k_F = (3.0*pi*pi*n(i1))**(1.0/3.0)
k_s = dsqrt( 4.0 * k_F / pi) k_s = sqrt( 4.0 * k_F / pi)
grad_n_squ = grad_n(i1,1)**2 + grad_n(i1,2)**2 + grad_n(i1,3)**2 grad_n_squ = grad_n(i1,1)**2 + grad_n(i1,2)**2 + grad_n(i1,3)**2
t = dsqrt(grad_n_squ)/(2.0*phi_zeta*k_s*n(i1)) t = sqrt(grad_n_squ)/(2.0*phi_zeta*k_s*n(i1))
AA = (beta_H/gamma_H)* & AA = (beta_H/gamma_H)* &
(1.0/(dexp(-1.0*eLDA_c/(gamma_H*phi_zeta**3))- 1.0 )) (1.0/(exp(-1.0*eLDA_c/(gamma_H*phi_zeta**3))- 1.0 ))
H = gamma_H*phi_zeta**3 * & H = gamma_H*phi_zeta**3 * &
dlog(1.0+ & log(1.0+ &
(beta_H/gamma_H)*(t**2)* & (beta_H/gamma_H)*(t**2)* &
(1.0+AA*t**2)/(1.0+AA*t**2+AA**2 * t**4)) (1.0+AA*t**2)/(1.0+AA*t**2+AA**2 * t**4))
! !
...@@ -289,7 +292,7 @@ ...@@ -289,7 +292,7 @@
! !
k_F = (3.0*pi*pi*n(i1))**(1.0/3.0) k_F = (3.0*pi*pi*n(i1))**(1.0/3.0)
! !
sqrt_grad_n = dsqrt(grad_n(i1,1)**2 + grad_n(i1,2)**2 + grad_n(i1,3)**2) sqrt_grad_n = sqrt(grad_n(i1,1)**2 + grad_n(i1,2)**2 + grad_n(i1,3)**2)
! !
ss=sqrt_grad_n/(2.0*k_F*n(i1)) ss=sqrt_grad_n/(2.0*k_F*n(i1))
! !
...@@ -1034,7 +1037,7 @@ ...@@ -1034,7 +1037,7 @@
beta_3*sqrt_r_s*r_s + beta_4*r_s*r_s) beta_3*sqrt_r_s*r_s + beta_4*r_s*r_s)
eLDA_x = (-3.0/(4.0*pi))*k_F eLDA_x = (-3.0/(4.0*pi))*k_F
eLDA_c = (-3.0/(4.0*pi))*LDA_dummy_1*dlog(1.0+1.0/LDA_dummy_2) eLDA_c = (-3.0/(4.0*pi))*LDA_dummy_1*log(1.0+1.0/LDA_dummy_2)
! !
! LDA correlation energy density is epsilon_c[n(r)]*n(r) ! LDA correlation energy density is epsilon_c[n(r)]*n(r)
! !
...@@ -1053,7 +1056,7 @@ ...@@ -1053,7 +1056,7 @@
enddo enddo
q_0(i) = q_cut*(1.0 - dexp(-sum_m)) q_0(i) = q_cut*(1.0 - exp(-sum_m))
enddo enddo
...@@ -1131,7 +1134,7 @@ ...@@ -1131,7 +1134,7 @@
ind = G_ind(i) ind = G_ind(i)
k = dsqrt(dble( G(i,1)**2 + G(i,2)**2 + G(i,3)**2)) k = sqrt(dble( G(i,1)**2 + G(i,2)**2 + G(i,3)**2))
call interpolate_kernel(k, phi_k) call interpolate_kernel(k, phi_k)
...@@ -1472,7 +1475,7 @@ ...@@ -1472,7 +1475,7 @@
beta_3*sqrt_r_s*r_s + beta_4*r_s*r_s) beta_3*sqrt_r_s*r_s + beta_4*r_s*r_s)
eLDA_x = (-3.0/(4.0*pi))*k_F eLDA_x = (-3.0/(4.0*pi))*k_F
eLDA_c = (-3.0/(4.0*pi))*LDA_dummy_1*dlog(1.0+1.0/LDA_dummy_2) eLDA_c = (-3.0/(4.0*pi))*LDA_dummy_1*log(1.0+1.0/LDA_dummy_2)
q = (1.0 + ( eLDA_c / eLDA_x) - (Zab/9.0)*grad_n_squ & q = (1.0 + ( eLDA_c / eLDA_x) - (Zab/9.0)*grad_n_squ &
/(4.0*(n(i)**2.0)*(k_F**2.0)))*k_F /(4.0*(n(i)**2.0)*(k_F**2.0)))*k_F
...@@ -1490,7 +1493,7 @@ ...@@ -1490,7 +1493,7 @@
enddo enddo
dq0_dq = dq0_dq * dexp(-sum_m) dq0_dq = dq0_dq * exp(-sum_m)
! here we calculate the derivatives needed for the potential later. i.e. we calculate ! here we calculate the derivatives needed for the potential later. i.e. we calculate
! n*dq0/dn and n*dq0/dgrad_n so that we do not have to divide by n which might be very ! n*dq0/dn and n*dq0/dgrad_n so that we do not have to divide by n which might be very
...@@ -1499,13 +1502,13 @@ ...@@ -1499,13 +1502,13 @@
dq0_dn(i) = dq0_dq * ( k_F/3.0 & dq0_dn(i) = dq0_dq * ( k_F/3.0 &
+ k_F*7.0/3.0*(Zab/9.0)*grad_n_squ &