Octopus
orbitalset.F90
Go to the documentation of this file.
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
22 use accel_oct_m
23 use batch_oct_m
25 use blas_oct_m
26 use comm_oct_m
27 use debug_oct_m
29 use global_oct_m
31 use ions_oct_m
34 use math_oct_m
35 use mesh_oct_m
42 use types_oct_m
44
45 implicit none
46
47 private
48
49 public :: &
64
65 type orbitalset_t
66 ! Components are public by default
67 integer :: nn, ll, ii
68 real(real64) :: jj
69 integer :: norbs
70 integer :: ndim
71 integer :: iatom
72 type(submesh_t) :: sphere
73 complex(real64), allocatable :: phase(:,:)
74 ! !< if the sphere cross the border of the box
75 real(real64) :: Ueff
76 real(real64) :: Ubar, Jbar
77 real(real64) :: alpha
78 integer :: nneighbors
79 ! !< interaction is considered
80 real(real64), allocatable :: V_ij(:,:)
81 real(real64), allocatable :: coulomb_IIJJ(:,:,:,:,:)
82 complex(real64), allocatable :: zcoulomb_IIJJ(:,:,:,:,:,:,:)
83 integer, allocatable:: map_os(:)
84 complex(real64), allocatable :: phase_shift(:,:)
85
86 real(real64) :: radius
87 class(species_t), pointer :: spec => null()
88 integer :: spec_index
89
90 real(real64), allocatable :: dorb(:,:,:)
91 complex(real64), allocatable :: zorb(:,:,:)
92 complex(real64), allocatable :: eorb_submesh(:,:,:,:)
93 ! !! (for isolated system with TD phase)
94 complex(real64), allocatable :: eorb_mesh(:,:,:,:)
95 ! !! (for periodic systems GS and TD)
96 integer :: ldorbs
97 type(accel_mem_t) :: dbuff_orb, zbuff_orb
98 type(accel_mem_t), allocatable :: buff_eorb (:)
99
100 logical :: use_submesh
101 logical :: allocated_on_mesh
102 ! This stores how the localized orbitals (without phase) were allocated.
103 ! This is different from "use_submesh", which dictate how the orbitals (with phase) will be applied.
104
105 type(poisson_t) :: poisson
106 contains
107 procedure :: set_species_index => orbitalset_set_species_index
108 end type orbitalset_t
109
110contains
111
112 subroutine orbitalset_init(this)
113 type(orbitalset_t), intent(inout) :: this
114
115 push_sub(orbitalset_init)
117 this%iatom = -1
118 this%nneighbors = 0
119 this%nn = 0
120 this%ll = 0
121 this%jj = m_one
122 this%ii = 0
123 this%iatom = 0
124 this%ndim = 1
125 this%spec_index = 0
126 this%norbs = 0
127
128 this%Ueff = m_zero
129 this%Ubar = m_zero
130 this%Jbar = m_zero
131 this%alpha = m_zero
132 this%radius = m_zero
133
134 pop_sub(orbitalset_init)
135 end subroutine orbitalset_init
136
137
138 subroutine orbitalset_end(this)
139 type(orbitalset_t), intent(inout) :: this
140
141 integer :: ik
142
143 push_sub(orbitalset_end)
144
145 safe_deallocate_a(this%phase)
146 safe_deallocate_a(this%dorb)
147 safe_deallocate_a(this%zorb)
148 safe_deallocate_a(this%eorb_submesh)
149 safe_deallocate_a(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 if(accel_is_enabled()) then
159 call accel_release_buffer(this%dbuff_orb)
160 call accel_release_buffer(this%zbuff_orb)
161 if(allocated(this%buff_eorb)) then
162 do ik = lbound(this%buff_eorb, dim=1), ubound(this%buff_eorb, dim=1)
163 call accel_release_buffer(this%buff_eorb(ik))
164 end do
165 safe_deallocate_a(this%buff_eorb)
166 end if
167 end if
169 pop_sub(orbitalset_end)
170 end subroutine orbitalset_end
172 subroutine orbitalset_set_jln(this, jj, ll, nn)
173 type(orbitalset_t), intent(inout) :: this
174 real(real64), intent(in) :: jj
175 integer, intent(in) :: ll, nn
179 this%jj = jj
180 this%ll = ll
181 this%nn = nn
184 end subroutine orbitalset_set_jln
188 subroutine orbitalset_update_phase(os, dim, kpt, kpoints, spin_polarized, vec_pot, vec_pot_var, kpt_max)
189 type(orbitalset_t), intent(inout) :: os
190 integer, intent(in) :: dim
191 type(distributed_t), intent(in) :: kpt
192 type(kpoints_t), intent(in) :: kpoints
193 logical, intent(in) :: spin_polarized
194 real(real64), optional, allocatable, intent(in) :: vec_pot(:)
195 real(real64), optional, allocatable, intent(in) :: vec_pot_var(:, :)
196 integer, optional, intent(in) :: kpt_max
197
198 integer :: ns, iq, is, ikpoint, im, idim, kpt_end
199 real(real64) :: kr, kpoint(1:dim), dx(1:dim)
200 integer :: iorb
201
203
204 ns = os%sphere%np
205
206 kpt_end = kpt%end
207 if (present(kpt_max)) kpt_end = min(kpt_max, kpt_end)
208
209 do iq = kpt%start, kpt_end
210 !This is durty but avoids to refer to states_get_kpoint_index
211 if (spin_polarized) then
212 ikpoint = 1 + (iq - 1)/2
213 else
214 ikpoint = iq
215 end if
216
217 ! if this fails, it probably means that sb is not compatible with std
218 assert(ikpoint <= kpoints_number(kpoints))
219
220 kpoint(1:dim) = kpoints%get_point(ikpoint)
221
222 !$omp parallel do private(dx, kr)
223 do is = 1, ns
224 ! this is only the correction to the global phase, that can
225 ! appear if the sphere crossed the boundary of the cell.
226 dx = os%sphere%rel_x(1:dim, is) - os%sphere%mesh%x(os%sphere%map(is), 1:dim) + os%sphere%center(1:dim)
227 kr = dot_product(kpoint, dx)
228 if (present(vec_pot)) then
229 if (allocated(vec_pot)) kr = kr + dot_product(vec_pot, dx)
230 end if
231
232 if (present(vec_pot_var)) then
233 if (allocated(vec_pot_var)) kr = kr + sum(vec_pot_var(1:dim, os%sphere%map(is)) &
234 *(os%sphere%rel_x(1:dim, is)+os%sphere%center))
235 end if
236
237 os%phase(is, iq) = exp(m_zi*kr)
238 end do
239
240 if (.not. os%use_submesh) then
241 ! We now compute the so-called Bloch sum of the localized orbitals
242 do idim = 1, os%ndim
243 do im = 1, os%norbs
244 os%eorb_mesh(:, im, idim, iq) = m_z0
245 do is = 1, ns
246 os%eorb_mesh(os%sphere%map(is),im,idim,iq) = os%eorb_mesh(os%sphere%map(is),im,idim,iq) &
247 + os%zorb(is, idim, im) * os%phase(is, iq)
248 end do
249 end do
250 end do
251 else ! In the case of the isolated system, we still use the submesh
252 !$omp parallel private(idim, is)
253 do im = 1, os%norbs
254 do idim = 1, os%ndim
255 !$omp do simd
256 do is = 1, ns
257 os%eorb_submesh(is,idim,im,iq) = os%zorb(is,idim,im)*os%phase(is, iq)
258 end do
259 !$omp end do simd
260 end do
261 end do
262 !$omp end parallel
263 end if
264
265 if(accel_is_enabled() .and. os%ndim == 1) then
266 if(os%use_submesh) then
267 do iorb = 1, os%norbs
268 call accel_write_buffer(os%buff_eorb(iq), ns, os%eorb_submesh(:, 1, iorb, iq), offset = (iorb - 1)*os%ldorbs)
269 end do
270 else
271 do iorb = 1, os%norbs
272 call accel_write_buffer(os%buff_eorb(iq), os%sphere%mesh%np, &
273 os%eorb_mesh(:, iorb, 1, iq), offset = (iorb - 1)*os%ldorbs)
274 end do
275 end if
276 end if
277 end do
278
280 end subroutine orbitalset_update_phase
281
282
284 subroutine orbitalset_update_phase_shift(os, dim, kpt, kpoints, spin_polarized, vec_pot, vec_pot_var, kpt_max)
285 type(orbitalset_t), intent(inout) :: os
286 integer, intent(in) :: dim
287 type(distributed_t), intent(in) :: kpt
288 type(kpoints_t), intent(in) :: kpoints
289 logical, intent(in) :: spin_polarized
290 real(real64), optional, allocatable, intent(in) :: vec_pot(:)
291 real(real64), optional, allocatable, intent(in) :: vec_pot_var(:, :)
292 integer, optional, intent(in) :: kpt_max
293
294 integer :: iq, ikpoint
295 real(real64) :: kr, kpoint(dim), dx(dim)
296 integer :: inn, kpt_end
297
299
300 kpt_end = kpt%end
301 if(present(kpt_max)) kpt_end = min(kpt_max, kpt_end)
302
303 do iq = kpt%start, kpt_end
304 !This is durty but avoids to refer to states_get_kpoint_index
305 if(spin_polarized) then
306 ikpoint = 1 + (iq - 1)/2
307 else
308 ikpoint = iq
309 end if
310
311 ! if this fails, it probably means that sb is not compatible with std
312 assert(ikpoint <= kpoints_number(kpoints))
313
314 kpoint(1:dim) = kpoints%get_point(ikpoint)
315
316 if (os%nneighbors > 0) then
317 do inn = 1, os%nneighbors
318 dx(1:dim) = os%V_ij(inn,1:dim)
319 kr = sum(kpoint(1:dim)*dx(1:dim))
320 if (present(vec_pot)) then
321 if (allocated(vec_pot)) kr = kr + sum(vec_pot(1:dim)*dx(1:dim))
322 end if
323
324 !At the moment the uniform vector potential is in vec_pot_var
325 if (present(vec_pot_var)) then
326 if (allocated(vec_pot_var)) kr = kr + sum(vec_pot_var(1:dim, 1)*dx(1:dim))
327 end if
328
329 !The sign is different as this is applied on the wavefunction and not the orbitals
330 os%phase_shift(inn, iq) = exp(-m_zi*kr)
331 end do
332 end if
333 end do
334
335
337 end subroutine orbitalset_update_phase_shift
338
339 ! ---------------------------------------------------------
345 subroutine orbitalset_set_species_index(os, ions)
346 class(orbitalset_t), intent(inout) :: os
347 type(ions_t), intent(in) :: ions
348
349 integer :: ja
350
352
353 os%spec_index = os%spec%get_index()
354 do ja = 1, ions%natoms
355 if(ions%atom(ja)%species == os%spec) then
356 os%spec_index = min(os%spec_index, ions%atom(ja)%species%get_index())
357 end if
358 end do
359
361 end subroutine orbitalset_set_species_index
362
363#include "undef.F90"
364#include "real.F90"
365#include "orbitalset_inc.F90"
366
367#include "undef.F90"
368#include "complex.F90"
369#include "orbitalset_inc.F90"
370
371end module orbitalset_oct_m
double exp(double __x) __attribute__((__nothrow__
subroutine, public accel_release_buffer(this, async)
Definition: accel.F90:919
pure logical function, public accel_is_enabled()
Definition: accel.F90:419
This module implements batches of mesh functions.
Definition: batch.F90:135
This module implements common operations on batches of mesh functions.
Definition: batch_ops.F90:118
This module contains interfaces for BLAS routines You should not use these routines directly....
Definition: blas.F90:120
real(real64), parameter, public m_zero
Definition: global.F90:191
complex(real64), parameter, public m_z0
Definition: global.F90:201
complex(real64), parameter, public m_zi
Definition: global.F90:205
real(real64), parameter, public m_one
Definition: global.F90:192
integer pure function, public kpoints_number(this)
Definition: kpoints.F90:1101
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:117
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:120
subroutine, public zorbitalset_get_position_matrix_elem(os, ndim, psib, idir, dot)
subroutine, public orbitalset_update_phase_shift(os, dim, kpt, kpoints, spin_polarized, vec_pot, vec_pot_var, kpt_max)
Build the phase shift for the intersite interaction.
Definition: orbitalset.F90:380
subroutine orbitalset_set_species_index(os, ions)
Set the species index for a given orbital set.
Definition: orbitalset.F90:441
subroutine, public orbitalset_init(this)
Definition: orbitalset.F90:208
subroutine, public orbitalset_end(this)
Definition: orbitalset.F90:234
subroutine, public dorbitalset_add_to_batch(os, ndim, psib, weight)
Definition: orbitalset.F90:833
subroutine, public dorbitalset_get_position_matrix_elem(os, ndim, psib, idir, dot)
subroutine, public zorbitalset_add_to_batch(os, ndim, psib, weight)
subroutine, public orbitalset_update_phase(os, dim, kpt, kpoints, spin_polarized, vec_pot, vec_pot_var, kpt_max)
Build the phase correction to the global phase in case the orbital crosses the border of the simulato...
Definition: orbitalset.F90:284
subroutine, public dorbitalset_get_coefficients(os, ndim, psi, ik, has_phase, dot, reduce)
Definition: orbitalset.F90:526
subroutine, public dorbitalset_get_coeff_batch(os, ndim, psib, dot, reduce)
Definition: orbitalset.F90:588
subroutine, public orbitalset_set_jln(this, jj, ll, nn)
Definition: orbitalset.F90:268
subroutine, public zorbitalset_get_coeff_batch(os, ndim, psib, dot, reduce)
subroutine, public zorbitalset_get_coefficients(os, ndim, psi, ik, has_phase, dot, reduce)
subroutine, public submesh_end(this)
Definition: submesh.F90:738
Distribution of N instances over mpi_grpsize processes, for the local rank mpi_grprank....