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_setup_hamiltonian
9      private
10      public :: setup_hamiltonian
11      CONTAINS
12
13      subroutine setup_hamiltonian( iscf )
14
15      USE siesta_options
16      use sparse_matrices, only: H_kin_1D, H_vkb_1D
17      use sparse_matrices, only: H_dftu_2D, H_so_2D
18      use sparse_matrices, only: listh, listhptr, numh, maxnh
19      use sparse_matrices, only: H, S, Hold
20      use sparse_matrices, only: Dscf, Escf, xijo
21      use class_dSpData1D,  only: val
22      use class_dSpData2D,  only: val
23
24      use siesta_geom
25      use atmfuncs, only: uion
26      use atomlist, only: no_u, iaorb, iphkb, qtot, indxuo, datm,
27     .                    lastkb, no_s, rmaxv, indxua, iphorb, lasto,
28     .                    rmaxo, no_l
29      use metaforce, only: lMetaForce, meta
30      use molecularmechanics, only : twobody
31      use dftu_specs,   only: switch_dftu     ! This variable determines whether
32                                              !   the subroutine to compute the
33                                              !   Hubbard terms should be called
34                                              !   or not
35      use m_dftu,       only: hubbard_term    ! Subroutine that compute the
36                                              !   Hubbard terms
37      use m_dhscf,      only: dhscf
38      use m_stress
39      use m_energies
40      use parallel, only: Node
41      use m_steps, only: istp
42      use m_ntm
43      use m_spin,         only: spin
44      use m_dipol
45      use alloc, only: re_alloc, de_alloc
46      use m_hsx, only: write_hsx
47      use sys, only: die, bye
48      use m_partial_charges, only: want_partial_charges
49      use files, only : filesOut_t    ! derived type for output file names
50      use m_rhog, only: rhog_in, rhog
51#ifdef MPI
52      use m_mpi_utils, only: globalize_sum
53#endif
54
55      implicit none
56      integer, intent(in) :: iscf
57      real(dp)            :: stressl(3,3)
58      real(dp), pointer   :: fal(:,:)   ! Local-node part of atomic F
59#ifdef MPI
60      real(dp)            :: buffer1
61#endif
62      integer             :: io, is, ispin
63      integer             :: ifa     ! Calc. forces?      0=>no, 1=>yes
64      integer             :: istr    ! Calc. stress?      0=>no, 1=>yes
65      integer             :: ihmat   ! Calc. hamiltonian? 0=>no, 1=>yes
66      real(dp)            :: g2max
67      type(filesOut_t)    :: filesOut  ! blank output file names
68      logical             :: use_rhog_in
69
70      real(dp), pointer   :: H_vkb(:), H_kin(:), H_dftu(:,:), H_so(:,:)
71
72!------------------------------------------------------------------------- BEGIN
73
74      call timer('setup_H',1)
75
76      ! Nullify pointers
77      nullify(fal)
78
79!$OMP parallel default(shared), private(ispin,io)
80
81!     Save present H matrix
82      do ispin = 1, spin%H
83!$OMP do
84         do io = 1,maxnh
85            Hold(io,ispin) = H(io,ispin)
86         enddo
87!$OMP end do nowait
88      enddo
89
90!$OMP single
91      H_kin => val(H_kin_1D)
92      H_vkb => val(H_vkb_1D)
93      if ( spin%SO ) then
94         ! Sadly some compilers (g95), does
95         ! not allow bounds for pointer assignments :(
96         H_so => val(H_so_2D)
97      end if
98!$OMP end single ! keep wait
99
100      ! We do not need to set the non-spinor components
101      ! For non-colinear they are set down below,
102      ! while for spin-orbit they are set to the H_so initial
103      ! spin-orbit.
104
105      do ispin = 1, spin%spinor
106!$OMP do
107          do io = 1,maxnh
108            H(io,ispin) = H_kin(io) + H_vkb(io)
109          end do
110!$OMP end do nowait
111      end do
112
113      if ( spin%NCol ) then
114
115         do ispin = 3 , spin%H
116!$OMP do
117            do io = 1, maxnh
118               H(io,ispin) = 0._dp
119            end do
120!$OMP end do nowait
121         end do
122
123      else if ( spin%SO ) then
124
125         do ispin = 3 , spin%H
126!$OMP do
127            do io = 1, maxnh
128               H(io,ispin) = H_so(io,ispin-2)
129            end do
130!$OMP end do nowait
131         end do
132
133      end if
134
135! ..................
136
137! Non-SCF part of total energy .......................................
138! Note that these will be "impure" for a mixed Dscf
139
140! If mixing the charge, Dscf is the previous step's DM_out. Since
141! the "scf" components of the energy are computed with the (mixed)
142! charge, this introduces an inconsistency. In this case the energies
143! coming out of this routine need to be corrected.
144!
145!$OMP single
146      Ekin = 0.0_dp
147      Enl = 0.0_dp
148!$OMP end single ! keep wait
149
150!$OMP do reduction(+:Ekin,Enl)
151      do io = 1,maxnh
152        do ispin = 1, spin%spinor
153          Ekin = Ekin + H_kin(io) * Dscf(io,ispin)
154          Enl  = Enl  + H_vkb(io) * Dscf(io,ispin)
155        end do
156      end do
157!$OMP end do nowait
158
159!$OMP single
160      Eso = 0._dp
161!$OMP end single
162      if ( spin%SO ) then
163!$OMP do reduction(+:Eso)
164         do io = 1, maxnh
165            Eso = Eso + H_so(io,1)*Dscf(io,7) + H_so(io,2)*Dscf(io,8)
166     .           + H_so(io,5)*Dscf(io,3) + H_so(io,6)*Dscf(io,4)
167     .           - H_so(io,3)*Dscf(io,5) - H_so(io,4)*Dscf(io,6)
168         end do
169!$OMP end do nowait
170      end if
171
172!$OMP end parallel
173
174#ifdef MPI
175      ! Global reduction of Ekin, Enl
176      call globalize_sum(Ekin,buffer1)
177      Ekin = buffer1
178      call globalize_sum(Enl,buffer1)
179      Enl = buffer1
180      if ( spin%SO ) then
181         ! Global reduction of Eso
182         call globalize_sum(Eso,buffer1)
183         Eso = buffer1
184      end if
185#endif
186
187!     Non-SCF part of total energy
188      call update_E0()
189
190! Hubbard term for LDA+U: energy, forces, stress and matrix elements ....
191      if( switch_dftu ) then
192        if ( spin%NCol ) then
193         call die('LDA+U cannot be used with non-colinear spin.')
194        end if
195        if ( spin%SO ) then
196         call die('LDA+U cannot be used with spin-orbit coupling.')
197        end if
198
199        call re_alloc( fal, 1, 3, 1, na_u, 'fal', 'setup_hamiltonian' )
200
201        H_dftu => val(H_dftu_2D)
202        call hubbard_term(scell, na_u, na_s, isa, xa, indxua,
203     .                    maxnh, maxnh, lasto, iphorb, no_u, no_l,
204     .                    numh, listhptr, listh, numh, listhptr, listh,
205     .                    spin%spinor, Dscf, Edftu, DEdftu, H_dftu,
206     .                    fal, stressl, H, iscf,
207     .                    matrix_elements_only=.true.)
208
209#ifdef MPI
210        ! Global reduction of energy terms
211        call globalize_sum(Edftu,buffer1)
212        Edftu = buffer1
213        ! DEdftu should not be globalized
214        ! as it is based on globalized occupations
215#endif
216        Edftu = Edftu + DEdftu
217
218        call de_alloc( fal, 'fal', 'setup_hamiltonian' )
219      endif
220! ..................
221
222
223! Add SCF contribution to energy and matrix elements ..................
224      g2max = g2cut
225
226      call re_alloc( fal, 1, 3, 1, na_u, 'fal', 'setup_hamiltonian' )
227
228      ifa  = 0
229      istr = 0
230      ihmat = 1
231      if ((hirshpop .or. voropop)
232     $     .and. partial_charges_at_every_scf_step) then
233         want_partial_charges = .true.
234      endif
235      use_rhog_in =  (mix_charge .and. iscf > 1)
236
237      call dhscf( spin%Grid, no_s, iaorb, iphorb, no_l,
238     .            no_u, na_u, na_s, isa, xa, indxua,
239     .            ntm, ifa, istr, ihmat, filesOut,
240     .            maxnh, numh, listhptr, listh, Dscf, Datm,
241     .            maxnh, H, Enaatm, Enascf, Uatm, Uscf, DUscf, DUext,
242     .            Exc, Dxc, dipol, stress, fal, stressl,
243     .            use_rhog_in)
244
245      ! This statement will apply to iscf = 1, for example, when
246      ! we do not use rhog_in. Rhog here is always the charge used to
247      ! build H, that is, rhog_in.
248      if (mix_charge) rhog_in = rhog
249
250      want_partial_charges = .false.
251      call de_alloc( fal, 'fal', 'setup_hamiltonian' )
252!  It is wasteful to write over and over H and S, as there are
253!  no different files.
254! Save Hamiltonian and overlap matrices ............................
255! Only in HSX format now.  Use Util/HSX/hsx2hs to generate an HS file
256
257      if (savehs .or. write_coop) then
258         call write_hsx( no_u, no_s, spin%H, indxuo,
259     &        maxnh, numh, listhptr, listh, H, S, qtot,
260     &        temp, xijo)
261      endif
262
263      call timer('setup_H',2)
264#ifdef SIESTA__PEXSI
265      if (node==0) call memory_snapshot("after setup_H")
266#endif
267
268      if ( h_setup_only ) then
269         call timer( 'all', 2 )  ! New call to close the tree
270         call timer( 'all', 3 )
271         call bye("H-Setup-Only requested")
272         STOP
273      endif
274
275!------------------------------------------------------------------------- END
276      END subroutine setup_hamiltonian
277      END module m_setup_hamiltonian
278