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