Octopus
orbitalset_utils.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#include "global.h"
18
21 use comm_oct_m
22 use debug_oct_m
25 use global_oct_m
27 use ions_oct_m
28 use, intrinsic :: iso_fortran_env
31 use loct_oct_m
32 use math_oct_m
33 use mesh_oct_m
35 use mpi_oct_m
40 use space_oct_m
43 use unit_oct_m
44
45 implicit none
46
47 private
48
49 public :: &
56
57 integer, public, parameter :: &
58 SM_POISSON_DIRECT = 0, &
59 sm_poisson_isf = 1, &
62
63
64contains
65
66
71 integer function orbitalset_utils_count(species, iselect) result(norb)
72 class(species_t), intent(in) :: species
73 integer, optional, intent(in) :: iselect
74
75 integer :: iorb, ii, ll, mm
76
77 norb = 0
78 do iorb = 1, species%get_niwfs()
79 call species%get_iwf_ilm(iorb, 1, ii, ll, mm)
80 if (present(iselect)) then
81 if (ii == iselect) norb = norb + 1
82 else
83 norb = max(norb, ii)
84 end if
85 end do
86 end function orbitalset_utils_count
87
88 subroutine orbitalset_init_intersite(this, namespace, space, ind, ions, der, psolver, os, nos, maxnorbs, &
89 rcut, kpt, has_phase, sm_poisson, basis_from_states, combine_j_orbitals)
90 type(orbitalset_t), intent(inout) :: this
91 type(namespace_t), intent(in) :: namespace
92 class(space_t), intent(in) :: space
93 integer, intent(in) :: ind
94 type(ions_t), intent(in) :: ions
95 type(derivatives_t), intent(in) :: der
96 type(poisson_t), intent(in) :: psolver
97 type(orbitalset_t), intent(inout) :: os(:)
98 integer, intent(in) :: nos, maxnorbs
99 real(real64), intent(in) :: rcut
100 type(distributed_t), intent(in) :: kpt
101 logical, intent(in) :: has_phase
102 integer, intent(in) :: sm_poisson
103 logical, intent(in) :: basis_from_states
104 logical, intent(in) :: combine_j_orbitals
105
106
107 type(lattice_iterator_t) :: latt_iter
108 real(real64) :: xat(space%dim), xi(space%dim)
109 real(real64) :: rr
110 integer :: inn, ist, jst
111 integer :: ip, ios, ip2, is1, is2
112 type(submesh_t) :: sm
113 real(real64), allocatable :: tmp(:), vv(:), nn(:)
114 complex(real64), allocatable :: ztmp(:), zvv(:,:), znn(:)
115 real(real64), allocatable :: dorb(:,:,:,:)
116 complex(real64), allocatable :: zorb(:,:,:,:)
117 real(real64), parameter :: TOL_INTERSITE = 1.e-5_real64
118 type(distributed_t) :: dist
119 integer :: sm_poisson_
120
122
123 call messages_print_with_emphasis(msg="Intersite Coulomb integrals", namespace=namespace)
124
125 sm_poisson_ = sm_poisson
126
127 this%nneighbors = 0
128 if (this%iatom /= -1) then
129 xat = ions%pos(:, this%iatom)
130 else
131 xat = this%sphere%center
132 end if
133
134 latt_iter = lattice_iterator_t(ions%latt, rcut)
135
136 !We first count first the number of neighboring atoms at a distance max rcut
137 do ios = 1, nos
138 do inn = 1, latt_iter%n_cells
139
140 xi = os(ios)%sphere%center(1:space%dim) + latt_iter%get(inn)
141 rr = norm2(xi - xat)
142
143 !This atom is too far
144 if( rr >rcut + tol_intersite ) cycle
145 !Intra atomic interaction
146 if( ios == ind .and. rr < tol_intersite) cycle
147
148 this%nneighbors = this%nneighbors +1
149 end do
150 end do
151
152 !The first three values are the position of the periodic copies
153 !and the zero value one is used to store the actual value of V_ij
154 safe_allocate(this%V_ij(1:this%nneighbors, 0:space%dim+1))
155 this%V_ij(1:this%nneighbors, 0:space%dim+1) = m_zero
156 safe_allocate(this%map_os(1:this%nneighbors))
157 this%map_os(1:this%nneighbors) = 0
158 if(has_phase) then
159 safe_allocate(this%phase_shift(1:this%nneighbors, kpt%start:kpt%end))
160 end if
161
162 this%nneighbors = 0
163 do ios = 1, nos
164 do inn = 1, latt_iter%n_cells
165 xi = os(ios)%sphere%center(1:space%dim) + latt_iter%get(inn)
166 rr = norm2(xi - xat)
167
168 if( rr > rcut + tol_intersite ) cycle
169 if( ios == ind .and. rr < tol_intersite) cycle
170
171 this%nneighbors = this%nneighbors +1
172
173 this%V_ij(this%nneighbors, 1:space%dim) = xi(1:space%dim) -os(ios)%sphere%center(1:space%dim)
174 this%V_ij(this%nneighbors, space%dim+1) = rr
175
176 this%map_os(this%nneighbors) = ios
177 end do
178 end do
179
180 write(message(1),'(a, i3, a)') 'Intersite interaction will be computed for ', this%nneighbors, ' neighboring atoms.'
181 call messages_info(1, namespace=namespace)
182
183
184 if (this%ndim == 1) then
185 safe_allocate(this%coulomb_IIJJ(1:this%norbs,1:this%norbs,1:maxnorbs,1:maxnorbs,1:this%nneighbors))
186 this%coulomb_IIJJ = m_zero
187 else
188 safe_allocate(this%zcoulomb_IIJJ(1:this%norbs,1:this%norbs,1:maxnorbs,1:maxnorbs, 1:this%ndim, 1:this%ndim, 1:this%nneighbors))
189 this%zcoulomb_IIJJ = m_zero
190 end if
191
192 if(this%nneighbors == 0) then
193 call messages_print_with_emphasis(namespace=namespace)
195 return
196 end if
197
198 call distributed_nullify(dist, this%nneighbors)
199#ifdef HAVE_MPI
200 if(.not. der%mesh%parallel_in_domains) then
201 call distributed_init(dist, this%nneighbors, mpi_comm_world, 'orbs')
202 end if
203#endif
204
205 do inn = dist%start, dist%end
206
207 ios = this%map_os(inn)
208
209 if(.not. basis_from_states) then
210 !Init a submesh from the union of two submeshes
211 call submesh_merge(sm, space, der%mesh, this%sphere, os(ios)%sphere, &
212 shift = this%V_ij(inn, 1:space%dim))
213
214 write(message(1),'(a, i3, a, f6.3, a, i7, a)') 'Neighbor ', inn, ' is located at ', &
215 this%V_ij(inn, space%dim+1), ' Bohr and has ', sm%np, ' grid points.'
216 call messages_info(1, namespace=namespace)
217
218 if (this%ndim == 1) then
219 safe_allocate(dorb(1:sm%np, 1:1, 1:max(this%norbs, os(ios)%norbs), 1:2))
220 else
221 safe_allocate(zorb(1:sm%np, 1:2, 1:max(this%norbs, os(ios)%norbs), 1:2))
222 end if
223
224 ! Get the orbitals from the first set of orbitals
225 if (this%ndim == 1) then
226 do ist = 1, this%norbs
227 call dget_orbital(this%spec, sm, this%ii, this%ll, this%jj, ist, 1, dorb(:, :, ist, 1), combine_j_orbitals)
228 end do
229 else
230 do ist = 1, this%norbs
231 call zget_orbital(this%spec, sm, this%ii, this%ll, this%jj, ist, 2, zorb(:, :, ist, 1), combine_j_orbitals)
232 end do
233 end if
234
235 call submesh_shift_center(sm, space, this%V_ij(inn, 1:space%dim)+os(ios)%sphere%center(1:space%dim))
236
237 ! Get the orbitals from the second set of orbitals
238 if (this%ndim == 1) then
239 do ist = 1, os(ios)%norbs
240 call dget_orbital(os(ios)%spec, sm, os(ios)%ii, os(ios)%ll, os(ios)%jj, ist, 1, dorb(:, :, ist, 2), combine_j_orbitals)
241 end do
242 else
243 do ist = 1, os(ios)%norbs
244 call zget_orbital(os(ios)%spec, sm, os(ios)%ii, os(ios)%ll, os(ios)%jj, ist, 2, zorb(:, :, ist, 2), combine_j_orbitals)
245 end do
246 end if
247
248 else
249 ! TODO: Replace datomic_orbital_get_submesh_safe by the logic above
250 assert(this%ndim == 1)
251 !Init a submesh from the union of two submeshes
252 call submesh_merge(sm, ions%space, der%mesh, this%sphere, os(ios)%sphere, &
253 shift = this%V_ij(inn, 1:ions%space%dim))
254
255 write(message(1),'(a, i3, a, f6.3, a, i5, a)') 'Neighbor ', inn, ' is located at ', &
256 this%V_ij(inn, ions%space%dim+1), ' Bohr and has ', sm%np, ' grid points.'
257 call messages_info(1, namespace=namespace)
258
259 safe_allocate(dorb(1:sm%np, 1:1, 1:max(this%norbs,os(ios)%norbs), 1:2))
260 dorb = m_zero
261
262 ! All the points of the first submesh are included in the union of the submeshes
263 do ist = 1, this%norbs
264 if(allocated(this%dorb)) then
265 dorb(1:this%sphere%np, 1, ist, 1) = this%dorb(:,1,ist)
266 else
267 dorb(1:this%sphere%np, 1, ist, 1) = real(this%zorb(:,1,ist), real64)
268 end if
269 end do
270
271 call submesh_shift_center(sm, ions%space, this%V_ij(inn, 1:ions%space%dim)+os(ios)%sphere%center(1:ions%space%dim))
272
273 ! TODO: This probably needs some optimization
274 ! However, this is only done at initialization time
275 do ist = 1, os(ios)%norbs
276 if(allocated(this%dorb)) then
277 !$omp parallel do private(ip)
278 do ip2 = 1, sm%np
279 do ip = 1, os(ios)%sphere%np
280 if(all(abs(sm%rel_x(1:ions%space%dim, ip2)-os(ios)%sphere%rel_x(1:ions%space%dim, ip)) < 1e-6_real64)) then
281 dorb(ip2, 1, ist, 2) = os(ios)%dorb(ip, 1, ist)
282 end if
283 end do
284 end do
285 else
286 !$omp parallel do private(ip)
287 do ip2 = 1, sm%np
288 do ip = 1, os(ios)%sphere%np
289 if(all(abs(sm%rel_x(1:ions%space%dim, ip2)-os(ios)%sphere%rel_x(1:ions%space%dim, ip)) < 1e-6_real64)) then
290 dorb(ip2, 1, ist, 2) = real(os(ios)%zorb(ip, 1, ist), real64)
291 end if
292 end do
293 end do
294 end if
295 end do
296
297 end if
298
299 select case (sm_poisson_)
300 case(sm_poisson_direct)
301 !Build information needed for the direct Poisson solver on the submesh
302 call submesh_build_global(sm, space)
303 call poisson_init_sm(this%poisson, namespace, space, psolver, der, sm, method = poisson_direct_sum)
304 case(sm_poisson_isf)
305 call poisson_init_sm(this%poisson, namespace, space, psolver, der, sm, method = poisson_isf)
307 call poisson_init_sm(this%poisson, namespace, space, psolver, der, sm, method = poisson_psolver)
308 case(sm_poisson_fft)
309 call poisson_init_sm(this%poisson, namespace, space, psolver, der, sm, method = poisson_fft, &
310 force_cmplx=(this%ndim==2))
311 end select
312
313 if (this%ndim == 1) then ! Real orbitals
314 safe_allocate(tmp(1:sm%np))
315 safe_allocate(nn(1:sm%np))
316 safe_allocate(vv(1:sm%np))
317
318 do ist = 1, this%norbs
319 !$omp parallel do
320 do ip = 1, sm%np
321 nn(ip) = dorb(ip, 1, ist, 1)*dorb(ip, 1, ist, 1)
322 end do
323 !$omp end parallel do
324
325 ! Here it is important to use a non-periodic poisson solver, e.g. the direct solver,
326 ! to not include contribution from periodic copies.
327 call dpoisson_solve_sm(this%poisson, namespace, sm, vv, nn)
328
329 do jst = 1, os(ios)%norbs
330
331 !$omp parallel do
332 do ip = 1, sm%np
333 tmp(ip) = vv(ip)*dorb(ip, 1, jst, 2)*dorb(ip, 1, jst, 2)
334 end do
335 !$omp end parallel do
336
337 this%coulomb_IIJJ(ist, ist, jst, jst, inn) = dsm_integrate(der%mesh, sm, tmp, reduce = .false.)
338 end do !jst
339 end do !ist
340
341 safe_deallocate_a(nn)
342 safe_deallocate_a(vv)
343 safe_deallocate_a(tmp)
344
345 else ! Complex orbitals
346 safe_allocate(ztmp(1:sm%np))
347 safe_allocate(znn(1:sm%np))
348 safe_allocate(zvv(1:sm%np, 1:this%ndim))
349
350 do ist = 1, this%norbs
351 do is1 = 1, this%ndim
352 !$omp parallel do
353 do ip = 1, sm%np
354 znn(ip) = conjg(zorb(ip, is1, ist, 1))*zorb(ip, is1, ist, 1)
355 end do
356 !$omp end parallel do
357
358 ! Here it is important to use a non-periodic poisson solver, e.g. the direct solver,
359 ! to not include contribution from periodic copies.
360 call zpoisson_solve_sm(this%poisson, namespace, sm, zvv(:,is1), znn)
361 end do !is1
362
363 do jst = 1, os(ios)%norbs
364
365 do is1 = 1, this%ndim
366 do is2 = 1, this%ndim
367 !$omp parallel do
368 do ip = 1, sm%np
369 ztmp(ip) = zvv(ip, is1)*conjg(zorb(ip, is2, jst, 2))*zorb(ip, is2, jst, 2)
370 end do
371 !$omp end parallel do
372
373 this%zcoulomb_IIJJ(ist, ist, jst, jst, is1, is2, inn) = zsm_integrate(der%mesh, sm, ztmp)
374 end do
375 end do
376 end do !jst
377 end do !ist
378
379 safe_deallocate_a(znn)
380 safe_deallocate_a(zvv)
381 safe_deallocate_a(ztmp)
382 end if
383
384 call poisson_end(this%poisson)
385
386
387 if (sm_poisson_ == sm_poisson_direct) call submesh_end_global(sm)
388 call submesh_end_cube_map(sm)
389 call submesh_end(sm)
390 safe_deallocate_a(dorb)
391 safe_deallocate_a(zorb)
392 end do !inn
393
394 if(der%mesh%parallel_in_domains .and. this%ndim == 1) then
395 call der%mesh%allreduce(this%coulomb_IIJJ)
396 end if
397
398 if(dist%parallel) then
399 if (this%ndim == 1) then
400 call comm_allreduce(dist%mpi_grp, this%coulomb_IIJJ)
401 else
402 do inn = 1, this%nneighbors
403 do is2 = 1, this%ndim
404 do is1 = 1, this%ndim
405 call comm_allreduce(dist%mpi_grp, this%zcoulomb_IIJJ(:,:,:,:, is1, is2, inn))
406 end do
407 end do
408 end do
409 end if
410 end if
411
412 call distributed_end(dist)
413
414 call messages_print_with_emphasis(namespace=namespace)
415
417 end subroutine orbitalset_init_intersite
418
419#include "undef.F90"
420#include "real.F90"
421#include "orbitalset_utils_inc.F90"
422
423#include "undef.F90"
424#include "complex.F90"
425#include "orbitalset_utils_inc.F90"
426
427end module orbitalset_utils_oct_m
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
subroutine, public distributed_end(this)
subroutine, public distributed_nullify(this, total)
subroutine, public distributed_init(this, total, comm, tag, scalapack_compat)
Distribute N instances across M processes of communicator comm
real(real64), parameter, public m_zero
Definition: global.F90:188
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:115
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:118
subroutine, public messages_print_with_emphasis(msg, iunit, namespace)
Definition: messages.F90:920
character(len=512), private msg
Definition: messages.F90:165
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:160
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:616
integer, parameter, public sm_poisson_psolver
subroutine zget_orbital(species, sm, i, l, j_in, iorb, ndim, orb, combine_j_orbitals)
Returns an orbital on the intersite submesh.
subroutine, public dorbitalset_utils_getorbitals(namespace, os, ions, mesh, use_mesh, normalize, index_shift)
integer, parameter, public sm_poisson_isf
subroutine, public zorbitalset_get_center_of_mass(os, space, mesh, latt)
subroutine, public zorbitalset_utils_getorbitals(namespace, os, ions, mesh, use_mesh, normalize, index_shift)
subroutine, public orbitalset_init_intersite(this, namespace, space, ind, ions, der, psolver, os, nos, maxnorbs, rcut, kpt, has_phase, sm_poisson, basis_from_states, combine_j_orbitals)
subroutine dget_orbital(species, sm, i, l, j_in, iorb, ndim, orb, combine_j_orbitals)
Returns an orbital on the intersite submesh.
subroutine, public dorbitalset_get_center_of_mass(os, space, mesh, latt)
integer, parameter, public sm_poisson_fft
integer function, public orbitalset_utils_count(species, iselect)
Count the number of orbital sets we have for a given atom.
subroutine, public zpoisson_solve_sm(this, namespace, sm, pot, rho, all_nodes)
Calculates the Poisson equation. Given the density returns the corresponding potential.
Definition: poisson.F90:2346
integer, parameter, public poisson_psolver
Definition: poisson.F90:189
integer, parameter, public poisson_fft
Definition: poisson.F90:189
subroutine, public dpoisson_solve_sm(this, namespace, sm, pot, rho, all_nodes)
Calculates the Poisson equation. Given the density returns the corresponding potential.
Definition: poisson.F90:2161
integer, parameter, public poisson_direct_sum
Definition: poisson.F90:189
subroutine, public poisson_init_sm(this, namespace, space, main, der, sm, method, force_cmplx)
Definition: poisson.F90:999
integer, parameter, public poisson_isf
Definition: poisson.F90:189
subroutine, public poisson_end(this)
Definition: poisson.F90:709
subroutine, public submesh_end_global(this)
Definition: submesh.F90:870
complex(real64) function, public zsm_integrate(mesh, sm, ff, reduce)
Definition: submesh.F90:1702
subroutine, public submesh_shift_center(this, space, newcenter)
Definition: submesh.F90:631
real(real64) function, public dsm_integrate(mesh, sm, ff, reduce)
Definition: submesh.F90:1165
subroutine, public submesh_merge(this, space, mesh, sm1, sm2, shift)
Definition: submesh.F90:565
subroutine, public submesh_end_cube_map(sm)
Definition: submesh.F90:1045
subroutine, public submesh_end(this)
Definition: submesh.F90:735
subroutine, public submesh_build_global(this, space)
Definition: submesh.F90:819
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:132
The following class implements a lattice iterator. It allows one to loop over all cells that are with...