Branch data 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
# # # # #
# # # ]
|
