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