Line data Source code
1 : !!****m* ABINIT/m_supercell_maker
2 : !! NAME
3 : !! m_supercell_maker
4 : !!
5 : !! FUNCTION
6 : !! This module define the supercell_maker file, which provide functions to help build
7 : !! potentials in supercell.
8 : !! Note: This module works at lower level than m_supercell. The purpose is not to get
9 : !! an supercell of crystal structure, but it provide the functions to help building such thing (but not limited to).
10 : !!
11 : !! TODO: update m_supercell based on this so that it can use non-diagonal supercell matrix.
12 : !! Then this need to be moved to level same as m_supercell.
13 : !!
14 : !! Datatypes:
15 : !! supercell_maker_t
16 : !!
17 : !! Subroutines:
18 : !!
19 : !!
20 : !! COPYRIGHT
21 : !! Copyright (C) 2001-2026 ABINIT group (hexu)
22 : !! This file is distributed under the terms of the
23 : !! GNU General Public License, see ~abinit/COPYING
24 : !! or http://www.gnu.org/copyleft/gpl.txt .
25 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
26 : !!
27 : !! SOURCE
28 :
29 :
30 : #if defined HAVE_CONFIG_H
31 : #include "config.h"
32 : #endif
33 : #include "abi_common.h"
34 :
35 : module m_supercell_maker
36 : use defs_basis
37 : use m_abicore
38 : use m_errors
39 : use m_xmpi
40 : use m_mathfuncs, only: binsearch_left_integerlist, rotate_by_angle_around_axis
41 : use m_matrix, only: matr3inv, mat33det
42 : use m_mpi_scheduler, only: init_mpi_info
43 : use m_supercell
44 : implicit none
45 : private
46 : !!***
47 : !-----------------------------------------------------------------------
48 : !> @brief A helper type to build supercell and potentials from primitive
49 : !> It keep the information of the supercell matrix, and the R-vectors (R1~Rn),
50 : !> so that the potentials in
51 : !> For a list of indices in primitive cell (e.g. (i, j, k)), it will become
52 : !> (i,R0), (j, R0), (k, R0), (i, R1), (i, R1), ... (i, Rn), (j,Rn), (k, Rn)
53 : !-----------------------------------------------------------------------
54 : type ,public :: supercell_maker_t
55 : integer :: scmat(3,3) ! supercell matrix
56 : real(dp) :: inv_scmat(3,3) ! inverse of supercell matrix
57 : integer :: ncells ! number of cells in supercell
58 : integer, allocatable :: rvecs(:,:) ! R vectors for cells in supercell. dim:(3, ncells)
59 : contains
60 : procedure :: initialize
61 : procedure :: finalize
62 : procedure :: to_red_sc
63 : procedure :: build_rvec
64 : procedure :: sc_cell
65 : procedure :: R_to_sc
66 :
67 : ! Translations: from one primitive cell element to ncell supercell elements
68 : ! with values unchanged (therefore just repeat ncells times)
69 : generic :: repeat => repeat_int1d, repeat_real1d, repeat_real2d, repeat_realmat
70 : procedure :: repeat_int1d
71 : procedure :: repeat_int1d_noalloc
72 : procedure :: repeat_real1d
73 : procedure :: repeat_real2d
74 : procedure :: repeat_realmat
75 :
76 : ! Translations: from one primitive cell element to ncell supercell elements
77 : ! with values changed (different from repeat)
78 : procedure :: trans_xred
79 : procedure :: trans_xcart
80 : procedure :: trans_i
81 : procedure :: trans_i_noalloc
82 : procedure :: trans_ilist
83 : procedure :: trans_j_and_Rj
84 : procedure :: trans_j_and_Rj_noalloc
85 : procedure :: trans_jlist_and_Rj
86 : procedure :: rvec_for_each
87 :
88 : ! Wave like quantities:
89 : procedure :: generate_wave_scalar ! s=s0*e^{ikR}
90 : procedure :: generate_spin_wave_vectorlist !vector with rotation angle= 2pi k * R
91 : end type supercell_maker_t
92 :
93 : public :: scmaker_unittest
94 : contains
95 :
96 : !-----------------------------------------------------------------------
97 : !> @brief initialize
98 : !> @param [in] sc_matrix: the supercell matrix which maps the primitive cell to sc
99 : !> it does not have to be diagonal.
100 : !-----------------------------------------------------------------------
101 6 : subroutine initialize(self, sc_matrix)
102 : class(supercell_maker_t), intent(inout) :: self
103 : integer, intent(in) :: sc_matrix(3, 3)
104 : real(dp) :: tmp(3,3)
105 : integer :: master, my_rank, comm, nproc, ierr
106 : logical :: iam_master
107 6 : call init_mpi_info(master, iam_master, my_rank, comm, nproc)
108 :
109 78 : self%scmat(:,:)=sc_matrix
110 6 : call xmpi_bcast(self%scmat, master, comm, ierr)
111 6 : self%ncells=abs(mat33det(self%scmat))
112 : ! call xmpi_bcast(self%ncells, master, comm, ierr)
113 18 : ABI_MALLOC(self%rvecs, (3, self%ncells))
114 : ! Why transpose?
115 78 : tmp(:,:)=transpose(self%scmat)
116 :
117 : ! save inverse of matrix so A^{-1} x can be easily calculated.
118 6 : call matr3inv(tmp, self%inv_scmat)
119 6 : call self%build_rvec()
120 : !call xmpi_bcast(self%inv_scmat, master, comm, ierr)
121 : !call xmpi_bcast(self%ncells, master, comm, ierr)
122 : !call xmpi_bcast(self%rvecs, master, comm, ierr)
123 6 : end subroutine initialize
124 :
125 : !-----------------------------------------------------------------------
126 : !> @brief Finalize
127 : !-----------------------------------------------------------------------
128 6 : subroutine finalize(self)
129 : class(supercell_maker_t), intent(inout) :: self
130 78 : self%scmat=0
131 6 : self%ncells=0
132 6 : if (allocated(self%rvecs)) then
133 6 : ABI_FREE(self%rvecs)
134 : end if
135 6 : end subroutine finalize
136 :
137 : !-----------------------------------------------------------------------
138 : !> @brief from reduced coordinate in primitive cell to that insupercell
139 : !> This does not use translation symmetry, so it is a ONE to ONE convert.
140 : !> @param [in] pos (reduced position in primitive cell)
141 : !> @param [out] ret reduce position in supercell
142 : !-----------------------------------------------------------------------
143 3619346 : function to_red_sc(self, pos) result(ret)
144 : class(supercell_maker_t), intent(inout) :: self
145 : real(dp), intent(in) :: pos(3)
146 : real(dp) :: ret(3)
147 61528882 : ret(:)=matmul(pos, self%inv_scmat)
148 3619346 : end function to_red_sc
149 :
150 : !-----------------------------------------------------------------------
151 : !> @brief build R vectors for supercells.
152 : !-----------------------------------------------------------------------
153 6 : subroutine build_rvec(self)
154 : class(supercell_maker_t), intent(inout) :: self
155 : real(dp):: scorners_newcell(8,3), corners(8,3), x(3), tmp(3)
156 : integer :: minr(3), maxr(3), ix, iy, iz, counter
157 : real(dp) :: eps=1d-9
158 6 : counter=0
159 : ! cornors
160 : scorners_newcell=transpose(reshape([0., 0., 0., 0., 0., 1., 0., 1., 0., &
161 : 0., 1., 1., 1., 0., 0., 1., 0., 1., &
162 6 : 1., 1., 0., 1., 1., 1.],[3,8]))
163 582 : corners = matmul(scorners_newcell, self%scmat)
164 : !rep=ceiling(maxval(corners, dim=1) - minval(corners, dim=1))
165 24 : minr=floor(minval(corners, dim=1))
166 24 : maxr=ceiling(maxval(corners, dim=1))
167 : ! NOTE: DO NOT CHANGE THE ORDER. It is used in the binary search.
168 48 : do ix = minr(1), maxr(1)
169 350 : do iy = minr(2), maxr(2)
170 2570 : do iz = minr(3), maxr(3)
171 8904 : x(:)=[ix, iy, iz]
172 8904 : tmp=self%to_red_sc(x)
173 : if ( (.not. any(tmp<=-1.0*eps)) &
174 14234 : .and. (.not. any(tmp>1.0-eps)) ) then
175 1440 : counter=counter+1
176 5760 : self%rvecs(:,counter)=nint(x)
177 : end if
178 : end do
179 : end do
180 : end do
181 6 : if (counter /= self%ncells ) then
182 0 : ABI_ERROR("Bug found. supercell_maker: Wrong number of supercell found in build_rvec. ")
183 : end if
184 6 : end subroutine build_rvec
185 :
186 : !-----------------------------------------------------------------------
187 : !> @brief cell parameters from unitcell to supercell
188 : !> @param [in] cell: cell parameter in primitive cell
189 : !> @param [out] sccell: supercell cell parameter
190 : !-----------------------------------------------------------------------
191 8 : function sc_cell(self, cell) result(sccell)
192 : class(supercell_maker_t), intent(inout) :: self
193 : real(dp), intent(in) :: cell(3,3)
194 : real(dp) :: sccell(3,3)
195 528 : sccell=matmul(self%scmat, cell)
196 8 : end function sc_cell
197 :
198 : !-----------------------------------------------------------------------
199 : !> @brief reduced position in primitive cell to supercell.
200 : !> It use translation symmetry. so (3,n) -> (3, n*ncells)
201 : !> Note that allocation will be done if not already allocated.
202 : !> @param [in] xred: reduced coordinates in primitive cell
203 : !> @param [out] sccell: reduced coordinates in supercell
204 : !-----------------------------------------------------------------------
205 0 : subroutine trans_xred(self, xred, scxred)
206 : class(supercell_maker_t), intent(inout) :: self
207 : real(dp), intent(in) :: xred(:,:)
208 : real(dp), allocatable, intent(inout) :: scxred(:, :)
209 : integer :: npos, icell, ipos, counter=0
210 0 : npos=size(xred, dim=2)
211 0 : if (.not. allocated(scxred)) then
212 0 : ABI_MALLOC(scxred, (3,size(xred, dim=2)*self%ncells))
213 : end if
214 0 : do icell = 1, self%ncells
215 0 : do ipos=1 , npos
216 0 : counter=counter+1
217 0 : scxred(:, counter) = self%to_red_sc(xred(:, ipos)+ self%rvecs(:, icell))
218 : end do
219 : end do
220 0 : end subroutine trans_xred
221 :
222 : !-----------------------------------------------------------------------
223 : !> @brief cartesion position in primitive cell to supercell.
224 : !> It use translation symmetry. so (3,n) -> (3, n*ncells)
225 : !> Note that allocation will be done if not already allocated.
226 : !> @param [in] primcell: cell parameter in primitive cell
227 : !> @param [in] xcart: cartesion coordinates in primitive cell
228 : !> @param [out] scxcart:cartesion coordinates in supercell
229 : !-----------------------------------------------------------------------
230 8 : subroutine trans_xcart(self, primcell, xcart, scxcart)
231 : class(supercell_maker_t), intent(inout) :: self
232 : real(dp), intent(in) :: xcart(:,:), primcell(3,3)
233 : real(dp), allocatable, intent(inout) :: scxcart(:,:)
234 : integer :: npos, icell, ipos, counter
235 8 : counter =0
236 8 : npos=size(xcart, dim=2)
237 8 : if (.not. allocated(scxcart)) then
238 18 : ABI_MALLOC(scxcart, (3,npos*self%ncells))
239 : end if
240 1880 : do icell = 1, self%ncells
241 10024 : do ipos=1 , npos
242 8144 : counter=counter+1
243 189184 : scxcart(:, counter) = matmul(self%rvecs(:, icell), primcell) + xcart(:,ipos)
244 : end do
245 : end do
246 8 : end subroutine trans_xcart
247 :
248 :
249 : !-----------------------------------------------------------------------
250 : !> @brief Find a R vector in primitive cell unit to a R vector in supercell unit
251 : !> Here is an example:
252 : !> A Rp for a 1*1*1 primitive cell. Rp=(0,0, 3)
253 : !> For a supercell of 1*1*2, Rp = (0, 0, 1) * supercell + (0, 0, 1) * primcell
254 : !> @param [in] R index using primitive cell parameter
255 : !> @param [out] R index using supercell parameter
256 : !> @param [out] index of cell INSIDE supercell in primitive cell.
257 : !-----------------------------------------------------------------------
258 3617120 : subroutine R_to_sc(self, R, R_sc, ind_sc)
259 : class(supercell_maker_t), intent(inout) :: self
260 : integer, intent(in) :: R(3)
261 : integer, intent(inout) :: ind_sc, R_sc(3)
262 : integer :: rprim(3)
263 25319840 : R_sc=floor(self%to_red_sc(R*1.0d0)+1.0d-8)
264 61491040 : rprim(:)= R-matmul(R_sc, self%scmat)
265 3617120 : ind_sc=binsearch_left_integerlist(self%rvecs, rprim)
266 3617120 : if (ind_sc==0) then
267 : !print *, "R: " , R
268 : !print *, "Rredsc:", self%to_red_sc(R*1.0d0)
269 : !print *, "R_sc: " , R_sc
270 : !print *, "rprim", rprim
271 0 : ABI_ERROR("Bug found. supercell_maker%R_to_sc: Cannot find rprim")
272 : end if
273 3617120 : end subroutine R_to_sc
274 :
275 :
276 :
277 : !-----------------------------------------------------------------------
278 : !> @brief generate a wave like quantity. A_sc(R) = A(0)* exp(ikR)
279 : !> The A_sc will be allocated. size : (ncells)
280 : !> @param [in] A: a scalar quantity.
281 : !> @param [in] kpoint: the wave vector
282 : !> @param [out] A_sc: indices in supercell.
283 : !-----------------------------------------------------------------------
284 0 : subroutine generate_wave_scalar(self, A, kpoint, A_sc)
285 : class(supercell_maker_t), intent(inout) :: self
286 : real(dp), intent(in) :: A
287 : real(dp), intent(in) :: kpoint(3)
288 : real(dp), allocatable, intent(inout) :: A_sc(:)
289 : integer :: i
290 0 : if(.not. allocated(A_sc)) then
291 0 : ABI_MALLOC(A_sc, (self%ncells))
292 : end if
293 0 : do i=1, self%ncells
294 : A_sc(i) = real(A*exp(cmplx(0.0,two_pi, kind=dp) * &
295 0 : &dot_product(kpoint, self%rvecs(:, i))), kind=dp)
296 : end do
297 0 : end subroutine generate_wave_scalar
298 :
299 : !-----------------------------------------------------------------------
300 : !> @brief generate a wave like quantity. A_sc(R) = A(0)* exp(ikR)
301 : !> The A_sc will be allocated. size : (nA*ncells)
302 : !> @param [in] A: a list scalar quantity in primitive cell.
303 : !> @param [in] kpoint: the wave vector
304 : !> @param [out] A_sc: indices in supercell.
305 : !-----------------------------------------------------------------------
306 : subroutine generate_wave_scalarlist(self, A, kpoint, A_sc)
307 : class(supercell_maker_t), intent(inout) :: self
308 : real(dp), intent(in) :: A(:)
309 : real(dp), intent(in) :: kpoint(3)
310 : real(dp), allocatable, intent(inout) :: A_sc(:)
311 : integer :: icell, iA, nA, counter
312 : nA=size(A)
313 : if(.not. allocated(A_sc)) then
314 : ABI_MALLOC(A_sc, (nA*self%ncells))
315 : end if
316 : counter=0
317 : do icell = 1, self%ncells
318 : do iA=1 , nA
319 : counter=counter+1
320 : A_sc(counter) = real(A(iA)*exp(cmplx(0.0,two_pi, kind=dp) * &
321 : &dot_product(kpoint, self%rvecs(:, icell))), kind=dp)
322 : end do
323 : end do
324 : end subroutine generate_wave_scalarlist
325 :
326 :
327 : !-----------------------------------------------------------------------
328 : !> @brief generate a wave like quantity. A_sc(R) = A(0)* exp(ikR)
329 : !> The A_sc will be allocated. size : (ndim, nA*ncells)
330 : !> NOTE: This generate a density wave. For a spin wave,
331 : !> use generate_spin_wave_vectorlist.
332 : !> @param [in] A: a vector quantity. dimension (xyz, iA)
333 : !> @param [in] kpoint: the wave vector
334 : !> @param [out] A_sc: indices in supercell.
335 : !-----------------------------------------------------------------------
336 : subroutine generate_wave_vectorlist(self, A, kpoint, A_sc)
337 : class(supercell_maker_t), intent(inout) :: self
338 : real(dp), intent(in) :: A(:, :)
339 : real(dp), intent(in) :: kpoint(3)
340 : real(dp), allocatable, intent(inout) :: A_sc(:, :)
341 : integer :: icell, iA, nA, counter, ndim
342 : ndim=size(A, 1)
343 : nA=size(A,2)
344 : if(.not. allocated(A_sc)) then
345 : ABI_MALLOC(A_sc, (ndim, nA*self%ncells))
346 : end if
347 : counter=0
348 : do icell = 1, self%ncells
349 : do iA=1 , nA
350 : counter=counter+1
351 : A_sc(:,counter) = real(A(:,iA)*exp(cmplx(0.0,two_pi, kind=dp) * &
352 : &dot_product(kpoint, self%rvecs(:, icell))), kind=dp)
353 : end do
354 : end do
355 : end subroutine generate_wave_vectorlist
356 :
357 :
358 : !-----------------------------------------------------------------------
359 : !> @brief generate a spin wave quantity.
360 : !> rotate by theta(R) = 2pi * k.dot.R, around axis.
361 : !> The A_sc will be allocated. size : (ndim, nA*ncells)
362 : !> @param [in] A: a vector quantity. dimension (xyz, iA)
363 : !> @param [in] kpoint: the wave vector
364 : !> @param [in] axis: the axis to rotate around
365 : !> @param [out] A_sc: indices in supercell.
366 : !-----------------------------------------------------------------------
367 2 : subroutine generate_spin_wave_vectorlist(self, A, kpoint, axis, A_sc)
368 : class(supercell_maker_t), intent(inout) :: self
369 : real(dp), intent(in) :: A(:, :) !(3, dimension of primitive cell)
370 : real(dp), intent(in) :: kpoint(3)
371 : real(dp), intent(in) :: axis(3)
372 : real(dp), allocatable, intent(inout) :: A_sc(:, :) !(3, dimension of supercell)
373 : integer :: icell, iA, nA, counter, ndim
374 : real(dp):: angle
375 2 : ndim=size(A, 1)
376 2 : nA=size(A,2)
377 2 : if(.not. allocated(A_sc)) then
378 0 : ABI_MALLOC(A_sc, (ndim, nA*self%ncells))
379 : end if
380 2 : counter=0
381 434 : do icell = 1, self%ncells
382 866 : do iA=1, nA
383 432 : counter=counter+1
384 1728 : angle=two_pi* DOT_PRODUCT(kpoint, self%rvecs(:, icell))
385 2160 : A_sc(:,counter) = rotate_by_angle_around_axis(angle=angle, axis=axis, vec=A(:, iA))
386 : end do
387 : end do
388 2 : end subroutine generate_spin_wave_vectorlist
389 :
390 :
391 : !-----------------------------------------------------------------------
392 : !> @brief Get the indices of all the supercell repeat of an index in primitive cell
393 : !> @param [in] nbasis: total number of indices in the primitive cell (eg. natom in a cell)
394 : !> Note that if the dimension xyz is also in the index, then use 3*natom.
395 : !> i_sc array will be allocated if not already done
396 : !> This function is mostly useful for building supercell potential with such format H_ij(R_j)
397 : !> @param [in] i: index in primitive cell
398 : !> @param [out] i_sc: indices in supercell.
399 : !-----------------------------------------------------------------------
400 8309 : subroutine trans_i(self, nbasis, i, i_sc)
401 : class(supercell_maker_t), intent(inout) :: self
402 : integer, intent(in) :: i, nbasis
403 : integer, allocatable, intent(inout) :: i_sc(:)
404 8309 : if(.not. allocated(i_sc)) then
405 16347 : ABI_MALLOC(i_sc, (self%ncells))
406 : end if
407 8309 : call self%trans_i_noalloc(nbasis, i, i_sc)
408 8309 : end subroutine trans_i
409 :
410 : !-----------------------------------------------------------------------
411 : !> @brief Get the indices of all the supercell repeat of an index in primitive cell
412 : !> @param [in] nbasis: total number of indices in the primitive cell (eg. natom in a cell)
413 : !> Note that if the dimension xyz is also in the index, then use 3*natom.
414 : !> i_sc array will NOT be allocated.
415 : !> This function is mostly useful for building supercell potential with such format H_ij(R_j)
416 : !> @param [in] i: index in primitive cell
417 : !> @param [out] i_sc: indices in supercell.
418 : !-----------------------------------------------------------------------
419 8309 : subroutine trans_i_noalloc(self, nbasis, i, i_sc)
420 : class(supercell_maker_t), intent(inout) :: self
421 : integer, intent(in) :: i, nbasis
422 : integer, intent(inout) :: i_sc(:)
423 : integer :: icell
424 1757397 : do icell =1, self%ncells
425 1757397 : i_sc(icell)=nbasis*(icell-1)+i
426 : end do
427 8309 : end subroutine trans_i_noalloc
428 :
429 :
430 : !-----------------------------------------------------------------------
431 : !> @brief Get the indices of all the supercell repeat of an LIST of indices in primitive cell
432 : !> @param [in] nbasis: total number of indices in the primitive cell (eg. natom in a cell)
433 : !> Note that if the dimension xyz is also in the index, then use 3*natom.
434 : !> i_sc array will be allocated if not already done.
435 : !> This function is mostly useful for building supercell potential with such format H_ij(R_j)
436 : !> @param [in] ilist: list of indices in primitive cell
437 : !> @param [out] ilist_sc: list of indices in supercell.
438 : !-----------------------------------------------------------------------
439 0 : subroutine trans_ilist(self, nbasis, ilist, ilist_sc)
440 : class(supercell_maker_t), intent(inout) :: self
441 : integer, intent(in) :: ilist(:), nbasis
442 : integer, allocatable, intent(inout) :: ilist_sc(:)
443 : integer :: i
444 0 : if(.not. allocated(ilist_sc)) then
445 0 : ABI_MALLOC(ilist_sc, (size(ilist)*self%ncells))
446 : end if
447 0 : do i =1, size(ilist)
448 0 : call trans_i_noalloc(self,nbasis, ilist(i), ilist_sc(self%ncells*(i-1)+1:self%ncells*i ))
449 : end do
450 0 : end subroutine trans_ilist
451 :
452 : !-----------------------------------------------------------------------
453 : !> @brief Get the indices and R vector for a index of j displaced by a R vector.
454 : !> This is often used to represent a vector outside the primitive cell.
455 : !> i.e. its position is (r_j + R_j) in primitive cell unit. Then it will have
456 : !> periodical repeats, which has index of j_sc inside supercell, and a R_sc vector to shift it.
457 : !> @param [in] nbasis: total number of indices in the primitive cell (eg. natom in a cell)
458 : !> Note that if the dimension xyz is also in the index, then use 3*natom.
459 : !> i_sc array will be allocated if not already done.
460 : !> This function is mostly useful for building supercell potential with such format H_ij(R_j)
461 : !> @param [in] j: index in primitive cell
462 : !> @param [in] Rj: R vector in primitive cell unit
463 : !> @param [out] j_sc: index in supercell
464 : !> @param [out] Rj_sc: R' vector in supercell unit.
465 : !-----------------------------------------------------------------------
466 15121 : subroutine trans_j_and_Rj(self, nbasis, j, Rj, j_sc, Rj_sc)
467 : class(supercell_maker_t), intent(inout) :: self
468 : integer, intent(in) :: j, Rj(3), nbasis
469 : integer, allocatable , intent(inout) :: j_sc(:), Rj_sc(:, :)
470 15121 : if(.not. allocated(j_sc)) then
471 36783 : ABI_MALLOC(j_sc, (self%ncells))
472 : endif
473 15121 : if(.not. allocated(Rj_sc)) then
474 36783 : ABI_MALLOC(Rj_sc, (3, self%ncells))
475 : endif
476 15121 : call self%trans_j_and_Rj_noalloc(nbasis, j, Rj, j_sc, Rj_sc)
477 15121 : end subroutine trans_j_and_Rj
478 :
479 : !-----------------------------------------------------------------------
480 : !> @brief same as trans_j_and_Rj, but it does not allocate data.
481 : !-----------------------------------------------------------------------
482 15121 : subroutine trans_j_and_Rj_noalloc(self, nbasis, j, Rj, j_sc, Rj_sc)
483 : class(supercell_maker_t), intent(inout) :: self
484 : integer, intent(in) :: j, Rj(3), nbasis
485 : integer, intent(inout) :: j_sc(:), Rj_sc(:, :)
486 : integer :: i,jj
487 3632241 : do i =1, self%ncells
488 14468480 : call self%R_to_sc(Rj + self%rvecs(:,i), Rj_sc(:,i), jj)
489 3632241 : j_sc(i)=nbasis*(jj-1)+j
490 : end do
491 15121 : end subroutine trans_j_and_Rj_noalloc
492 :
493 : !-----------------------------------------------------------------------
494 : !> @brief same as trans_j_and_Rj, but it loop over a list of j with the
495 : !> same Rj
496 : !-----------------------------------------------------------------------
497 0 : subroutine trans_jlist_and_Rj(self, nbasis, jlist, Rj, ind_sc, R_sc)
498 : class(supercell_maker_t), intent(inout) :: self
499 : integer, intent(in) :: jlist(:), Rj(3), nbasis
500 : integer, allocatable, intent(inout) :: ind_sc(:), R_sc(:,:)
501 : integer :: i,jj, counter, indj
502 0 : if (.not. allocated(ind_sc)) then
503 0 : ABI_MALLOC(ind_sc, (self%ncells*size(jlist)) )
504 : endif
505 0 : if (.not. allocated(R_sc)) then
506 0 : ABI_MALLOC(R_sc, (3, self%ncells))
507 : endif
508 0 : counter=0
509 0 : do i =1, self%ncells
510 0 : call self%R_to_sc(Rj + self%rvecs(:,i), R_sc(:,i), jj)
511 0 : do indj=1, size(jlist)
512 0 : counter=counter+1
513 0 : ind_sc(counter)=nbasis*(jj-1)+jlist(indj)
514 : end do
515 : end do
516 0 : end subroutine trans_jlist_and_Rj
517 :
518 : !-----------------------------------------------------------------------
519 : !> @brief same as trans_j_and_Rj, but it loop over a list of j with the
520 : !> same Rj
521 : !-----------------------------------------------------------------------
522 : subroutine trans_jlist_and_Rj_noalloc(self, nbasis, jlist, Rj, ind_sc, R_sc)
523 : class(supercell_maker_t), intent(inout) :: self
524 : integer, optional, intent(in) :: jlist(:)
525 : integer, intent(in) :: Rj(3), nbasis
526 : integer, intent(inout) :: ind_sc(:), R_sc(:,:)
527 : integer :: i,jj, counter, indj, n
528 : counter=0
529 : if (present(jlist)) then
530 : n=size(jlist)
531 : else
532 : n=nbasis
533 : end if
534 : do i =1, self%ncells
535 : call self%R_to_sc(Rj + self%rvecs(:,i), R_sc(:,i), jj)
536 : do indj=1, n
537 : counter=counter+1
538 : ind_sc(counter)=nbasis*(jj-1)+jlist(indj)
539 : end do
540 : end do
541 : end subroutine trans_jlist_and_Rj_noalloc
542 :
543 :
544 :
545 : !-----------------------------------------------------------------------
546 : !> @brief same as trans_j_and_Rj, but it loop over a list of j with the
547 : !> same Rj
548 : !-----------------------------------------------------------------------
549 : subroutine trans_ijR(self, nbasis_i, nbasis_j, ilist, jlist, Rj, i_sc, j_sc, R_sc)
550 : class(supercell_maker_t), intent(inout) :: self
551 : integer, intent(in) ::ilist(:), jlist(:), Rj(3), nbasis_i, nbasis_j
552 : integer, allocatable, intent(inout) :: i_sc(:), j_sc(:), R_sc(:,:)
553 : integer :: icell,jj, counter, indj
554 : if (.not. allocated(i_sc)) then
555 : ABI_MALLOC(i_sc, (self%ncells*size(jlist)) )
556 : endif
557 : if (.not. allocated(j_sc)) then
558 : ABI_MALLOC(j_sc, (self%ncells*size(jlist)) )
559 : endif
560 : if (.not. allocated(R_sc)) then
561 : ABI_MALLOC(R_sc, (3, self%ncells))
562 : endif
563 : counter=0
564 : do icell =1, self%ncells
565 : call self%R_to_sc(Rj + self%rvecs(:,icell), R_sc(:,icell), jj)
566 : do indj=1, size(jlist)
567 : counter=counter+1
568 : i_sc(counter)=nbasis_i*(icell-1)+ilist(indj)
569 : j_sc(counter)=nbasis_j*(jj-1)+jlist(indj)
570 : end do
571 : end do
572 : end subroutine trans_ijR
573 :
574 :
575 :
576 : !-----------------------------------------------------------------------
577 : !> @brief repeat a quantity (which is a scalar for each index)
578 : !> memory will be allocated if not already done
579 : !> @param [in] a: the quantity to be repeated (loop over index i in primcell)
580 : !> @param [out] ret: the repeat. also a 1D matrix. (loop over index in supercell)
581 : !-----------------------------------------------------------------------
582 9 : subroutine repeat_int1d(self, a, ret)
583 : class(supercell_maker_t), intent(inout) :: self
584 : integer, intent(in) :: a(:)
585 : integer, allocatable :: ret(:)
586 : integer :: n, i
587 9 : n=size(a)
588 9 : if (.not. allocated(ret)) then
589 21 : ABI_MALLOC(ret, (n*self%ncells))
590 : end if
591 2393 : do i =1, self%ncells
592 11561 : ret((i-1)*n+1: i*n) = a(:)
593 : end do
594 9 : end subroutine repeat_int1d
595 :
596 : !-----------------------------------------------------------------------
597 : !> @brief repeat a quantity (which is a scalar for each index)
598 : !> memory will NOT be allocated.
599 : !> @param [in] a: the quantity to be repeated (loop over index i in primcell)
600 : !> @param [out] ret: the repeat. also a 1D matrix. (loop over index in supercell)
601 : !-----------------------------------------------------------------------
602 0 : subroutine repeat_int1d_noalloc(self, a, ret)
603 : class(supercell_maker_t), intent(inout) :: self
604 : integer, intent(in) :: a(:)
605 : integer :: ret(:)
606 : integer :: n, i
607 0 : n=size(a)
608 0 : do i =1, self%ncells
609 0 : ret((i-1)*n+1: i*n) = a(:)
610 : end do
611 0 : end subroutine repeat_int1d_noalloc
612 :
613 : !-----------------------------------------------------------------------
614 : !> @brief repeat a quantity (which is a 1D matrix for each index)
615 : !> memory will be allocated if not already done
616 : !> @param [in] a: the quantity to be repeated (loop over index i in primcell)
617 : !> @param [out] ret: the repeat. also a 1D matrix. (loop over index in supercell)
618 : !-----------------------------------------------------------------------
619 13 : subroutine repeat_real1d(self, a, ret)
620 : class(supercell_maker_t), intent(inout) :: self
621 : real(dp), intent(in) :: a(:)
622 : real(dp), allocatable:: ret(:)
623 : integer :: n, i
624 13 : n=size(a)
625 13 : if (.not. allocated(ret)) then
626 21 : ABI_MALLOC(ret, (n*self%ncells))
627 : end if
628 3261 : do i =1, self%ncells
629 13293 : ret((i-1)*n+1: i*n) = a(:)
630 : end do
631 13 : end subroutine repeat_real1d
632 :
633 : !-----------------------------------------------------------------------
634 : !> @brief repeat a quantity (which is a 1D matrix for each index)
635 : !> memory will be allocated if not already done
636 : !> @param [in] a: the quantity to be repeated ( 2nd dim loop over index i in primcell)
637 : !> @param [out] ret: the repeat. also a 2D matrix. (2nd dim loop over index in supercell)
638 : !-----------------------------------------------------------------------
639 0 : subroutine repeat_real2d(self, a, ret)
640 : class(supercell_maker_t), intent(inout) :: self
641 : real(dp), intent(in) :: a(:,:)
642 : real(dp), allocatable:: ret(:, :)
643 : integer :: n1, n2, i
644 0 : n1=size(a, dim=1)
645 0 : n2=size(a,dim=2)
646 0 : if (.not. allocated(ret)) then
647 0 : ABI_MALLOC(ret, (n1, n2*self%ncells))
648 : end if
649 0 : do i =1, self%ncells
650 0 : ret(:,(i-1)*n2+1: i*n2) = a(:,:)
651 : end do
652 0 : end subroutine repeat_real2d
653 :
654 : !-----------------------------------------------------------------------
655 : !> @brief repeat a quantity (which is a 2D matrix for each index)
656 : !> memory will be allocated if not already done
657 : !> @param [in] a: the quantity to be repeated ( 3rd dim loop over index i in primcell)
658 : !> @param [out] ret: the repeat. also a 2D matrix. (3rd dim loop over index in supercell)
659 : !-----------------------------------------------------------------------
660 0 : subroutine repeat_realmat(self, a, ret)
661 : class(supercell_maker_t), intent(inout) :: self
662 : real(dp), intent(in) :: a(:,:,:)
663 : real(dp), allocatable, intent(inout) :: ret(:,:,:)
664 : integer :: n, i
665 0 : n=size(a, 3)
666 0 : if (.not. allocated(ret)) then
667 0 : ABI_MALLOC(ret, (size(a, dim=1), size(a, dim=2), n*self%ncells))
668 : end if
669 0 : do i=1, self%ncells
670 0 : ret(:,:,(i-1)*n+1: i*n)=a(:, :, :)
671 : end do
672 0 : end subroutine repeat_realmat
673 :
674 : !-----------------------------------------------------------------------
675 : !> @brief a list of R vectors for all the indices in supercell
676 : !> memory will be allocated if not already done.
677 : !> @param [in] nbasis: number of indices in primitivecell
678 : !> @param [ret] A R vector for each index in supercell mat(3, nbasis*ncells)
679 : !-----------------------------------------------------------------------
680 3 : subroutine rvec_for_each(self, nbasis, ret)
681 : class(supercell_maker_t), intent(inout) :: self
682 : integer, intent(in) :: nbasis
683 : integer, allocatable, intent(inout) :: ret(:,:)
684 : integer :: icell, ibasis
685 3 : if (.not. allocated(ret)) then
686 3 : ABI_MALLOC(ret, (3, nbasis*self%ncells))
687 : end if
688 947 : do icell=1, self%ncells
689 2403 : do ibasis=1, nbasis
690 6768 : ret(:, (icell-1)*nbasis+ibasis)= self%rvecs(:, icell)
691 : end do
692 : end do
693 3 : end subroutine rvec_for_each
694 :
695 : !============================Unit test=====================================
696 :
697 0 : function test1() result(err)
698 0 : type(supercell_maker_t) :: maker
699 : integer :: err
700 : integer :: scmat(3,3)
701 0 : integer, allocatable :: ind_sc(:)
702 : integer :: R(3), R_sc(3), ind_sc2
703 0 : integer, allocatable :: j_sc(:), R_sc3(:, :)
704 0 : scmat=reshape([2,0, 0, 0,2, 0, 0,0,2], [3,3])
705 0 : call maker%initialize(scmat)
706 0 : err=0
707 0 : call maker%trans_i(nbasis=3, i=1, i_sc=ind_sc)
708 0 : R=[0,0,0]
709 0 : call maker%R_to_sc(R, R_sc, ind_sc2)
710 0 : if (.not.( (all(R_sc==[0,0,0])) .and. (ind_sc2==1 ))) then
711 0 : ABI_ERROR("R_to_sc is wrong")
712 : end if
713 :
714 0 : R=[0,0,3]
715 0 : call maker%trans_j_and_Rj(nbasis=3, j=1, Rj=R, j_sc=j_sc, Rj_sc=R_sc3 )
716 0 : if (.not. ( all(j_sc(1:3)==[4,1,10]) &
717 : .and. all(R_sc3(:,1)==[0,0,1]) &
718 : .and. all(R_sc3(:,3)==[0,0,1]) ) ) then
719 0 : ABI_ERROR("Wrong trans_j_and_Rj")
720 0 : err=1
721 : end if
722 0 : ABI_FREE(ind_sc)
723 0 : ABI_FREE(j_sc)
724 0 : ABI_FREE(R_sc3)
725 0 : end function test1
726 :
727 0 : function test2() result(err)
728 0 : type(supercell_maker_t) :: maker
729 : integer :: err
730 : integer :: scmat(3,3)
731 : real(dp) :: sccell(3,3)
732 0 : real(dp), allocatable :: scxred(:, :)
733 0 : integer, allocatable :: rep1(:)
734 0 : real(dp), allocatable :: rep2(:)
735 0 : scmat=transpose(reshape([1,2, 3, 4,5,6,7,8,6], [3,3]))
736 : !scmat=reshape([1,2, 3, 4,5,6,7,8,6], [3,3])
737 0 : call maker%initialize(scmat)
738 :
739 0 : err=0
740 :
741 : ! test build_rvecs
742 0 : if (.not. all(maker%rvecs==reshape([0, 0, 0, 2, 3, 4, 3, 4, 4, &
743 : 3, 4, 5, 4, 5, 5, 5, 6, 5, &
744 : 5, 6, 6, 6, 7, 6, 8, 10, 10], [3, 8]))) then
745 0 : ABI_ERROR("Wrong Rvecs found !")
746 0 : err=1
747 : end if
748 :
749 : ! test sc_cell
750 0 : sccell(:,:)=maker%sc_cell(transpose(reshape([1.d0,2.d0,3.d0,4.d0,5.d0,6.d0,3.d0,2.d0,3.d0], [3,3])))
751 0 : if (.not. all(abs(sccell-transpose(reshape([18,18, 24, 42, 45, 60, 57, 66, 87], [3,3])))<1e-4)) then
752 0 : ABI_ERROR("Wrong cell paramter found !")
753 0 : err=1
754 : end if
755 :
756 :
757 : ! test trans_xred
758 0 : call maker%trans_xred(reshape([0.1d0, 0.2d0, 0.3d0, 0.1d0, 0.4d0, 0.6d0], [3, 2]), scxred)
759 0 : if (.not. all(abs(scxred(:,18)-[1.0666666666666684,0.53333333333333277,0.70000000000000040])<1e-4)) then
760 0 : ABI_ERROR("Wrong sc_xred !")
761 0 : err=1
762 : end if
763 :
764 : ! test repeat
765 0 : call maker%repeat([1,2], rep1)
766 :
767 0 : call maker%repeat([1.0d0,2.0d0], rep2)
768 :
769 0 : end function test2
770 :
771 0 : subroutine scmaker_unittest()
772 0 : if (test1()/=0) ABI_ERROR("Supercell maker: Test1 Failed")
773 0 : if (test2()/=0) ABI_ERROR("Supercell maker: Test2 Failed")
774 0 : end subroutine scmaker_unittest
775 :
776 :
777 : !===================================================================================
778 : ! The functions below are deprecated!!!
779 : ! ! They are used with the m_supercell module.
780 : !
781 : ! ! R (in term of primitive cell) to R_sc(in term of supercell) + R_prim
782 : ! subroutine find_R_PBC(scell, R, R_sc, R_prim)
783 : ! type(supercell_type) , intent(in):: scell
784 : ! integer, intent(in):: R(3)
785 : ! integer, intent(out):: R_sc(3), R_prim(3)
786 : ! real(dp) :: R_sc_d(3), sc_mat(3,3)
787 :
788 : ! integer:: ipriv(3), info
789 : ! !call dgesv( n, nrhs, a, lda, ipiv, b, ldb, info )
790 : ! sc_mat(:,:)=scell%rlatt
791 : ! R_sc_d(:)=R(:)
792 : ! call dgesv(3, 1, sc_mat, 3, ipriv, R_sc_d, 3, info)
793 : ! if ( info/=0 ) then
794 : ! ABI_ERROR("Failed to find R_sc")
795 : ! end if
796 :
797 : ! ! if only diagonal of rlatt works.
798 : ! !R_sc_d(1)= real(R(1))/real(scell%rlatt(1,1))
799 : ! !R_sc_d(2)= real(R(2))/real(scell%rlatt(2,2))
800 : ! !R_sc_d(3)= real(R(3))/real(scell%rlatt(3,3))
801 : ! ! TODO hexu: R_prim should be non-negative, which is assumed in m_supercell.
802 : ! ! But should we make it more general?
803 : ! R_sc(1)=floor(R_sc_d(1))
804 : ! R_sc(2)=floor(R_sc_d(2))
805 : ! R_sc(3)=floor(R_sc_d(3))
806 : ! R_prim(1)=(R(1)-R_sc(1)*scell%rlatt(1,1))
807 : ! R_prim(2)=(R(2)-R_sc(2)*scell%rlatt(2,2))
808 : ! R_prim(3)=(R(3)-R_sc(3)*scell%rlatt(3,3))
809 : ! end subroutine find_R_PBC
810 :
811 : ! ! TODO hexu: move this to m_supercell?
812 : ! ! find the spercelll atom index from index of atom in primitive cell and R vector
813 : ! function find_supercell_index(scell, iatom_prim, rvec) result(iatom_supercell)
814 : ! type(supercell_type) , intent(in):: scell
815 : ! integer, intent(in) :: iatom_prim, rvec(3)
816 : ! integer :: iatom_supercell
817 : ! integer :: i
818 : ! iatom_supercell=-1
819 : ! do i=1, scell%natom, 1
820 : ! if ( scell%atom_indexing(i) == iatom_prim .and. &
821 : ! all(scell%uc_indexing(:,i)==rvec) ) then
822 : ! iatom_supercell=i
823 : ! return
824 : ! end if
825 : ! end do
826 : ! ABI_ERROR("BUG found. supercell_maker%find_supercell_index cannot find iatom_prim, rvec pair in supercell")
827 : ! end function find_supercell_index
828 :
829 : ! !i0, j0+R0 shifted by R to i1=i0+0+R->periodic, j1=j0+R0+R->periodic
830 : ! subroutine find_supercell_ijR(scell, i0, j0, R0, R, i1, j1, R1, R_sc)
831 :
832 : ! type(supercell_type) , intent(in):: scell
833 : ! integer, intent(in) :: i0, j0, R0(3), R(3)
834 : ! integer, intent(out) :: i1, j1, R1(3), R_sc(3)
835 : ! i1=find_supercell_index(scell,i0,R)
836 : ! call find_R_PBC(scell,R0+R,R_sc,R1)
837 : ! j1=find_supercell_index(scell, j0, R1)
838 : ! end subroutine find_supercell_ijR
839 :
840 :
841 0 : end module m_supercell_maker
|
