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_final_H_f_stress 9 private 10 public :: final_H_f_stress 11 CONTAINS 12 13 subroutine final_H_f_stress( istep, iscf , SCFconverged ) 14 use units, only: eV 15 16 use files, only : slabel 17#ifdef NCDF_4 18 use siesta_options, only: write_cdf 19#endif 20 use siesta_options, only: fixspin 21 use siesta_options, only: g2cut, savehs, temp, nmove 22 use siesta_options, only: recompute_H_after_scf 23 use sparse_matrices, only: numh, listh, listhptr 24 use sparse_matrices, only: H, S, Dscf, Escf, maxnh, xijo 25 use sparse_matrices, only: H_dftu_2D 26 use class_dSpData1D, only: val 27 use class_dSpData2D, only: val 28 29 use siesta_geom 30 use atomlist, only: no_u, iaorb, iphkb, qtot, indxuo, datm, 31 & lastkb, no_s, rmaxv, indxua, iphorb, lasto, 32 & rmaxo, no_l 33 use metaforce, only: lMetaForce, meta 34 use molecularmechanics, only : twobody 35 use m_nlefsm, only: nlefsm 36 use m_overfsm, only: overfsm 37 use m_kinefsm, only: kinefsm 38 use m_naefs, only: naefs 39 use m_dnaefs, only: dnaefs 40 use m_grdsam, only: grdsam 41 use dftu_specs, only: switch_dftu ! This variable determines whether 42 ! the subroutine to compute the 43 ! Hubbard terms should be called 44 ! or not 45 use m_dftu, only: hubbard_term ! Subroutine that compute the 46 ! Hubbard terms 47 use m_stress 48 use m_energies 49 use m_ntm 50 use m_spin, only: spin 51 use spinorbit, only: spinorb 52 use m_dipol 53 use m_forces, only: fa 54 use alloc, only: re_alloc, de_alloc 55 use m_hsx, only: write_hsx 56 use sys, only: die 57 use fdf 58#ifdef MPI 59 use m_mpi_utils, only: globalize_sum 60#endif 61 use parallel, only: IOnode 62#ifdef CDF 63#ifdef NCDF_4 64 use dictionary 65 use m_ncdf_siesta 66#endif 67#endif 68 use siesta_options, only : idyn, ia1, write_tshs_history 69 use sparse_matrices, only: H_2D, S_1D 70 use files, only : label_length 71 use m_ts_options, only : TS_HS_save 72 use m_ts_kpoints, only : ts_Gamma, ts_kscell, ts_kdispl 73 use m_ts_io, only : ts_write_TSHS,fname_TSHS, FC_index 74 75#ifdef FINAL_CHECK_HS 76 use m_compute_max_diff, only: compute_max_diff 77#endif 78 79 implicit none 80 81 ! MD-step, SCF-step 82 integer, intent(in) :: istep, iscf 83 logical, intent(in) :: SCFconverged 84 real(dp) :: stressl(3,3) 85 real(dp), pointer :: fal(:,:) ! Local-node part of atomic F 86 integer :: ifa ! Calculate forces? 0=>no, 1=>yes 87 integer :: istr ! Calculate stress? 0=>no, 1=>yes 88 real(dp) :: g2max 89 ! To avoid overwriting the current Hamiltonian and S 90 real(dp), pointer :: H_tmp(:,:) => null() 91 real(dp), pointer :: S_tmp(:) => null() 92 real(dp), pointer :: H_dftu(:,:) 93#ifdef FINAL_CHECK_HS 94 real(dp) :: diff_H 95#endif 96#ifdef MPI 97 real(dp) :: stresstmp(3,3) 98 real(dp), pointer :: fatmp(:,:) 99 real(dp) :: buffer1 100#endif 101 character(len=label_length+13) :: fname 102 integer :: io_istep, io_ia1 103 logical :: not_using_auxcell 104#ifdef CDF 105#ifdef NCDF_4 106 type(dictionary_t) :: dic_save 107#endif 108#endif 109!------------------------------------------------------------------------- BEGIN 110 111 not_using_auxcell = (no_s == no_u) 112 113! Initialize Hamiltonian ........................................ 114 call re_alloc(H_tmp, 1,maxnh, 1, spin%H, 115 $ 'H_tmp','final_H_f_stress') 116 117 H_tmp(:,:) = 0.0_dp 118 119! Initialize forces and stress ................... 120 nullify(fal) 121 call re_alloc( fal, 1, 3, 1, na_u, 'fal', 'final_H_f_stress' ) 122 123 fa(1:3,1:na_u) = 0.0_dp 124 fal(1:3,1:na_u) = 0.0_dp 125 stress(1:3,1:3) = 0.0_dp 126 stressl(1:3,1:3) = 0.0_dp 127 128 129! Neutral-atom: energy, forces and stress ............................ 130 call naefs( na_u, na_s, scell, xa, indxua, rmaxv, 131 & isa, Ena, fa, stress, forces_and_stress=.true. ) 132 call dnaefs( na_u, na_s, scell, xa, indxua, rmaxv, 133 & isa, DEna, fa, stress, forces_and_stress=.true. ) 134 Ena = Ena + DEna 135 136! Kinetic: energy, forces, stress and matrix elements ................. 137 call kinefsm( na_u, na_s, no_s, scell, xa, indxua, rmaxo, 138 & maxnh, maxnh, lasto, iphorb, isa, 139 & numh, listhptr, listh, numh, listhptr, listh, 140 & spin%spinor, Dscf, Ekin, fal, stressl, H_tmp, 141 & matrix_elements_only=.false. ) 142#ifdef MPI 143! Global reduction of energy terms 144 call globalize_sum(Ekin,buffer1) 145 Ekin = buffer1 146#endif 147! .................. 148 149! Non-local-pseudop: energy, forces, stress and matrix elements ....... 150 call nlefsm( scell, na_u, na_s, isa, xa, indxua, 151 & maxnh, maxnh, lasto, lastkb, iphorb, iphKB, 152 & numh, listhptr, listh, numh, listhptr, listh, 153 & spin%spinor, Dscf, Enl, fal, stressl, H_tmp, 154 & matrix_elements_only=.false. ) 155 156#ifdef MPI 157! Global reduction of energy terms 158 call globalize_sum(Enl,buffer1) 159 Enl = buffer1 160#endif 161 162! If in the future the spin-orbit routine is able to compute 163! forces and stresses, then "last" will be needed. If we are not 164! computing forces and stresses, calling it in the first iteration 165! should be enough 166! 167 if ( spin%SO ) then 168 call spinorb(no_u,no_l,iaorb,iphorb,isa,indxuo, 169 & maxnh,numh,listhptr,listh,Dscf,H_tmp(:,3:),Eso) 170 else 171 Eso = 0._dp 172 endif 173 174! Non-SCF part of total energy 175 call update_E0() 176 177! Hubbard term for LDA+U: energy, forces, stress and matrix elements .... 178 if( switch_dftu ) then 179 H_dftu => val(H_dftu_2D) 180 call hubbard_term(scell, na_u, na_s, isa, xa, indxua, 181 . maxnh, maxnh, lasto, iphorb, no_u, no_l, 182 . numh, listhptr, listh, numh, listhptr, listh, 183 . spin%spinor, Dscf, Edftu, DEdftu, H_dftu, 184 . fal, stressl, H_tmp, iscf, 185 . matrix_elements_only=.false.) 186#ifdef MPI 187 ! Global reduction of energy terms 188 call globalize_sum(Edftu,buffer1) 189 Edftu = buffer1 190 ! DEdftu should not be globalized as it 191 ! is based on globalized occupations 192#endif 193 Edftu = Edftu + DEdftu 194 endif 195! .................. 196 197 198! Non-local-pseudop: energy, forces, stress and matrix elements 199! Add SCF contribution to energy and matrix elements 200 g2max = g2cut 201! Last call to dhscf and grid-cell sampling if requested 202 ifa = 1 203 istr = 1 204 205 ! This will call dhscf with the final DM coming out of 206 ! the scf cycle 207 208 call grdsam( spin%Grid, no_s, iaorb, iphorb, 209 & no_l, no_u, na_u, na_s, isa, xa, indxua, 210 & ucell, ntm, ifa, istr, maxnh, 211 & maxnh, numh, listhptr, listh, Dscf, Datm, H_tmp, 212 & Enaatm, Enascf, Uatm, Uscf, DUscf, DUext, 213 & Exc, Dxc, dipol, stress, fal, stressl ) 214 215 ! Orthonormalization forces 216 call overfsm( na_u, na_s, no_l, no_s, scell, xa, indxua, rmaxo, 217 & lasto, iphorb, isa, maxnh, numh, listhptr, listh, 218 & spin, Escf, fal, stressl ) 219 220! Metadynamics forces 221 if (lMetaForce) then 222 call meta(xa,na_u,ucell,Emeta,fa,stress,.true.,.true.) 223 endif 224 225! Add on force field contribution if required 226 call twobody( na_u, xa, isa, ucell, Emm, 227 & ifa=1, fa=fa, istr=1, stress=stress ) 228 229#ifdef MPI 230! Global reduction of forces and stresses 231 nullify(fatmp) 232 call re_alloc( fatmp, 1, 3, 1, na_u, 'fatmp', 'final_H_f_stress' ) 233 call globalize_sum(stressl(1:3,1:3),stresstmp(1:3,1:3)) 234 call globalize_sum( fal(1:3,1:na_u), fatmp(1:3,1:na_u) ) 235 stress(1:3,1:3) = stress(1:3,1:3) + stresstmp(1:3,1:3) 236 fa(1:3,1:na_u) = fa(1:3,1:na_u) + fatmp(1:3,1:na_u) 237 call de_alloc( fatmp, 'fatmp', 'final_H_f_stress' ) 238#else 239 stress(1:3,1:3) = stress(1:3,1:3) + stressl(1:3,1:3) 240 fa(1:3,1:na_u) = fa(1:3,1:na_u) + fal(1:3,1:na_u) 241#endif 242 243 ! If backward compatibility is requested, recover 244 ! the old behavior with respect to H 245 246#ifdef FINAL_CHECK_HS 247 ! We also print-out the max-difference between the 248 ! final H after SCF and the final computed H used for 249 ! the forces. 250 ! This value should correspond to dHmax 251 call compute_max_diff(H, H_tmp, diff_H) 252 if ( IONode ) then 253 ! Print out the final differences 254 write(*,"(a,f10.6)") ":!: Final H_scf - H_force (eV) : ", 255 & diff_H/eV 256 end if 257#endif 258 259 if (recompute_H_after_scf) then 260 if (ionode) then 261 write(6,"(a)") ":!: Updating H after scf cycle" // 262 $ " as requested by compat. option" 263 endif 264 H = H_tmp 265 endif 266 267 call de_alloc( fal, 'fal', 'final_H_f_stress' ) 268 call de_alloc( H_tmp, 'H_tmp', 'final_H_f_stress' ) 269 270 ! Determine whether anything should be saved 271 ! If not, return immediately 272 if ( .not. SCFconverged ) return 273 274! Save Hamiltonian and overlap matrices ............................ 275! Only in HSX format now. Use Util/HSX/hsx2hs to generate an HS file 276! 277! Note that we save the Hamiltonian coming out of the scf cycle, 278! not the one computed (from DM_out) in this routine. 279! This call could be moved to a more appropriate place 280! 281 if (savehs) then 282 call write_hsx( no_u, no_s, spin%H, indxuo, 283 & maxnh, numh, listhptr, listh, H, S, Qtot, 284 & Temp, xijo) 285 endif 286#ifdef CDF 287#ifdef NCDF_4 288 if ( write_cdf ) then 289 if ( idyn == 6 ) then 290 call cdf_save_fc(trim(slabel)//'.nc', istep) 291 if ( istep == 0 ) then 292 dic_save = ('fa'.kv.1)//('stress'.kv.1)//('Ef'.kv.1) 293 if ( savehs .or. TS_HS_save ) 294 & dic_save = dic_save//('H'.kv.1) 295 call cdf_save_state(trim(slabel)//'.nc',dic_save) 296 call delete(dic_save) 297 end if 298 else if ( idyn /=0 .or. nmove /= 0 ) then 299 !call cdf_save_md(trim(slabel)//'.nc') 300 else 301 dic_save = ('fa'.kv.1)//('stress'.kv.1)//('Ef'.kv.1) 302 if ( savehs .or. TS_HS_save ) dic_save = dic_save //('H'.kv.1) 303 call cdf_save_state(trim(slabel)//'.nc',dic_save) 304 call delete(dic_save) 305 end if 306 end if 307#endif 308#endif 309 310 if ( fixspin .and. (write_tshs_history .or. TS_HS_save) ) then 311 ! For fixed spin calculations we shift the Hamiltonian according 312 ! to the first spin such that we have a "common" E_F, after writing 313 ! we shift back again. 314 H_tmp => val(H_2D) 315 S_tmp => val(S_1D) 316 call daxpy(size(S_tmp),Efs(1)-Efs(2),S_tmp(1),1,H_tmp(1,2),1) 317 Ef = Efs(1) 318 end if 319 320 if ( write_tshs_history ) then 321 ! This is "pure" MD and we only write consecutive numbers 322 ! Together with this you cannot also save FC 323 fname = fname_TSHS(slabel, istep = istep ) 324 call ts_write_tshs(fname, .false., not_using_auxcell, ts_Gamma, 325 & ucell, nsc, isc_off, na_u, no_s, spin%H, 326 & ts_kscell, ts_kdispl, 327 & xa, lasto, 328 & H_2D, S_1D, indxuo, 329 & Ef, Qtot, Temp, istep, 0) 330 else if ( TS_HS_save ) then 331 if ( idyn == 6 ) then 332 ! Correct the FC indices such that 333 ! the TSHS format contains, atom displaced, index of 334 ! displacement: 335 ! 1 = -x 336 ! 2 = +x 337 ! 3 = -y 338 ! 4 = +y 339 ! 5 = -z 340 ! 6 = +z 341 call FC_index(istep,ia1,io_istep,io_ia1) 342 fname = fname_TSHS(slabel,istep = io_istep, ia1 = io_ia1) 343 call ts_write_tshs(fname, .false., 344 & not_using_auxcell, ts_Gamma, 345 & ucell, nsc, isc_off, na_u, no_s, spin%H, 346 & ts_kscell, ts_kdispl, 347 & xa, lasto, 348 & H_2D, S_1D, indxuo, 349 & Ef, Qtot, Temp, io_istep, io_ia1) 350 else 351 fname = fname_TSHS(slabel) 352 call ts_write_tshs(fname, .false., 353 & not_using_auxcell, ts_Gamma, 354 & ucell, nsc, isc_off, na_u, no_s, spin%H, 355 & ts_kscell, ts_kdispl, 356 & xa, lasto, 357 & H_2D, S_1D, indxuo, 358 & Ef, Qtot, Temp, 0, 0) 359 end if 360 endif 361 362 if ( fixspin .and. (write_tshs_history .or. TS_HS_save) ) then 363 ! Shift back 364 call daxpy(size(S_tmp),Efs(2)-Efs(1),S_tmp(1),1,H_tmp(1,2),1) 365 end if 366 367!----------------------------------------------------------------------- END 368 END subroutine final_H_f_stress 369 END module m_final_H_f_stress 370