1!! Copyright (C) 2016 N. Tancogne-Dejean 2!! 3!! This program is free software; you can redistribute it and/or modify 4!! it under the terms of the GNU General Public License as published by 5!! the Free Software Foundation; either version 2, or (at your option) 6!! any later version. 7!! 8!! This program is distributed in the hope that it will be useful, 9!! but WITHOUT ANY WARRANTY; without even the implied warranty of 10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11!! GNU General Public License for more details. 12!! 13!! You should have received a copy of the GNU General Public License 14!! along with this program; if not, write to the Free Software 15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 16!! 02110-1301, USA. 17!! 18 19#include "global.h" 20 21module orbitalset_oct_m 22 use batch_oct_m 23 use batch_ops_oct_m 24 use blas_oct_m 25 use comm_oct_m 26 use distributed_oct_m 27 use global_oct_m 28 use hardware_oct_m 29 use kpoints_oct_m 30 use lalg_basic_oct_m 31 use mesh_oct_m 32 use mesh_function_oct_m 33 use messages_oct_m 34 use periodic_copy_oct_m 35 use poisson_oct_m 36 use profiling_oct_m 37 use simul_box_oct_m 38 use species_oct_m 39 use submesh_oct_m 40 use wfs_elec_oct_m 41 42 implicit none 43 44 private 45 46 public :: & 47 orbitalset_t, & 48 orbitalset_nullify, & 49 orbitalset_init, & 50 orbitalset_end, & 51 orbitalset_update_phase, & 52 dorbitalset_get_coefficients, & 53 zorbitalset_get_coefficients, & 54 dorbitalset_get_coeff_batch, & 55 zorbitalset_get_coeff_batch, & 56 dorbitalset_add_to_batch, & 57 zorbitalset_add_to_batch, & 58 dorbitalset_add_to_psi, & 59 zorbitalset_add_to_psi, & 60 orbitalset_set_jln 61 62 type orbitalset_t 63 ! Components are public by default 64 integer :: nn, ll, ii 65 FLOAT :: jj 66 integer :: norbs 67 integer :: ndim 68 integer :: iatom 69 type(submesh_t) :: sphere !> The submesh of the orbital 70 CMPLX, pointer :: phase(:,:) !> Correction to the global phase 71 !> if the sphere cross the border of the box 72 FLOAT :: Ueff !> The effective U of the simplified rotational invariant form 73 FLOAT :: Ubar, Jbar 74 FLOAT :: alpha !> A potential used to constrained occupations, as defined in PRB 71, 035105 (2005) 75 integer :: nneighbors !> Number of neighbouring atoms on which the intersite 76 !> interaction is considered 77 FLOAT, allocatable :: V_ij(:,:) !> The list of intersite interaction parameters 78 FLOAT, allocatable :: coulomb_IIJJ(:,:,:,:,:) !> Coulomb integrales with neighboring atoms 79 integer, allocatable:: map_os(:) 80 CMPLX, allocatable :: phase_shift(:,:) 81 82 FLOAT :: radius 83 type(species_t), pointer :: spec 84 85 FLOAT, pointer :: dorb(:,:,:) !> The orbital, if real, on the submesh 86 CMPLX, pointer :: zorb(:,:,:) !> The orbital, if complex, on the submesh 87 CMPLX, pointer :: eorb_submesh(:,:,:,:) !> Orbitals with its phase factor, on the submesh (for isolated system with TD phase) 88 CMPLX, pointer :: eorb_mesh(:,:,:,:) !> Orbitals with its phase factor, on the mesh (for periodic systems GS and TD) 89 90 logical :: submeshforperiodic !> Do we use or not submeshes for the orbitals 91 92 type(poisson_t) :: poisson !> For computing the Coulomb integrals 93 end type orbitalset_t 94 95contains 96 97 subroutine orbitalset_nullify(this) 98 type(orbitalset_t), intent(inout) :: this 99 100 PUSH_SUB(orbitalset_nullify) 101 102 nullify(this%phase) 103 nullify(this%spec) 104 nullify(this%dorb) 105 nullify(this%zorb) 106 nullify(this%eorb_submesh) 107 nullify(this%eorb_mesh) 108 109 call submesh_null(this%sphere) 110 call orbitalset_init(this) 111 112 POP_SUB(orbitalset_nullify) 113 114 end subroutine orbitalset_nullify 115 116 subroutine orbitalset_init(this) 117 type(orbitalset_t), intent(inout) :: this 118 119 PUSH_SUB(orbitalset_init) 120 121 this%iatom = -1 122 this%nneighbors = 0 123 this%nn = 0 124 this%ll = 0 125 this%jj = M_ONE 126 this%ii = 0 127 this%iatom = 0 128 this%ndim = 1 129 130 this%Ueff = M_ZERO 131 this%Ubar = M_ZERO 132 this%Jbar = M_ZERO 133 this%alpha = M_ZERO 134 this%radius = M_ZERO 135 136 POP_SUB(orbitalset_init) 137 end subroutine orbitalset_init 138 139 140 subroutine orbitalset_end(this) 141 type(orbitalset_t), intent(inout) :: this 142 143 PUSH_SUB(orbitalset_end) 144 145 SAFE_DEALLOCATE_P(this%phase) 146 SAFE_DEALLOCATE_P(this%dorb) 147 SAFE_DEALLOCATE_P(this%zorb) 148 SAFE_DEALLOCATE_P(this%eorb_submesh) 149 SAFE_DEALLOCATE_P(this%eorb_mesh) 150 nullify(this%spec) 151 call submesh_end(this%sphere) 152 153 SAFE_DEALLOCATE_A(this%V_ij) 154 SAFE_DEALLOCATE_A(this%coulomb_IIJJ) 155 SAFE_DEALLOCATE_A(this%map_os) 156 SAFE_DEALLOCATE_A(this%phase_shift) 157 158 POP_SUB(orbitalset_end) 159 end subroutine orbitalset_end 160 161 subroutine orbitalset_set_jln(this, jj, ll, nn) 162 type(orbitalset_t), intent(inout) :: this 163 FLOAT, intent(in) :: jj 164 integer, intent(in) :: ll, nn 165 166 PUSH_SUB(orbitalset_set_jln) 167 168 this%jj = jj 169 this%ll = ll 170 this%nn = nn 171 172 POP_SUB(orbitalset_set_jln) 173 end subroutine orbitalset_set_jln 174 175 176 !> Build the phase correction to the global phase in case the orbital crosses the border of the simulaton box 177 subroutine orbitalset_update_phase(os, sb, kpt, spin_polarized, vec_pot, vec_pot_var, kpt_max) 178 type(orbitalset_t), intent(inout) :: os 179 type(simul_box_t), intent(in) :: sb 180 type(distributed_t), intent(in) :: kpt 181 logical, intent(in) :: spin_polarized 182 FLOAT, optional, allocatable, intent(in) :: vec_pot(:) !< (sb%dim) 183 FLOAT, optional, allocatable, intent(in) :: vec_pot_var(:, :) !< (1:sb%dim, 1:ns) 184 integer, optional, intent(in) :: kpt_max 185 186 integer :: ns, iq, is, ikpoint, im, idim 187 FLOAT :: kr, kpoint(1:MAX_DIM), dx(1:MAX_DIM) 188 integer :: ndim, inn, kpt_end 189 190 PUSH_SUB(orbitalset_update_phase) 191 192 ns = os%sphere%np 193 ndim = sb%dim 194 195 kpt_end = kpt%end 196 if(present(kpt_max)) kpt_end = min(kpt_max, kpt_end) 197 198 do iq = kpt%start, kpt_end 199 !This is durty but avoids to refer to states_get_kpoint_index 200 if(spin_polarized) then 201 ikpoint = 1 + (iq - 1)/2 202 else 203 ikpoint = iq 204 end if 205 206 ! if this fails, it probably means that sb is not compatible with std 207 ASSERT(ikpoint <= kpoints_number(sb%kpoints)) 208 209 kpoint = M_ZERO 210 kpoint(1:ndim) = kpoints_get_point(sb%kpoints, ikpoint) 211 212 do is = 1, ns 213 ! this is only the correction to the global phase, that can 214 ! appear if the sphere crossed the boundary of the cell. 215 dx(1:ndim) = os%sphere%x(is, 1:ndim) - os%sphere%mesh%x(os%sphere%map(is), 1:ndim) + os%sphere%center(1:ndim) 216 kr = sum(kpoint(1:ndim)*dx(1:ndim)) 217 if(present(vec_pot)) then 218 if(allocated(vec_pot)) kr = kr + sum(vec_pot(1:ndim)*dx(1:ndim)) 219 end if 220 221 if(present(vec_pot_var)) then 222 if(allocated(vec_pot_var)) kr = kr + sum(vec_pot_var(1:ndim, os%sphere%map(is))*os%sphere%x(is, 1:ndim)) 223 end if 224 225 os%phase(is, iq) = exp(M_zI*kr) 226 end do 227 228 if(simul_box_is_periodic(sb) .and. .not. os%submeshforperiodic) then 229 !We now compute the so-called Bloch sum of the localized orbitals 230 os%eorb_mesh(:,:,:,iq) = M_Z0 231 do idim = 1, os%ndim 232 do im = 1, os%norbs 233 do is = 1, ns 234 os%eorb_mesh(os%sphere%map(is),im,idim,iq) = os%eorb_mesh(os%sphere%map(is),im,idim,iq) & 235 + os%zorb(is,idim,im)*os%phase(is, iq) 236 end do 237 end do 238 end do 239 else !In the case of the isolated system, we still use the submesh 240 do im = 1, os%norbs 241 do idim = 1, os%ndim 242 do is = 1, ns 243 os%eorb_submesh(is,idim,im,iq) = os%zorb(is,idim,im)*os%phase(is, iq) 244 end do 245 end do 246 end do 247 endif 248 249 250 if(os%nneighbors > 0) then 251 do inn = 1, os%nneighbors 252 dx(1:ndim) = os%V_ij(inn,1:ndim) 253 kr = sum(kpoint(1:ndim)*dx(1:ndim)) 254 if(present(vec_pot)) then 255 if(allocated(vec_pot)) kr = kr + sum(vec_pot(1:ndim)*dx(1:ndim)) 256 end if 257 258 !At the moment the uniform vector potential is in vec_pot_var 259 if(present(vec_pot_var)) then 260 if(allocated(vec_pot_var)) kr = kr + sum(vec_pot_var(1:ndim, 1)*dx(1:ndim)) 261 end if 262 263 !The sign is different as this is applied on the wavefunction and not the orbitals 264 os%phase_shift(inn, iq) = exp(-M_zI*kr) 265 end do 266 end if 267 end do 268 269 270 POP_SUB(orbitalset_update_phase) 271 end subroutine orbitalset_update_phase 272 273#include "undef.F90" 274#include "real.F90" 275#include "orbitalset_inc.F90" 276 277#include "undef.F90" 278#include "complex.F90" 279#include "orbitalset_inc.F90" 280 281end module orbitalset_oct_m 282