LCOV - code coverage report
Current view: top level - cdn - cdnovlp.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 184 255 72.2 %
Date: 2019-09-08 04:53:50 Functions: 6 6 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             :       MODULE m_cdnovlp
       8             :       USE m_juDFT
       9             :       IMPLICIT NONE
      10             :       PRIVATE
      11             :       PUBLIC :: cdnovlp 
      12             :       CONTAINS
      13         172 :         SUBROUTINE cdnovlp(mpi,&
      14             :              &                   sphhar,stars,atoms,sym,&
      15             :              &                   DIMENSION,vacuum,cell,&
      16             :              &                   input,oneD,l_st,&
      17         172 :              &                   jspin,rh,&
      18         172 :              &                   qpw,rhtxy,rho,rht)
      19             :           !*****************************************************************
      20             :           !     calculates the overlapping core tail density and adds
      21             :           !     its contribution to the corresponging components of
      22             :           !     valence density.      
      23             :           !
      24             :           !     OLD VERSION:
      25             :           !     The previous version to calculate the overlap of the
      26             :           !     core tails was done in real space:
      27             :           !     A three dimensional real space grid was set up. On each
      28             :           !     of these grid points the charge density of all atoms inside
      29             :           !     the unit cell and neigboring unit cells was calculated.
      30             :           !     This calculation required a lagrange fit since the core
      31             :           !     charge is on a radial mesh. The same was done in the vacuum
      32             :           !     region, except on each z-plane we worked on a two dimension
      33             :           !     grid. After the charge density was generated on a equidistant
      34             :           !     grid the charge density was FFTed into G-space. The set up
      35             :           !     of the charge density on the real space grid is rather time
      36             :           !     consuming. 3-loops are required for the 3D grid
      37             :           !                1-loop over all atoms in the unit cells
      38             :           !                Larange interpolation
      39             :           !                3D and 2D FFTs
      40             :           !     In order to save time the non-spherical contributions inside
      41             :           !     the sphere had been ignored. It turns out that the later
      42             :           !     approximation is pure in the context of force calculations.
      43             :           !     
      44             :           !     PRESENT VERSION:
      45             :           !     The present version is written from scratch. It is based on the
      46             :           !     idea that the FFT of an overlap of spherically symmetric
      47             :           !     charges can be expressed by the product of
      48             :           !
      49             :           !     sum_natype{ F(G,ntype) * sum_atom(atype) {S(\vec{G},atom)}}
      50             :           !
      51             :           !     of form factor F and structure factor S. The Form factor
      52             :           !     depends only G while the structure factor depends on \vec{G}
      53             :           !     and can build up in G-space. F of a gaussian chargedensity can
      54             :           !     be calculated analytically.
      55             : 
      56             :           !     The core-tails to the vacuum region are described by an
      57             :           !     exponentially decaying function into the vacuum:
      58             : 
      59             :           !     rho(r_||,z,ivac)= sum_n{ rho(n,ivac) * exp(-kappa*(z-z_v))
      60             :           !                                          * exp(iG(n)r_||) }
      61             : 
      62             :           !     And the plane waves are expanded into lattice harmonics
      63             :           !     up to a l_cutoff. Tests of the accuracy inside the sphere
      64             :           !     have shown that a reduction of the l_cutoff inside the 
      65             :           !     in order to save time leads to sizable errors and should 
      66             :           !     be omitted.
      67             : 
      68             :           !     rho_L(r) =  4 \pi i^l \sum_{g =|= 0}  \rho_int(g) r_i^{2} \times
      69             :           !                              j_{l} (gr_i) \exp{ig\xi_i} Y^*_{lm} (g)
      70             : 
      71             :           !     Tests have shown that the present version is about 10 times
      72             :           !     faster than the previous one also all nonspherical terms are
      73             :           !     included up to l=8 and the previous one included only l=0.
      74             :           !     For l=16 it is still faster by factor 2.
      75             : 
      76             :           !     coded                  Stefan Bl"ugel, IFF Nov. 1997
      77             :           !     tested                 RObert Abt    , IFF Dez. 1997
      78             :           !*****************************************************************
      79             :           !
      80             :           USE m_constants
      81             :           USE m_qpwtonmt
      82             :           USE m_cylbes
      83             :           USE m_dcylbs
      84             :           USE m_od_cylbes
      85             :           USE m_diflgr
      86             :           USE m_types
      87             : #ifdef CPP_MPI
      88             :           USE m_mpi_bc_st
      89             : #endif
      90             :           !
      91             :           !     .. Parameters ..
      92             :           TYPE(t_mpi),INTENT(IN)     :: mpi
      93             :           TYPE(t_sphhar),INTENT(IN)   :: sphhar
      94             :           TYPE(t_atoms),INTENT(IN)    :: atoms
      95             :           TYPE(t_stars),INTENT(IN)    :: stars
      96             :           TYPE(t_cell),INTENT(IN)     :: cell
      97             :           TYPE(t_sym),INTENT(IN)      :: sym
      98             :           TYPE(t_oneD),INTENT(IN)     :: oneD
      99             :           TYPE(t_dimension),INTENT(IN)::DIMENSION
     100             :           TYPE(t_vacuum),INTENT(in):: vacuum
     101             :           TYPE(t_input),INTENT(in)::input
     102             : 
     103             :           INTEGER    method1, method2
     104             :           PARAMETER (method1 = 2, method2 = 1)
     105             :           !     ..
     106             :           !     .. Scalar Arguments ..
     107             : 
     108             :           INTEGER,INTENT (IN) :: jspin       
     109             :           LOGICAL,INTENT (IN) :: l_st
     110             :           !     ..
     111             :           !     .. Array Arguments ..
     112             :           COMPLEX,INTENT (INOUT) :: qpw(stars%ng3,input%jspins)
     113             :           COMPLEX,INTENT (INOUT) :: rhtxy(vacuum%nmzxyd,oneD%odi%n2d-1,2,input%jspins)
     114             :           REAL,   INTENT (INOUT) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
     115             :           REAL,   INTENT (INOUT) :: rht(vacuum%nmzd,2,input%jspins)
     116             :           REAL,   INTENT (INOUT) :: rh(DIMENSION%msh,atoms%ntype)
     117             :           !     ..
     118             :           !     .. Local Scalars ..
     119             :           COMPLEX czero,carg,VALUE,slope,c_ph
     120             :           REAL    dif,dxx,g,gz,dtildh,&
     121             :                &        rkappa,sign,signz,tol_14,z,zero,zvac,&
     122             :                &        g2,phi,gamma,qq
     123             :           INTEGER ig3,imz,ivac,j,j1,k,kz,k1,k2,l_cutoff,m0,&
     124             :                &        n,nz,nrz,nzvac,&
     125             :                &        irec2,irec3,irec1,m,gzi
     126             :           !     ..
     127             :           !     .. Local Arrays ..
     128         172 :           COMPLEX, ALLOCATABLE :: qpwc(:)
     129         344 :           REAL    acoff(atoms%ntype),alpha(atoms%ntype),rho_out(2)
     130         344 :           REAL    rat(DIMENSION%msh,atoms%ntype)
     131         344 :           INTEGER mshc(atoms%ntype)
     132         344 :           REAL    fJ(-oneD%odi%M:oneD%odi%M),dfJ(-oneD%odi%M:oneD%odi%M)
     133             :           !     ..
     134             :           DATA  czero /(0.0,0.0)/, zero /0.0/, tol_14 /1.0e-10/!-14
     135             : #ifdef CPP_MPI
     136             :       EXTERNAL MPI_BCAST
     137             :       INTEGER ierr
     138             : #include "cpp_double.h"
     139             :       INCLUDE "mpif.h"
     140             : #endif
     141             :           !
     142             :           !----> Abbreviation
     143             :           !
     144             :           !     l_cutoff : cuts off the l-expansion of the nonspherical charge
     145             :           !                density due to the core-tails of neigboring atoms
     146             :           !     mshc     : maximal radial meshpoint for which the radial coretail
     147             :           !                density is larger than tol_14
     148             :           !     method1  : two different ways to calculate the derivative of the
     149             :           !                charge density at the sphere boundary.
     150             :           !                (1) use subroutine diflgr based on lagrange interpol.
     151             :           !                (2) use two point formular in real space, 
     152             :           !                    see notes of SB.
     153             :           !                Tests have shown that (2) is more accurate.
     154             :           !     method2  : two different integration routines to calculate form
     155             :           !                factor of coretails outside the sphere.
     156             :           !                (1) use subroutine intgrz to integrate the tails from
     157             :           !                    outside to inside.
     158             :           !                (2) use subroutine intgr3 to integrate the tails from
     159             :           !                    muffin-tin radius to outside and include correction
     160             :           !                    for start up.
     161             :           !                Tests have shown that (1) is more accurate.
     162             :           !           
     163             :           !
     164             : 
     165         172 :           ALLOCATE (qpwc(stars%ng3))
     166             :           !
     167             :           !----> prepare local array to store pw-expansion of pseudo core charge
     168             :           !
     169      949352 :           DO k = 1 , stars%ng3    
     170      474762 :              qpwc(k) = czero
     171             :           ENDDO
     172             :           ! 
     173             :           !----> (1) set up radial mesh beyond muffin-tin radius
     174             :           !      (2) cut_off core tails from noise 
     175             :           !
     176             : #ifdef CPP_MPI
     177         172 :           CALL MPI_BCAST(rh,DIMENSION%msh*atoms%ntype,CPP_MPI_REAL,0,mpi%mpi_comm,ierr)
     178             : #endif
     179         556 :           nloop: DO  n = 1 , atoms%ntype
     180         556 :               IF ((atoms%ncst(n).GT.0).OR.l_st) THEN
     181      260552 :                    DO  j = 1 , atoms%jri(n)
     182      260552 :                       rat(j,n) = atoms%rmsh(j,n)
     183             :                    ENDDO
     184         384 :                    dxx = EXP(atoms%dx(n))
     185       76584 :                    DO j = atoms%jri(n) + 1 , DIMENSION%msh
     186       76584 :                       rat(j,n) = rat(j-1,n)*dxx
     187             :                    ENDDO
     188       77352 :                    DO j = atoms%jri(n) - 1 , DIMENSION%msh
     189       77352 :                       rh(j,n) = rh(j,n)/ (fpi_const*rat(j,n)*rat(j,n))
     190             :                    ENDDO
     191      108508 :                    DO j = DIMENSION%msh , atoms%jri(n) , -1
     192       54446 :                       IF ( rh(j,n) .GT. tol_14 ) THEN
     193         312 :                          mshc(n) = j
     194         312 :                          CYCLE nloop
     195             :                       END IF
     196             :                    ENDDO
     197          72 :                    mshc(n) = atoms%jri(n)
     198             :               ENDIF
     199             :           ENDDO nloop
     200             :           !
     201             :           !-----> the core density inside the spheres is replaced by a
     202             :           !       gaussian-like pseudo density : n(r) = acoff*exp(-alpha*r*r)
     203             :           !       acoff and alpha determined to obtain a continous and 
     204             :           !       differentiable density at the sphere boundary.
     205             :           !       IF mshc = jri  either core tail too small or no core (i.e. H)
     206             :           !
     207         940 :           DO  n = 1,atoms%ntype
     208         556 :               IF ((mshc(n).GT.atoms%jri(n)).AND.((atoms%ncst(n).GT.0).OR.l_st)) THEN
     209             : 
     210         312 :                    j1 = atoms%jri(n) - 1
     211             :                    IF ( method1 .EQ. 1) THEN
     212             :                       dif = diflgr(rat(j1,n),rh(j1,n))
     213             :                       WRITE (6,FMT=8000) n,rh(atoms%jri(n),n),dif
     214             :                       alpha(n) = -0.5 * dif / ( rh(atoms%jri(n),n)*atoms%rmt(n) )
     215             :                    ELSEIF ( method1 .EQ. 2) THEN
     216         312 :                       alpha(n) = LOG( rh(j1,n) / rh(atoms%jri(n),n) )
     217             :                       alpha(n) = alpha(n)&
     218         312 :                            &                   / ( atoms%rmt(n)*atoms%rmt(n)*( 1.0-EXP( -2.0*atoms%dx(n) ) ) )
     219             :                    ELSE
     220             :                       WRITE (6,'('' error in choice of method1 in cdnovlp '')')
     221             :                       CALL juDFT_error("error in choice of method1 in cdnovlp"&
     222             :                            &              ,calledby ="cdnovlp")
     223             :                    ENDIF
     224         312 :                    acoff(n) = rh(atoms%jri(n),n) * EXP( alpha(n)*atoms%rmt(n)*atoms%rmt(n) )
     225             :                    !WRITE (6,FMT=8010) alpha(n),acoff(n)
     226      206144 :                    DO j = 1,atoms%jri(n) - 1
     227      206144 :                       rh(j,n) = acoff(n) * EXP( -alpha(n)*rat(j,n)**2 )
     228             :                    ENDDO
     229             : 
     230             :                 ELSE
     231          72 :                    alpha(n) = 0.0
     232             :               ENDIF
     233             :           ENDDO
     234             :           !
     235             :           IF (mpi%irank ==0) THEN
     236             : 8000         FORMAT (/,10x,'core density and its first derivative ',&
     237             :                   &                 'at sph. bound. for atom type',&
     238             :                   &             i2,' is',3x,2e15.7)
     239             : 8010         FORMAT (/,10x,'alpha=',f10.5,5x,'acoff=',f10.5)
     240             :           END IF          
     241             :           !
     242             :           !=====> calculate the fourier transform of the core-pseudocharge
     243             : 
     244             :           CALL ft_of_CorePseudocharge(mpi,DIMENSION,atoms,mshc,alpha,tol_14,rh, &
     245         172 :                           acoff,stars,method2,rat,cell,oneD,sym,qpwc)
     246             : 
     247      474762 :           DO k = 1 , stars%ng3    
     248      474762 :               qpw(k,jspin) = qpw(k,jspin) + qpwc(k) 
     249             :           ENDDO
     250             : 
     251         172 :           IF (mpi%irank ==0) THEN
     252             :              !
     253             :              !=====> calculate core-tails to the vacuum region                
     254             :              !       Coretails expanded in exponentially decaying functions.
     255             :              !       Describe vacuum by: 
     256             :              !       rho(r_||,z,ivac)= sum_n{ rho(n,ivac) * exp(-kappa*(z-z_v)) 
     257             :              !                                           * exp(iG(n)r_||) }
     258          86 :              IF (input%film .AND. .NOT.oneD%odi%d1) THEN
     259             :                 !+gu
     260          15 :                 dtildh = 0.5 * tpi_const / cell%bmat(3,3)
     261          15 :                 IF (vacuum%nvac.EQ.1) THEN
     262          15 :                    rho_out(1) = qpwc(1)*cell%z1
     263       53916 :                    DO k = 2,stars%ng3
     264       53916 :                       IF ((stars%kv3(1,k).EQ.0).AND.(stars%kv3(2,k).EQ.0)) THEN
     265         315 :                          nz = 1
     266         315 :                          IF (sym%invs.OR.sym%zrfs) nz = 2
     267         315 :                          g = stars%kv3(3,k) * cell%bmat(3,3)
     268         315 :                          rho_out(1) = rho_out(1) + nz*qpwc(k)*SIN(g*cell%z1)/g
     269             :                       ENDIF
     270             :                    ENDDO
     271          15 :                    rho_out(1) =  qpwc(1) * dtildh - rho_out(1)
     272             :                 ELSE
     273           0 :                    DO ivac = 1, vacuum%nvac
     274             :                       carg = CMPLX(0.0,0.0)
     275           0 :                       DO k = 2,stars%ng3
     276           0 :                          IF ((stars%kv3(1,k).EQ.0).AND.(stars%kv3(2,k).EQ.0)) THEN
     277           0 :                             g = stars%kv3(3,k) * cell%bmat(3,3) * (3. - 2.*ivac)
     278           0 :                             carg = carg -qpwc(k)*(EXP(ImagUnit*g*dtildh)-EXP(ImagUnit*g*cell%z1))/g
     279             :                          ENDIF
     280             :                       ENDDO
     281           0 :                       rho_out(ivac) = qpwc(1) * ( dtildh-cell%z1 ) - AIMAG(carg)
     282             :                    ENDDO
     283             :                 ENDIF
     284             :                 !-gu
     285             :                 !        nzvac = min(50,nmz)
     286          15 :                 m0 = -stars%mx3
     287          15 :                 IF (sym%zrfs) m0 = 0
     288             :                 !
     289             :                 !---> loop over 2D stars
     290             :                 !
     291        4090 :                 DO 280 k = 1,stars%ng2
     292        4075 :                    k1 = stars%kv2(1,k)
     293        4075 :                    k2 = stars%kv2(2,k)
     294        8150 :                    DO 270 ivac = 1,vacuum%nvac
     295        4075 :                       VALUE = czero
     296        4075 :                       slope = czero
     297        4075 :                       sign = 3. - 2.*ivac
     298             :                       !
     299             :                       ! ---> sum over gz-stars
     300      143090 :                       DO 250 kz = m0,stars%mx3
     301      139015 :                          ig3 = stars%ig(k1,k2,kz)
     302      139015 :                          c_ph = stars%rgphs(k1,k2,kz) ! phase factor for invs=T & zrfs=F
     303             :                          !        ----> use only stars within the g_max sphere (oct.97 shz)
     304      139015 :                          IF (ig3.NE.0) THEN
     305       89609 :                             nz = 1
     306       89609 :                             IF (sym%zrfs) nz = stars%nstr(ig3)/stars%nstr2(k)
     307       89609 :                             gz = kz*cell%bmat(3,3)
     308      193366 :                             DO 240 nrz = 1,nz
     309      103757 :                                signz = 3. - 2.*nrz
     310      103757 :                                carg = ImagUnit*sign*signz*gz
     311      103757 :                                VALUE = VALUE + c_ph*qpwc(ig3)* EXP(carg*cell%z1)
     312      103757 :                                slope = slope + c_ph*carg*qpwc(ig3)* EXP(carg*cell%z1)
     313       89609 : 240                         ENDDO
     314             :                          END IF
     315        4075 : 250                   ENDDO
     316             :                       ! roa work-around
     317        4075 :                       IF (  ABS(REAL(VALUE)).GT.zero ) THEN
     318             :                          ! roa work-around
     319             :                          ! gb works also around
     320        4075 :                          rkappa = - REAL( slope/VALUE )
     321        4075 :                          IF (k.EQ.1) rkappa = VALUE/rho_out(ivac)
     322             :                          !               rkappa = - sign * real( slope/value )
     323        4075 :                          IF (rkappa.LE.zero) rkappa=MIN(rkappa,-tol_14)
     324        4075 :                          IF (rkappa.GT.zero) rkappa=MAX(rkappa,tol_14)
     325             :                          ! gb works also around
     326        4075 :                          zvac   = - LOG( tol_14/cabs(VALUE) ) / rkappa
     327        4075 :                          zvac   = MIN (2.*vacuum%nmz,abs(zvac)) ! avoid int-overflow in next line
     328        4075 :                          nzvac  = INT( zvac/vacuum%delz ) + 1
     329             :                          !               IF ( rkappa.GT.zero .AND. real(value).GT.zero ) THEN
     330        4075 :                          IF ( rkappa.GT.zero ) THEN
     331        2264 :                             z = 0. 
     332        2264 :                             IF ( k.EQ.1 ) THEN
     333        1231 :                                DO imz = 1 , MIN( nzvac,vacuum%nmz )
     334             :                                   rht(imz,ivac,jspin) = &
     335        1216 :                                        &                  rht(imz,ivac,jspin) + VALUE*EXP(-rkappa*z)
     336        1216 :                                   z = z + vacuum%delz
     337          15 : 220                            ENDDO
     338             :                             ELSE
     339       53769 :                                DO imz = 1 , MIN( nzvac,vacuum%nmzxy )
     340             :                                   rhtxy(imz,k-1,ivac,jspin) = &
     341       51520 :                                        &                  rhtxy(imz,k-1,ivac,jspin) + VALUE*EXP(-rkappa*z)
     342       51520 :                                   z = z + vacuum%delz
     343        2249 : 230                            ENDDO
     344             :                             END IF
     345             :                          END IF
     346             :                          ! roa work-around
     347             :                       END IF
     348             :                       ! roa work-around
     349        4075 : 270                ENDDO
     350          15 : 280             ENDDO
     351          71 :              ELSEIF (oneD%odi%d1) THEN
     352             :                 !-odim
     353             :                 !--->  rho_out is the charge lost between the interstitial and the
     354             :                 !--->  rectangular unit cell
     355             : 
     356           0 :                 rho_out(1) = cell%vol*qpwc(1)
     357             : 
     358           0 :                 DO k = 2,stars%ng3
     359           0 :                    IF (stars%kv3(3,k).EQ.0) THEN
     360           0 :                       irec3 = stars%ig(stars%kv3(1,k),stars%kv3(2,k),stars%kv3(3,k))
     361           0 :                       IF (irec3.NE.0) THEN
     362           0 :                          irec2 = stars%ig2(irec3)
     363           0 :                          IF (irec2.NE.1) THEN
     364           0 :                             g2 = stars%sk2(irec2)
     365           0 :                             CALL cylbes(oneD%odi%M,g2*cell%z1,fJ)
     366             :                             rho_out(1) = rho_out(1) +&
     367           0 :                                  &                    qpwc(k)*2.*cell%vol*fJ(1)/(cell%z1*g2)
     368             :                          END IF
     369             :                       END IF
     370             :                    END IF
     371             :                 END DO
     372             : 
     373           0 :                 rho_out(1) = cell%bmat(3,3)*(qpwc(1)*cell%omtil - rho_out(1))/(tpi_const*tpi_const)
     374             :                 !     then we are constructing our radial components of the vacuum
     375             :                 !     charge density so that the so that they are continuous in the
     376             :                 !     value and slope on the boundary
     377             :                 !     value = sum{gpar}[qpw(gz,gpar)exp{-i*m*phi(gpar))J_m(gpar*z1)]
     378             :                 !     slope = .......the same...................*gpar*dJ_m(gpar*z1)]
     379             :                 !     for determining the rht we need only continuity, because the
     380             :                 !     second condition is the charge neutrality of the system
     381             :                 !                         Y.Mokrousov
     382           0 :                 DO irec1 = 1,oneD%odi%nq2
     383           0 :                    VALUE = czero
     384           0 :                    slope = czero
     385           0 :                    gzi = oneD%odi%kv(1,irec1)
     386           0 :                    m = oneD%odi%kv(2,irec1)
     387           0 :                    DO irec2 = 1,stars%ng2
     388           0 :                       irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),gzi)
     389           0 :                       IF (irec3.NE.0) THEN
     390           0 :                          g2 = stars%sk2(irec2)
     391           0 :                          phi = stars%phi2(irec2)
     392           0 :                          CALL cylbes(oneD%odi%M,g2*cell%z1,fJ)
     393           0 :                          CALL dcylbs(oneD%odi%M,g2*cell%z1,fJ,dfJ)
     394             :                          VALUE = VALUE + (ImagUnit**m)*qpwc(irec3)*&
     395           0 :                               &                 EXP(CMPLX(0.,-m*phi))*fJ(m)
     396             :                          slope = slope + (ImagUnit**m)*g2*qpwc(irec3)*&
     397           0 :                               &                 EXP(CMPLX(0.,-m*phi))*dfJ(m)
     398             :                       END IF
     399             :                    END DO
     400             : 
     401           0 :                    IF (ABS(REAL(VALUE)).GT.zero) THEN
     402           0 :                       IF (irec1.EQ.1) THEN
     403           0 :                          qq = REAL(VALUE/rho_out(1))
     404           0 :                          rkappa = (cell%z1*qq+SQRT(cell%z1*cell%z1*qq*qq + 4*qq))/2.
     405           0 :                          gamma = rkappa
     406             :                       ELSE
     407             :                          rkappa = gamma
     408             :                          !                  rkappa = -real(slope/value)
     409             :                       END IF
     410           0 :                       IF (rkappa.LE.zero) rkappa=MIN(rkappa,-tol_14)
     411           0 :                       IF (rkappa.GT.zero) rkappa=MAX(rkappa,tol_14)
     412           0 :                       zvac=-LOG(tol_14/cabs(VALUE))/rkappa
     413           0 :                       nzvac=INT(zvac/vacuum%delz)+1
     414           0 :                       IF (rkappa.GT.zero .AND. REAL(VALUE).GT.zero) THEN
     415           0 :                          z = 0.0
     416           0 :                          IF (irec1.EQ.1) THEN
     417           0 :                             DO imz = 1,MIN(nzvac,vacuum%nmz)
     418             :                                rht(imz,vacuum%nvac,jspin)=rht(imz,vacuum%nvac,jspin)+&
     419           0 :                                     &                       VALUE*EXP(-rkappa*z)
     420           0 :                                z = z + vacuum%delz
     421             :                             END DO
     422             :                          ELSE
     423           0 :                             DO imz = 1,MIN(nzvac,vacuum%nmzxy)
     424             :                                rhtxy(imz,irec1-1,vacuum%nvac,jspin)=&
     425             :                                     &                          rhtxy(imz,irec1-1,vacuum%nvac,jspin)+&
     426           0 :                                     &                          VALUE*EXP(-rkappa*z)
     427           0 :                                z = z + vacuum%delz
     428             :                             END DO
     429             :                          END IF
     430             :                       END IF
     431             :                    END IF
     432             :                 END DO
     433             :                 !+odim
     434             : 
     435             :              END IF
     436             :              !
     437             :              !=====> update density inside the spheres                        
     438             :              !
     439             :              ! ----> (1) subtract on-site contribution, because 
     440             :              !           they are contained in the plane wave part 
     441             :              !
     442         278 :              DO n = 1,atoms%ntype
     443         278 :                 IF  ((mshc(n).GT.atoms%jri(n)).AND.((atoms%ncst(n).GT.0).OR.l_st)) THEN
     444      206300 :                    DO j = 1,atoms%jri(n)
     445             :                       rho(j,0,n,jspin) = rho(j,0,n,jspin)&
     446      103228 :                            &                          - sfp_const*rat(j,n)*rat(j,n)*rh(j,n)
     447             :                    ENDDO
     448             :                 ENDIF
     449             :              ENDDO
     450             : !
     451             :              ! ----> (2) add the plane wave contribution of (core tails + on-site 
     452             :              !           contribution) to the m.t. density, include full nonspherical 
     453             :              !           components
     454             :              !
     455             :           ENDIF ! mpi%irank ==0
     456         172 :           l_cutoff=input%coretail_lmax
     457             : #ifdef CPP_MPI
     458         172 :           IF ( mpi%isize > 1 ) CALL mpi_bc_st(mpi,stars,qpwc)
     459             : #endif
     460             : 
     461             :           CALL qpw_to_nmt(&
     462             :                &                sphhar,atoms,stars,&
     463             :                &                sym,cell,oneD,mpi,&
     464             :                &                jspin,l_cutoff,qpwc,&
     465         172 :                &                rho)
     466             : 
     467             : #ifdef CPP_MPI
     468         172 :           IF ( mpi%isize > 1) CALL mpi_col_st(mpi,atoms,sphhar,rho(1,0,1,jspin))
     469             : #endif
     470             : 
     471         172 :           DEALLOCATE (qpwc)
     472             : 
     473         172 :         END SUBROUTINE cdnovlp
     474             : 
     475             : !***********************************************************************
     476             : !     INTERNAL SUBROUTINES
     477             : !***********************************************************************
     478             : 
     479         344 :       subroutine ft_of_CorePseudocharge(mpi,DIMENSION,atoms,mshc,alpha,&
     480         344 :             tol_14,rh,acoff,stars,method2,rat,cell,oneD,sym,qpwc)
     481             : 
     482             :       !=====> calculate the fourier transform of the core-pseudocharge
     483             :       !
     484             :       !     qpw(\vec{G}) = Sum_{n} [ F(G,n) * Sum_{atm{n}} S(\vec{G},atm) ]
     485             :       !                  n = atom_type
     486             :       !                  F = Formfactor = F_in_sphere + F_outsphere
     487             :       !                  S = Structure factor
     488             : 
     489             :       USE m_types
     490             : 
     491             :       type(t_mpi)      ,intent(in) :: mpi
     492             :       type(t_dimension),intent(in) :: DIMENSION
     493             :       type(t_atoms)    ,intent(in) :: atoms
     494             :       integer          ,intent(in) :: mshc(atoms%ntype)
     495             :       real             ,intent(in) :: alpha(atoms%ntype), tol_14
     496             :       real             ,intent(in) :: rh(DIMENSION%msh,atoms%ntype)
     497             :       real             ,intent(in) :: acoff(atoms%ntype)
     498             :       type(t_stars)    ,intent(in) :: stars
     499             :       integer          ,intent(in) :: method2
     500             :       real             ,intent(in) :: rat(DIMENSION%msh,atoms%ntype)
     501             :       type(t_cell)     ,intent(in) :: cell
     502             :       type(t_oneD)     ,intent(in) :: oneD
     503             :       type(t_sym)      ,intent(in) :: sym
     504             :       complex         ,intent(out) :: qpwc(stars%ng3)
     505             : 
     506             : !     ..Local variables
     507             :       integer nat1, n, n_out_p, k
     508             :       complex czero
     509             : 
     510             : !     ..Local arrays
     511         344 :       real :: qf(stars%ng3)
     512         344 :       complex qpwc_at(stars%ng3)
     513             : #ifdef CPP_MPI
     514             :       external mpi_bcast
     515         344 :       complex :: qpwc_loc(stars%ng3)
     516             :       integer :: ierr
     517             : #include "cpp_double.h"
     518             :       include "mpif.h"
     519             : #endif
     520             : 
     521         172 :       czero = (0.0,0.0)
     522             : #ifdef CPP_MPI
     523      474762 :       DO k = 1 , stars%ng3
     524         172 :          qpwc_loc(k) = czero
     525             :       ENDDO
     526             : #endif
     527      949352 :       DO k = 1 , stars%ng3    
     528      474762 :           qpwc(k) = czero
     529             :       ENDDO
     530             : 
     531             :       !
     532             :       !*****> start loop over the atom type
     533             :       !
     534         364 :       DO  n = 1 + mpi%irank, atoms%ntype, mpi%isize
     535         192 :           IF ( ( mshc(n) .GT. atoms%jri(n) ).AND.&
     536         172 :               &        ( alpha(n) .GT. tol_14 ) )    THEN
     537             :                    
     538         156 :               n_out_p = mshc(n)-atoms%jri(n)+1
     539             :               
     540             :               ! (1) Form factor for each atom type
     541             :              
     542             :               CALL FormFactor_forAtomType(DIMENSION,method2,n_out_p,&
     543             :                                  atoms%rmt(n),atoms%jri(n),atoms%dx(n),mshc(n),rat(:,n), &
     544         156 :                                  rh(:,n),alpha(n),stars,cell,acoff(n),qf)
     545             : 
     546             :               ! (2) structure constant for each atom of atom type
     547             :             
     548         156 :               nat1 = 1
     549         156 :               IF (n>1) THEN
     550         324 :                  DO k = 1, n-1
     551         203 :                     nat1 = nat1 + atoms%neq(k)
     552             :                  END DO
     553             :               END IF       
     554             :               CALL StructureConst_forAtom(nat1,stars,oneD,sym,&
     555             :                                  atoms%neq(n),atoms%nat,atoms%taual,&
     556         156 :                                  cell,qf,qpwc_at)
     557             : #ifdef CPP_MPI
     558      331655 :               DO k = 1, stars%ng3
     559      331655 :                  qpwc_loc(k) = qpwc_loc(k)  + qpwc_at(k)
     560             :               END DO
     561             : #else
     562             :               DO k = 1 , stars%ng3    
     563             :                  qpwc(k) = qpwc(k) + qpwc_at(k)
     564             :               END DO
     565             : #endif
     566             : 
     567             :           END IF
     568             :        ENDDO
     569             : #ifdef CPP_MPI
     570             :        CALL mpi_allreduce(qpwc_loc,qpwc,stars%ng3,CPP_MPI_COMPLEX,mpi_sum, &
     571         172 :                mpi%mpi_comm,ierr)
     572             : #endif
     573             : 
     574         172 :       end subroutine ft_of_CorePseudocharge
     575             : 
     576             : 
     577             : !----------------------------------------------------------------------
     578         156 :       subroutine StructureConst_forAtom(nat1,stars,oneD,sym,&
     579         156 :                           neq,natd,taual,cell,qf,qpwc_at)
     580             : !       calculates structure constant for each atom of atom type
     581             : 
     582             :       USE m_types
     583             :       USE m_spgrot
     584             :       USE m_constants
     585             :       USE m_od_chirot
     586             : 
     587             :        integer          ,intent(in) :: nat1
     588             :        type(t_stars)    ,intent(in) :: stars
     589             :        type(t_oneD)     ,intent(in) :: oneD
     590             :        type(t_sym)      ,intent(in) :: sym
     591             :        integer          ,intent(in) :: neq,natd
     592             :        real             ,intent(in) :: taual(3,natd)
     593             :        type(t_cell)     ,intent(in) :: cell
     594             :        real             ,intent(in) :: qf(stars%ng3)
     595             :        complex         ,intent(out) :: qpwc_at(stars%ng3)
     596             : 
     597             : !      ..Local variables
     598             :       integer k, nat2, nat, j
     599             :       real x
     600             :       complex sf, czero
     601             : 
     602             : !      ..Local arrays
     603         156 :        integer kr(3,sym%nop)
     604         156 :        real    kro(3,oneD%ods%nop)
     605         156 :        complex phas(sym%nop)
     606         156 :        complex phaso(oneD%ods%nop)
     607             : 
     608         156 :       czero = (0.0,0.0)
     609      331655 :       DO k = 1 , stars%ng3    
     610      331655 :           qpwc_at(k) = czero
     611             :       ENDDO
     612             : 
     613             :       !
     614             :       !    first G=0
     615             :       !
     616         156 :       k=1
     617         156 :       qpwc_at(k)      = qpwc_at(k)      + neq * qf(k)
     618             : 
     619             :       !
     620             :       !    then  G>0
     621             :       !
     622             : !$OMP PARALLEL DO DEFAULT(none) &
     623             : !$OMP SHARED(stars,oneD,sym,neq,natd,nat1,taual,cell,qf,qpwc_at) &
     624             : !$OMP FIRSTPRIVATE(czero) &
     625             : !$OMP PRIVATE(k,kr,phas,nat2,nat,sf,j,x,kro,phaso)
     626             :       DO  k = 2,stars%ng3
     627         312 :           IF (.NOT.oneD%odi%d1) THEN
     628             :               CALL spgrot(&
     629             :                    &                       sym%nop,sym%symor,sym%mrot,sym%tau,sym%invtab,&
     630             :                    &                       stars%kv3(:,k),&
     631      331343 :                    &                       kr,phas)
     632             :               !
     633             :               ! ----> start loop over equivalent atoms
     634             :               !
     635      331343 :                nat2 = nat1 + neq - 1
     636      708841 :                DO  nat = nat1,nat2
     637      377498 :                    sf = czero
     638     2675542 :                    DO  j = 1,sym%nop
     639             :                        x = -tpi_const* ( kr(1,j) * taual(1,nat)&
     640             :                            &          + kr(2,j) * taual(2,nat)&
     641     2298044 :                            &          + kr(3,j) * taual(3,nat))
     642             :                        !gb      sf = sf + CMPLX(COS(x),SIN(x))*phas(j)
     643     2675542 :                        sf = sf + CMPLX(COS(x),SIN(x))*conjg(phas(j))
     644             :                    ENDDO
     645      377498 :                    sf = sf / REAL( sym%nop )
     646      708841 :                    qpwc_at(k)      = qpwc_at(k)      + sf * qf(k)
     647             :                ENDDO
     648             :           ELSE
     649             :                !-odim
     650           0 :                CALL od_chirot(oneD%odi,oneD%ods,cell%bmat,stars%kv3(:,k),kro,phaso)
     651           0 :                nat2 = nat1 + neq - 1
     652           0 :                DO  nat = nat1,nat2
     653             :                    !                  sf = cmplx(1.,0.)
     654           0 :                    sf = czero
     655           0 :                    DO  j = 1,oneD%ods%nop
     656             :                        x = -tpi_const* ( kro(1,j)*taual(1,nat)&
     657             :                            &          + kro(2,j)*taual(2,nat)&
     658           0 :                            &          + kro(3,j)*taual(3,nat))
     659           0 :                        sf = sf + CMPLX(COS(x),SIN(x))*phaso(j)
     660             :                    ENDDO
     661           0 :                    sf = sf / REAL( oneD%ods%nop )
     662           0 :                    qpwc_at(k)      = qpwc_at(k)      + sf * qf(k)
     663             :                ENDDO
     664             :                !+odim
     665             :           END IF
     666             :       ENDDO
     667             : !$OMP END PARALLEL DO
     668         156 :       end subroutine StructureConst_forAtom
     669             : 
     670             : !----------------------------------------------------------------------
     671         156 :       subroutine FormFactor_forAtomType(DIMENSION,method2,n_out_p,&
     672         156 :                        rmt,jri,dx,mshc,rat,&
     673         156 :                        rh,alpha,stars,cell,acoff,qf)
     674             : 
     675             :       USE m_types
     676             :       USE m_constants
     677             :       USE m_rcerf
     678             :       USE m_intgr, ONLY : intgr3,intgz0
     679             : 
     680             :       type(t_dimension),intent(in) :: DIMENSION
     681             :       integer          ,intent(in) :: method2, n_out_p
     682             :       real             ,intent(in) :: rmt
     683             :       integer          ,intent(in) :: jri
     684             :       real             ,intent(in) :: dx
     685             :       integer          ,intent(in) :: mshc
     686             :       real             ,intent(in) :: rat(DIMENSION%msh)
     687             :       real             ,intent(in) :: rh(DIMENSION%msh)
     688             :       real             ,intent(in) :: alpha
     689             :       type(t_stars)    ,intent(in) :: stars
     690             :       type(t_cell)     ,intent(in) :: cell
     691             :       real             ,intent(in) :: acoff
     692             :       real            ,intent(out) :: qf(stars%ng3)
     693             : 
     694             : 
     695             : !     ..Local variables
     696             :       real f11, f12, ar, g, ai, qfin, qfout, gr, a4, alpha3, zero
     697             :       integer k, ir, j
     698             :       logical tail
     699             : 
     700             : !     ..Local arrays
     701         156 :       real rhohelp(DIMENSION%msh)
     702             : 
     703         156 :       zero = 0.0
     704      331655 :       DO k = 1,stars%ng3
     705      331655 :         qf(k) = 0.0
     706             :       END DO
     707             : 
     708         156 :       tail = .FALSE.
     709         156 :       f11 = tpi_const * rmt * rh(jri) / alpha
     710         156 :       f12 = acoff * ( pi_const/alpha ) *SQRT(pi_const/alpha)
     711         156 :       ar  = SQRT( alpha ) * rmt 
     712             : 
     713             : !$OMP PARALLEL DO DEFAULT(none) & 
     714             : !$OMP SHARED(stars,f11,f12,ar,method2,n_out_p,jri,rat,rh,dx,tail) &
     715             : !$OMP SHARED(alpha,cell,mshc,rmt,qf) &
     716             : !$OMP FIRSTPRIVATE(zero) &
     717             : !$OMP PRIVATE(k,g,ai,qfin,ir,j,rhohelp,qfout,gr,a4,alpha3)
     718             :       DO  k = 1,stars%ng3
     719      331499 :           g = stars%sk3(k)
     720             :           !    first G=0
     721      331499 :           IF ( k.EQ.1 ) THEN
     722         156 :               ai = zero
     723             :               !
     724             :               ! ---->     calculate form factor inside the mt-sphere
     725             :               !           (use analytic integration of gaussian)
     726             :               !
     727         156 :               qfin = - f11 + f12 * rcerf(ar,ai)
     728             :               !
     729             :               ! ---->     calculate form factor outside the mt-sphere
     730             :               !           (do numerical integration of tails)
     731             :               !
     732         156 :               IF ( method2 .EQ. 1) THEN
     733             :  
     734       11417 :                   DO ir = 1 , n_out_p
     735       11261 :                      j  = jri+ir-1
     736             :                      rhohelp(mshc+1-j) =  rat(j) * rat(j) &
     737       11417 :                            &                       * rat(j) *  rh(j)
     738             :                   END DO
     739         156 :                   CALL intgz0(rhohelp,dx,n_out_p,qfout,tail)
     740             :  
     741             :               ELSE
     742             :  
     743           0 :                   DO ir = 1 , n_out_p
     744           0 :                      j  = jri+ir-1
     745           0 :                      rhohelp(ir) = rat(j) * rat(j) * rh(j)
     746             :                   END DO
     747             :                   CALL intgr3(rhohelp,rat(jri),dx,&
     748           0 :                            &                        n_out_p,qfout)
     749             :                   !--->             have to remove the small r-correction from intgr3
     750           0 :                   qfout=qfout-rmt*rhohelp(1)
     751             :  
     752             :               END IF
     753             :  
     754         156 :               qfout = fpi_const * qfout
     755             :               !
     756             :           ELSE 
     757             :               !    then  G>0
     758      331343 :               ai = 0.5*g/SQRT(alpha)
     759      331343 :               gr = g*rmt
     760      331343 :               a4 = 0.25/alpha
     761             :               !
     762             :               ! ---->     calculate form factor inside the mt-sphere
     763             :               !           (use analytic integration of gaussian)
     764             :               !
     765             :               qfin = - f11 * SIN(gr)/gr &
     766      331343 :                            &                + f12 * rcerf(ar,ai) * EXP(-a4*g*g) 
     767             :               !
     768             :               ! ---->     calculate form factor outside the mt-sphere
     769             :               !           (do numerical integration of tails)
     770             : 
     771      331343 :               IF ( method2 .EQ. 1) THEN
     772             : 
     773    23812655 :                   DO ir = 1 , n_out_p 
     774    23481312 :                       j  = jri+ir-1
     775             :                       rhohelp(mshc-jri+2-ir) =  rat(j)*rat(j) &
     776    23812655 :                                     &                                     * rh(j) * SIN( g*rat(j) )
     777             :                   END DO
     778             :                   !
     779             :                   !--->       note we use here the integration routine for vacuum. Because 
     780             :                   !           the vacuum integration is made for an inwards integration 
     781             :                   !           from outside to inside. Outside the starting value will be 
     782             :                   !           nearly zero since the core density is small. if the tail 
     783             :                   !           correction (tail=.true.) is included for the integrals, the 
     784             :                   !           integrand is from infinity inward. This might not be 
     785             :                   !           necessary. Further the integration routine is made for 
     786             :                   !           equidistant meshpoints, therefore the term r(i) of
     787             :                   !           dr/di = dx*r(i) is included in rhohelp
     788             : 
     789             : 
     790      331343 :                   CALL intgz0(rhohelp,dx,n_out_p,qfout,tail)
     791             : 
     792             :               ELSE
     793             : 
     794           0 :                   DO ir = 1 , n_out_p
     795           0 :                       j  = jri+ir-1
     796           0 :                       rhohelp(ir) = rat(j) * rh(j) * SIN(g*rat(j))
     797             :                   END DO
     798             :                   CALL intgr3(rhohelp,rat(jri),dx,&
     799           0 :                               &                        n_out_p,qfout)
     800             :                   !--->             have to remove the small r-correction from intgr3
     801             :                   !roa...correction.from.intgr3.......................
     802           0 :                   IF (rhohelp(1)*rhohelp(2).GT.zero) THEN
     803           0 :                       alpha3 = 1.0 + LOG(rhohelp(2)/rhohelp(1))/dx
     804           0 :                       IF (alpha3.GT.zero)&
     805           0 :                               &                 qfout = qfout - rat(jri)*rhohelp(1)/alpha3
     806             :                       ENDIF
     807             :                   !roa...end.correction...............................
     808             : 
     809             : 
     810             :                   END IF
     811             : 
     812      331343 :                   qfout = fpi_const * qfout / g
     813             :                   !
     814             :           END IF
     815             :           !
     816      331811 :           qf(k)    = (qfin + qfout)/cell%omtil
     817             :       ENDDO
     818             : !$OMP END PARALLEL DO
     819         156 :       end subroutine FormFactor_forAtomType
     820             : 
     821             :       END MODULE m_cdnovlp
     822             : 

Generated by: LCOV version 1.13