Branch data Line data Source code
1 : : !!****m* ABINIT/m_tdep_sym
2 : : !! NAME
3 : : !! m_tdep_sym
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains the TDEP Symmetries data type, which holds the
7 : : !! symmetry operations of the crystal (in reduced and cartesian
8 : : !! coordinates), and the routines that find the symmetry operations
9 : : !! relating atoms, pairs, triplets and quadruplets of atoms of the supercell.
10 : : !!
11 : : !! COPYRIGHT
12 : : !! Copyright (C) 2011-2026 ABINIT group (GA,FB,JB)
13 : : !! This file is distributed under the terms of the
14 : : !! GNU General Public License, see ~abinit/COPYING
15 : : !! or http://www.gnu.org/copyleft/gpl.txt .
16 : : !!
17 : : !! SOURCE
18 : :
19 : : #if defined HAVE_CONFIG_H
20 : : #include "config.h"
21 : : #endif
22 : :
23 : : #include "abi_common.h"
24 : :
25 : : module m_tdep_sym
26 : :
27 : : use defs_basis
28 : : use m_abicore
29 : : use m_errors
30 : : use m_xmpi
31 : : use m_matrix, only : mati3inv
32 : : use m_symtk, only : symatm
33 : : use m_symfind, only : symfind, symanal, symlatt
34 : : use m_tdep_dataset, only : atdep_dataset_type, MPI_enreg_type
35 : : use m_tdep_latt, only : Lattice_type, tdep_make_inbox
36 : : use m_tdep_supercell, only : tdep_Supercell_type
37 : :
38 : : implicit none
39 : : !!***
40 : :
41 : : !----------------------------------------------------------------------
42 : :
43 : : !!****t* m_tdep_sym/Symmetries_type
44 : : !! NAME
45 : : !! Symmetries_type
46 : : !!
47 : : !! FUNCTION
48 : : !! Holds the symmetry operations of the crystal, in reduced and cartesian
49 : : !! coordinates, and the action of each of them on the atoms of the supercell.
50 : : !!
51 : : !! SOURCE
52 : :
53 : : type,public :: Symmetries_type
54 : :
55 : : integer :: msym
56 : : ! Maximum number of symmetry operations
57 : :
58 : : integer :: nptsym
59 : : ! Number of point symmetry operations of the Bravais lattice
60 : :
61 : : integer :: nsym
62 : : ! Number of symmetry operations of the crystal
63 : :
64 : : integer :: spgroup
65 : : ! Space group number
66 : :
67 : : integer, allocatable :: ptsymrel(:,:,:)
68 : : ! ptsymrel(3,3,msym)
69 : : ! Point symmetry operations of the Bravais lattice, in reduced coordinates
70 : :
71 : : integer, allocatable :: symrec(:,:,:)
72 : : ! symrec(3,3,nsym)
73 : : ! Symmetry operations in reciprocal space, in reduced coordinates
74 : : ! (inverse of the transpose of symrel)
75 : :
76 : : integer, allocatable :: symrel(:,:,:)
77 : : ! symrel(3,3,nsym)
78 : : ! Symmetry operations in real space, in reduced coordinates
79 : :
80 : : integer, allocatable :: symafm(:)
81 : : ! symafm(nsym)
82 : : ! Ferromagnetic (1) or antiferromagnetic (-1) character of each symmetry operation
83 : :
84 : : integer, allocatable :: indsym(:,:,:)
85 : : ! indsym(4,nsym,natom)
86 : : ! For the symmetry operation isym applied to the atom iatom of the supercell,
87 : : ! indsym(4,isym,iatom) gives the equivalent atom in the reference unitcell,
88 : : ! and indsym(1:3,isym,iatom) the corresponding lattice vector
89 : : ! (see tdep_SearchS_1at).
90 : :
91 : : double precision, allocatable :: S_ref(:,:,:,:)
92 : : ! S_ref(3,3,nsym,2)
93 : : ! Symmetry operations in cartesian coordinates (only S_ref(:,:,:,1) is used)
94 : :
95 : : double precision, allocatable :: S_inv(:,:,:,:)
96 : : ! S_inv(3,3,nsym,2)
97 : : ! Inverse (i.e. transpose) of S_ref (only S_inv(:,:,:,1) is used)
98 : :
99 : : double precision, allocatable :: tnons(:,:)
100 : : ! tnons(3,nsym)
101 : : ! Non-symmorphic translations, in reduced coordinates
102 : :
103 : : double precision, allocatable :: xred_zero(:,:)
104 : : ! xred_zero(3,natom_unitcell)
105 : : ! Copy of Invar%xred_unitcell (atoms anchored so atom 1 sits at the
106 : : ! origin (0,0,0), wrapped into [0,1); see tdep_read_input). This
107 : : ! anchor-on-atom-1 convention is not local to symmetry-finding: it is
108 : : ! also the convention used for the reference unitcell throughout the
109 : : ! rest of atdep (Crystal and Supercell objects are built from this same
110 : : ! array, see tdep_init_crystal).
111 : :
112 : : contains
113 : :
114 : : procedure :: init => tdep_make_sym
115 : : procedure :: free => tdep_destroy_sym
116 : :
117 : : end type Symmetries_type
118 : : !!***
119 : :
120 : : public :: tdep_SearchS_1at
121 : : public :: tdep_SearchS_2at
122 : : public :: tdep_SearchS_3at
123 : : public :: tdep_SearchS_4at
124 : : public :: tdep_destroy_sym
125 : : public :: tdep_calc_indsym2
126 : :
127 : : contains
128 : :
129 : : !=====================================================================================================
130 : :
131 : : !!****f* m_tdep_sym/tdep_make_sym
132 : : !! NAME
133 : : !! tdep_make_sym
134 : : !!
135 : : !! FUNCTION
136 : : !! Initialize the Symmetries object: find the point symmetry operations of
137 : : !! the Bravais lattice (symlatt), then the symmetry operations of the
138 : : !! crystal and their non-symmorphic translations (symfind), compute the
139 : : !! corresponding operations in reciprocal space (symrec) and the space group
140 : : !! (symanal). Finally, transform the symmetry operations into cartesian
141 : : !! coordinates (S_ref), compute their inverse (S_inv), check that they are
142 : : !! orthogonal, and write them in the _sym.dat file.
143 : : !!
144 : : !! INPUTS
145 : : !! MPIdata = Info on MPI parallelism.
146 : : !! Lattice = Lattice object describing the ideal structure.
147 : : !!
148 : : !! OUTPUT
149 : : !! Sym = Symmetries object, initialized (except Sym%indsym, which is
150 : : !! computed in tdep_SearchS_1at).
151 : : !!
152 : : !! SIDE EFFECTS
153 : : !! Invar = Input object containing the input variables.
154 : : !! Invar%bravais is used by symlatt.
155 : : !!
156 : : !! SOURCE
157 : :
158 : 46 : subroutine tdep_make_sym(Sym,Invar,MPIdata,Lattice)
159 : :
160 : : class(Symmetries_type),intent(out) :: Sym
161 : : type(atdep_dataset_type),intent(inout) :: Invar
162 : : type(MPI_enreg_type), intent(in) :: MPIdata
163 : : type(Lattice_type),intent(in) :: Lattice
164 : :
165 : : integer :: nspden,use_inversion,chkprim
166 : : integer :: ptgroupma,isym,ii,jj
167 : : integer :: bravais(11)
168 : 46 : integer, allocatable :: symrel(:,:,:),symafm(:)
169 : : double precision :: genafm(3)
170 : 46 : double precision, allocatable :: spinat(:,:)
171 : 46 : double precision, allocatable :: tnons(:,:)
172 [ + - ]: 46 : double precision, allocatable :: tmp1(:,:)
173 : :
174 : :
175 : : ! Compute all the symetries coming from the bravais lattice
176 : : ! The routine used is symlatt (from Abinit code)
177 : : ! GA: FIXME msym needs to be an input variable
178 : 46 : Sym%msym=1000 !msym needs to be very large due to non-primitive cell calculations
179 [ - + - + : 598046 : ABI_MALLOC(Sym%ptsymrel,(3,3,Sym%msym)) ; Sym%ptsymrel(:,:,:)=0
+ + + + +
+ ]
180 : 46 : call symlatt(Invar%bravais,std_out,Sym%msym,Sym%nptsym,Sym%ptsymrel,Lattice%rprimdt,tol8)
181 : 46 : write(ab_out,'(a,1x,11(i4,1x))')' bravais=',Invar%bravais(:)
182 : 46 : Sym%nsym=Sym%nptsym
183 : :
184 : : ! Invar%xred_unitcell is already anchored on atom 1 (shifted and wrapped into
185 : : ! [0,1), see tdep_read_input): symlatt/symfind below need atom 1 at the
186 : : ! origin to find the space group, and this same Sym%xred_zero array is
187 : : ! reused, unmodified, as the reference unitcell (Crystal%xred, and from
188 : : ! there SC%xred_unitcell) for the rest of the code (see tdep_init_crystal).
189 [ - + + - : 138 : ABI_MALLOC(Sym%xred_zero,(3,Invar%natom_unitcell))
+ - + - -
+ - + -
+ ]
190 [ + + + + ]: 486 : Sym%xred_zero(:,:) = Invar%xred_unitcell(:,:)
191 : :
192 : : ! Calculation of the (non-symmorphic) translations
193 : 46 : nspden=1
194 [ - + + - : 532 : ABI_MALLOC(spinat,(3,Invar%natom_unitcell)); spinat(:,:)=0.d0
- + - + +
+ + + ]
195 : 46 : use_inversion=1
196 [ - + + - : 598138 : ABI_MALLOC(symrel ,(3,3,Sym%msym)) ; symrel (:,:,:)=0
+ - + - -
+ - + + +
+ + + + ]
197 [ - + + - : 184138 : ABI_MALLOC(tnons ,(3,Sym%msym)) ; tnons (:,:)=0.d0
+ - + - -
+ - + + +
+ + ]
198 [ - + + - : 46138 : ABI_MALLOC(symafm ,(Sym%msym)) ; symafm (:)=1
+ - + - -
+ - + +
+ ]
199 : : call symfind(Lattice%gprimd,Sym%msym,Invar%natom_unitcell,Sym%nptsym,nspden,Sym%nsym,&
200 : 46 : & 0,Sym%ptsymrel,spinat,symafm,symrel,tnons,tol8,Invar%typat_unitcell,use_inversion,Sym%xred_zero)
201 [ - + ]: 46 : ABI_FREE(spinat)
202 [ - + + - : 22160 : ABI_MALLOC(Sym%symrel,(3,3,Sym%nsym)) ; Sym%symrel(:,:,:)=0
+ - + - -
+ - + - +
+ + + + +
+ ]
203 [ - + + - : 6914 : ABI_MALLOC(Sym%tnons ,(3,Sym%nsym)) ; Sym%tnons (:,:)=0.d0
+ - + - -
+ - + - +
+ + + + ]
204 [ - + + - : 1832 : ABI_MALLOC(Sym%symafm ,(Sym%nsym)) ; Sym%symafm (:)=1
+ - + - -
+ - + - +
+ + ]
205 [ + + + + : 22068 : Sym%symrel(:,:,:)=symrel(:,:,1:Sym%nsym)
+ + ]
206 [ + + + + ]: 6822 : Sym%tnons (:,:)=tnons (:,1:Sym%nsym)
207 [ + + ]: 1740 : Sym%symafm (:)=symafm (1:Sym%nsym)
208 [ - + ]: 46 : ABI_FREE(symrel)
209 [ - + ]: 46 : ABI_FREE(tnons)
210 [ - + ]: 46 : ABI_FREE(symafm)
211 [ + + ]: 46 : if (Sym%nptsym.ne.Sym%nsym) then
212 : 3 : write(std_out,'(a,i4,a,i4)') 'WARNING: nsym=',Sym%nsym,' is not equal to nptsym=',Sym%nptsym
213 : 3 : write(std_out,'(a)') ' Symrel='
214 [ + + ]: 293 : do isym=1,Sym%nsym
215 : 290 : write(std_out,*) ' For sym=',isym
216 : 290 : write(std_out,*) Sym%symrel(1,:,isym)
217 : 290 : write(std_out,*) Sym%symrel(2,:,isym)
218 : 293 : write(std_out,*) Sym%symrel(3,:,isym)
219 : : end do
220 : 3 : write(std_out,'(a)') ' Ptsymrel='
221 [ + + ]: 111 : do isym=1,Sym%nptsym
222 : 108 : write(std_out,*) ' For sym=',isym
223 : 108 : write(std_out,*) Sym%ptsymrel(1,:,isym)
224 : 108 : write(std_out,*) Sym%ptsymrel(2,:,isym)
225 : 111 : write(std_out,*) Sym%ptsymrel(3,:,isym)
226 : : end do
227 : : end if
228 : :
229 : : ! Calculation of symrec
230 [ - + + - : 22160 : ABI_MALLOC(Sym%symrec,(3,3,Sym%nsym)); Sym%symrec(:,:,:)=0
+ - + - -
+ - + - +
+ + + + +
+ ]
231 [ + + ]: 1740 : do isym=1,Sym%nsym
232 : 1740 : call mati3inv(Sym%symrel(:,:,isym),Sym%symrec(:,:,isym))
233 : : end do
234 : : ! Calculation of spgroup
235 : 46 : chkprim=0
236 : 46 : genafm(:)=0.d0
237 : : ptgroupma=0
238 : 46 : Sym%spgroup=0
239 : 46 : call symanal(bravais,chkprim,genafm,Sym%msym,Sym%nsym,ptgroupma,Lattice%rprimdt,Sym%spgroup,Sym%symafm,Sym%symrel,Sym%tnons,tol8)
240 : :
241 : 46 : write(ab_out,'(a)') ' Writing symmetries in file sym.dat'
242 : :
243 : : ! Transform S_ref in cartesian coordinates
244 [ - + + - : 22206 : ABI_MALLOC(Sym%S_ref,(3,3,Sym%nsym,2)) ; Sym%S_ref(:,:,:,1)=real(Sym%symrel(:,:,1:Sym%nsym))
+ - + - +
- - + - +
- + + + +
+ + + ]
245 [ - + + - : 22114 : ABI_MALLOC(Sym%S_inv,(3,3,Sym%nsym,2)) ; Sym%S_inv(:,:,:,1)=zero
- + - + -
+ + + + +
+ + ]
246 [ - + + + : 598 : ABI_MALLOC(tmp1,(3,3)); tmp1(:,:)=0.d0
+ + ]
247 [ + - ]: 46 : if (MPIdata%iam_master) open(unit=75,file=trim(Invar%output_prefix)//'_sym.dat')
248 [ + + ]: 1740 : do isym=1,Sym%nsym
249 [ + - ]: 1694 : if (MPIdata%iam_master) then
250 : 1694 : write(75,*) ' '
251 : 1694 : write(75,*) 'For isym=',isym
252 : 1694 : write(75,*) 'In reduced coordinates:'
253 : 1694 : write(75,'(3(f5.2,1x))') Sym%S_ref(1,1,isym,1),Sym%S_ref(1,2,isym,1),Sym%S_ref(1,3,isym,1)
254 : 1694 : write(75,'(3(f5.2,1x))') Sym%S_ref(2,1,isym,1),Sym%S_ref(2,2,isym,1),Sym%S_ref(2,3,isym,1)
255 : 1694 : write(75,'(3(f5.2,1x))') Sym%S_ref(3,1,isym,1),Sym%S_ref(3,2,isym,1),Sym%S_ref(3,3,isym,1)
256 : : end if
257 : 1694 : call DGEMM('N','N',3,3,3,1.d0,Lattice%rprimdt,3,Sym%S_ref(:,:,isym,1),3,0.d0,tmp1,3)
258 : 1694 : call DGEMM('N','N',3,3,3,1.d0,tmp1,3,Lattice%rprimdtm1,3,0.d0,Sym%S_ref(:,:,isym,1),3)
259 [ + + ]: 6776 : do ii=1,3
260 [ + + ]: 22022 : do jj=1,3
261 [ + + ]: 20328 : if (abs(Sym%S_ref(ii,jj,isym,1)).lt.tol8) Sym%S_ref(ii,jj,isym,1)=zero
262 : : end do
263 : : end do
264 : : ! Inversion of S_ref (which is equivalent to the transposition)
265 [ + + ]: 6776 : do ii=1,3
266 [ + + ]: 22022 : do jj=1,3
267 : 20328 : Sym%S_inv(ii,jj,isym,1)=Sym%S_ref(jj,ii,isym,1)
268 : : end do
269 : : end do
270 [ + - ]: 1740 : if (MPIdata%iam_master) then
271 : 1694 : write(75,*) 'In cartesian coordinates:'
272 : 1694 : write(75,'(3(f5.2,1x))') Sym%S_ref(1,1,isym,1),Sym%S_ref(1,2,isym,1),Sym%S_ref(1,3,isym,1)
273 : 1694 : write(75,'(3(f5.2,1x))') Sym%S_ref(2,1,isym,1),Sym%S_ref(2,2,isym,1),Sym%S_ref(2,3,isym,1)
274 : 1694 : write(75,'(3(f5.2,1x))') Sym%S_ref(3,1,isym,1),Sym%S_ref(3,2,isym,1),Sym%S_ref(3,3,isym,1)
275 : : end if
276 : : end do
277 [ + - ]: 46 : if (MPIdata%iam_master) close(75)
278 : :
279 : : ! We verify that S.S^T = Id
280 [ + + + + ]: 598 : tmp1(:,:)=zero
281 [ + + ]: 1740 : do isym=1,Sym%nsym
282 : 1694 : call DGEMM('T','N',3,3,3,1.d0,Sym%S_ref(:,:,isym,1),3,Sym%S_ref(:,:,isym,1),3,0.d0,tmp1,3)
283 [ + + ]: 6822 : do ii=1,3
284 [ + + ]: 22022 : do jj=1,3
285 [ + + + - : 20328 : if ((ii/=jj.and.abs(tmp1(ii,jj)).gt.tol8).or.(ii==jj.and.abs(tmp1(ii,jj)-1.d0).gt.tol8)) then
+ + - + ]
286 : 0 : write(ab_out,'(a)') ' STOP : the matrix is not orthogonal : '
287 : 0 : write(ab_out,'(a,1x,i3,1x,i3,1x,f15.10)') ' ii, jj, sym(ii,jj)=',ii,jj,tmp1(ii,jj)
288 : 0 : ABI_ERROR('The matrix is not orthogonal')
289 : : end if
290 : : end do
291 : : end do
292 : : end do
293 : :
294 [ - + ]: 46 : ABI_FREE(tmp1)
295 : :
296 : 46 : end subroutine tdep_make_sym
297 : : !!***
298 : :
299 : : !=====================================================================================================
300 : :
301 : : !!****f* m_tdep_sym/tdep_SearchS_1at
302 : : !! NAME
303 : : !! tdep_SearchS_1at
304 : : !!
305 : : !! FUNCTION
306 : : !! Find, for every atom of the supercell SC and every symmetry operation,
307 : : !! the equivalent atom of the reference unitcell and the lattice vector
308 : : !! connecting them, and store the result in Sym%indsym.
309 : : !!
310 : : !! Concretely, this routine, in order:
311 : : !! 1) Calls symatm to compute indsym for the natom_unitcell atoms of the
312 : : !! reference unitcell only (Sym%indsym(:,:,1:SC%natom_unitcell)), from
313 : : !! the symmetry operations (Sym%symrec/tnons) and the "zeroed"
314 : : !! unitcell positions Sym%xred_zero (computed once in tdep_make_sym).
315 : : !! 2) Extends indsym to every atom of the supercell: the unitcell atom
316 : : !! equivalent to a supercell atom iatom is trivially
317 : : !! mod(iatom-1,natom_unitcell)+1 (SC's natural ordering repeats the
318 : : !! unitcell motif atom-by-atom, cell-by-cell), and the lattice vector
319 : : !! connecting the two is recovered by comparing the ideal position of
320 : : !! iatom against that of its unitcell-equivalent site (the "motif"),
321 : : !! both expressed in unitcell-reduced coordinates.
322 : : !! 3) For every pair of atoms (iatom,jatom) and every symmetry operation,
323 : : !! combines the above into indsym2, an 8-component vector describing
324 : : !! how the operation maps both atoms at once (the first four
325 : : !! components as in indsym for iatom, the last four for jatom).
326 : : !! This quantity is only used, in this routine, for the optional debug
327 : : !! output controlled by the internal debug parameter: it is not
328 : : !! returned nor stored in Sym. Elsewhere in the code, the equivalent per-pair, per-symmetry
329 : : !! quantity is instead computed on the fly by tdep_calc_indsym2.
330 : : !!
331 : : !! INPUTS
332 : : !! Invar = Input object containing the input variables.
333 : : !! MPIdata = Info on MPI parallelism.
334 : : !! SC = Supercell object describing the ideal structure. Its atoms provide
335 : : !! the reference ordering used throughout the rest of the code, so
336 : : !! this routine must be called after tdep_match_supercell.
337 : : !!
338 : : !! OUTPUT
339 : : !!
340 : : !! SIDE EFFECTS
341 : : !! Sym = Symmetries object. Sym%indsym(4,nsym,natom) is allocated and
342 : : !! filled: for symmetry operation isym applied to atom iatom,
343 : : !! indsym(4,isym,iatom) gives the equivalent atom iat_sym in the
344 : : !! reference unitcell, and indsym(1:3,isym,iatom) gives the lattice
345 : : !! vector R_0 such that
346 : : !! R^{-1} (xred(:,iatom)-tau) = xred(:,iat_sym) + R_0 .
347 : : !!
348 : : !! SOURCE
349 : :
350 : 46 : subroutine tdep_SearchS_1at(Invar,MPIdata,Sym,SC)
351 : :
352 : : type(atdep_dataset_type), intent(in) :: Invar
353 : : type(Symmetries_type), intent(inout) :: Sym
354 : : type(MPI_enreg_type), intent(in) :: MPIdata
355 : : type(tdep_Supercell_type), intent(in) :: SC
356 : :
357 : : logical :: debug
358 : : integer :: isym,jatom,mu
359 : : integer :: iatom_unitcell,jatom_unitcell,iatom
360 : 92 : integer :: vecti(3),vectj(3),vectsym(4,Sym%nsym)
361 : 46 : integer, allocatable :: indsym2(:,:,:,:)
362 : : double precision :: temp3(3)
363 : : double precision :: xred_unit_i(3),xred_unit_j(3)
364 : 92 : double precision :: xred_unit_motif(3,SC%natom_unitcell)
365 : : double precision :: origin(3)
366 : 46 : double precision, allocatable :: xred_ideal(:,:)
367 : :
368 : 46 : debug = (Invar%prtvol > 1)
369 : :
370 : 46 : write(ab_out,*) ' '
371 : 46 : write(ab_out,*) '#############################################################################'
372 : 46 : write(ab_out,*) '################ Find symmetry operations in the supercell ##################'
373 : 46 : write(ab_out,*) '#############################################################################'
374 : :
375 : : ! write(ab_out,'(a)') 'Begining SearchS'
376 : : ! TODO: Initialize the local variables at 0
377 : :
378 : : ! Positions are anchored on atom 1 (as Sym%xred_zero is, see tdep_make_sym)
379 : : ! before wrapping, since wrapping does not commute with a shift of origin.
380 : 46 : origin = SC%get_xred_scaled(1)
381 : :
382 : : ! Cache the wrapped ideal position of every atom: it is used
383 : : ! repeatedly below, in a loop over all pairs of atoms.
384 [ - + + - : 138 : ABI_MALLOC(xred_ideal,(3,SC%natom))
+ - + - -
+ - + ]
385 [ + + ]: 7324 : do iatom=1,SC%natom
386 [ + + ]: 29158 : xred_ideal(:,iatom) = SC%wrap_inbox(SC%get_xred_scaled(iatom) - origin)
387 : : end do
388 : :
389 : : ! Compute the indsym fundamental quantity
390 : : ! === Obtain a list of rotated atoms ===
391 : : ! $ R^{-1} (xred(:,iat)-\tau) = xred(:,iat_sym) + R_0 $
392 : : ! * indsym(4, isym,iat) gives iat_sym in the original unit cell.
393 : : ! * indsym(1:3,isym,iat) gives the lattice vector $R_0$.
394 : 46 : write(ab_out,'(a)') ' Search the matrix transformation going from (k) to (i)...'
395 [ - + + - : 1281382 : ABI_MALLOC(Sym%indsym,(4,Sym%nsym,SC%natom)); Sym%indsym(:,:,:)=zero
- + + - +
- + - - +
- + - + +
+ + + +
+ ]
396 : : call symatm(Sym%indsym(:,:,1:SC%natom_unitcell),SC%natom_unitcell,Sym%nsym,&
397 : 46 : & Sym%symrec,Sym%tnons,tol8,SC%typat_unitcell,Sym%xred_zero)
398 : :
399 : : ! Store the positions of the atoms in the motif, in unitcell-reduced coordinates
400 [ + + ]: 156 : do iatom=1,SC%natom_unitcell
401 : 156 : xred_unit_motif(:,iatom) = SC%scaled_to_unit(xred_ideal(:,iatom))
402 : : end do
403 : : ! Write the Indsym of the atoms included in the (reference) unitcell (i.e.: the motif)
404 [ - + - - ]: 46 : if (debug.and.MPIdata%iam_master) then
405 : 0 : write(ab_out,'(a)') ' Writing file Indsym-unitcell.dat'
406 : 0 : open(unit=40,file=trim(Invar%output_prefix)//'_Indsym-unitcell.dat')
407 [ # # ]: 0 : do iatom=1,SC%natom_unitcell
408 : 0 : write(40,*) '=========================================='
409 : 0 : write(40,'(a,i4,a,3(f10.5,1x))') 'For iatom=',iatom,' with xred (supercell)=',xred_ideal(:,iatom)
410 [ # # ]: 0 : do isym=1,Sym%nsym
411 : 0 : write(40,'(a,i2,a,i4,a,3(i4,1x),a,i2,a,3(f10.5,1x))') ' indsym(isym=',isym,',',iatom,')=',&
412 : 0 : & Sym%indsym(1:3,isym,iatom),'|iat=',Sym%indsym(4,isym,iatom),'| with tnons=',Sym%tnons(:,isym)
413 : : end do
414 : : end do
415 : 0 : close(40)
416 : : end if
417 : :
418 : : ! Search the matrix transformation going from (k,l) to (i,j)
419 : 46 : write(ab_out,'(a)') ' Search the matrix transformation going from (k,l) to (i,j)...'
420 [ - + + - : 431985816 : ABI_MALLOC(indsym2,(8,Sym%nsym,SC%natom,SC%natom)); indsym2(:,:,:,:)=0
- + + - -
+ + - + -
+ - - + -
+ + + + +
+ + + + ]
421 [ + + ]: 7324 : do iatom=1,SC%natom
422 : : ! For a single iatom
423 : 7278 : xred_unit_i = SC%scaled_to_unit(xred_ideal(:,iatom))
424 : :
425 : : ! GA: TODO This operation could be refactored using SC's logic,
426 : : ! simply storing icell.
427 : : ! and work with integers for arrays like.
428 : 7278 : iatom_unitcell=mod(iatom-1,SC%natom_unitcell)+1
429 [ + + ]: 29112 : vecti(:)=nint(xred_unit_i(:)-xred_unit_motif(:,iatom_unitcell))
430 [ + + ]: 262108 : do isym=1,Sym%nsym
431 [ + + + + ]: 72617024 : vectsym(:,:)=0
432 [ + + ]: 1019136 : do mu=1,3 ! Apply inverse transformation to original coordinates. Note transpose of symrec.
433 : 1019136 : vectsym(mu,isym) = Sym%symrec(1,mu,isym)*vecti(1)+Sym%symrec(2,mu,isym)*vecti(2)+Sym%symrec(3,mu,isym)*vecti(3)
434 : : end do
435 [ + + ]: 1273920 : Sym%indsym(1:4,isym,iatom)=Sym%indsym(1:4,isym,iatom_unitcell)+vectsym(1:4,isym)
436 [ + + ]: 48087278 : do jatom=1,SC%natom
437 [ + + ]: 239380864 : indsym2(1:4,isym,iatom,jatom)=Sym%indsym(1:4,isym,iatom_unitcell)+vectsym(1:4,isym)
438 : : end do
439 : : end do
440 : : end do
441 [ - + - - ]: 46 : if (debug.and.MPIdata%iam_master) then
442 : 0 : write(ab_out,'(a)') ' Writing file Indsym-supercell.dat'
443 : 0 : open(unit=40,file=trim(Invar%output_prefix)//'_Indsym-supercell.dat')
444 [ # # ]: 0 : do iatom=1,SC%natom
445 : 0 : write(40,*) '=========================================='
446 : 0 : write(40,'(a,i4,a,3(f10.5,1x))') 'For iatom=',iatom,' with xred (supercell)=',xred_ideal(:,iatom)
447 [ # # ]: 0 : do isym=1,Sym%nsym
448 : 0 : write(40,'(a,i2,a,i4,a,3(i4,1x),a,i2,a,3(f10.5,1x))') ' indsym(isym=',isym,',',iatom,')=',&
449 : 0 : & Sym%indsym(1:3,isym,iatom),'|iat=',Sym%indsym(4,isym,iatom),'| with tnons=',Sym%tnons(:,isym)
450 : : end do
451 : 0 : write(40,'(a,i4)') ' '
452 : : end do
453 : 0 : close(40)
454 : : end if
455 : :
456 : : ! For a couple of (iatom,jatom). The (iatom,jatom) vector depends on the position of iatom (due to PBC)
457 [ + + ]: 7324 : do iatom=1,SC%natom
458 [ + + ]: 1558688 : do jatom=1,SC%natom
459 [ + + ]: 6205456 : temp3(:)=xred_ideal(:,jatom)-xred_ideal(:,iatom)
460 : 1551364 : call tdep_make_inbox(temp3,1,1d-3,temp3)
461 [ + + ]: 6205456 : temp3(:)=xred_ideal(:,iatom)+temp3(:)
462 : 1551364 : xred_unit_j = SC%scaled_to_unit(temp3)
463 : 1551364 : jatom_unitcell=mod(jatom-1,SC%natom_unitcell)+1
464 [ + + ]: 6205456 : vectj(:)=nint(xred_unit_j(:)-xred_unit_motif(:,jatom_unitcell))
465 [ + + ]: 49383858 : do isym=1,Sym%nsym
466 [ + + + + ]: 11929973056 : vectsym(:,:)=0
467 [ + + ]: 191300864 : do mu=1,3 ! Apply inverse transformation to original coordinates. Note transpose of symrec.
468 : 191300864 : vectsym(mu,isym) = Sym%symrec(1,mu,isym)*vectj(1)+Sym%symrec(2,mu,isym)*vectj(2)+Sym%symrec(3,mu,isym)*vectj(3)
469 : : end do
470 [ + + ]: 240677444 : indsym2(5:8,isym,iatom,jatom)=Sym%indsym(1:4,isym,jatom_unitcell)+vectsym(1:4,isym)
471 : : end do
472 : : end do
473 : : end do
474 [ - + - - ]: 46 : if (debug.and.MPIdata%iam_master) then
475 : 0 : write(ab_out,'(a)') ' Writing file Indsym-2atoms.dat'
476 : 0 : open(unit=40,file=trim(Invar%output_prefix)//'_Indsym-2atoms.dat')
477 [ # # ]: 0 : do iatom=1,SC%natom
478 : 0 : write(40,*) '=========================================='
479 : 0 : write(40,'(a,i4,a,3(f10.5,1x))') 'For iatom=',iatom,' with xred (supercell)=',xred_ideal(:,iatom)
480 [ # # ]: 0 : do jatom=1,SC%natom
481 : 0 : write(40,*) ' =========================================='
482 : 0 : write(40,'(a,i4,a,3(f10.5,1x))') ' For jatom=',jatom,' with xred (supercell)=',xred_ideal(:,jatom)
483 [ # # ]: 0 : do isym=1,Sym%nsym
484 : 0 : write(40,'(a,i2,a,i4,a,i4,a,3(i4,1x),a,i2,a,3(i4,1x),a,i2,a)') ' indsym2(isym=',isym,',',iatom,',',jatom,')=',&
485 : 0 : & indsym2(1:3,isym,iatom,jatom),&
486 : 0 : & '|iat=',indsym2(4,isym,iatom,jatom),'|',indsym2(5:7,isym,iatom,jatom),'|iat=',indsym2(8,isym,iatom,jatom),'|'
487 : : end do
488 : : end do
489 : 0 : write(40,'(a,i4)') ' '
490 : : end do
491 : 0 : close(40)
492 : : end if
493 [ - + ]: 46 : ABI_FREE(indsym2)
494 [ - + ]: 46 : ABI_FREE(xred_ideal)
495 : :
496 : 46 : end subroutine tdep_SearchS_1at
497 : : !!***
498 : :
499 : : !=====================================================================================================
500 : :
501 : : !!****f* m_tdep_sym/tdep_SearchS_2at
502 : : !! NAME
503 : : !! tdep_SearchS_2at
504 : : !!
505 : : !! FUNCTION
506 : : !! Search the symmetry operation that transforms the pair of atoms (e,f)
507 : : !! into the reference pair of atoms (i,j), i.e. such that e and f are
508 : : !! respectively equivalent to i and j, and that the lattice vectors
509 : : !! connecting them are consistent (same bond, modulo a lattice translation).
510 : : !! The first symmetry operation found is retained.
511 : : !! Only the direct transformation is searched here: the reversed bond
512 : : !! (j,i) is handled by the caller (see tdep_init_shell2at).
513 : : !!
514 : : !! INPUTS
515 : : !! SC = Supercell object describing the ideal structure.
516 : : !! iatom, jatom = Atoms of the reference pair (iatom in the unitcell).
517 : : !! eatom, fatom = Atoms of the pair to be transformed.
518 : : !! Sym = Symmetries object describing all the symmetry operations of the crystal.
519 : : !!
520 : : !! OUTPUT
521 : : !!
522 : : !! SIDE EFFECTS
523 : : !! Isym2at(natom,natom,2) = If a symmetry operation isym is found,
524 : : !! Isym2at(eatom,fatom,1)=isym and
525 : : !! Isym2at(eatom,fatom,2)=1 (direct transformation).
526 : : !! Unchanged otherwise.
527 : : !!
528 : : !! SOURCE
529 : :
530 : 330346 : subroutine tdep_SearchS_2at(SC,iatom,jatom,eatom,fatom,Isym2at,Sym)
531 : :
532 : : integer, intent(in) :: iatom,jatom,eatom,fatom
533 : : type(tdep_Supercell_type),intent(in) :: SC
534 : : type(Symmetries_type),intent(in) :: Sym
535 : : integer, intent(inout) :: Isym2at(SC%natom,SC%natom,2)
536 : :
537 : : integer :: isym,ee,ff,ii
538 : : integer :: iatom_unitcell,jatom_unitcell
539 : : integer :: vecti(3),vectj(3),latte(3),lattf(3),indsym2(8)
540 : : double precision :: xred_i(3),xred_j(3),temp(3)
541 : : double precision :: xred_unit_i(3),xred_unit_j(3)
542 : 660692 : double precision :: xred_unit_motif(3,SC%natom_unitcell)
543 : : double precision :: origin(3)
544 : : logical :: ok
545 : :
546 : : ! Search the couple of atoms (e,f) obtained through the transformation (R+t)
547 : : ! starting from (i,j)
548 : : ! In that case, note that two cases are possible:
549 : : ! - The bond is just transformed, so indsym(4,e)=i and indsym(4,f)=j
550 : : ! - The bond is transformed and reversed, so indsym(4,e)=j and indsym(4,f)=i
551 : :
552 : : ! Positions are anchored on atom 1 (as Sym%xred_zero is, see tdep_make_sym)
553 : : ! before wrapping, since wrapping does not commute with a shift of origin.
554 : 330346 : origin = SC%get_xred_scaled(1)
555 : :
556 : : ! Store the positions of the atoms in the motif, in unitcell-reduced coordinates
557 [ + + ]: 1002474 : do ii=1,SC%natom_unitcell
558 [ + + ]: 3018858 : xred_unit_motif(:,ii) = SC%scaled_to_unit(SC%wrap_inbox(SC%get_xred_scaled(ii) - origin))
559 : : end do
560 : :
561 : : ! Search the atom equivalent to iatom in the (reference) unitcell
562 : : ! Note that in the (present) particular case : iatom_unitcell=iatom and vecti(:)=zero
563 [ + + ]: 1321384 : xred_i = SC%wrap_inbox(SC%get_xred_scaled(iatom) - origin)
564 : 330346 : xred_unit_i = SC%scaled_to_unit(xred_i)
565 : 330346 : iatom_unitcell=mod(iatom-1,SC%natom_unitcell)+1
566 [ + + ]: 1321384 : vecti(:)=nint(xred_unit_i(:)-xred_unit_motif(:,iatom_unitcell))
567 : :
568 : : ! Search the atom equivalent to jatom in the (reference) unitcell.
569 : : ! jatom can be outside the box (according to the value of iatom)
570 : : ! so we use the inbox procedure to put the distance within [-0.5,0.5[
571 [ + + ]: 1321384 : xred_j = SC%wrap_inbox(SC%get_xred_scaled(jatom) - origin)
572 [ + + ]: 1321384 : temp(:)=xred_j(:)-xred_i(:)
573 : 330346 : call tdep_make_inbox(temp,1,1d-3,temp)
574 [ + + ]: 1321384 : temp(:)=xred_i(:)+temp(:)
575 : 330346 : xred_unit_j = SC%scaled_to_unit(temp)
576 : 330346 : jatom_unitcell=mod(jatom-1,SC%natom_unitcell)+1
577 [ + + ]: 1321384 : vectj(:)=nint(xred_unit_j(:)-xred_unit_motif(:,jatom_unitcell))
578 : :
579 : : ! To understand the meaning of "latt", see SearchS_1at
580 : 330346 : ok=.false.
581 [ + + ]: 5722998 : do isym=1,Sym%nsym
582 : 5666528 : indsym2(:)=0
583 : 5666528 : call tdep_calc_indsym2(SC,eatom,fatom,indsym2,isym,Sym)
584 : 5666528 : ee=indsym2(4)
585 : 5666528 : ff=indsym2(8)
586 [ + + ]: 22666112 : latte(:)=indsym2(1:3)
587 [ + + ]: 22666112 : lattf(:)=indsym2(5:7)
588 [ + + ]: 5666528 : if (ee==iatom_unitcell.and.ff==jatom_unitcell) then
589 [ + + + + ]: 11854000 : if (sum(abs((vecti(:)-latte(:))-(vectj(:)-lattf(:)))).lt.tol8) then
590 : 273876 : Isym2at(eatom,fatom,1)=isym
591 : 273876 : Isym2at(eatom,fatom,2)=1
592 : : ok=.true.
593 : : end if
594 : : end if
595 : 56470 : if (ok) exit
596 : : end do
597 : :
598 : 330346 : end subroutine tdep_SearchS_2at
599 : : !!***
600 : :
601 : : !=====================================================================================================
602 : :
603 : : !!****f* m_tdep_sym/tdep_SearchS_3at
604 : : !! NAME
605 : : !! tdep_SearchS_3at
606 : : !!
607 : : !! FUNCTION
608 : : !! Search the symmetry operation that transforms the triplet of atoms
609 : : !! (e,f,g) into the reference triplet of atoms (i,j,k), i.e. such that
610 : : !! e, f and g are respectively equivalent to i, j and k, and that the
611 : : !! lattice vectors connecting each pair of atoms are consistent.
612 : : !! The first symmetry operation found is retained.
613 : : !! The permutations of the triplet are handled by the caller
614 : : !! (see tdep_init_shell3at).
615 : : !!
616 : : !! INPUTS
617 : : !! SC = Supercell object describing the ideal structure.
618 : : !! iatom, jatom, katom = Atoms of the reference triplet (iatom in the unitcell).
619 : : !! eatom, fatom, gatom = Atoms of the triplet to be transformed.
620 : : !! Sym = Symmetries object describing all the symmetry operations of the crystal.
621 : : !!
622 : : !! OUTPUT
623 : : !!
624 : : !! SIDE EFFECTS
625 : : !! Isym3at(2) = If a symmetry operation isym is found, Isym3at(1)=isym
626 : : !! and Isym3at(2)=1. Unchanged otherwise.
627 : : !!
628 : : !! SOURCE
629 : :
630 : 43178 : subroutine tdep_SearchS_3at(SC,iatom,jatom,katom,eatom,fatom,gatom,Isym3at,Sym)
631 : :
632 : : integer, intent(in) :: iatom,jatom,katom,eatom,fatom,gatom
633 : : type(tdep_Supercell_type),intent(in) :: SC
634 : : type(Symmetries_type),intent(in) :: Sym
635 : : integer, intent(inout) :: Isym3at(2)
636 : :
637 : : integer :: isym,ee,ff,gg,ii
638 : : integer :: iatom_unitcell,jatom_unitcell,katom_unitcell
639 : : integer :: vecti(3),vectj(3),vectk(3),indsym2(8)
640 : : integer :: lattef(3),lattfg(3),lattge(3),lattfe(3),lattgf(3),latteg(3)
641 : : double precision :: xred_i(3),xred_j(3),xred_k(3),temp(3)
642 : : double precision :: xred_unit_i(3),xred_unit_j(3),xred_unit_k(3)
643 : 86356 : double precision :: xred_unit_motif(3,SC%natom_unitcell)
644 : : double precision :: origin(3)
645 : : logical :: ok
646 : :
647 : : ! Search the couple of atoms (e,f) obtained through the transformation (R+t)
648 : : ! starting from (i,j)
649 : : ! In that case, note that two cases are possible:
650 : : ! - The bond is just transformed, so indsym(4,e)=i and indsym(4,f)=j
651 : : ! - The bond is transformed and reversed, so indsym(4,e)=j and indsym(4,f)=i
652 : :
653 : : ! Positions are anchored on atom 1 (as Sym%xred_zero is, see tdep_make_sym)
654 : : ! before wrapping, since wrapping does not commute with a shift of origin.
655 : 43178 : origin = SC%get_xred_scaled(1)
656 : :
657 : : ! Store the positions of the atoms in the motif, in unitcell-reduced coordinates
658 [ + + ]: 143022 : do ii=1,SC%natom_unitcell
659 [ + + ]: 442554 : xred_unit_motif(:,ii) = SC%scaled_to_unit(SC%wrap_inbox(SC%get_xred_scaled(ii) - origin))
660 : : end do
661 : :
662 : : ! Search the atom equivalent to iatom in the (reference) unitcell
663 : : ! Note that in the (present) particular case : iatom_unitcell=iatom and vecti(:)=zero
664 [ + + ]: 172712 : xred_i = SC%wrap_inbox(SC%get_xred_scaled(iatom) - origin)
665 : 43178 : xred_unit_i = SC%scaled_to_unit(xred_i)
666 : 43178 : iatom_unitcell=mod(iatom-1,SC%natom_unitcell)+1
667 [ + + ]: 172712 : vecti(:)=nint(xred_unit_i(:)-xred_unit_motif(:,iatom_unitcell))
668 : :
669 : : ! Search the atom equivalent to jatom in the (reference) unitcell.
670 : : ! jatom can be outside the box (according to the value of iatom)
671 : : ! so we use the inbox procedure to put the distance within [-0.5,0.5[
672 [ + + ]: 172712 : xred_j = SC%wrap_inbox(SC%get_xred_scaled(jatom) - origin)
673 [ + + ]: 172712 : temp(:)=xred_j(:)-xred_i(:)
674 : 43178 : call tdep_make_inbox(temp,1,1d-3,temp)
675 [ + + ]: 172712 : temp(:)=xred_i(:)+temp(:)
676 : 43178 : xred_unit_j = SC%scaled_to_unit(temp)
677 : 43178 : jatom_unitcell=mod(jatom-1,SC%natom_unitcell)+1
678 [ + + ]: 172712 : vectj(:)=nint(xred_unit_j(:)-xred_unit_motif(:,jatom_unitcell))
679 : :
680 : : ! Search the atom equivalent to katom in the (reference) unitcell.
681 : : ! katom can be outside the box (according to the value of iatom)
682 : : ! so we use the inbox procedure to put the distance within [-0.5,0.5[
683 [ + + ]: 172712 : xred_k = SC%wrap_inbox(SC%get_xred_scaled(katom) - origin)
684 [ + + ]: 172712 : temp(:)=xred_k(:)-xred_i(:)
685 : 43178 : call tdep_make_inbox(temp,1,1d-3,temp)
686 [ + + ]: 172712 : temp(:)=xred_i(:)+temp(:)
687 : 43178 : xred_unit_k = SC%scaled_to_unit(temp)
688 : 43178 : katom_unitcell=mod(katom-1,SC%natom_unitcell)+1
689 [ + + ]: 172712 : vectk(:)=nint(xred_unit_k(:)-xred_unit_motif(:,katom_unitcell))
690 : :
691 : : ! To understand the meaning of "latt", see SearchS_1at
692 : 43178 : ok=.false.
693 [ + + ]: 975613 : do isym=1,Sym%nsym
694 : : ! TODO : A CHECKER!!!!!!!!!!!!!!!!!
695 : 959861 : indsym2(:)=zero
696 : 959861 : call tdep_calc_indsym2(SC,eatom,fatom,indsym2,isym,Sym)
697 : 959861 : ee=indsym2(4)
698 [ + + ]: 3839444 : lattef(:)=indsym2(1:3)
699 [ + + ]: 3839444 : lattfe(:)=indsym2(5:7)
700 : :
701 : 959861 : indsym2(:)=zero
702 : 959861 : call tdep_calc_indsym2(SC,fatom,gatom,indsym2,isym,Sym)
703 : 959861 : ff=indsym2(4)
704 [ + + ]: 3839444 : lattfg(:)=indsym2(1:3)
705 [ + + ]: 3839444 : lattgf(:)=indsym2(5:7)
706 : :
707 : 959861 : indsym2(:)=zero
708 : 959861 : call tdep_calc_indsym2(SC,gatom,eatom,indsym2,isym,Sym)
709 : 959861 : gg=indsym2(4)
710 [ + + ]: 3839444 : lattge(:)=indsym2(1:3)
711 [ + + ]: 3839444 : latteg(:)=indsym2(5:7)
712 [ + + + + ]: 959861 : if (ee==iatom_unitcell.and.ff==jatom_unitcell.and.gg==katom_unitcell) then
713 : : if ((sum(abs((vecti(:)-lattef(:))-(vectj(:)-lattfe(:)))).lt.tol8).and.&
714 [ + + + + : 1963720 : & (sum(abs((vectj(:)-lattfg(:))-(vectk(:)-lattgf(:)))).lt.tol8).and.&
+ + + + +
- ]
715 : : & (sum(abs((vectk(:)-lattge(:))-(vecti(:)-latteg(:)))).lt.tol8)) then
716 : 27426 : Isym3at(1)=isym
717 : 27426 : Isym3at(2)=1
718 : : ok=.true.
719 : : end if
720 : : end if
721 : 15752 : if (ok) exit
722 : : end do
723 : :
724 : 43178 : end subroutine tdep_SearchS_3at
725 : : !!***
726 : :
727 : : !=====================================================================================================
728 : :
729 : : !!****f* m_tdep_sym/tdep_SearchS_4at
730 : : !! NAME
731 : : !! tdep_SearchS_4at
732 : : !!
733 : : !! FUNCTION
734 : : !! Search the symmetry operation that transforms the quadruplet of atoms
735 : : !! (e,f,g,h) into the reference quadruplet of atoms (i,j,k,l), i.e. such
736 : : !! that e, f, g and h are respectively equivalent to i, j, k and l, and
737 : : !! that the lattice vectors connecting each pair of atoms are consistent.
738 : : !! The first symmetry operation found is retained.
739 : : !! The permutations of the quadruplet are handled by the caller
740 : : !! (see tdep_init_shell4at).
741 : : !!
742 : : !! INPUTS
743 : : !! SC = Supercell object describing the ideal structure.
744 : : !! iatom, jatom, katom, latom = Atoms of the reference quadruplet
745 : : !! (iatom in the unitcell).
746 : : !! eatom, fatom, gatom, hatom = Atoms of the quadruplet to be transformed.
747 : : !! Sym = Symmetries object describing all the symmetry operations of the crystal.
748 : : !!
749 : : !! OUTPUT
750 : : !!
751 : : !! SIDE EFFECTS
752 : : !! Isym4at(2) = If a symmetry operation isym is found, Isym4at(1)=isym
753 : : !! and Isym4at(2)=1. Unchanged otherwise.
754 : : !!
755 : : !! SOURCE
756 : :
757 : 68388 : subroutine tdep_SearchS_4at(SC,iatom,jatom,katom,latom,eatom,fatom,gatom,hatom,Isym4at,Sym)
758 : :
759 : : integer, intent(in) :: iatom,jatom,katom,latom,eatom,fatom,gatom,hatom
760 : : type(tdep_Supercell_type),intent(in) :: SC
761 : : type(Symmetries_type),intent(in) :: Sym
762 : : integer, intent(inout) :: Isym4at(2)
763 : :
764 : : integer :: isym,ee,ff,gg,hh,ii
765 : : integer :: iatom_unitcell,jatom_unitcell,katom_unitcell,latom_unitcell
766 : : integer :: vecti(3),vectj(3),vectk(3),vectl(3),indsym2(8)
767 : : integer :: lattef(3),lattfe(3)
768 : : integer :: latteg(3),lattge(3)
769 : : integer :: latteh(3),latthe(3)
770 : : integer :: lattfg(3),lattgf(3)
771 : : integer :: lattfh(3),latthf(3)
772 : : integer :: lattgh(3),latthg(3)
773 : : double precision :: xred_i(3),xred_j(3),xred_k(3),xred_l(3),temp(3)
774 : : double precision :: xred_unit_i(3),xred_unit_j(3),xred_unit_k(3),xred_unit_l(3)
775 : 136776 : double precision :: xred_unit_motif(3,SC%natom_unitcell)
776 : : double precision :: origin(3)
777 : : logical :: ok
778 : :
779 : : ! Search the couple of atoms (e,f) obtained through the transformation (R+t)
780 : : ! starting from (i,j)
781 : : ! In that case, note that two cases are possible:
782 : : ! - The bond is just transformed, so indsym(4,e)=i and indsym(4,f)=j
783 : : ! - The bond is transformed and reversed, so indsym(4,e)=j and indsym(4,f)=i
784 : :
785 : : ! Positions are anchored on atom 1 (as Sym%xred_zero is, see tdep_make_sym)
786 : : ! before wrapping, since wrapping does not commute with a shift of origin.
787 : 68388 : origin = SC%get_xred_scaled(1)
788 : :
789 : : ! Store the positions of the atoms in the motif, in unitcell-reduced coordinates
790 [ + + ]: 230744 : do ii=1,SC%natom_unitcell
791 [ + + ]: 717812 : xred_unit_motif(:,ii) = SC%scaled_to_unit(SC%wrap_inbox(SC%get_xred_scaled(ii) - origin))
792 : : end do
793 : :
794 : : ! Search the atom equivalent to iatom in the (reference) unitcell
795 : : ! Note that in the (present) particular case : iatom_unitcell=iatom and vecti(:)=zero
796 [ + + ]: 273552 : xred_i = SC%wrap_inbox(SC%get_xred_scaled(iatom) - origin)
797 : 68388 : xred_unit_i = SC%scaled_to_unit(xred_i)
798 : 68388 : iatom_unitcell=mod(iatom-1,SC%natom_unitcell)+1
799 [ + + ]: 273552 : vecti(:)=nint(xred_unit_i(:)-xred_unit_motif(:,iatom_unitcell))
800 : :
801 : : ! Search the atom equivalent to jatom in the (reference) unitcell.
802 : : ! jatom can be outside the box (according to the value of iatom)
803 : : ! so we use the inbox procedure to put the distance within [-0.5,0.5[
804 [ + + ]: 273552 : xred_j = SC%wrap_inbox(SC%get_xred_scaled(jatom) - origin)
805 [ + + ]: 273552 : temp(:)=xred_j(:)-xred_i(:)
806 : 68388 : call tdep_make_inbox(temp,1,1d-3,temp)
807 [ + + ]: 273552 : temp(:)=xred_i(:)+temp(:)
808 : 68388 : xred_unit_j = SC%scaled_to_unit(temp)
809 : 68388 : jatom_unitcell=mod(jatom-1,SC%natom_unitcell)+1
810 [ + + ]: 273552 : vectj(:)=nint(xred_unit_j(:)-xred_unit_motif(:,jatom_unitcell))
811 : :
812 : : ! Search the atom equivalent to katom in the (reference) unitcell.
813 : : ! katom can be outside the box (according to the value of iatom)
814 : : ! so we use the inbox procedure to put the distance within [-0.5,0.5[
815 [ + + ]: 273552 : xred_k = SC%wrap_inbox(SC%get_xred_scaled(katom) - origin)
816 [ + + ]: 273552 : temp(:)=xred_k(:)-xred_i(:)
817 : 68388 : call tdep_make_inbox(temp,1,1d-3,temp)
818 [ + + ]: 273552 : temp(:)=xred_i(:)+temp(:)
819 : 68388 : xred_unit_k = SC%scaled_to_unit(temp)
820 : 68388 : katom_unitcell=mod(katom-1,SC%natom_unitcell)+1
821 [ + + ]: 273552 : vectk(:)=nint(xred_unit_k(:)-xred_unit_motif(:,katom_unitcell))
822 : :
823 : : ! Search the atom equivalent to latom in the (reference) unitcell.
824 : : ! latom can be outside the box (according to the value of iatom)
825 : : ! so we use the inbox procedure to put the distance within [-0.5,0.5[
826 [ + + ]: 273552 : xred_l = SC%wrap_inbox(SC%get_xred_scaled(latom) - origin)
827 [ + + ]: 273552 : temp(:)=xred_l(:)-xred_i(:)
828 : 68388 : call tdep_make_inbox(temp,1,1d-3,temp)
829 [ + + ]: 273552 : temp(:)=xred_i(:)+temp(:)
830 : 68388 : xred_unit_l = SC%scaled_to_unit(temp)
831 : 68388 : latom_unitcell=mod(latom-1,SC%natom_unitcell)+1
832 [ + + ]: 273552 : vectl(:)=nint(xred_unit_l(:)-xred_unit_motif(:,latom_unitcell))
833 : :
834 : : ! To understand the meaning of "latt", see SearchS_1at
835 : 68388 : ok=.false.
836 [ + + ]: 2086763 : do isym=1,Sym%nsym
837 : : ! TODO : A CHECKER!!!!!!!!!!!!!!!!!
838 : 2050751 : indsym2(:)=zero
839 : 2050751 : call tdep_calc_indsym2(SC,eatom,fatom,indsym2,isym,Sym)
840 : 2050751 : ee=indsym2(4)
841 : 2050751 : ff=indsym2(8)
842 [ + + ]: 8203004 : lattef(:)=indsym2(1:3)
843 [ + + ]: 8203004 : lattfe(:)=indsym2(5:7)
844 : :
845 : 2050751 : indsym2(:)=zero
846 : 2050751 : call tdep_calc_indsym2(SC,eatom,gatom,indsym2,isym,Sym)
847 : 2050751 : gg=indsym2(8)
848 [ + + ]: 8203004 : latteg(:)=indsym2(1:3)
849 [ + + ]: 8203004 : lattge(:)=indsym2(5:7)
850 : :
851 : 2050751 : indsym2(:)=zero
852 : 2050751 : call tdep_calc_indsym2(SC,eatom,hatom,indsym2,isym,Sym)
853 : 2050751 : hh=indsym2(8)
854 [ + + ]: 8203004 : latteh(:)=indsym2(1:3)
855 [ + + ]: 8203004 : latthe(:)=indsym2(5:7)
856 : :
857 : 2050751 : indsym2(:)=zero
858 : 2050751 : call tdep_calc_indsym2(SC,fatom,gatom,indsym2,isym,Sym)
859 [ + + ]: 8203004 : lattfg(:)=indsym2(1:3)
860 [ + + ]: 8203004 : lattgf(:)=indsym2(5:7)
861 : :
862 : 2050751 : indsym2(:)=zero
863 : 2050751 : call tdep_calc_indsym2(SC,fatom,hatom,indsym2,isym,Sym)
864 [ + + ]: 8203004 : lattfh(:)=indsym2(1:3)
865 [ + + ]: 8203004 : latthf(:)=indsym2(5:7)
866 : :
867 : 2050751 : indsym2(:)=zero
868 : 2050751 : call tdep_calc_indsym2(SC,gatom,hatom,indsym2,isym,Sym)
869 [ + + ]: 8203004 : lattgh(:)=indsym2(1:3)
870 [ + + ]: 8203004 : latthg(:)=indsym2(5:7)
871 : :
872 [ + + + + ]: 2050751 : if (ee==iatom_unitcell.and.ff==jatom_unitcell.and.gg==katom_unitcell.and.hh==latom_unitcell) then
873 : : if ((sum(abs((vecti(:)-lattef(:))-(vectj(:)-lattfe(:)))).lt.tol8).and.&
874 : : & (sum(abs((vecti(:)-latteg(:))-(vectk(:)-lattge(:)))).lt.tol8).and.&
875 : : & (sum(abs((vecti(:)-latteh(:))-(vectl(:)-latthe(:)))).lt.tol8).and.&
876 : : & (sum(abs((vectj(:)-lattfg(:))-(vectk(:)-lattgf(:)))).lt.tol8).and.&
877 [ + + + + : 5374606 : & (sum(abs((vectj(:)-lattfh(:))-(vectl(:)-latthf(:)))).lt.tol8).and.&
+ + + + +
+ + + + +
+ + + - ]
878 : : & (sum(abs((vectk(:)-lattgh(:))-(vectl(:)-latthg(:)))).lt.tol8)) then
879 : 32376 : Isym4at(1)=isym
880 : 32376 : Isym4at(2)=1
881 : : ok=.true.
882 : : end if
883 : : end if
884 : 36012 : if (ok) exit
885 : : end do
886 : :
887 : 68388 : end subroutine tdep_SearchS_4at
888 : : !!***
889 : :
890 : : !=====================================================================================================
891 : :
892 : : !!****f* m_tdep_sym/tdep_calc_indsym2
893 : : !! NAME
894 : : !! tdep_calc_indsym2
895 : : !!
896 : : !! FUNCTION
897 : : !! Compute, for a pair of atoms (iatom,jatom) of the supercell and a given
898 : : !! symmetry operation, the equivalent atoms in the reference unitcell and
899 : : !! the corresponding lattice vectors (as in Sym%indsym, see
900 : : !! tdep_SearchS_1at). The position of jatom is taken as its periodic image
901 : : !! closest to iatom, so that the (iatom,jatom) bond is preserved.
902 : : !!
903 : : !! INPUTS
904 : : !! SC = Supercell object describing the ideal structure.
905 : : !! iatom, jatom = Atoms of the pair.
906 : : !! isym = Index of the symmetry operation.
907 : : !! Sym = Symmetries object. Sym%indsym must already be computed
908 : : !! (see tdep_SearchS_1at).
909 : : !!
910 : : !! OUTPUT
911 : : !! indsym2(8) = indsym2(1:3) and indsym2(4) are the lattice vector and the
912 : : !! equivalent unitcell atom for iatom; indsym2(5:7) and
913 : : !! indsym2(8) are the same quantities for jatom.
914 : : !!
915 : : !! SOURCE
916 : :
917 : 20850617 : subroutine tdep_calc_indsym2(SC,iatom,jatom,indsym2,isym,Sym)
918 : :
919 : : integer, intent(in) :: iatom,jatom
920 : : type(tdep_Supercell_type),intent(in) :: SC
921 : : type(Symmetries_type),intent(in) :: Sym
922 : : integer, intent(out) :: indsym2(8)
923 : :
924 : : integer :: isym,mu,ii
925 : : integer :: iatom_unitcell,jatom_unitcell
926 : 41701234 : integer :: vecti(3),vectj(3),vectsym(4,Sym%nsym)
927 : : double precision :: xred_i(3),xred_j(3),temp3(3)
928 : : double precision :: xred_unit_i(3),xred_unit_j(3)
929 : 41701234 : double precision :: xred_unit_motif(3,SC%natom_unitcell)
930 : : double precision :: origin(3)
931 : :
932 : : ! Positions are anchored on atom 1 (as Sym%xred_zero is, see tdep_make_sym)
933 : : ! before wrapping, since wrapping does not commute with a shift of origin.
934 : 20850617 : origin = SC%get_xred_scaled(1)
935 : :
936 : : ! Store the positions of the atoms in the motif, in unitcell-reduced coordinates
937 [ + + ]: 69177220 : do ii=1,SC%natom_unitcell
938 [ + + ]: 214157029 : xred_unit_motif(:,ii) = SC%scaled_to_unit(SC%wrap_inbox(SC%get_xred_scaled(ii) - origin))
939 : : end do
940 : :
941 : : ! For a single iatom
942 [ + + ]: 83402468 : xred_i = SC%wrap_inbox(SC%get_xred_scaled(iatom) - origin)
943 : 20850617 : xred_unit_i = SC%scaled_to_unit(xred_i)
944 : 20850617 : iatom_unitcell=mod(iatom-1,SC%natom_unitcell)+1
945 [ + + ]: 83402468 : vecti(:)=nint(xred_unit_i(:)-xred_unit_motif(:,iatom_unitcell))
946 [ + + + + ]: 4850947337 : vectsym(:,:)=0
947 [ + + ]: 83402468 : do mu=1,3 ! Apply inverse transformation to original coordinates. Note transpose of symrec.
948 : 83402468 : vectsym(mu,isym) = Sym%symrec(1,mu,isym)*vecti(1)+Sym%symrec(2,mu,isym)*vecti(2)+Sym%symrec(3,mu,isym)*vecti(3)
949 : : end do
950 [ + + ]: 104253085 : indsym2(1:4)=Sym%indsym(1:4,isym,iatom_unitcell)+vectsym(1:4,isym)
951 : :
952 : : ! For a couple of (iatom,jatom). The (iatom,jatom) vector depends on the position of iatom (due to PBC)
953 [ + + ]: 83402468 : xred_j = SC%wrap_inbox(SC%get_xred_scaled(jatom) - origin)
954 [ + + ]: 83402468 : temp3(:)=xred_j(:)-xred_i(:)
955 : 20850617 : call tdep_make_inbox(temp3,1,1d-3,temp3)
956 [ + + ]: 83402468 : temp3(:)=xred_i(:)+temp3(:)
957 : 20850617 : xred_unit_j = SC%scaled_to_unit(temp3)
958 : 20850617 : jatom_unitcell=mod(jatom-1,SC%natom_unitcell)+1
959 [ + + ]: 83402468 : vectj(:)=nint(xred_unit_j(:)-xred_unit_motif(:,jatom_unitcell))
960 [ + + + + ]: 4850947337 : vectsym(:,:)=0
961 [ + + ]: 83402468 : do mu=1,3 ! Apply inverse transformation to original coordinates. Note transpose of symrec.
962 : 83402468 : vectsym(mu,isym) = Sym%symrec(1,mu,isym)*vectj(1)+Sym%symrec(2,mu,isym)*vectj(2)+Sym%symrec(3,mu,isym)*vectj(3)
963 : : end do
964 [ + + ]: 104253085 : indsym2(5:8)=Sym%indsym(1:4,isym,jatom_unitcell)+vectsym(1:4,isym)
965 : :
966 : 20850617 : end subroutine tdep_calc_indsym2
967 : : !!***
968 : :
969 : : !=====================================================================================================
970 : :
971 : : !!****f* m_tdep_sym/tdep_destroy_sym
972 : : !! NAME
973 : : !! tdep_destroy_sym
974 : : !!
975 : : !! FUNCTION
976 : : !! Free the memory allocated in the Symmetries object.
977 : : !!
978 : : !! SIDE EFFECTS
979 : : !! Sym = Symmetries object. All its allocated arrays are deallocated.
980 : : !!
981 : : !! SOURCE
982 : :
983 : 46 : subroutine tdep_destroy_sym(Sym)
984 : :
985 : : class(Symmetries_type),intent(inout) :: Sym
986 : :
987 [ - + ]: 46 : ABI_FREE(Sym%ptsymrel)
988 [ - + ]: 46 : ABI_FREE(Sym%S_ref)
989 [ - + ]: 46 : ABI_FREE(Sym%S_inv)
990 [ - + ]: 46 : ABI_FREE(Sym%xred_zero)
991 [ - + ]: 46 : ABI_FREE(Sym%tnons)
992 [ - + ]: 46 : ABI_FREE(Sym%symafm)
993 [ - + ]: 46 : ABI_FREE(Sym%symrec)
994 [ - + ]: 46 : ABI_FREE(Sym%symrel)
995 [ - + ]: 46 : ABI_FREE(Sym%indsym)
996 : :
997 : 46 : end subroutine tdep_destroy_sym
998 : : !!***
999 : :
1000 : : !=====================================================================================================
1001 : :
1002 [ - + - - : 138 : end module m_tdep_sym
+ + - + +
- - + - +
+ - - + -
+ - + - +
- + - + -
+ ]
1003 : : !!***
|
