Refactoring of vmtxc.F and vmtxcg.F started. Not correct yet..

types/types_xcpot.F90

@@ -39,6 +39,7 @@ MODULE m_types_xcpot
      REAL,ALLOCATABLE :: gggrd(:),grgru(:),grgrd(:)
      !These are the contracted Gradients used in libxc
      REAL,ALLOCATABLE :: sigma(:,:)
+     REAL,ALLOCATABLE :: gr(:,:,:)
   END TYPE t_gradients
 
 CONTAINS

types/types_xcpot_libxc.F90

@@ -172,6 +172,8 @@ CONTAINS
     TYPE(t_gradients),INTENT(OUT):: grad
     !For libxc we only need the sigma array...
     ALLOCATE(grad%sigma(MERGE(1,3,jspins==1),ngrid))
+    ALLOCATE(grad%gr(3,ngrid,jspins))
+    
   END SUBROUTINE xcpot_alloc_gradients
 
 

vgen/CMakeLists.txt

@@ -19,6 +19,7 @@ set(fleur_F90 ${fleur_F90}
 vgen/b_field.F90
 vgen/convol.f90
 vgen/grdrsis.f90
+vgen/mt_tofrom_grid.F90
 vgen/int_nv.F90
 vgen/lhglptg.f90
 vgen/lhglpts.f90
@@ -36,8 +37,7 @@ vgen/visxc.f90
 vgen/visxcg.f90
 vgen/vmatgen.f90
 vgen/vmts.F90
-vgen/vmtxc.f90
-vgen/vmtxcg.F90
+vgen/vmt_xc.F90
 vgen/vvac.f90
 vgen/vvacis.f90
 vgen/vvacxc.f90

vgen/mt_tofrom_grid.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_mt_tofrom_grid
+  USE m_types
+  PRIVATE
+  REAL,PARAMETER    :: d_15 = 1.e-15
+  INTEGER,PARAMETER :: ndvgrd=6 ! this should be consistent across GGA derivative routines
+  INTEGER           :: nsp
+  REAL, ALLOCATABLE :: ylh(:,:,:),ylht(:,:,:),ylhtt(:,:,:)
+  REAL, ALLOCATABLE :: ylhf(:,:,:),ylhff(:,:,:),ylhtf(:,:,:)
+  REAL, ALLOCATABLE :: wt(:),rx(:,:),thet(:)
+  PUBLIC :: init_mt_grid,mt_to_grid,mt_from_grid,finish_mt_grid
+CONTAINS
+  SUBROUTINE init_mt_grid(nsp,jspins,atoms,sphhar,xcpot,sym,l_grad)
+    USE m_gaussp
+    USE m_lhglptg
+    USE m_lhglpts
+    IMPLICIT NONE
+    INTEGER,INTENT(IN)          :: nsp,jspins
+    TYPE(t_atoms),INTENT(IN)    :: atoms
+    TYPE(t_sphhar),INTENT(IN)   :: sphhar
+    CLASS(t_xcpot),INTENT(IN)   :: xcpot
+    TYPE(t_sym),INTENT(IN)      :: sym
+    LOGICAL,INTENT(IN)          :: l_grad
+    
+    ! generate nspd points on a sherical shell with radius 1.0
+    ! angular mesh equidistant in phi,
+    ! theta are zeros of the legendre polynomials
+    ALLOCATE(wt(nsp),rx(3,nsp),thet(nsp))
+    CALL gaussp(atoms%lmaxd, rx,wt)
+    ! generate the lattice harmonics on the angular mesh
+    ALLOCATE ( ylh(nsp,0:sphhar%nlhd,sphhar%ntypsd))
+    IF (l_grad) ALLOCATE(ylht,ylhtt,ylhf,ylhff,ylhtf,MOLD=ylh )
+
+    IF (l_grad) THEN
+       CALL lhglptg(sphhar,atoms,rx,nsp,xcpot,sym,&
+            ylh,thet,ylht,ylhtt,ylhf,ylhff,ylhtf)
+    ELSE
+       CALL lhglpts( sphhar,atoms, rx,nsp, sym, ylh)
+    END IF
+  END SUBROUTINE init_mt_grid
+  
+  SUBROUTINE mt_to_grid(atoms,sphhar,den,jspins,n,l_grad,ch,grad)
+    USE m_grdchlh
+    USE m_mkgylm
+    IMPLICIT NONE
+    TYPE(t_atoms),INTENT(IN)    :: atoms
+    TYPE(t_sphhar),INTENT(IN)   :: sphhar
+    TYPE(t_potden),INTENT(IN)   :: den
+    INTEGER,INTENT(IN)          :: n,jspins
+    LOGICAL,INTENT(IN)          :: l_grad
+    REAL,INTENT(OUT)               :: ch(:,:)
+    TYPE(t_gradients),INTENT(INOUT)::grad
+    
+    REAL, ALLOCATABLE :: chlh(:,:,:),chlhdr(:,:,:),chlhdrr(:,:,:)
+    REAL, ALLOCATABLE :: chdr(:,:),chdt(:,:),chdf(:,:)
+    REAL, ALLOCATABLE :: chdrr(:,:),chdtt(:,:),chdff(:,:),chdtf(:,:)
+    REAL, ALLOCATABLE :: chdrt(:,:),chdrf(:,:)
+    INTEGER:: nd,lh,js,jr,kt,k
+
+    nd = atoms%ntypsy(SUM(atoms%neq(:n-1))+1)
+
+    ALLOCATE ( chlh(atoms%jmtd,0:sphhar%nlhd,jspins))
+
+    IF (l_grad)  THEN
+       ALLOCATE(chdr(nsp,jspins),chdt(nsp,jspins),chdf(nsp,jspins),chdrr(nsp,jspins),&
+            chdtt(nsp,jspins),chdff(nsp,jspins),chdtf(nsp,jspins),chdrt(nsp,jspins),&
+            chdrf(nsp,jspins) )
+       ALLOCATE (chlhdr(atoms%jmtd,0:sphhar%nlhd,jspins))
+       ALLOCATE (chlhdrr(atoms%jmtd,0:sphhar%nlhd,jspins))
+    ENDIF
+    
+    DO lh = 0,sphhar%nlh(nd)
+
+       !         calculates gradients of radial charge densities of l=> 0.
+       !         rho*ylh/r**2 is charge density. chlh=rho/r**2.
+       !         charge density=sum(chlh*ylh).
+       !         chlhdr=d(chlh)/dr, chlhdrr=dd(chlh)/drr.
+       
+       DO js = 1,jspins      
+          DO jr = 1,atoms%jri(n)
+             chlh(jr,lh,js) = den%mt(jr,lh,n,js)/(atoms%rmsh(jr,n)*atoms%rmsh(jr,n))
+          ENDDO
+          IF (l_grad) CALL grdchlh(1,1,atoms%jri(n),atoms%dx(n),atoms%rmsh(1,n),&
+               chlh(1,lh,js),ndvgrd, chlhdr(1,lh,js),chlhdrr(1,lh,js))
+          
+       ENDDO ! js
+    ENDDO   ! lh
+    kt=0
+    DO jr = 1,atoms%jri(n)
+       !         following are at points on jr-th sphere.
+       ch(:,:)    = 0.0     ! charge density (on extended grid for all jr)
+       !  generate the densities on an angular mesh
+       DO js = 1,jspins
+          DO lh = 0,sphhar%nlh(nd)
+             DO k = 1,nsp
+                ch(kt+k,js) = ch(kt+k,js) + ylh(k,lh,nd)*chlh(jr,lh,js)
+             ENDDO
+          ENDDO
+       ENDDO
+       IF (l_grad) THEN
+          chdr(:,:)  = 0.0     ! d(ch)/dr
+          chdt(:,:)  = 0.0     ! d(ch)/dtheta
+          chdf(:,:)  = 0.0     ! d(ch)/dfai
+          chdrr(:,:) = 0.0     ! dd(ch)/drr
+          chdtt(:,:) = 0.0     ! dd(ch)/dtt
+          chdff(:,:) = 0.0     ! dd(ch)/dff
+          chdtf(:,:) = 0.0     ! dd(ch)/dtf
+          chdrt(:,:) = 0.0     ! d(d(ch)/dr)dt
+          chdrf(:,:) = 0.0     ! d(d(ch)/dr)df
+          !  generate the derivatives on an angular mesh
+          DO js = 1,jspins
+             DO lh = 0,sphhar%nlh(nd)
+                ! 
+                DO k = 1,nsp
+                   chdr(k,js) =chdr(k,js)+ ylh(k,lh,nd)*chlhdr(jr,lh,js)
+                   chdrr(k,js)=chdrr(k,js)+ylh(k,lh,nd)*chlhdrr(jr,lh,js)
+                ENDDO
+                
+                DO k = 1,nsp
+                   chdrt(k,js)=chdrt(k,js)+ylht(k,lh,nd)*chlhdr(jr,lh,js)
+                   chdrf(k,js)=chdrf(k,js)+ylhf(k,lh,nd)*chlhdr(jr,lh,js)
+                   chdt(k,js) =chdt(k,js) +ylht(k,lh,nd)*chlh(jr,lh,js)
+                   chdf(k,js) =chdf(k,js) +ylhf(k,lh,nd)*chlh(jr,lh,js)
+                   chdtt(k,js)=chdtt(k,js)+ylhtt(k,lh,nd)*chlh(jr,lh,js)
+                   chdff(k,js)=chdff(k,js)+ylhff(k,lh,nd)*chlh(jr,lh,js)
+                   chdtf(k,js)=chdtf(k,js)+ylhtf(k,lh,nd)*chlh(jr,lh,js)
+                ENDDO
+             ENDDO ! lh
+          ENDDO   ! js
+          
+       
+          CALL mkgylm(jspins,atoms%rmsh(jr,n),thet,nsp,&
+               ch(kt+1:,:),chdr,chdt,chdf,chdrr,chdtt,chdff,chdtf,chdrt,chdrf,grad,kt)
+       ENDIF
+       kt=kt+nsp
+    END DO
+    !Set charge to minimum value
+    ch(:,:)=MAX(ch(:,:),d_15)
+    
+  END SUBROUTINE mt_to_grid
+
+
+  SUBROUTINE mt_from_grid(atoms,sphhar,nsp,n,jspins,v_in,vr)
+    IMPLICIT NONE
+    TYPE(t_atoms),INTENT(IN) :: atoms
+    TYPE(t_sphhar),INTENT(IN):: sphhar
+    INTEGER,INTENT(IN)       :: nsp,jspins,n
+    REAL,INTENT(IN)          :: v_in(:,:)
+    REAL,INTENT(INOUT)       :: vr(:,0:,:,:)
+    
+    REAL    :: vpot(nsp),vlh
+    INTEGER :: js,kt,lh,jr,nd
+    nd=atoms%ntypsy(SUM(atoms%neq(:n-1))+1)
+    DO js = 1,jspins
+       !
+       kt=0
+       DO jr=1,atoms%jri(n)
+          vpot=v_in(kt+1:kt+nsp,js)*wt(:)!  multiplicate v_in with the weights of the k-points
+          
+          DO lh = 0,sphhar%nlh(nd)
+             !
+             ! --->        determine the corresponding potential number
+             !c            through gauss integration
+             !
+             vlh=dot_PRODUCT(vpot(:),ylh(:nsp,lh,nd))
+             vr(jr,lh,n,js) = vr(jr,lh,n,js) + vlh
+          ENDDO ! lh
+       ENDDO   ! jr
+       kt=kt+nsp
+    ENDDO
+
+  END SUBROUTINE mt_from_grid
+    
+  SUBROUTINE finish_mt_grid()
+    DEALLOCATE(ylh,ylht,ylhtt)
+    DEALLOCATE(ylhf,ylhff,ylhtf)
+    DEALLOCATE(wt,rx,thet)
+  END SUBROUTINE finish_mt_grid
+
+END MODULE m_mt_tofrom_grid

vgen/vgen_xcpot.F90

@@ -19,8 +19,7 @@ CONTAINS
     !     ***********************************************************
     USE m_constants
     USE m_intnv
-    USE m_vmtxcg
-    USE m_vmtxc
+    USE m_vmt_xc
     USE m_vvacxc
     USE m_vvacxcg
     USE m_visxc
@@ -169,13 +168,9 @@ CONTAINS
 #ifdef CPP_MPI
     CALL MPI_BCAST(den%mt,atoms%jmtd*(1+sphhar%nlhd)*atoms%ntype*dimension%jspd,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
 #endif
-    IF (xcpot%is_gga()) THEN
-       CALL vmtxcg(dimension,mpi,sphhar,atoms, den,xcpot,input,sym,&
-            obsolete, vTot,vx,exc)
-    ELSE
-       CALL vmtxc(DIMENSION,sphhar,atoms, den,xcpot,input,sym, vTot,exc,vx)
-    ENDIF
-
+    CALL vmt_xc(DIMENSION,mpi,sphhar,atoms, den,xcpot,input,sym,&
+         obsolete, vTot,vx,exc)
+    
 
     !
     ! add MT EXX potential to vr

vgen/vmt_xc.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_vmt_xc
+  !.....------------------------------------------------------------------
+  !     Calculate the GGA xc-potential in the MT-spheres
+  !.....------------------------------------------------------------------
+  !     instead of vmtxcor.f: the different exchange-correlation
+  !     potentials defined through the key icorr are called through
+  !     the driver subroutine vxcallg.f, subroutines vectorized
+  !     ** r.pentcheva 22.01.96
+  !     *********************************************************
+  !     angular mesh calculated on speacial gauss-legendre points
+  !     in order to use orthogonality of lattice harmonics and
+  !     avoid a least square fit
+  !     ** r.pentcheva 04.03.96
+  !     *********************************************************
+  !     MPI and OpenMP parallelization
+  !             U.Alekseeva, February 2017
+  !     *********************************************************
+
+CONTAINS
+  SUBROUTINE vmt_xc(DIMENSION,mpi,sphhar,atoms,&
+       den,xcpot,input,sym, obsolete,vxc,vx,exc)
+
+#include"cpp_double.h"
+   
+    USE m_mt_tofrom_grid
+    USE m_types_xcpot_inbuild
+    USE m_types
+    IMPLICIT NONE
+
+    CLASS(t_xcpot),INTENT(IN)      :: xcpot
+    TYPE(t_dimension),INTENT(IN)   :: dimension
+    TYPE(t_mpi),INTENT(IN)         :: mpi
+    TYPE(t_obsolete),INTENT(IN)    :: obsolete
+    TYPE(t_input),INTENT(IN)       :: input
+    TYPE(t_sym),INTENT(IN)         :: sym
+    TYPE(t_sphhar),INTENT(IN)      :: sphhar
+    TYPE(t_atoms),INTENT(IN)       :: atoms
+    TYPE(t_potden),INTENT(IN)      :: den
+    TYPE(t_potden),INTENT(INOUT)   :: vxc,vx,exc
+#ifdef CPP_MPI
+    include "mpif.h"
+#endif
+    !     ..
+    !     .. Local Scalars ..
+    TYPE(t_gradients)     :: grad
+    TYPE(t_xcpot_inbuild) :: xcpot_tmp
+    REAL, ALLOCATABLE     :: ch(:,:)
+    INTEGER               :: n,nsp,nt,jr
+    REAL                  :: divi
+   
+    !     ..
+ 
+    !locals for mpi
+    integer :: ierr
+    integer:: n_start,n_stride
+    REAL:: v_x((atoms%lmaxd+1+MOD(atoms%lmaxd+1,2))*(2*atoms%lmaxd+1),input%jspins)
+    REAL:: v_xc((atoms%lmaxd+1+MOD(atoms%lmaxd+1,2))*(2*atoms%lmaxd+1),input%jspins)
+    REAL:: e_xc((atoms%lmaxd+1+MOD(atoms%lmaxd+1,2))*(2*atoms%lmaxd+1),1)
+    REAL:: xcl(DIMENSION%nspd,DIMENSION%jspd)
+    LOGICAL :: lda_atom(atoms%ntype)
+    !.....------------------------------------------------------------------
+    lda_atom=.FALSE.
+    SELECT TYPE(xcpot)
+    TYPE IS(t_xcpot_inbuild)
+       lda_atom=xcpot%lda_atom
+       IF (ANY(lda_atom)) THEN
+          IF((.NOT.xcpot%is_name("pw91"))) &
+               CALL judft_warn("Using locally LDA only possible with pw91 functional")
+          CALL xcpot_tmp%init("l91",.FALSE.,atoms%ntype)
+       ENDIF
+    END SELECT
+   
+    nsp=(atoms%lmaxd+1+MOD(atoms%lmaxd+1,2))*(2*atoms%lmaxd+1)
+    ALLOCATE(ch(nsp*atoms%jmtd,input%jspins))
+    IF (xcpot%is_gga()) CALL xcpot%alloc_gradients(SIZE(ch,1),input%jspins,grad)
+    
+    CALL init_mt_grid(nsp,input%jspins,atoms,sphhar,xcpot,sym,xcpot%is_gga())
+  
+
+#ifdef CPP_MPI
+    n_start=mpi%irank+1
+    n_stride=mpi%isize
+#else
+    n_start=1
+    n_stride=1
+#endif
+
+    DO n = n_start,atoms%ntype,n_stride
+       CALL mt_to_grid(atoms,sphhar,den,input%jspins,n,xcpot%is_gga(),ch,grad)
+       !
+       !         calculate the ex.-cor. potential
+       CALL xcpot%get_vxc(input%jspins,ch(:nsp,:),v_xc,v_x,grad)
+       IF (lda_atom(n)) THEN
+          ! Use local part of pw91 for this atom
+          CALL xcpot_tmp%get_vxc(input%jspins,ch(:nsp,:),xcl,v_x,grad)
+          !Mix the potentials
+          divi = 1.0 / (atoms%rmsh(atoms%jri(n),n) - atoms%rmsh(1,n))
+          nt=0
+          DO jr=1,atoms%jri(n)
+             v_xc(nt+1:nt+nsp,:) = ( xcl(nt+1:nt+nsp,:) * ( atoms%rmsh(atoms%jri(n),n) - atoms%rmsh(jr,n) ) +&
+                  v_xc(nt+1:nt+nsp,:) * ( atoms%rmsh(jr,n) - atoms%rmsh(1,n) ) ) * divi
+             nt=nt+nsp
+          ENDDO
+       ENDIF
+       
+       CALL mt_from_grid(atoms,sphhar,nsp,n,input%jspins,v_xc,vxc%mt)
+       CALL mt_from_grid(atoms,sphhar,nsp,n,input%jspins,v_x,vx%mt)
+       
+       IF (ALLOCATED(exc%mt)) THEN
+          !
+          !           calculate the ex.-cor energy density
+          !
+          CALL xcpot%get_exc(input%jspins,ch(:nsp,:),e_xc(:,1),grad)
+          IF (lda_atom(n)) THEN
+             ! Use local part of pw91 for this atom
+             CALL xcpot_tmp%get_exc(input%jspins,ch(:nsp,:),xcl(:,1),grad)
+             !Mix the potentials
+             nt=0
+             DO jr=1,atoms%jri(n)
+                e_xc(nt+1:nt+nsp,1) = ( xcl(nt+1:nt+nsp,1) * ( atoms%rmsh(atoms%jri(n),n) - atoms%rmsh(jr,n) ) +&
+                     e_xc(nt+1:nt+nsp,1) * ( atoms%rmsh(jr,n) - atoms%rmsh(1,n) ) ) * divi
+                nt=nt+nsp
+             END DO
+          ENDIF
+          CALL mt_from_grid(atoms,sphhar,nsp,n,input%jspins,e_xc,exc%mt)
+       ENDIF
+    ENDDO
+    
+#ifdef CPP_MPI
+    CALL MPI_ALLREDUCE(vxr_local,vx%mt,atoms%jmtd*(1+sphhar%nlhd)*atoms%ntype*DIMENSION%jspd,CPP_MPI_REAL,MPI_SUM,mpi%mpi_comm,ierr)     !ToDo:CPP_MPI_REAL?
+    !using vxr_local as a temporal buffer
+    CALL MPI_ALLREDUCE(vr_local,vxr_local,atoms%jmtd*(1+sphhar%nlhd)*atoms%ntype*DIMENSION%jspd,CPP_MPI_REAL,MPI_SUM,mpi%mpi_comm,ierr)    
+    vxc%mt = vxc%mt + vxr_local
+    CALL MPI_ALLREDUCE(excr_local,exc%mt(:,:,:,1),atoms%jmtd*(1+sphhar%nlhd)*atoms%ntype,CPP_MPI_REAL,MPI_SUM,mpi%mpi_comm,ierr)    
+#endif
+    !
+    
+    RETURN
+  END SUBROUTINE vmt_xc
+END MODULE m_vmt_xc

vgen/vmtxcg.F90

@@ -26,7 +26,7 @@ CONTAINS
        den,xcpot,input,sym, obsolete,vxc,vx,exc)
 
 #include"cpp_double.h"
-    USE m_lhglptg
+   
     USE m_grdchlh
     USE m_mkgylm
     USE m_gaussp
@@ -236,26 +236,9 @@ CONTAINS
           CALL xcpot%alloc_gradients(nsp,input%jspins,grad)
           CALL mkgylm(input%jspins,atoms%rmsh(jr,n),thet,nsp,DIMENSION%nspd,DIMENSION%jspd,ch,chdr,&
                chdt,chdf,chdrr,chdtt,chdff,chdtf,chdrt,chdrf,grad)
-          
-          !
-          !         rhmnm: rho_minimum_muffin-tin..
-
-          rhmnm=10.e+10
-
-          DO js=1,input%jspins
-             DO i=1,nsp
-                ch(i,js) = max(ch(i,js),d_15)
-                rhmnm = min(rhmnm,ch(i,js))
-             ENDDO
-          ENDDO
-
-
-          IF (rhmnm.LT.obsolete%chng) THEN
-             WRITE (6,'(/'' rhmn.lt.obsolete%chng in vmtxc. rhmn,obsolete%chng='',&
-                  &        2d9.2)') rhmnm,obsolete%chng
-             !             CALL juDFT_error("vmtxcg: rhmn.lt.chng",calledby="vmtxcg")
-          ENDIF
 
+          !Set charge to minimum value
+          ch(:,:)=MAX(ch(:,:),d_15)
           !
           !         calculate the ex.-cor. potential