1!
2! Copyright (C) 1996-2016	The SIESTA group
3!  This file is distributed under the terms of the
4!  GNU General Public License: see COPYING in the top directory
5!  or http://www.gnu.org/copyleft/gpl.txt.
6! See Docs/Contributors.txt for a list of contributors.
7!
8      subroutine madelung(cell, shape, charnet, Emad)
9
10c *******************************************************************
11c Finds Madelung correction for charged systems
12c (Makov & Payne, PRB 51, 4014 (1995))
13c (Values of Madelung constant:
14C  Coldwell-Horsfall & Maradudin, J. Math. Phys. 1, 395 (1960))
15c
16c Written by P. Ordejon, March 1999.
17c ********* INPUT ***************************************************
18c real*8       cell(3,3): Lattice vectors (in Bohr)
19c character*10 shape    : Cell shape (atom,molecule,chain,slab,bulk)
20c real*8       charnet  : Net charge of the system
21c ********* OUTPUT **************************************************
22c real*8       Emad     : Madelung correction energy (in Ry)
23c *******************************************************************
24C
25C  Modules
26C
27      use precision
28      use parallel,   only : Node
29      use sys, only: die
30
31      implicit none
32
33      real(dp)
34     .  cell(3,3), charnet, Emad
35
36      character
37     .  shape*10
38
39C Local variables.........
40
41      real(dp)
42     .  alpha, lv
43
44      character
45     .  ctype*4
46
47C If system is charged, find if energy correction terms can be applied
48C (cluster or molecule, with SC, FCC or BCC cell) .....................
49      if (shape .ne. 'molecule' .and. shape .ne. 'atom') then
50        if (Node.eq.0) then
51          write(6,'(2a)')
52     .      'madelung: WARNING: Charged system, but not'
53     .      ,' an atom or molecule.'
54          write(6,'(2a)')
55     .      'madelung: WARNING: Energy correction terms'
56     .      ,' can not be applied. '
57          write(6,'(2a)')
58     .      'madelung: WARNING: Continue only if you really know'
59     .      ,' what you are doing.'
60        endif
61        Emad = 0.0_dp
62        return
63      else
64        call typecell(cell,ctype,lv)
65        if (ctype .eq. 'none') then
66          if (Node.eq.0) then
67            write(6,'(2a)')
68     .        'madelung: WARNING: Charged system, but not'
69     .        ,' SC, FCC or BCC cell.'
70            write(6,'(2a)')
71     .        'madelung: WARNING: Energy correction terms'
72     .        ,' can not be applied. '
73            write(6,'(2a)')
74     .        'madelung: WARNING: Continue only if you really know'
75     .        ,' what you are doing.'
76          endif
77          Emad = 0.0_dp
78          return
79        else if (ctype .eq. 'sc') then
80          alpha = 5.6745947_dp
81          if (Node.eq.0) then
82            write(6,'(a)')
83     .        'madelung: Charged system with SC cell'
84          endif
85        else if (ctype .eq. 'bcc') then
86          alpha = 7.278473_dp
87          if (Node.eq.0) then
88            write(6,'(a)')
89     .        'madelung: Charged system with BCC cell'
90          endif
91        else if (ctype .eq. 'fcc') then
92          alpha = 9.1697536_dp
93          if (Node.eq.0) then
94            write(6,'(a)')
95     .        'madelung: Charged system with FCC cell'
96          endif
97        else
98          call die('madelung: ERROR: Wrong type of cell')
99        endif
100      endif
101
102      Emad = charnet**2 * alpha / (2.0_dp * lv)
103C ...................
104      return
105      end
106
107