Branch data Line data Source code
1 : : !!****m* ABINIT/m_ifc_supercell
2 : : !! NAME
3 : : !! m_ifc_supercell
4 : : !!
5 : : !! FUNCTION
6 : : !! Module for defining a supercell for the purpose of IFC
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2026-2026 ABINIT group (GA)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public Licence, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
14 : : !!
15 : : !! SOURCE
16 : :
17 : :
18 : : #if defined HAVE_CONFIG_H
19 : : #include "config.h"
20 : : #endif
21 : :
22 : : #include "abi_common.h"
23 : :
24 : : module m_ifc_supercell
25 : :
26 : : use defs_basis
27 : : use m_io_tools
28 : : use m_abicore
29 : : use m_errors
30 : : use m_xmpi
31 : : use m_numeric_tools, only : l2norm, wrap2_pmhalf
32 : : use m_matrix, only : matr3inv, mati3det, mati3HNF
33 : : use m_atomdata, only : znucl2symbol
34 : : use m_crystal, only : crystal_t
35 : :
36 : : implicit none
37 : :
38 : : private
39 : :
40 : : public :: ifc_supercell_type
41 : : !!***
42 : :
43 : : !----------------------------------------------------------------------
44 : :
45 : : !!****t* m_ifc_supercell/ifc_supercell_type
46 : : !! NAME
47 : : !! ifc_supercell_type
48 : : !!
49 : : !! FUNCTION
50 : : !! Supercell defined by an integer matrix (S).
51 : : !!
52 : : !! A supercell holds N unitcell and N lattice vectors R, with
53 : : !! N = det(S).
54 : : !! Different supercell matrices are equivalent if they have the same Hermite normal form (H).
55 : : !!
56 : : !! The Hermite normal form allows to enumerate a set of canonical lattice vectors (Rcan),
57 : : !! which form a compact box, in reduced coordinates.
58 : : !!
59 : : !! Any lattice vector R can be decomposed as
60 : : !! R = (H * t) + Rcan
61 : : !! where t(3) is an integer vector. This decomposition defines the modulo operation
62 : : !! with respect to H (modH) as
63 : : !! Rcan = modH(R)
64 : : !! Two R vectors are equivalent if they correspond to the same canonical vector.
65 : : !!
66 : : !! A supercell shift (R_sh) is a lattice vector whose canonical coordinates are zero.
67 : : !!
68 : : !! Since the canonical lattice vectors form a compact box,
69 : : !! there exist a natural order to enumerate the lattice vectors,
70 : : !! and therefore there is a natural ordering for the atoms of the supercell.
71 : : !!
72 : : !! The supercell defines three different coordinate systems within a supercell:
73 : : !! cartesian coordinates, unitcell reduced coordinates, supercell reduced coordinates.
74 : : !!
75 : : !! GA: Should eventually replace m_dynmat/make_bigbox and m_supercell
76 : : !!
77 : : !! SOURCE
78 : :
79 : : type ifc_supercell_type
80 : :
81 : : integer :: natom
82 : : ! Number of atoms in the supercell.
83 : :
84 : : integer :: ntypat
85 : : ! Number of atom types.
86 : :
87 : : integer :: natom_unitcell
88 : : ! Number of atoms in the unitcell
89 : :
90 : : integer :: Ncells
91 : : ! Number of unit cells within the supercell (determinant of the matrix).
92 : :
93 : : integer :: N1, N2, N3
94 : : ! Limits for canonical coordinates. N1*N2*N3=Ncells
95 : :
96 : : integer :: matrix(3,3)
97 : : ! Defines the supercell. Allows for non-diagonal supercells.
98 : :
99 : : real(dp) :: matrixT(3,3)
100 : : ! Transposed supercell matrix.
101 : :
102 : : real(dp) :: matrixinvT(3,3)
103 : : ! Inverse and Transposed supercell matrix.
104 : :
105 : : integer :: HH(3,3)
106 : : ! Hermite normal form of the supercell matrix.
107 : :
108 : : real(dp) :: rprimd(3,3)
109 : : ! Primitive vectors of the unit cell, in Bohr. Columns are the lattice vectors.
110 : :
111 : : real(dp) :: rprimdt(3,3)
112 : : ! Transposed primitive vectors of the unit cell, in Bohr.
113 : :
114 : : real(dp) :: rprimd_inv(3,3)
115 : : ! Inverse primitive vectors of the unit cell.
116 : :
117 : : real(dp) :: rprimd_scaled(3,3)
118 : : real(dp) :: rprimdt_scaled(3,3)
119 : : ! Primitive vectors of the supercell, in Bohr.
120 : :
121 : : real(dp) :: rprimd_scaled_inv(3,3)
122 : : ! Inverse primitive vectors of the supercell.
123 : :
124 : :
125 : : real(dp) :: acell(3)
126 : : ! Cell parameters of the unit cell, in Bohr.
127 : : ! DEPRECATED
128 : : ! Avoid using acell internally.
129 : :
130 : : real(dp), allocatable :: znucl(:)
131 : : ! znucl(ntypat)
132 : : ! Atomic numbers of each atom type.
133 : :
134 : : real(dp), allocatable :: xred_unitcell(:,:)
135 : : ! xred_unitcell(3,natom_unitcell)
136 : : ! Reduced coordinates of atoms expressed in the unitcell coordinate system.
137 : :
138 : : real(dp), allocatable :: xcart_unit(:,:)
139 : : ! xcart_unit(3,natom_unitcell)
140 : : ! Cartesian coordinates of xred_unitcell, folded into ]-1/2,1/2] and expressed in Bohr.
141 : : ! Used by calc_ws_wghatm's minimum-image search, which only scans a few neighboring
142 : : ! supercells and therefore needs atoms already close to the origin.
143 : :
144 : : integer, allocatable :: typat_unitcell(:)
145 : : ! typat_unitcell(natom_unitcell)
146 : : ! Atom type in the unitcell.
147 : :
148 : : real(dp), allocatable :: R_red(:,:)
149 : : ! R_red(3, ncells)
150 : : ! Reduced coordinate of the lattice vectors of the unitcell
151 : : ! within the supercell.
152 : : ! This represents a particular choice for the supercell definition.
153 : :
154 : : integer, allocatable :: R_sh(:,:)
155 : : ! R_sh(3, natom)
156 : : ! An arbitrary supercell shift to apply to each atom.
157 : : ! Different atoms of the same unitcell may have different shifts
158 : : ! but the shift preserve canonical coordinates.
159 : :
160 : : real(dp), allocatable :: dij(:,:,:)
161 : : ! dij(natom, natom, 4)
162 : : ! Distance (dij(:,:,1)) and cartesian vector components (dij(:,:,2:4))
163 : : ! between every pair of atoms in the supercell,
164 : : ! using the minimum image convention.
165 : :
166 : : contains
167 : :
168 : : procedure :: init => supercell_init
169 : : ! Construct the object from the input definition.
170 : :
171 : : procedure :: free => supercell_free
172 : : ! Construct the object from the input definition.
173 : :
174 : : ! --------------------------
175 : :
176 : : procedure :: print => supercell_print
177 : : ! Write out supercell information.
178 : :
179 : : procedure :: write_xyz => supercell_write_xyz
180 : : ! Write out supercell in cartesian coordinates.
181 : :
182 : : procedure :: write_cif => supercell_write_cif
183 : : ! Write out supercell to a CIF file.
184 : :
185 : : ! --------------------------
186 : :
187 : : procedure :: modH => supercell_modH
188 : : ! Take the modulo with the superlattice.
189 : :
190 : : procedure :: is_mod0 => supercell_is_mod0
191 : : ! Is this vector modulo H equals to zero.
192 : :
193 : : procedure :: get_icell => supercell_get_icell
194 : : ! Get linear index of lattice vector.
195 : :
196 : : procedure :: get_Rcan => supercell_get_Rcan
197 : : !Get R_can from linear index.
198 : :
199 : : procedure :: get_iat_supercell => supercell_get_iat_supercell
200 : : ! Get supercell atom index from cell index and unitcell index.
201 : :
202 : : procedure :: get_iat_unitcell => supercell_get_iat_unitcell
203 : : ! Get cell index and unitcell index from supercell atom index.
204 : :
205 : : ! --------------------------
206 : :
207 : : procedure :: make_R_box => supercell_make_R_box
208 : : ! Create the list of lattice vectors.
209 : :
210 : : procedure :: get_R_at => supercell_get_R_at
211 : : ! Get the R vector corresponding to an atom.
212 : :
213 : : procedure :: set_R_at => supercell_set_R_at
214 : : ! Set the R vector of an atom.
215 : :
216 : : ! --------------------------
217 : : ! Coordinate transformations
218 : :
219 : : procedure :: cart_to_red => supercell_cart_to_red_unit
220 : : ! Transform from cartesian to reduced in the unitcell system of coordinate.
221 : :
222 : : procedure :: red_to_cart => supercell_red_unit_to_cart
223 : : ! Transform from reduced in the unitcell system of coordinate to cartesian.
224 : :
225 : : procedure :: cart_to_red_scaled => supercell_cart_to_red_scaled
226 : : ! Transform from cartesian to reduced in the supercell system of coordinate.
227 : :
228 : : procedure :: unit_to_scaled => supercell_unit_to_scaled
229 : : ! Transform from unit cell coordinates to scaled coordinates.
230 : :
231 : : procedure :: scaled_to_unit => supercell_scaled_to_unit
232 : : ! Transform from scaled coordinates to unit cell coordinates.
233 : :
234 : : procedure :: red_scaled_to_cart => supercell_red_scaled_to_cart
235 : : ! Transform from reduced in the supercell system of coordinate to cartesian.
236 : :
237 : : ! ------------------------
238 : :
239 : : procedure :: wrap_inbox => supercell_wrap_inbox
240 : : ! Wrap a vector in supercell-scaled reduced coordinates into [-0.5,0.5).
241 : :
242 : : procedure :: get_inbox_wrap => supercell_get_inbox_wrap
243 : : ! Find R shift that brings a R vector to [-0.5,0.5) in scaled coordinates.
244 : :
245 : : procedure :: make_box_centered => supercell_make_box_centered
246 : : ! Set all R_at so that scaled reduced coordinates fold into [-0.5,0.5).
247 : :
248 : : ! ------------------------
249 : :
250 : : procedure :: get_xred => supercell_get_xred_atom
251 : : ! Return the position of an atom in unitcell reduced coordinates.
252 : :
253 : : procedure :: get_xred_scaled => supercell_get_xred_scaled_atom
254 : : ! Return the position of an atom in supercell reduced coordinates.
255 : :
256 : : procedure :: get_xcart => supercell_get_xcart_atom
257 : : ! Return the position of an atom in cartesian coordinates.
258 : :
259 : : ! ------------------------
260 : :
261 : : procedure :: find_nearest_atom => supercell_find_nearest_atom
262 : : ! Get lattice vector from linear index.
263 : :
264 : : procedure :: atoms_distance => supercell_atoms_distance
265 : : ! Compute the distance between two atoms in the supercell.
266 : :
267 : : procedure :: calc_dij => supercell_calc_dij
268 : : ! Compute the distance and cartesian vector between all pairs of atoms
269 : : ! in the supercell, and store the result in SC%dij. Optional: call this
270 : : ! only if SC%dij is actually needed.
271 : :
272 : : procedure :: calc_wghatm => supercell_calc_wghatm
273 : : ! Compute the multiplicity weight of every (unitcell atom pair, canonical
274 : : ! lattice vector) triplet, generalizing the historical Bravais-lattice
275 : : ! dependent Wigner-Seitz weighting (wght9 in m_dynmat) to an arbitrary
276 : : ! supercell matrix.
277 : :
278 : : procedure :: calc_ws_wghatm => supercell_calc_ws_wghatm
279 : : ! Explicit list of the real-space lattice vectors lying in the Wigner-Seitz cell of the
280 : : ! supercell lattice for each atom pair, with 1/N weights for N tied images.
281 : :
282 : : procedure :: shortest_vector_norm => supercell_shortest_vector_norm
283 : : ! Norm of the shortest nonzero superlattice vector (generalizes the
284 : : ! "r_inscribed_sphere" quantity computed inline in wght9, m_dynmat).
285 : :
286 : : end type ifc_supercell_type
287 : : !!***
288 : :
289 : : contains
290 : : !!***
291 : :
292 : : !----------------------------------------------------------------------
293 : :
294 : : !----------------------------------------------------------------------
295 : :
296 : : !!****f* m_ifc_supercell/supercell_init
297 : : !! NAME
298 : : !! supercell_init
299 : : !!
300 : : !! FUNCTION
301 : : !! Initialize the supercell object from a Crystal object describing the
302 : : !! unitcell, and the supercell matrix.
303 : : !!
304 : : !! INPUTS
305 : : !! Crystal=Crystal object describing the unitcell.
306 : : !! multiplicity(3,3)=supercell matrix (rounded to the nearest integers)
307 : : !! acell(3)=optional, cell parameters of the unitcell, in Bohr. Defaults to [one,one,one]
308 : : !! (see the type's own note: acell is deprecated internally, kept only for callers that
309 : : !! still want to read it back from SC%acell).
310 : : !!
311 : : !! SIDE EFFECTS
312 : : !! SC=supercell object, allocated and filled on output
313 : : !!
314 : : !! SOURCE
315 : :
316 : 310 : subroutine supercell_init(SC,Crystal,multiplicity,acell)
317 : :
318 : : class(ifc_supercell_type),intent(inout) :: SC
319 : : type(crystal_t), intent(in) :: Crystal
320 : : real(dp), intent(in) :: multiplicity(3,3) ! TODO Should be integer
321 : : real(dp), intent(in), optional :: acell(3)
322 : :
323 : : integer :: natom_unitcell
324 : : integer :: iatom,ii,jj
325 : : real(dp) :: rmatrix(3,3)
326 : 620 : real(dp) :: xred_unit(3,Crystal%natom), shift_red(3,Crystal%natom)
327 : :
328 : 310 : natom_unitcell = Crystal%natom
329 : 310 : SC%natom_unitcell = natom_unitcell
330 : 310 : SC%ntypat = Crystal%ntypat
331 : :
332 [ - + + - : 1411 : ABI_CALLOC(SC%znucl, (SC%ntypat))
+ - + - -
+ - + - +
+ + ]
333 [ + + ]: 791 : SC%znucl(:) = Crystal%znucl(:)
334 : :
335 [ + + ]: 1240 : SC%acell = [one, one, one]
336 [ + + + + ]: 448 : if (present(acell)) SC%acell = acell
337 : :
338 [ + + + + ]: 4030 : SC%rprimd = Crystal%rprimd
339 [ + + + + : 7750 : SC%rprimdt = transpose(SC%rprimd)
+ + + + ]
340 : 310 : call matr3inv(SC%rprimdt, SC%rprimd_inv)
341 : :
342 [ - + + - : 4254 : ABI_CALLOC(SC%xred_unitcell,(3,natom_unitcell))
+ - + - -
+ - + - +
+ + + + ]
343 [ - + + - : 1761 : ABI_CALLOC(SC%typat_unitcell,(natom_unitcell))
+ - + - -
+ - + - +
+ + ]
344 [ + + + + ]: 3634 : SC%xred_unitcell(:,:) = Crystal%xred(:,:)
345 [ + + ]: 1141 : SC%typat_unitcell(:) = Crystal%typat(:)
346 : :
347 : : ! Cartesian positions folded into ]-1/2,1/2], needed by calc_ws_wghatm.
348 [ - + + - : 3944 : ABI_CALLOC(SC%xcart_unit,(3,natom_unitcell))
- + - + -
+ + + +
+ ]
349 [ + + + + ]: 3634 : call wrap2_pmhalf(Crystal%xred, xred_unit, shift_red)
350 [ + + ]: 1141 : do iatom=1,natom_unitcell
351 [ + + + + : 13606 : SC%xcart_unit(:,iatom) = matmul(SC%rprimd, xred_unit(:,iatom))
+ + + - ]
352 : : end do
353 : :
354 : : ! Round multiplicity to nearest integers
355 [ + + ]: 1240 : do ii=1,3
356 [ + + ]: 4030 : do jj=1,3
357 : 3720 : SC%matrix(ii, jj) = NINT(multiplicity(ii,jj))
358 : : end do
359 : : end do
360 : :
361 : : ! Compute supercell vectors
362 [ + + + + : 19840 : SC%rprimdt_scaled = matmul(SC%matrix, SC%rprimdt)
+ + + + +
+ + + +
+ ]
363 [ + + + + : 7750 : SC%rprimd_scaled = transpose(SC%rprimdt_scaled)
+ + + + ]
364 : 310 : call matr3inv(SC%rprimdt_scaled, SC%rprimd_scaled_inv)
365 : :
366 : : ! Compute number of primitive cells in the supercell
367 : 310 : call mati3det(SC%matrix, SC%Ncells)
368 : :
369 : : ! Compute inverse matrix
370 [ + + + + ]: 4030 : rmatrix = REAL(SC%matrix, KIND=dp)
371 : 310 : call matr3inv(rmatrix, SC%matrixinvT)
372 [ + + + + ]: 4030 : SC%matrixT = transpose(rmatrix)
373 : :
374 : : ! Compute Hermite normal form
375 : 310 : call mati3HNF(SC%matrix, SC%HH)
376 : :
377 : : ! Save limits
378 : 310 : SC%N1 = SC%HH(1,1); SC%N2 = SC%HH(2,2); SC%N3 = SC%HH(3,3);
379 : 310 : SC%natom = SC%Ncells * SC%natom_unitcell
380 : :
381 : : ! ------------------------------------
382 : :
383 : : ! Construct supercell box containing all different canonical vectors.
384 [ - + + - : 930 : ABI_MALLOC(SC%R_red, (3,SC%Ncells))
+ - + - -
+ - + -
+ ]
385 : 310 : call SC%make_R_box()
386 : :
387 : : ! Finally, set up R_sh so that atoms wrap in centered unit cell.
388 [ - + + - : 330590 : ABI_CALLOC(SC%R_sh, (3,SC%natom))
+ - + - -
+ - + - +
+ + + + ]
389 : 310 : call SC%make_box_centered()
390 : :
391 : 310 : end subroutine supercell_init
392 : : !!***
393 : :
394 : : !----------------------------------------------------------------------
395 : :
396 : : !!****f* m_ifc_supercell/supercell_free
397 : : !! NAME
398 : : !! supercell_free
399 : : !!
400 : : !! FUNCTION
401 : : !! Free the memory allocated in the supercell object.
402 : : !!
403 : : !! SIDE EFFECTS
404 : : !! SC=supercell object, deallocated on output
405 : : !!
406 : : !! SOURCE
407 : :
408 : 1563 : subroutine supercell_free(SC)
409 : :
410 : : class(ifc_supercell_type),intent(inout) :: SC
411 : :
412 [ + + ]: 1563 : ABI_SFREE(SC%xred_unitcell)
413 [ + + ]: 1563 : ABI_SFREE(SC%xcart_unit)
414 [ + + ]: 1563 : ABI_SFREE(SC%typat_unitcell)
415 [ + + ]: 1563 : ABI_SFREE(SC%znucl)
416 [ + + ]: 1563 : ABI_SFREE(SC%R_red)
417 [ + + ]: 1563 : ABI_SFREE(SC%R_sh)
418 [ + + ]: 1563 : ABI_SFREE(SC%dij)
419 : :
420 : 1563 : end subroutine supercell_free
421 : : !!***
422 : :
423 : : !----------------------------------------------------------------------
424 : :
425 : : !!****f* m_ifc_supercell/supercell_print
426 : : !! NAME
427 : : !! supercell_print
428 : : !!
429 : : !! FUNCTION
430 : : !! Write out supercell information to a file.
431 : : !!
432 : : !! INPUTS
433 : : !! iout=output unit number
434 : : !!
435 : : !! SOURCE
436 : :
437 : 46 : subroutine supercell_print(SC,iout)
438 : :
439 : : class(ifc_supercell_type),intent(in) :: SC
440 : : integer, intent(in) :: iout
441 : :
442 : : logical :: debug = .false.
443 : : integer :: ii,jj,icell
444 : : integer :: RR(3)
445 : : integer :: Rcan(3)
446 : :
447 : 46 : write(iout,'(a)') '===================================================='
448 : 46 : write(iout,'(a)') '--------------- Supercell definition ---------------'
449 : 46 : write(iout,'(a)') '----------------------------------------------------'
450 : 46 : write(iout,'(a)') 'Unit cell'
451 : 46 : write(iout,'(a)') '---------'
452 : 46 : write(iout,'(a)') 'Primitive vectors of the unit cell'
453 [ + + ]: 184 : do ii=1,3
454 [ + - + + ]: 598 : write(iout,'(24x,3(f15.10,1x))') (SC%rprimdt(ii,jj),jj=1,3)
455 : : end do
456 : 46 : write(iout,'(a)') 'Supercell vectors '
457 [ + + ]: 184 : do ii=1,3
458 [ + - + + ]: 598 : write(iout,'(24x,3(f15.10,1x))') (SC%rprimdt_scaled(ii,jj),jj=1,3)
459 : : end do
460 : : !write(iout,'(a)') 'natom_unitcell ',SC%natom_unitcell
461 : :
462 : 46 : write(iout,'(a)') 'Supercell matrix'
463 [ + + ]: 184 : do ii=1,3
464 : 184 : write(iout,'(24x,3i6)') SC%matrix(ii,:)
465 : : end do
466 : :
467 : 46 : write(iout,'(a)') 'Hermite Normal Form'
468 [ + + ]: 184 : do ii=1,3
469 : 184 : write(iout,'(24x,3i6)') SC%HH(ii,:)
470 : : end do
471 : :
472 : 46 : write(iout,'(a)') 'Number of cells'
473 : 46 : write(iout,'(24x,i6)') SC%Ncells
474 : :
475 : : !write(iout,'(a)') 'Lattice vectors in reduced coordinates w.r.t. primitive cell'
476 : : !do icell=1,SC%Ncells
477 : : ! RR(:) = SC%R_red(:,icell)
478 : : ! write(iout,'(8x,i8,8x,3i6)') icell, (RR(ii), ii=1,3)
479 : : !end do
480 : :
481 [ - + ]: 46 : if (debug) then
482 : :
483 : 0 : write(std_out,'(a)') '------------------------------'
484 : 0 : write(std_out,'(a)') 'Lattice vectors in reduced coordinates and canonical coordinates'
485 [ # # ]: 0 : do icell=1,SC%Ncells
486 : 0 : write(std_out, '(a)') '---------'
487 [ # # ]: 0 : RR(:) = SC%R_red(:,icell)
488 [ # # ]: 0 : Rcan = SC%modH(RR)
489 : 0 : write(std_out, *) 'icell =',icell
490 : 0 : write(std_out, *) 'R_red =',RR
491 : 0 : write(std_out, *) 'Rcan =',Rcan
492 : 0 : jj = SC%get_icell(RR)
493 : 0 : write(std_out, *) 'jcell =',jj
494 : 0 : write(std_out, '(a)') '---------'
495 : : end do
496 : 0 : write(std_out, '(a)') '------------------------------'
497 : : end if
498 : :
499 : 46 : write(iout,'(a)') '----------------------------------------------------'
500 : :
501 : 46 : end subroutine supercell_print
502 : : !!***
503 : :
504 : : !----------------------------------------------------------------------
505 : :
506 : : !!****f* m_ifc_supercell/supercell_make_R_box
507 : : !! NAME
508 : : !! supercell_make_R_box
509 : : !!
510 : : !! FUNCTION
511 : : !! Build the list of lattice vectors R_red forming a compact box
512 : : !! around the origin, one per canonical coordinate.
513 : : !!
514 : : !! SIDE EFFECTS
515 : : !! SC=supercell object, R_red is filled on output
516 : : !!
517 : : !! SOURCE
518 : :
519 : 310 : subroutine supercell_make_R_box(SC)
520 : :
521 : : class(ifc_supercell_type),intent(inout) :: SC
522 : :
523 : : integer :: icell
524 : : integer :: ii,jj,kk
525 : : integer :: di,dj,dk
526 : : integer :: Rred(3)
527 : : integer :: R_try(3)
528 : : integer :: T_try(3)
529 : : integer :: Ht(3,3)
530 : : real(dp) :: Rcan(3)
531 : : real(dp) :: Rcart(3), Rcan_cart(3)
532 : : real(dp) :: H1(3),H2(3),H3(3)
533 : : real(dp) :: R1(3),R2(3),R3(3)
534 : : real(dp) :: Rnorm, Rnorm_sm
535 : :
536 [ + + + + ]: 4030 : Ht = transpose(SC%HH)
537 [ + + ]: 1240 : H1 = REAL(SC%HH(1,:))
538 [ + + ]: 1240 : H2 = REAL(SC%HH(2,:))
539 [ + + ]: 1240 : H3 = REAL(SC%HH(3,:))
540 : 310 : R1 = SC%red_to_cart(H1)
541 : 310 : R2 = SC%red_to_cart(H2)
542 : 310 : R3 = SC%red_to_cart(H3)
543 : :
544 : :
545 [ + + + + ]: 130126 : SC%R_red(:,:) = zero
546 : :
547 : 310 : icell = 0
548 [ + + ]: 1699 : do kk=0,SC%N3-1
549 [ + + ]: 9323 : do jj=0,SC%N2-1
550 [ + + ]: 41467 : do ii=0,SC%N1-1
551 : :
552 : 32454 : icell = icell+1
553 [ + + ]: 32454 : if (icell == 1) cycle ! Origin
554 : :
555 [ + + ]: 128576 : Rcan(:) = (/ii, jj, kk/)
556 : :
557 [ + + - + ]: 128576 : Rcan_cart = SC%red_to_cart(real(Rcan,kind=dp))
558 : :
559 : 32144 : Rnorm_sm = greatest_real
560 : 32144 : Rred = zero
561 : :
562 [ + + ]: 128576 : do di = 1,-1,-1
563 [ + + ]: 417872 : do dj = 1,-1,-1
564 [ + + ]: 1253616 : do dk = 1,-1,-1
565 [ + + ]: 3471552 : T_try(:) = (/di, dj, dk/)
566 [ + + + + : 13886208 : R_try(:) = Rcan(:) + matmul(Ht, T_try)
+ + ]
567 [ + + - + : 6943104 : Rcart = SC%red_to_cart(real(R_try,kind=dp))
+ + ]
568 : :
569 : 867888 : Rnorm = l2norm(Rcart)
570 : :
571 : : ! The tol12 margin discriminates symmetry-equivalent R vectors,
572 [ + + ]: 1157184 : if (Rnorm < Rnorm_sm - tol12) then
573 : 245959 : Rnorm_sm = Rnorm
574 : 245959 : Rred = R_try
575 : : end if
576 : :
577 : : end do
578 : : end do
579 : : end do
580 : :
581 [ + + ]: 136200 : SC%R_red(:,icell) = Rred(:)
582 : :
583 : : end do
584 : : end do
585 : : end do
586 : :
587 : 310 : end subroutine supercell_make_R_box
588 : : !!***
589 : :
590 : : !----------------------------------------------------------------------
591 : :
592 : : !!****f* m_ifc_supercell/supercell_write_xyz
593 : : !! NAME
594 : : !! supercell_write_xyz
595 : : !!
596 : : !! FUNCTION
597 : : !! Write the ideal supercell atomic positions to an XYZ file, in cartesian coordinates.
598 : : !!
599 : : !! INPUTS
600 : : !! prefix=prefix of the output file name
601 : : !!
602 : : !! SOURCE
603 : :
604 : 0 : subroutine supercell_write_xyz(SC,prefix)
605 : : class(ifc_supercell_type),intent(in) :: SC
606 : : character(len=fnlen), intent(in) :: prefix
607 : :
608 : : integer :: unt
609 : : integer :: iatom,icell,iat_unit
610 : : real(dp) :: xcart(3)
611 : :
612 : 0 : unt = get_unit()
613 : 0 : open(unit=unt,file=trim(prefix)//'_xcart_ideal.xyz')
614 : 0 : write(unt,'(a1,1x,a8,4x,a6,4x,3(a5,30x))') '#', 'xcart', 'iatom', 'icell', 'iat_unitcell'
615 : 0 : write(unt,'(a)')''
616 : :
617 [ # # ]: 0 : do iatom=1,SC%natom
618 : :
619 [ # # ]: 0 : xcart = SC%get_xcart(iatom)
620 : 0 : call SC%get_iat_unitcell(iatom, icell, iat_unit)
621 : :
622 : 0 : write(unt,'(2x,3(f10.6,1x),4x,3(i6,2x))') xcart, iatom, icell, iat_unit
623 : 0 : write(unt,'(a)')''
624 : :
625 : : end do
626 : :
627 : 0 : close(unt)
628 : :
629 : 0 : end subroutine supercell_write_xyz
630 : : !!***
631 : :
632 : : !----------------------------------------------------------------------
633 : :
634 : : !!****f* m_ifc_supercell/supercell_write_cif
635 : : !! NAME
636 : : !! supercell_write_cif
637 : : !!
638 : : !! FUNCTION
639 : : !! Write the ideal supercell (including R_sh shifts) to a CIF file.
640 : : !! znucl(ntypat) is used to resolve atomic numbers into element symbols.
641 : : !!
642 : : !! INPUTS
643 : : !! prefix=prefix of the output file name
644 : : !!
645 : : !! SOURCE
646 : :
647 : 46 : subroutine supercell_write_cif(SC,prefix)
648 : : class(ifc_supercell_type),intent(in) :: SC
649 : : character(len=fnlen), intent(in) :: prefix
650 : :
651 : : integer :: unt
652 : : integer :: iatom,icell,iat_unit
653 : : real(dp) :: aa,bb,cc,alpha,beta,gamma_,ZZ
654 : : real(dp) :: xred_unit(3),xred_scaled(3)
655 : : character(len=2) :: symbol
656 : : character(len=16) :: label
657 : :
658 : 46 : aa = l2norm(SC%rprimdt_scaled(1,:))
659 : 46 : bb = l2norm(SC%rprimdt_scaled(2,:))
660 : 46 : cc = l2norm(SC%rprimdt_scaled(3,:))
661 : :
662 [ + + ]: 184 : alpha = ACOS(dot_product(SC%rprimdt_scaled(2,:),SC%rprimdt_scaled(3,:)) / (bb*cc)) * 180.0_dp / pi
663 [ + + ]: 184 : beta = ACOS(dot_product(SC%rprimdt_scaled(1,:),SC%rprimdt_scaled(3,:)) / (aa*cc)) * 180.0_dp / pi
664 [ + + ]: 184 : gamma_= ACOS(dot_product(SC%rprimdt_scaled(1,:),SC%rprimdt_scaled(2,:)) / (aa*bb)) * 180.0_dp / pi
665 : :
666 : 46 : unt = get_unit()
667 : 46 : open(unit=unt,file=trim(prefix)//'_supercell.cif')
668 : :
669 : 46 : write(unt,'(a)') 'data_supercell'
670 : 46 : write(unt,'(a,f14.8)') '_cell_length_a ', aa*Bohr_Ang
671 : 46 : write(unt,'(a,f14.8)') '_cell_length_b ', bb*Bohr_Ang
672 : 46 : write(unt,'(a,f14.8)') '_cell_length_c ', cc*Bohr_Ang
673 : 46 : write(unt,'(a,f14.8)') '_cell_angle_alpha ', alpha
674 : 46 : write(unt,'(a,f14.8)') '_cell_angle_beta ', beta
675 : 46 : write(unt,'(a,f14.8)') '_cell_angle_gamma ', gamma_
676 : 46 : write(unt,'(a)') "_symmetry_space_group_name_H-M 'P 1'"
677 : 46 : write(unt,'(a)') '_symmetry_Int_Tables_number 1'
678 : 46 : write(unt,'(a)') 'loop_'
679 : 46 : write(unt,'(a)') '_atom_site_label'
680 : 46 : write(unt,'(a)') '_atom_site_type_symbol'
681 : 46 : write(unt,'(a)') '_atom_site_fract_x'
682 : 46 : write(unt,'(a)') '_atom_site_fract_y'
683 : 46 : write(unt,'(a)') '_atom_site_fract_z'
684 : :
685 [ + + ]: 7324 : do iatom=1,SC%natom
686 : :
687 : 7278 : call SC%get_iat_unitcell(iatom, icell, iat_unit)
688 : :
689 : 7278 : ZZ = SC%znucl(SC%typat_unitcell(iat_unit))
690 : 7278 : symbol = adjustl(znucl2symbol(ZZ))
691 : :
692 : 7278 : write(label,'(a,i0)') trim(symbol), iatom
693 : :
694 [ + + ]: 29112 : xred_unit = SC%get_xred(iatom)
695 [ + + ]: 29112 : xred_scaled = SC%unit_to_scaled(xred_unit)
696 : :
697 : 14602 : write(unt,'(a8,2x,a2,2x,3(f12.8,1x))') label, symbol, xred_scaled
698 : :
699 : : end do
700 : :
701 : 46 : close(unt)
702 : :
703 : 46 : end subroutine supercell_write_cif
704 : : !!***
705 : :
706 : : !----------------------------------------------------------------------
707 : :
708 : : !!****f* m_ifc_supercell/supercell_modH
709 : : !! NAME
710 : : !! supercell_modH
711 : : !!
712 : : !! FUNCTION
713 : : !! Take the modulo of a lattice vector RR with respect to the supercell
714 : : !! matrix, achieving the unique decomposition RR = (H * t) + Rcan,
715 : : !! where t(3) is an integer vector.
716 : : !!
717 : : !! INPUTS
718 : : !! RR(3)=lattice vector, in reduced coordinates
719 : : !!
720 : : !! OUTPUT
721 : : !! Rcan(3)=canonical lattice vector, RR modulo H
722 : : !!
723 : : !! SOURCE
724 : :
725 [ - + ]: 96971 : function supercell_modH(SC,RR) result(Rcan)
726 : :
727 : : class(ifc_supercell_type),intent(in) :: SC
728 : : integer, intent(in) :: RR(3)
729 : : integer :: Rcan(3)
730 : :
731 : : integer :: tt(3)
732 : : integer :: H11,H22,H33,H12,H13,H23
733 : :
734 : 96971 : H11=SC%HH(1,1); H22=SC%HH(2,2); H33=SC%HH(3,3);
735 : 96971 : H12=SC%HH(1,2); H13=SC%HH(1,3); H23=SC%HH(2,3);
736 : :
737 : 96971 : Rcan(1) = modulo(RR(1), H11)
738 : 96971 : tt(1) = (RR(1) - Rcan(1)) / H11
739 : :
740 : 96971 : Rcan(2) = modulo(RR(2) - H12*tt(1), H22)
741 : 96971 : tt(2) = (RR(2) - H12*tt(1) - Rcan(2)) / H22
742 : :
743 : 96971 : Rcan(3) = modulo(RR(3) - H13*tt(1) - H23*tt(2), H33)
744 : 96971 : tt(3) = (RR(3) - H13*tt(1) - H23*tt(2) - Rcan(3)) / H33
745 : :
746 : 96971 : end function supercell_modH
747 : : !!***
748 : :
749 : : !----------------------------------------------------------------------
750 : :
751 : : !!****f* m_ifc_supercell/supercell_is_mod0
752 : : !! NAME
753 : : !! supercell_is_mod0
754 : : !!
755 : : !! FUNCTION
756 : : !! Check whether a lattice vector is equivalent to zero modulo the supercell matrix.
757 : : !!
758 : : !! INPUTS
759 : : !! RR(3)=lattice vector, in reduced coordinates
760 : : !!
761 : : !! OUTPUT
762 : : !! is_mod0=.true. if modH(RR) is zero
763 : : !!
764 : : !! SOURCE
765 : :
766 : 82415 : function supercell_is_mod0(SC,RR) result(is_mod0)
767 : : class(ifc_supercell_type),intent(in) :: SC
768 : : integer,intent(in) :: RR(3)
769 : : logical :: is_mod0
770 : : integer :: ii
771 : : integer :: Rmod(3)
772 : :
773 [ + + ]: 329660 : Rmod = SC%modH(RR)
774 : 82415 : is_mod0 = .true.
775 [ + + ]: 329660 : do ii=1,3
776 [ - + ]: 329660 : if (Rmod(ii).ne.0) is_mod0 =.false.
777 : : end do
778 : :
779 : 82415 : end function supercell_is_mod0
780 : : !!***
781 : :
782 : : !----------------------------------------------------------------------
783 : :
784 : : !!****f* m_ifc_supercell/supercell_get_icell
785 : : !! NAME
786 : : !! supercell_get_icell
787 : : !!
788 : : !! FUNCTION
789 : : !! Get the linear cell index corresponding to a lattice vector.
790 : : !!
791 : : !! INPUTS
792 : : !! RR(3)=lattice vector, in reduced coordinates
793 : : !!
794 : : !! OUTPUT
795 : : !! icell=linear index of the cell containing RR
796 : : !!
797 : : !! SOURCE
798 : :
799 : 7278 : function supercell_get_icell(SC,RR) result(icell)
800 : : class(ifc_supercell_type),intent(in) :: SC
801 : : integer, intent(in) :: RR(3)
802 : : integer :: icell
803 : : integer :: Rcan(3)
804 [ + + ]: 29112 : Rcan = SC%modH(RR)
805 : 7278 : icell = 1 + Rcan(1) + Rcan(2) * SC%N1 + Rcan(3) * SC%N2 * SC%N1
806 : 7278 : end function supercell_get_icell
807 : : !!***
808 : :
809 : : !----------------------------------------------------------------------
810 : :
811 : : !!****f* m_ifc_supercell/supercell_get_Rcan
812 : : !! NAME
813 : : !! supercell_get_Rcan
814 : : !!
815 : : !! FUNCTION
816 : : !! Get the canonical lattice vector corresponding to a linear cell index.
817 : : !!
818 : : !! INPUTS
819 : : !! icell=linear index of the cell
820 : : !!
821 : : !! OUTPUT
822 : : !! Rcan(3)=canonical lattice vector
823 : : !!
824 : : !! SOURCE
825 : :
826 [ # # ]: 0 : function supercell_get_Rcan(SC,icell) result(Rcan)
827 : : class(ifc_supercell_type),intent(in) :: SC
828 : : integer, intent(in) :: icell
829 : : integer :: Rcan(3)
830 : : integer :: ii,jj,kk
831 : : integer :: idx
832 : 0 : idx = icell - 1
833 : 0 : ii = mod(idx, SC%N1)
834 : 0 : jj = mod((idx-ii)/SC%N1, SC%N2)
835 : 0 : kk = ((idx-ii)/SC%N1 - jj) / SC%N2 ! i.e. remainder after removing ii and jj
836 [ # # ]: 0 : Rcan(:) = (/ii, jj, kk/)
837 : 0 : end function supercell_get_Rcan
838 : : !!***
839 : :
840 : : !----------------------------------------------------------------------
841 : :
842 : : !!****f* m_ifc_supercell/supercell_get_R_at
843 : : !! NAME
844 : : !! supercell_get_R_at
845 : : !!
846 : : !! FUNCTION
847 : : !! Get the lattice vector R associated with an atom of the supercell.
848 : : !!
849 : : !! INPUTS
850 : : !! iat_super=index of the atom in the supercell
851 : : !!
852 : : !! OUTPUT
853 : : !! R_at(3)=lattice vector of the atom
854 : : !!
855 : : !! SOURCE
856 : :
857 [ - + ]: 89693 : function supercell_get_R_at(SC, iat_super) result(R_at)
858 : : class(ifc_supercell_type),intent(in) :: SC
859 : : integer, intent(in) :: iat_super
860 : : integer :: icell, iat_unit
861 : : integer :: R_at(3)
862 : 89693 : call SC%get_iat_unitcell(iat_super, icell, iat_unit)
863 [ + + ]: 358772 : R_at = SC%R_red(:,icell) + SC%R_sh(:,iat_super)
864 : 89693 : end function supercell_get_R_at
865 : : !!***
866 : :
867 : : !----------------------------------------------------------------------
868 : :
869 : : !!****f* m_ifc_supercell/supercell_get_iat_supercell
870 : : !! NAME
871 : : !! supercell_get_iat_supercell
872 : : !!
873 : : !! FUNCTION
874 : : !! Get the supercell atom index from the cell index and the unitcell atom index.
875 : : !!
876 : : !! INPUTS
877 : : !! icell=linear index of the cell
878 : : !! iat_unit=index of the atom in the unit cell
879 : : !!
880 : : !! OUTPUT
881 : : !! iat_super=index of the atom in the supercell
882 : : !!
883 : : !! SOURCE
884 : :
885 : 30280 : function supercell_get_iat_supercell(SC,icell,iat_unit) result(iat_super)
886 : : class(ifc_supercell_type),intent(in) :: SC
887 : : integer, intent(in) :: icell
888 : : integer, intent(in) :: iat_unit
889 : : integer :: iat_super
890 : :
891 : 30280 : iat_super = iat_unit + (icell - 1) * SC%natom_unitcell
892 : :
893 : 30280 : end function supercell_get_iat_supercell
894 : : !!***
895 : :
896 : : !----------------------------------------------------------------------
897 : :
898 : : !!****f* m_ifc_supercell/supercell_get_iat_unitcell
899 : : !! NAME
900 : : !! supercell_get_iat_unitcell
901 : : !!
902 : : !! FUNCTION
903 : : !! Get the cell index and the unitcell atom index from the supercell atom index.
904 : : !!
905 : : !! INPUTS
906 : : !! iat_super=index of the atom in the supercell
907 : : !!
908 : : !! OUTPUT
909 : : !! icell=linear index of the cell
910 : : !! iat_unit=index of the atom in the unit cell
911 : : !!
912 : : !! SOURCE
913 : :
914 : 116748567 : subroutine supercell_get_iat_unitcell(SC,iat_super,icell,iat_unit)
915 : : class(ifc_supercell_type),intent(in) :: SC
916 : : integer, intent(in) :: iat_super
917 : : integer, intent(out) :: icell
918 : : integer, intent(out) :: iat_unit
919 : :
920 : 116748567 : iat_unit = 1 + modulo(iat_super - 1, SC%natom_unitcell)
921 : 116748567 : icell = 1 + (iat_super - iat_unit) / SC%natom_unitcell
922 : :
923 : 116748567 : end subroutine supercell_get_iat_unitcell
924 : : !!***
925 : :
926 : : !----------------------------------------------------------------------
927 : :
928 : : !!****f* m_ifc_supercell/supercell_red_unit_to_cart
929 : : !! NAME
930 : : !! supercell_red_unit_to_cart
931 : : !!
932 : : !! FUNCTION
933 : : !! Transform from unitcell reduced coordinates to cartesian coordinates.
934 : : !!
935 : : !! INPUTS
936 : : !! xred(3)=position, in unitcell reduced coordinates
937 : : !!
938 : : !! OUTPUT
939 : : !! xcart(3)=position, in cartesian coordinates
940 : : !!
941 : : !! SOURCE
942 : :
943 [ - + ]: 1311998 : function supercell_red_unit_to_cart(SC,xred) result(xcart)
944 : : class(ifc_supercell_type),intent(in) :: SC
945 : : real(dp), intent(in) :: xred(3)
946 : : real(dp) :: xcart(3)
947 [ + + + + : 22303966 : xcart = matmul(SC%rprimd, xred)
+ + ]
948 : 1311998 : end function supercell_red_unit_to_cart
949 : : !!***
950 : :
951 : : !----------------------------------------------------------------------
952 : :
953 : : !!****f* m_ifc_supercell/supercell_cart_to_red_unit
954 : : !! NAME
955 : : !! supercell_cart_to_red_unit
956 : : !!
957 : : !! FUNCTION
958 : : !! Transform from cartesian coordinates to unitcell reduced coordinates.
959 : : !!
960 : : !! INPUTS
961 : : !! xcart(3)=position, in cartesian coordinates
962 : : !!
963 : : !! OUTPUT
964 : : !! xred(3)=position, in unitcell reduced coordinates
965 : : !!
966 : : !! SOURCE
967 : :
968 [ - + ]: 7278 : function supercell_cart_to_red_unit(SC,xcart) result(xred)
969 : : class(ifc_supercell_type),intent(in) :: SC
970 : : real(dp), intent(in) :: xcart(3)
971 : : real(dp) :: xred(3)
972 [ + + + + : 123726 : xred = matmul(SC%rprimd_inv, xcart)
+ + ]
973 : 7278 : end function supercell_cart_to_red_unit
974 : : !!***
975 : :
976 : : !----------------------------------------------------------------------
977 : :
978 : : !!****f* m_ifc_supercell/supercell_cart_to_red_scaled
979 : : !! NAME
980 : : !! supercell_cart_to_red_scaled
981 : : !!
982 : : !! FUNCTION
983 : : !! Transform from cartesian coordinates to supercell reduced coordinates.
984 : : !!
985 : : !! INPUTS
986 : : !! xcart(3)=position, in cartesian coordinates
987 : : !!
988 : : !! OUTPUT
989 : : !! xred(3)=position, in supercell reduced coordinates
990 : : !!
991 : : !! SOURCE
992 : :
993 [ # # ]: 0 : function supercell_cart_to_red_scaled(SC,xcart) result(xred)
994 : : class(ifc_supercell_type),intent(in) :: SC
995 : : real(dp), intent(in) :: xcart(3)
996 : : real(dp) :: xred(3)
997 [ # # # # : 0 : xred = matmul(SC%rprimd_scaled_inv, xcart)
# # ]
998 : 0 : end function supercell_cart_to_red_scaled
999 : : !!***
1000 : :
1001 : : !----------------------------------------------------------------------
1002 : :
1003 : : !!****f* m_ifc_supercell/supercell_unit_to_scaled
1004 : : !! NAME
1005 : : !! supercell_unit_to_scaled
1006 : : !!
1007 : : !! FUNCTION
1008 : : !! Transform from unitcell reduced coordinates to supercell reduced coordinates.
1009 : : !!
1010 : : !! INPUTS
1011 : : !! xred_unit(3)=position, in unitcell reduced coordinates
1012 : : !!
1013 : : !! OUTPUT
1014 : : !! xred_scaled(3)=position, in supercell reduced coordinates
1015 : : !!
1016 : : !! SOURCE
1017 : :
1018 [ - + ]: 116563191 : function supercell_unit_to_scaled(SC,xred_unit) result(xred_scaled)
1019 : : class(ifc_supercell_type),intent(in) :: SC
1020 : : real(dp), intent(in) :: xred_unit(3)
1021 : : real(dp) :: xred_scaled(3)
1022 [ + + + + : 1981574247 : xred_scaled = matmul(SC%matrixinvT, xred_unit)
+ + ]
1023 : 116563191 : end function supercell_unit_to_scaled
1024 : : !!***
1025 : :
1026 : : !----------------------------------------------------------------------
1027 : :
1028 : : !!****f* m_ifc_supercell/supercell_scaled_to_unit
1029 : : !! NAME
1030 : : !! supercell_scaled_to_unit
1031 : : !!
1032 : : !! FUNCTION
1033 : : !! Transform from supercell reduced coordinates to unitcell reduced coordinates.
1034 : : !!
1035 : : !! INPUTS
1036 : : !! xred_scaled(3)=position, in supercell reduced coordinates
1037 : : !!
1038 : : !! OUTPUT
1039 : : !! xred_unit(3)=position, in unitcell reduced coordinates
1040 : : !!
1041 : : !! SOURCE
1042 : :
1043 [ - + ]: 93687906 : function supercell_scaled_to_unit(SC,xred_scaled) result(xred_unit)
1044 : : class(ifc_supercell_type),intent(in) :: SC
1045 : : real(dp), intent(in) :: xred_scaled(3)
1046 : : real(dp) :: xred_unit(3)
1047 [ + + + + : 1592694402 : xred_unit = matmul(SC%matrixT, xred_scaled)
+ + ]
1048 : 93687906 : end function supercell_scaled_to_unit
1049 : : !!***
1050 : :
1051 : : !----------------------------------------------------------------------
1052 : :
1053 : : !!****f* m_ifc_supercell/supercell_red_scaled_to_cart
1054 : : !! NAME
1055 : : !! supercell_red_scaled_to_cart
1056 : : !!
1057 : : !! FUNCTION
1058 : : !! Transform from supercell reduced coordinates to cartesian coordinates.
1059 : : !!
1060 : : !! INPUTS
1061 : : !! xred(3)=position, in supercell reduced coordinates
1062 : : !!
1063 : : !! OUTPUT
1064 : : !! xcart(3)=position, in cartesian coordinates
1065 : : !!
1066 : : !! SOURCE
1067 : :
1068 [ - + ]: 1728670 : function supercell_red_scaled_to_cart(SC,xred) result(xcart)
1069 : : class(ifc_supercell_type),intent(in) :: SC
1070 : : real(dp), intent(in) :: xred(3)
1071 : : real(dp) :: xcart(3)
1072 [ + + + + : 29387390 : xcart = matmul(SC%rprimd_scaled, xred)
+ + ]
1073 : 1728670 : end function supercell_red_scaled_to_cart
1074 : : !!***
1075 : :
1076 : : !----------------------------------------------------------------------
1077 : :
1078 : : !!****f* m_ifc_supercell/supercell_wrap_inbox
1079 : : !! NAME
1080 : : !! supercell_wrap_inbox
1081 : : !!
1082 : : !! FUNCTION
1083 : : !! Wrap a vector expressed in supercell reduced coordinates into [-0.5,0.5),
1084 : : !! component-wise. Mirrors ifc_make_inbox; matches the minimum-image
1085 : : !! convention used in m_ifc.
1086 : : !!
1087 : : !! INPUTS
1088 : : !! xred_scaled(3)=position, in supercell reduced coordinates
1089 : : !!
1090 : : !! OUTPUT
1091 : : !! xred_wrapped(3)=position, wrapped into [-0.5,0.5)
1092 : : !!
1093 : : !! SOURCE
1094 : :
1095 [ - + ]: 93674278 : function supercell_wrap_inbox(SC,xred_scaled) result(xred_wrapped)
1096 : : class(ifc_supercell_type),intent(in) :: SC
1097 : : real(dp), intent(in) :: xred_scaled(3)
1098 : : real(dp) :: xred_wrapped(3)
1099 : :
1100 : : integer :: ii,jj
1101 : :
1102 : 93674278 : ABI_UNUSED_A(SC)
1103 [ + + ]: 374697112 : xred_wrapped = xred_scaled
1104 [ + + ]: 374697112 : do ii=1,3
1105 [ + + ]: 281022834 : if (xred_wrapped(ii) < zero) then
1106 : 85649960 : jj = int(xred_wrapped(ii) - half + tol8)
1107 : : else
1108 : 195372874 : jj = int(xred_wrapped(ii) + half + tol8)
1109 : : end if
1110 : 374697112 : xred_wrapped(ii) = xred_wrapped(ii) - real(jj, kind=dp)
1111 : : end do
1112 : :
1113 : 93674278 : end function supercell_wrap_inbox
1114 : : !!***
1115 : :
1116 : : !----------------------------------------------------------------------
1117 : :
1118 : : !!****f* m_ifc_supercell/supercell_get_inbox_wrap
1119 : : !! NAME
1120 : : !! supercell_get_inbox_wrap
1121 : : !!
1122 : : !! FUNCTION
1123 : : !! Find the lattice vector shift R_shift such that R_at + R_shift lies
1124 : : !! inside the [-0.5,0.5) box, in supercell reduced coordinates.
1125 : : !!
1126 : : !! INPUTS
1127 : : !! iat_super=index of the atom in the supercell
1128 : : !!
1129 : : !! OUTPUT
1130 : : !! R_shift(3)=lattice vector shift, in unitcell reduced coordinates
1131 : : !! R_shift_scaled(3)=lattice vector shift, in supercell reduced coordinates
1132 : : !!
1133 : : !! SOURCE
1134 : :
1135 : 82415 : subroutine supercell_get_inbox_wrap(SC,iat_super,R_shift,R_shift_scaled)
1136 : :
1137 : : class(ifc_supercell_type),intent(in) :: SC
1138 : : integer, intent(in) :: iat_super
1139 : : integer, intent(out) :: R_shift(3)
1140 : : real(dp), intent(out) :: R_shift_scaled(3)
1141 : : real(dp) :: xred_unit(3), xred_scaled(3), xred_wrap(3)
1142 : :
1143 : : ! Fold into [-0.5,0.5) and determine the (integer) shift that was added..
1144 [ + + ]: 329660 : xred_unit = SC%get_xred(iat_super)
1145 [ + + ]: 329660 : xred_scaled = SC%unit_to_scaled(xred_unit)
1146 [ + + ]: 329660 : xred_wrap = SC%wrap_inbox(xred_scaled)
1147 : :
1148 [ + + ]: 329660 : R_shift_scaled = xred_wrap - xred_scaled
1149 [ + + ]: 329660 : R_shift = NINT(SC%scaled_to_unit(R_shift_scaled))
1150 : :
1151 : 82415 : end subroutine supercell_get_inbox_wrap
1152 : : !!***
1153 : :
1154 : : !----------------------------------------------------------------------
1155 : :
1156 : : !!****f* m_ifc_supercell/supercell_make_box_centered
1157 : : !! NAME
1158 : : !! supercell_make_box_centered
1159 : : !!
1160 : : !! FUNCTION
1161 : : !! Set the R_sh shift of every atom so that its supercell reduced
1162 : : !! coordinates fold into [-0.5,0.5).
1163 : : !!
1164 : : !! SIDE EFFECTS
1165 : : !! SC=supercell object, R_sh is updated on output
1166 : : !!
1167 : : !! SOURCE
1168 : :
1169 : 310 : subroutine supercell_make_box_centered(SC)
1170 : :
1171 : : class(ifc_supercell_type),intent(inout) :: SC
1172 : : integer :: iat
1173 : : integer :: R_at(3), R_shift(3)
1174 : : real(dp) :: R_shift_scaled(3)
1175 : :
1176 [ + + ]: 82725 : do iat=1,SC%natom
1177 : 82415 : call SC%get_inbox_wrap(iat, R_shift, R_shift_scaled)
1178 [ + + ]: 329660 : R_at = SC%get_R_at(iat)
1179 [ + + ]: 329970 : call SC%set_R_at(iat, R_at+R_shift)
1180 : : end do
1181 : :
1182 : 310 : end subroutine supercell_make_box_centered
1183 : : !!***
1184 : :
1185 : : !----------------------------------------------------------------------
1186 : :
1187 : : !!****f* m_ifc_supercell/supercell_get_xred_atom
1188 : : !! NAME
1189 : : !! supercell_get_xred_atom
1190 : : !!
1191 : : !! FUNCTION
1192 : : !! Get the position of an atom, in unitcell reduced coordinates.
1193 : : !!
1194 : : !! INPUTS
1195 : : !! iat_super=index of the atom in the supercell
1196 : : !!
1197 : : !! OUTPUT
1198 : : !! xred(3)=position, in unitcell reduced coordinates
1199 : : !!
1200 : : !! SOURCE
1201 : :
1202 [ - + ]: 116561857 : function supercell_get_xred_atom(SC,iat_super) result(xred)
1203 : : class(ifc_supercell_type),intent(in) :: SC
1204 : : integer, intent(in) :: iat_super
1205 : : integer :: icell, iat_unit
1206 : : real(dp) :: xred(3)
1207 : :
1208 : 116561857 : call SC%get_iat_unitcell(iat_super, icell, iat_unit)
1209 [ + + ]: 466247428 : xred(:) = SC%R_red(:,icell) + SC%R_sh(:,iat_super) + SC%xred_unitcell(:,iat_unit)
1210 : :
1211 : 116561857 : end function supercell_get_xred_atom
1212 : : !!***
1213 : :
1214 : : !----------------------------------------------------------------------
1215 : :
1216 : : !!****f* m_ifc_supercell/supercell_get_xred_scaled_atom
1217 : : !! NAME
1218 : : !! supercell_get_xred_scaled_atom
1219 : : !!
1220 : : !! FUNCTION
1221 : : !! Get the position of an atom, in supercell reduced coordinates.
1222 : : !! If wrap is given and true, the result is wrapped into the centered
1223 : : !! box [-0.5,0.5).
1224 : : !!
1225 : : !! INPUTS
1226 : : !! iat_super=index of the atom in the supercell
1227 : : !! wrap=optional, if true, wrap the result into [-0.5,0.5). Defaults to
1228 : : !! false (no wrapping) when absent.
1229 : : !!
1230 : : !! OUTPUT
1231 : : !! xred_scaled(3)=position, in supercell reduced coordinates
1232 : : !!
1233 : : !! SOURCE
1234 : :
1235 [ - + ]: 113340398 : function supercell_get_xred_scaled_atom(SC,iat_super,wrap) result(xred_scaled)
1236 : : class(ifc_supercell_type),intent(in) :: SC
1237 : : integer, intent(in) :: iat_super
1238 : : logical, intent(in), optional :: wrap
1239 : : real(dp) :: xred_scaled(3)
1240 : : real(dp) :: xred(3)
1241 : :
1242 [ + + ]: 453361592 : xred = SC%get_xred(iat_super)
1243 [ + + ]: 453361592 : xred_scaled = SC%unit_to_scaled(xred)
1244 : :
1245 [ - + ]: 113340398 : if (present(wrap)) then
1246 [ # # # # ]: 0 : if (wrap) xred_scaled = SC%wrap_inbox(xred_scaled)
1247 : : end if
1248 : :
1249 : 113340398 : end function supercell_get_xred_scaled_atom
1250 : : !!***
1251 : :
1252 : : !----------------------------------------------------------------------
1253 : :
1254 : : !!****f* m_ifc_supercell/supercell_get_xcart_atom
1255 : : !! NAME
1256 : : !! supercell_get_xcart_atom
1257 : : !!
1258 : : !! FUNCTION
1259 : : !! Get the position of an atom, in cartesian coordinates.
1260 : : !!
1261 : : !! INPUTS
1262 : : !! iat_super=index of the atom in the supercell
1263 : : !!
1264 : : !! OUTPUT
1265 : : !! xcart(3)=position, in cartesian coordinates
1266 : : !!
1267 : : !! SOURCE
1268 : :
1269 [ - + ]: 6036 : function supercell_get_xcart_atom(SC,iat_super) result(xcart)
1270 : : class(ifc_supercell_type),intent(in) :: SC
1271 : : integer, intent(in) :: iat_super
1272 : : real(dp) :: xcart(3)
1273 : : real(dp) :: xred(3)
1274 : :
1275 [ + + ]: 24144 : xred = SC%get_xred(iat_super)
1276 [ + + ]: 24144 : xcart = SC%red_to_cart(xred)
1277 : :
1278 : 6036 : end function supercell_get_xcart_atom
1279 : : !!***
1280 : :
1281 : : !----------------------------------------------------------------------
1282 : :
1283 : : !!****f* m_ifc_supercell/supercell_set_R_at
1284 : : !! NAME
1285 : : !! supercell_set_R_at
1286 : : !!
1287 : : !! FUNCTION
1288 : : !! Set the lattice vector R of an atom, by adjusting its R_sh shift.
1289 : : !!
1290 : : !! INPUTS
1291 : : !! iat_super=index of the atom in the supercell
1292 : : !! R_at(3)=lattice vector to assign to the atom
1293 : : !!
1294 : : !! SIDE EFFECTS
1295 : : !! SC=supercell object, R_sh is updated on output
1296 : : !!
1297 : : !! SOURCE
1298 : :
1299 : 82415 : subroutine supercell_set_R_at(SC,iat_super,R_at)
1300 : : class(ifc_supercell_type),intent(inout) :: SC
1301 : : integer, intent(in) :: iat_super
1302 : : integer, intent(in) :: R_at(3)
1303 : : integer :: icell, iat_unit
1304 : : integer :: R_red(3)
1305 : : integer :: R_sh(3)
1306 : :
1307 : 82415 : call SC%get_iat_unitcell(iat_super, icell, iat_unit)
1308 : :
1309 [ + + ]: 329660 : R_red = SC%R_red(:,icell)
1310 [ + + ]: 329660 : R_sh = R_at - R_red
1311 : :
1312 [ - + ]: 82415 : if (.not. SC%is_mod0(R_sh)) then
1313 : 0 : ABI_ERROR('Invalid R_at assigned.')
1314 : : end if
1315 : :
1316 [ + + ]: 329660 : SC%R_sh(:,iat_super) = R_sh(:)
1317 : :
1318 : 82415 : end subroutine supercell_set_R_at
1319 : : !!***
1320 : :
1321 : : !----------------------------------------------------------------------
1322 : :
1323 : : !!****f* m_ifc_supercell/supercell_find_nearest_atom
1324 : : !! NAME
1325 : : !! supercell_find_nearest_atom
1326 : : !!
1327 : : !! FUNCTION
1328 : : !! Search for an atom in the supercell, and find the matching lattice vector and
1329 : : !! unitcell atom.
1330 : : !!
1331 : : !! INPUTS
1332 : : !! xcart(3)=cartesian coordinates of the atom we are searching
1333 : : !! typat=type of the atom we are searching
1334 : : !!
1335 : : !! OUTPUT
1336 : : !! iat_super=index of the matching atom in the supercell
1337 : : !! R_at(3)=lattice vector of the matching atom
1338 : : !! ucart(3)=displacement vector of xcart with respect to the position
1339 : : !! of the matching atom in the supercell
1340 : : !!
1341 : : !! SOURCE
1342 : :
1343 : 7278 : subroutine supercell_find_nearest_atom(SC,xcart,typat,iat_super,R_at,ucart)
1344 : : class(ifc_supercell_type),intent(in) :: SC
1345 : : real(dp), intent(in) :: xcart(3)
1346 : : integer, intent(in) :: typat
1347 : : integer, intent(out) :: iat_super
1348 : : integer, intent(out) :: R_at(3)
1349 : : real(dp), intent(out) :: ucart(3)
1350 : :
1351 : : logical :: found
1352 : : integer :: icell
1353 : : integer :: iat_unit
1354 : : integer :: iat
1355 : : integer :: ii,jj,kk
1356 : : integer :: RR(3)
1357 : : integer :: Rcan(3)
1358 : : integer :: Rshift_found(3)
1359 : : integer :: Rshift_try(3)
1360 : : real(dp) :: xred(3),xred_try(3)
1361 : : real(dp) :: xred_at_unitcell(3)
1362 : : real(dp) :: xred_origin_cell(3),xcart_origin_cell(3)
1363 : : real(dp) :: xcart_try(3), ucart_try(3)
1364 : : real(dp) :: xcart_match(3)
1365 : : real(dp) :: unorm_try, unorm
1366 : :
1367 [ + + ]: 29112 : xred = SC%cart_to_red(xcart)
1368 : :
1369 [ + + ]: 29112 : do ii=1,3
1370 : 29112 : RR(ii) = NINT(xred(ii))
1371 : : end do
1372 : :
1373 [ + + ]: 29112 : xred_origin_cell = xred - RR
1374 [ + + ]: 29112 : xcart_origin_cell = SC%red_to_cart(xred_origin_cell)
1375 : :
1376 : : ! Loop over atoms of the unitcell
1377 : 7278 : found = .false.
1378 : 7278 : iat_unit = zero
1379 : : xcart_match = zero
1380 : 7278 : Rshift_found = zero
1381 : 7278 : ucart = zero
1382 : 7278 : unorm = greatest_real
1383 : :
1384 [ + + ]: 22956 : do iat=1,SC%natom_unitcell
1385 : :
1386 [ + + ]: 15678 : if (SC%typat_unitcell(iat) .ne. typat) cycle
1387 : :
1388 [ + + ]: 54072 : xred_at_unitcell(:) = SC%xred_unitcell(:,iat)
1389 : :
1390 [ + + ]: 61350 : do ii=-1,1
1391 [ + + ]: 177894 : do jj=-1,1
1392 [ + + ]: 527202 : do kk=-1,1
1393 : :
1394 [ + + ]: 1459944 : Rshift_try(:) = (/ii, jj, kk/)
1395 [ + + ]: 1459944 : xred_try = Rshift_try + xred_at_unitcell
1396 [ + + ]: 1459944 : xcart_try = SC%red_to_cart(xred_try)
1397 : :
1398 [ + + ]: 1459944 : ucart_try(:) = xcart_try(:) - xcart_origin_cell(:)
1399 : 364986 : unorm_try = l2norm(ucart_try)
1400 : :
1401 : : ! Store the best match
1402 [ + + ]: 486648 : if (unorm_try < unorm) then
1403 : 42336 : found = .true.
1404 : : xcart_match = xcart_try
1405 : 42336 : unorm = unorm_try
1406 : 42336 : ucart = ucart_try
1407 : 42336 : Rshift_found = Rshift_try ! As it turns out, this is zero in test aifc_t01.
1408 : 42336 : iat_unit = iat
1409 : : end if
1410 : :
1411 : : end do
1412 : : end do
1413 : : end do
1414 : :
1415 : : end do
1416 : :
1417 [ - + ]: 7278 : if (.not. found) then
1418 : 0 : ABI_ERROR('Atom not found.')
1419 : : end if
1420 : :
1421 : : ! Shift RR with the matching Rshift to recover the unitcell lattice vector
1422 [ + + ]: 29112 : R_at = RR + Rshift_found
1423 : :
1424 [ + + ]: 29112 : Rcan = SC%modH(R_at)
1425 : 7278 : icell = SC%get_icell(Rcan)
1426 : 7278 : iat_super = SC%get_iat_supercell(icell, iat_unit)
1427 : :
1428 : 7278 : end subroutine supercell_find_nearest_atom
1429 : : !!***
1430 : :
1431 : : !----------------------------------------------------------------------
1432 : :
1433 : : !!****f* m_ifc_supercell/supercell_atoms_distance
1434 : : !! NAME
1435 : : !! supercell_atoms_distance
1436 : : !!
1437 : : !! FUNCTION
1438 : : !! Compute the cartesian vector and distance between two atoms,
1439 : : !! using the minimum image convention.
1440 : : !!
1441 : : !! INPUTS
1442 : : !! iat_super=index of the first atom in the supercell
1443 : : !! jat_super=index of the second atom in the supercell
1444 : : !!
1445 : : !! OUTPUT
1446 : : !! xij(3)=cartesian vector from atom iat_super to atom jat_super
1447 : : !! dist=distance between the two atoms
1448 : : !!
1449 : : !! SOURCE
1450 : :
1451 : 1551364 : subroutine supercell_atoms_distance(SC, iat_super, jat_super, xij, dist)
1452 : : class(ifc_supercell_type),intent(in) :: SC
1453 : : integer, intent(in) :: iat_super, jat_super
1454 : : real(dp), intent(out) :: xij(3)
1455 : : real(dp), intent(out) :: dist
1456 : :
1457 : : real(dp) :: xred_i(3), xred_j(3)
1458 : : real(dp) :: xred_ij_scaled(3)
1459 : :
1460 [ - + + + ]: 7756820 : xred_i = SC%unit_to_scaled(SC%get_xred(iat_super))
1461 [ - + + + ]: 7756820 : xred_j = SC%unit_to_scaled(SC%get_xred(jat_super))
1462 : :
1463 [ + + ]: 6205456 : xred_ij_scaled = xred_j - xred_i
1464 [ + + ]: 6205456 : xred_ij_scaled = SC%wrap_inbox(xred_ij_scaled)
1465 : :
1466 [ + + ]: 6205456 : xij = SC%red_scaled_to_cart(xred_ij_scaled)
1467 : 1551364 : dist = l2norm(xij)
1468 : :
1469 : 1551364 : end subroutine supercell_atoms_distance
1470 : : !!***
1471 : :
1472 : : !----------------------------------------------------------------------
1473 : :
1474 : : !!****f* m_ifc_supercell/supercell_calc_dij
1475 : : !! NAME
1476 : : !! supercell_calc_dij
1477 : : !!
1478 : : !! FUNCTION
1479 : : !! Compute the cartesian vector and distance between all pairs of atoms
1480 : : !! in the supercell, using the minimum image convention, and store the
1481 : : !! result in SC%dij. This array is not needed by all callers, so it is
1482 : : !! not computed as part of supercell_init: call this routine explicitly
1483 : : !! whenever SC%dij is required.
1484 : : !!
1485 : : !! SIDE EFFECTS
1486 : : !! SC=supercell object, dij is allocated (if needed) and filled on output
1487 : : !!
1488 : : !! SOURCE
1489 : :
1490 : 0 : subroutine supercell_calc_dij(SC)
1491 : :
1492 : : class(ifc_supercell_type),intent(inout) :: SC
1493 : :
1494 : : integer :: iatom,jatom
1495 : : real(dp) :: xij(3)
1496 : : real(dp) :: dist
1497 : :
1498 [ # # ]: 0 : if (.not. allocated(SC%dij)) then
1499 [ # # # # : 0 : ABI_MALLOC(SC%dij, (SC%natom,SC%natom,4))
# # # # #
# # # # #
# # # # #
# ]
1500 : : end if
1501 : :
1502 [ # # ]: 0 : do iatom=1,SC%natom
1503 [ # # ]: 0 : do jatom=1,SC%natom
1504 : 0 : call SC%atoms_distance(iatom, jatom, xij, dist)
1505 : 0 : SC%dij(iatom,jatom,1) = dist
1506 [ # # ]: 0 : SC%dij(iatom,jatom,2:4) = xij(1:3)
1507 : : end do
1508 : : end do
1509 : :
1510 : 0 : end subroutine supercell_calc_dij
1511 : : !!***
1512 : :
1513 : : !----------------------------------------------------------------------
1514 : :
1515 : : !!****f* m_ifc_supercell/supercell_calc_wghatm
1516 : : !! NAME
1517 : : !! supercell_calc_wghatm
1518 : : !!
1519 : : !! FUNCTION
1520 : : !! Compute the multiplicity weight of every (unitcell atom pair, canonical
1521 : : !! lattice vector) triplet needed to downfold a fine q-mesh (the
1522 : : !! Born-von Karman supercell described by SC) into real-space interatomic
1523 : : !! force constants.
1524 : : !!
1525 : : !! The historical implementation (wght9 in m_dynmat) works with an
1526 : : !! EXPLICIT "big box" of real-space lattice vectors that can contain
1527 : : !! several distinct entries representing the very same periodic image
1528 : : !! (e.g. both R=+1 and R=-1 along some direction, when they are equidistant
1529 : : !! from a given atom pair): each such duplicate entry gets a fractional
1530 : : !! weight 1/N, so that summing over all N of them (as ftifc_r2q, dynmat_dq
1531 : : !! and asrif9's ASR sum all do) reconstructs the correct total.
1532 : : !!
1533 : : !! SC%R_red, by contrast, is already a DEDUPLICATED enumeration: exactly
1534 : : !! one canonical representative per coset of the superlattice (that is the
1535 : : !! whole point of building it from the Hermite Normal Form). Any two
1536 : : !! candidate images that are equally close to a given atom pair and
1537 : : !! related by a superlattice translation are therefore, by construction,
1538 : : !! the SAME canonical cell here -- there is no second, duplicate slot for
1539 : : !! wght9-style fractional weight to distribute over. Applying a 1/N
1540 : : !! multiplicity weight on top of this deduplicated enumeration would
1541 : : !! under-count reconstruction sums (ASR in particular) by exactly that
1542 : : !! factor for atom pairs sitting on a supercell boundary -- this was
1543 : : !! tried and confirmed wrong against the anaddb test suite for an
1544 : : !! extreme (ngqpt=2x2x2) edge case where such boundary coincidences are
1545 : : !! common. The correct weight for this one-slot-per-coset representation
1546 : : !! is simply 1: atmfrc(icell), as produced by the ordinary (unweighted)
1547 : : !! Fourier transform, already IS the full aliased coefficient needed to
1548 : : !! exactly reconstruct D(q) at every original mesh q-point.
1549 : : !!
1550 : : !! (Note: unlike wght9's explicit big box, this uniform weighting does not
1551 : : !! attempt to symmetrize the real-space representation over tied
1552 : : !! candidates for BETTER INTERPOLATION AWAY FROM the original mesh -- only
1553 : : !! one arbitrary tied representative is kept, chosen by SC%R_red's own
1554 : : !! minimum-distance search. This only matters for atom pairs that sit
1555 : : !! exactly on a supercell boundary, and only affects interpolated values
1556 : : !! at q-points off the original mesh.)
1557 : : !!
1558 : : !! OUTPUT
1559 : : !! wghatm(natom_unitcell,natom_unitcell,Ncells)=weight of each triplet
1560 : : !! (uniformly 1).
1561 : : !!
1562 : : !! SOURCE
1563 : :
1564 : 63 : subroutine supercell_calc_wghatm(SC, wghatm)
1565 : :
1566 : : class(ifc_supercell_type),intent(in) :: SC
1567 : : real(dp),allocatable,intent(out) :: wghatm(:,:,:)
1568 : :
1569 [ - + + - : 315 : ABI_MALLOC(wghatm, (SC%natom_unitcell,SC%natom_unitcell,SC%Ncells))
- + + - -
+ + - + -
+ - - + -
+ - + ]
1570 [ + + + + : 153002 : wghatm = one
+ + ]
1571 : :
1572 : 63 : end subroutine supercell_calc_wghatm
1573 : : !!***
1574 : :
1575 : : !----------------------------------------------------------------------
1576 : :
1577 : : !!****f* m_ifc_supercell/supercell_calc_ws_wghatm
1578 : : !! NAME
1579 : : !! supercell_calc_ws_wghatm
1580 : : !!
1581 : : !! FUNCTION
1582 : : !! Build the explicit "big box" of real-space lattice vectors R needed to Fourier transform
1583 : : !! dynamical matrices from the q-mesh dual to the supercell, together with the Wigner-Seitz
1584 : : !! weights, for an arbitrary supercell matrix (generalizes wght9/make_bigbox of m_dynmat).
1585 : : !!
1586 : : !! For every atom pair (i,j) and every coset of the superlattice, the images R of that
1587 : : !! coset are ranked by the distance |R + xcart_unit(j) - xcart_unit(i)|. Only the closest
1588 : : !! ones (all of them, when several are equidistant) are kept, each with weight 1/N for N
1589 : : !! tied images. For each atom pair the weights of all the R sum to Ncells, and summing over
1590 : : !! the tied images (as ftifc_r2q, dynmat_dq and asrif9 do) reconstructs exactly the aliased
1591 : : !! coefficient at every q-point of the mesh, while making the interpolation symmetric
1592 : : !! between the tied images.
1593 : : !!
1594 : : !! INPUTS
1595 : : !! toldist=tolerance (squared distances are compared within 2*toldist).
1596 : : !!
1597 : : !! OUTPUT
1598 : : !! nrpt=number of R vectors kept (those with a nonzero weight for at least one pair).
1599 : : !! R_red(3,nrpt)=integer coordinates of the R vectors, in units of the unit cell vectors.
1600 : : !! wghatm(natom_unitcell,natom_unitcell,nrpt)=weights.
1601 : : !!
1602 : : !! SOURCE
1603 : :
1604 : 201 : subroutine supercell_calc_ws_wghatm(SC, toldist, nrpt, R_red, wghatm)
1605 : :
1606 : : class(ifc_supercell_type),intent(in) :: SC
1607 : : real(dp),intent(in) :: toldist
1608 : : integer,intent(out) :: nrpt
1609 : : integer,allocatable,intent(out) :: R_red(:,:)
1610 : : real(dp),allocatable,intent(out) :: wghatm(:,:,:)
1611 : :
1612 : : integer,parameter :: nmax=2, ntr=(2*nmax+1)**3
1613 : : integer :: n1,n2,n3,it,ic,ia,ib,ipass,nties,ii
1614 : : integer :: Tint(3,ntr), Rp(3), lo(3), hi(3)
1615 : 201 : integer,allocatable :: idx(:,:,:)
1616 : : real(dp) :: Tcart(3,ntr), vec0(3), d2(ntr), dmin
1617 : :
1618 : : ! Superlattice translations (integer coordinates of the unit cell): rows of the supercell matrix.
1619 : 201 : it = 0
1620 [ + + ]: 1206 : do n3=-nmax,nmax
1621 [ + + ]: 6231 : do n2=-nmax,nmax
1622 [ + + ]: 31155 : do n1=-nmax,nmax
1623 : 25125 : it = it + 1
1624 [ + + ]: 100500 : Tint(:,it) = n1*SC%matrix(1,:) + n2*SC%matrix(2,:) + n3*SC%matrix(3,:)
1625 [ + + + + : 532650 : Tcart(:,it) = matmul(SC%rprimd, real(Tint(:,it), kind=dp))
+ + + + ]
1626 : : end do
1627 : : end do
1628 : : end do
1629 : :
1630 [ + + ]: 804 : do ii=1,3
1631 [ + - - + : 115521 : lo(ii) = minval(SC%R_red(ii,:)) + minval(Tint(ii,:))
- - + + +
+ ]
1632 [ + - - + : 116325 : hi(ii) = maxval(SC%R_red(ii,:)) + maxval(Tint(ii,:))
- - + + +
+ ]
1633 : : end do
1634 [ - + + - : 2557762 : ABI_CALLOC(idx, (lo(1):hi(1),lo(2):hi(2),lo(3):hi(3)))
- + + - -
+ + - + -
+ - - + +
- + + + +
+ + ]
1635 : :
1636 [ + + ]: 603 : do ipass=1,2
1637 [ + + ]: 402 : if (ipass == 2) then
1638 : : ! Number the R vectors in lexicographic order (first coordinate slowest), the order in which
1639 : : ! the historical big box (bigbx9) enumerated them: input files that refer to R by its index
1640 : : ! (e.g. a-TDEP's ifc_in.dat) rely on it.
1641 [ + + + + : 2556757 : nrpt = count(idx /= 0)
+ + + + ]
1642 [ - + + - : 603 : ABI_MALLOC(R_red, (3,nrpt))
+ - + - -
+ - + -
+ ]
1643 [ - + + - : 419862 : ABI_CALLOC(wghatm, (SC%natom_unitcell,SC%natom_unitcell,nrpt))
- + + - -
+ + - + -
+ - - + -
+ - + + +
+ + + + ]
1644 : 201 : nrpt = 0
1645 [ + + ]: 3783 : do n1=lo(1),hi(1)
1646 [ + + ]: 100759 : do n2=lo(2),hi(2)
1647 [ + + ]: 2575478 : do n3=lo(3),hi(3)
1648 [ + + ]: 2474920 : if (idx(n1,n2,n3) == 0) cycle
1649 : 23775 : nrpt = nrpt + 1
1650 : 23775 : idx(n1,n2,n3) = nrpt
1651 [ + + ]: 2643221 : R_red(:,nrpt) = [n1,n2,n3]
1652 : : end do
1653 : : end do
1654 : : end do
1655 : : end if
1656 [ + + ]: 26563 : do ic=1,SC%Ncells
1657 [ + + ]: 112072 : do ia=1,SC%natom_unitcell
1658 [ + + ]: 467620 : do ib=1,SC%natom_unitcell
1659 [ + + + + : 6763050 : vec0 = matmul(SC%rprimd, real(SC%R_red(:,ic), kind=dp)) + SC%xcart_unit(:,ib) - SC%xcart_unit(:,ia)
+ + + - +
+ ]
1660 [ + + ]: 44849700 : do it=1,ntr
1661 [ + + ]: 178330950 : d2(it) = sum((vec0 + Tcart(:,it))**2)
1662 : : end do
1663 [ - + + - : 45205650 : dmin = minval(d2)
+ + ]
1664 [ + + + + ]: 44849700 : nties = count(d2 - dmin <= two*toldist)
1665 [ + + ]: 44935410 : do it=1,ntr
1666 [ + + ]: 44493750 : if (d2(it) - dmin > two*toldist) cycle
1667 [ + + ]: 1739024 : Rp = SC%R_red(:,ic) + Tint(:,it)
1668 [ + + ]: 790706 : if (ipass == 1) then
1669 : 217378 : idx(Rp(1),Rp(2),Rp(3)) = 1
1670 : : else
1671 : 217378 : wghatm(ia,ib,idx(Rp(1),Rp(2),Rp(3))) = one / real(nties, kind=dp)
1672 : : end if
1673 : : end do
1674 : : end do
1675 : : end do
1676 : : end do
1677 : : end do
1678 : :
1679 [ - + ]: 201 : ABI_FREE(idx)
1680 : :
1681 : 201 : end subroutine supercell_calc_ws_wghatm
1682 : : !!***
1683 : :
1684 : : !----------------------------------------------------------------------
1685 : :
1686 : : !!****f* m_ifc_supercell/supercell_shortest_vector_norm
1687 : : !! NAME
1688 : : !! supercell_shortest_vector_norm
1689 : : !!
1690 : : !! FUNCTION
1691 : : !! Norm of the shortest nonzero superlattice vector, i.e. the smallest
1692 : : !! cartesian distance between two periodic images of the origin under the
1693 : : !! supercell's superlattice. Generalizes the "r_inscribed_sphere" quantity
1694 : : !! computed inline in wght9 (m_dynmat.F90) for a diagonal ngqpt*rprimd
1695 : : !! superlattice to an arbitrary supercell matrix.
1696 : : !!
1697 : : !! SOURCE
1698 : :
1699 : 264 : function supercell_shortest_vector_norm(SC) result(vnorm)
1700 : :
1701 : : class(ifc_supercell_type),intent(in) :: SC
1702 : : real(dp) :: vnorm
1703 : :
1704 : : integer :: di,dj,dk
1705 : : integer :: Ht(3,3), T_try(3)
1706 : : real(dp) :: cand_cart(3), this_norm
1707 : :
1708 [ + + + + ]: 3432 : Ht = transpose(SC%HH)
1709 : : vnorm = greatest_real
1710 : :
1711 [ + + ]: 1584 : do di=-2,2
1712 [ + + ]: 8184 : do dj=-2,2
1713 [ + + ]: 40920 : do dk=-2,2
1714 [ + + ]: 33000 : if (di==0 .and. dj==0 .and. dk==0) cycle
1715 [ + + ]: 130944 : T_try(:) = (/di, dj, dk/)
1716 [ + + + + : 654720 : cand_cart(:) = SC%red_to_cart(real(matmul(Ht, T_try), kind=dp))
+ + - + +
+ ]
1717 : 32736 : this_norm = l2norm(cand_cart)
1718 : 39336 : vnorm = min(vnorm, this_norm)
1719 : : end do
1720 : : end do
1721 : : end do
1722 : :
1723 : 264 : end function supercell_shortest_vector_norm
1724 : : !!***
1725 : :
1726 : : !----------------------------------------------------------------------
1727 : :
1728 : : ! Find the list of R vectors belonging to the shell of lattice vector R.
1729 : :
1730 : : ! subroutine supercell_find_R_shell(SC, RR, N_shell, R_shell, dist)
1731 : : ! class(ifc_supercell_type),intent(in) :: SC
1732 : : ! integer, intent(in) :: RR(3)
1733 : : ! real(dp), intent(out) :: dist
1734 : : ! integer, intent(out) :: N_shell
1735 : : ! integer, allocatable, intent(out) :: R_shell(:,:) ! R_shell(N_shell, 3)
1736 : : !
1737 : : ! end subroutine supercell_find_R_shell
1738 : :
1739 : : ! Find the list of R vectors belonging to the shell of lattice vector R.
1740 : : !subroutine supercell_find_atom_shell(SC, iat, N_shell, R_shell, dist)
1741 : : ! class(ifc_supercell_type),intent(in) :: SC
1742 : : ! integer, intent(in) :: iat
1743 : : ! real(dp), intent(out) :: dist
1744 : : ! integer, intent(out) :: N_shell
1745 : : ! integer, allocatable, intent(out) :: R_shell(:,:)
1746 : : !
1747 : : !end subroutine supercell_find_atom_shell
1748 : :
1749 : : !----------------------------------------------------------------------
1750 : :
1751 [ # # # # : 355950 : end module m_ifc_supercell
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
|
