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      module m_compute_rhog
9      private
10      public :: compute_rhog
11      CONTAINS
12
13      subroutine compute_rhog(DM)
14
15!      USE siesta_options
16      use sparse_matrices, only: listh, listhptr, numh, maxnh
17      use sparse_matrices, only: H
18      use siesta_geom
19      use siesta_options, only: g2cut
20      use atomlist, only: no_u, iaorb, iphkb, qtot, indxuo, datm,
21     .                    lastkb, no_s, rmaxv, indxua, iphorb, lasto,
22     .                    rmaxo, no_l
23
24      use m_dhscf,      only: dhscf
25      use m_stress
26      use m_energies
27      use m_ntm
28      use m_spin,         only: nspin, h_spin_dim
29      use m_dipol
30      use alloc, only: re_alloc, de_alloc
31      use files, only : filesOut_t    ! derived type for output file names
32
33      implicit none
34
35      real(dp), intent(in):: DM(maxnh,h_spin_dim)
36
37      real(dp)            :: stressl(3,3)
38      real(dp), pointer   :: fal(:,:)   ! Local-node part of atomic F
39      integer             :: io, is, ispin
40      integer             :: ifa     ! Calc. forces?      0=>no, 1=>yes
41      integer             :: istr    ! Calc. stress?      0=>no, 1=>yes
42      integer             :: ihmat   ! Calc. hamiltonian? 0=>no, 1=>yes
43      real(dp)            :: g2max
44      type(filesOut_t)    :: filesOut  ! blank output file names
45
46!------------------------------------------------------------------------- BEGIN
47
48      call timer('compute_rhog',1)
49
50! Add SCF contribution to energy and matrix elements ..................
51      g2max = g2cut
52
53      nullify(fal)
54      call re_alloc( fal, 1, 3, 1, na_u, 'fal', 'setup_hamiltonian' )
55
56      ifa  = 0
57      istr = 0
58      ihmat = 0
59
60      ! Will compute rhog if mix_charge is .true., and
61      ! then return
62      call dhscf( nspin, no_s, iaorb, iphorb, no_l,
63     .            no_u, na_u, na_s, isa, xa, indxua,
64     .            ntm, ifa, istr, ihmat, filesOut,
65     .            maxnh, numh, listhptr, listh, DM, Datm,
66     .            maxnh, H, Enaatm, Enascf, Uatm, Uscf, DUscf, DUext,
67     .            Exc, Dxc, dipol, stress, fal, stressl,
68     .            charge_density_only=.true.)
69
70      call de_alloc( fal, 'fal', 'setup_hamiltonian' )
71
72      call timer('compute_rhog',2)
73
74
75!------------------------------------------------------------------------- END
76      END subroutine compute_rhog
77      END module m_compute_rhog
78