Branch data Line data Source code
1 : : !!****m* ABINIT/m_tdep_shell
2 : : !! NAME
3 : : !! m_tdep_shell
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains the TDEP Shell data type, which holds the shells
7 : : !! of interactions (sets of 1, 2, 3 or 4 atoms that are equivalent by
8 : : !! symmetry) and, for each shell, the projector onto the independent
9 : : !! IFC coefficients.
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_shell
26 : :
27 : : use defs_basis
28 : : use m_errors
29 : : use m_abicore
30 : : use m_xmpi
31 : : use m_io_tools
32 : : use m_numeric_tools, only : uniformrandom
33 : : use m_sort, only : sort_dp
34 : : use m_tdep_dataset, only : atdep_dataset_type, MPI_enreg_type
35 : : use m_tdep_sym, only : Symmetries_type, tdep_SearchS_2at, tdep_SearchS_3at, tdep_SearchS_4at
36 : : use m_tdep_supercell, only : tdep_Supercell_type
37 : : !!***
38 : :
39 : : !----------------------------------------------------------------------
40 : :
41 : : !!****t* m_tdep_shell/List_of_neighbours
42 : : !! NAME
43 : : !! List_of_neighbours
44 : : !!
45 : : !! FUNCTION
46 : : !! For a given atom (i) and a given shell, list of the interactions
47 : : !! (i,j), (i,j,k) or (i,j,k,l) belonging to this shell, together with
48 : : !! the symmetry operation (and the permutation) that transforms the
49 : : !! reference interaction of the shell into each of them.
50 : : !!
51 : : !! SOURCE
52 : :
53 : : type List_of_neighbours
54 : : integer :: n_interactions
55 : : integer, allocatable :: atomj_in_shell(:)
56 : : integer, allocatable :: atomk_in_shell(:)
57 : : integer, allocatable :: atoml_in_shell(:)
58 : : integer, allocatable :: sym_in_shell(:)
59 : : integer, allocatable :: transpose_in_shell(:)
60 : : end type List_of_neighbours
61 : : !!***
62 : :
63 : : !----------------------------------------------------------------------
64 : :
65 : : !!****t* m_tdep_shell/Shell_type
66 : : !! NAME
67 : : !! Shell_type
68 : : !!
69 : : !! FUNCTION
70 : : !! Holds the shells of interactions at a given order, the reference
71 : : !! atoms of each shell, the list of interactions of each atom in each
72 : : !! shell, and the projectors onto the independent IFC coefficients.
73 : : !!
74 : : !! SOURCE
75 : :
76 : : type Shell_type
77 : :
78 : : integer :: order
79 : : ! Order of interaction (1, 2, 3, or 4)
80 : :
81 : : integer :: norder
82 : : ! Dimension of IFC matrix at that order (3**order)
83 : :
84 : : integer :: natom
85 : : ! Number of atoms
86 : :
87 : : integer :: ntotcoeff
88 : : ! Total number of coefficients
89 : :
90 : : integer :: nshell
91 : : ! Number of shells
92 : :
93 : : integer, allocatable :: ncoeff(:)
94 : : ! ncoeff(nshell)
95 : : ! Number of coefficients in each shell
96 : :
97 : : integer, allocatable :: ncoeff_prev(:)
98 : : ! ncoeff_prev(nshell)
99 : :
100 : : integer, allocatable :: ishell_self(:)
101 : : ! ncoeff_prev(natom_unitcell)
102 : :
103 : : integer, allocatable :: iatref(:)
104 : : ! iatref(nshell)
105 : :
106 : : integer, allocatable :: jatref(:)
107 : : ! jatref(nshell)
108 : :
109 : : integer, allocatable :: katref(:)
110 : : ! katref(nshell)
111 : :
112 : : integer, allocatable :: latref(:)
113 : : ! latref(nshell)
114 : :
115 : : double precision, allocatable :: proj(:,:,:)
116 : : ! proj(norder, norder, nshell)
117 : : ! Projector onto the subset of non-zero coefficients
118 : :
119 : : type(List_of_neighbours),allocatable :: neighbours(:,:)
120 : : ! neighbours(natom, nshell)
121 : :
122 : : contains
123 : :
124 : : procedure :: init => tdep_init_shell
125 : : procedure :: free => tdep_free_shell
126 : :
127 : : end type Shell_type
128 : : !!***
129 : :
130 : : public :: tdep_init_shell1at
131 : : public :: tdep_init_shell2at
132 : : public :: tdep_init_shell3at
133 : : public :: tdep_init_shell4at
134 : : public :: tdep_destroy_shell
135 : : public :: tdep_calc_nbcoeff
136 : :
137 : : contains
138 : :
139 : : !====================================================================================================
140 : :
141 : : !!****f* ABINIT/m_tdep_shell/tdep_init_shell
142 : : !! NAME
143 : : !! tdep_init_shell
144 : : !!
145 : : !! FUNCTION
146 : : !! Build the shells of interactions at a given order, by calling
147 : : !! tdep_init_shell1at, tdep_init_shell2at, tdep_init_shell3at or
148 : : !! tdep_init_shell4at.
149 : : !!
150 : : !! INPUTS
151 : : !! Invar = Input object containing the input variables.
152 : : !! MPIdata = Info on MPI parallelism.
153 : : !! SC = Supercell object describing the ideal structure.
154 : : !! order = Order of the interactions (1, 2, 3 or 4).
155 : : !!
156 : : !! OUTPUT
157 : : !!
158 : : !! SIDE EFFECTS
159 : : !! Shell = Shell object, built at the requested order.
160 : : !! Sym = Symmetries object describing all the symmetry operations of the crystal.
161 : : !!
162 : : !! SOURCE
163 : :
164 : 117 : subroutine tdep_init_shell(Shell,Invar,MPIdata,SC,Sym,order)
165 : : class(Shell_type), intent(inout) :: Shell
166 : : type(atdep_dataset_type),intent(in) :: Invar
167 : : type(MPI_enreg_type), intent(in) :: MPIdata
168 : : type(tdep_Supercell_type),intent(in) :: SC
169 : : type(Symmetries_type),intent(inout) :: Sym
170 : : integer, intent(in) :: order
171 : :
172 : 163 : SELECT CASE (order)
173 : : CASE (1)
174 : 46 : call tdep_init_shell1at(Invar,MPIdata,Shell,SC,Sym)
175 : : CASE (2)
176 : 46 : call tdep_init_shell2at(Invar,MPIdata,Shell,SC,Sym)
177 : : CASE (3)
178 : 16 : call tdep_init_shell3at(Invar,MPIdata,Shell,SC,Sym)
179 : : CASE (4)
180 : 9 : call tdep_init_shell4at(Invar,MPIdata,Shell,SC,Sym)
181 : : CASE DEFAULT
182 [ + + + + : 117 : ABI_ERROR('Wrong value for order')
- ]
183 : : END SELECT
184 : :
185 : 117 : end subroutine tdep_init_shell
186 : : !!***
187 : :
188 : : !====================================================================================================
189 : :
190 : : !!****f* ABINIT/m_tdep_shell/tdep_init_shell1at
191 : : !! NAME
192 : : !! tdep_init_shell1at
193 : : !!
194 : : !! FUNCTION
195 : : !! Build the shells of first-order interactions: identify the sets of
196 : : !! atoms that are equivalent by symmetry, store for each atom
197 : : !! and each shell the list of interactions and the symmetry operations
198 : : !! relating them to the reference interaction of the shell, and compute
199 : : !! (see tdep_calc_nbcoeff) the number of independent IFC coefficients of
200 : : !! each shell together with the corresponding projectors.
201 : : !! The details are written in the output file and in the
202 : : !! _nbcoeff-phi1.dat file.
203 : : !! Note that the number of independent coefficients is currently set to
204 : : !! zero for every shell at the first order (and the corresponding
205 : : !! interactions are removed).
206 : : !!
207 : : !! INPUTS
208 : : !! Invar = Input object containing the input variables.
209 : : !! MPIdata = Info on MPI parallelism.
210 : : !! SC = Supercell object, used for the interatomic distances (SC%dij).
211 : : !!
212 : : !! OUTPUT
213 : : !! Shell1at = Shell object for first-order interactions.
214 : : !!
215 : : !! SIDE EFFECTS
216 : : !! Sym = Symmetries object describing all the symmetry operations of the crystal.
217 : : !!
218 : : !! SOURCE
219 : :
220 : 46 : subroutine tdep_init_shell1at(Invar,MPIdata,Shell1at,SC,Sym)
221 : :
222 : : type(Shell_type),intent(inout) :: Shell1at
223 : : type(atdep_dataset_type),intent(in) :: Invar
224 : : type(MPI_enreg_type), intent(in) :: MPIdata
225 : : type(tdep_Supercell_type),intent(in) :: SC
226 : : type(Symmetries_type),intent(inout) :: Sym
227 : :
228 : : integer :: ishell,iatcell,iatom,eatom,iatref,isym
229 : : integer :: natom,natom_unitcell,counter,ncoeff,ncoeff_prev
230 : : integer :: norder,order,nshell_max,nshell
231 : 46 : integer, allocatable :: ref1at(:,:),Isym1at(:,:)
232 : :
233 : 46 : natom = Invar%natom
234 : 46 : natom_unitcell = Invar%natom_unitcell
235 : 46 : nshell_max = Invar%nshell_max
236 : 46 : order = 1
237 : 46 : norder = 3
238 : 46 : Shell1at%order = order
239 : 46 : Shell1at%norder = norder
240 : 46 : Shell1at%natom = natom
241 : :
242 : 46 : write(ab_out,*) ' '
243 : 46 : write(ab_out,*) '#############################################################################'
244 : 46 : write(ab_out,*) '####### FIRST ORDER : find the number of coefficients #######################'
245 : 46 : write(ab_out,*) '#############################################################################'
246 : :
247 : :
248 : : ! - Identify the shells
249 : : ! - Store the index of the atoms included in each shell
250 : : ! - Store the reference atoms for each shell
251 : : ! - Compute the symetry operation between the reference atom and another one
252 : 46 : write(ab_out,*) ' Build the ref1at and Isym1at tables...'
253 [ - + + - : 14786 : ABI_MALLOC(ref1at ,(natom,2)) ; ref1at (:,:)=zero
+ - + - -
+ - + + +
+ + ]
254 [ - + + - : 7462 : ABI_MALLOC(Isym1at,(natom,1)) ; Isym1at(:,:)=zero
+ - + - -
+ - + + +
+ + ]
255 : 46 : ishell=0
256 [ + + ]: 156 : do iatcell=1,natom_unitcell
257 [ + + ]: 110 : if (ref1at(iatcell,1).ne.0) cycle
258 : 65 : ishell=ishell+1
259 [ + + ]: 9805 : do eatom=1,natom
260 [ + + ]: 9804 : if (ref1at(eatom,1).eq.0) then
261 [ + + ]: 52766 : do isym=1,Sym%nsym
262 [ + + ]: 52766 : if (Sym%indsym(4,isym,eatom).eq.iatcell) then
263 : 7278 : Isym1at(eatom,1)=isym
264 : 7278 : ref1at(eatom,1)=iatcell
265 : 7278 : ref1at(eatom,2)=ishell
266 : : ! GA: This polutes the output. Write this in a file if needed.
267 : : !if (Invar%debug)
268 : : ! write(ab_out,'(a,1x,2(i4,1x),a,i4)') 'For:',iatcell,eatom,' direct transformation with isym=',Isym1at(eatom,1)
269 : : !end if
270 : 7278 : exit
271 : : end if
272 : : end do !isym
273 : : end if !already treated
274 : : end do !eatom
275 : : end do !iatcell
276 : 46 : Shell1at%nshell = ishell
277 : 46 : nshell = Shell1at%nshell
278 [ - + ]: 46 : if (nshell.gt.nshell_max) then
279 : 0 : write(ab_out,*) ' STOP : The maximum number of shells allowed by the code is:',nshell_max
280 : 0 : write(ab_out,*) ' In the present calculation, the number of shells is:',nshell
281 : 0 : write(ab_out,*) ' Action: increase nshell_max'
282 : 0 : ABI_ERROR('The maximum number of shells allowed by the code is reached')
283 : : end if
284 : :
285 : :
286 : : ! Store all the previous quantities in a better way than in ref1at (without using too memory).
287 : 46 : write(ab_out,*) ' Build the Shell1at datatype...'
288 [ - + + - : 268 : ABI_MALLOC(Shell1at%neighbours,(1,nshell))
+ - + - -
+ - + - +
+ + + + -
+ ]
289 [ - + + - : 203 : ABI_CALLOC(Shell1at%iatref, (nshell))
+ - + - -
+ - + - +
+ + ]
290 [ + + ]: 111 : do ishell=1,Shell1at%nshell
291 : : counter=0
292 [ + + ]: 9759 : do iatom=1,natom
293 [ + + ]: 9759 : if (ref1at(iatom,2).eq.ishell) counter=counter+1
294 : : end do
295 : 65 : Shell1at%neighbours(1,ishell)%n_interactions=counter
296 [ + - ]: 65 : if (counter.eq.0) then
297 : : cycle
298 : : end if
299 [ - + + - : 7473 : ABI_CALLOC(Shell1at%neighbours(1,ishell)%atomj_in_shell,(counter))
+ - - + -
+ - + +
+ ]
300 [ - + + - : 7408 : ABI_CALLOC(Shell1at%neighbours(1,ishell)%sym_in_shell,(counter))
- + - + -
+ + + ]
301 : : counter=0
302 [ + + ]: 9805 : do iatom=1,natom
303 [ + + ]: 9759 : if (ref1at(iatom,2).eq.ishell) then
304 : 7278 : counter=counter+1
305 : 7278 : Shell1at%neighbours(1,ishell)%atomj_in_shell(counter)=iatom
306 : 7278 : Shell1at%iatref(ishell)=ref1at(iatom,1)
307 : 7278 : Shell1at%neighbours(1,ishell)%sym_in_shell(counter)=Isym1at(iatom,1)
308 : : end if
309 : : end do
310 : : end do
311 [ - + ]: 46 : ABI_FREE(ref1at)
312 [ - + ]: 46 : ABI_FREE(Isym1at)
313 : :
314 : : ! Find the number of coefficients of the (3x3) Phi2 for a given shell
315 [ - + + - : 983 : ABI_CALLOC(Shell1at%proj, (norder,norder,nshell))
+ - + - -
+ - + - +
+ + + + +
+ ]
316 [ - + + - : 157 : ABI_CALLOC(Shell1at%ncoeff, (nshell))
- + - + -
+ + + ]
317 [ - + + - : 157 : ABI_CALLOC(Shell1at%ncoeff_prev, (nshell))
- + - + -
+ + + ]
318 : 46 : write(ab_out,*) ' Number of shells=',nshell
319 : 46 : write(ab_out,*) '============================================================================'
320 [ + - ]: 46 : if (MPIdata%iam_master) open(unit=16,file=trim(Invar%output_prefix)//'_nbcoeff-phi1.dat')
321 : 46 : ncoeff_prev=0
322 [ + + ]: 111 : do ishell=1,nshell
323 : 65 : ncoeff=0
324 : 65 : iatref=Shell1at%iatref(ishell)
325 : 65 : write(ab_out,*) 'Shell number:',ishell
326 : 65 : write(ab_out,'(a,i5,a)') ' For atom',iatref,':'
327 : 65 : call tdep_calc_nbcoeff(SC,MPIdata,iatref,ishell,1,1,1,ncoeff,norder,nshell,order,Shell1at%proj,Sym)
328 : 65 : ncoeff=0
329 : : if (ncoeff.eq.0) then
330 : 65 : Shell1at%neighbours(1,ishell)%n_interactions=0
331 [ + - ]: 65 : if(allocated(Shell1at%neighbours(1,ishell)%atomj_in_shell)) then
332 : 65 : ABI_FREE(Shell1at%neighbours(1,ishell)%atomj_in_shell)
333 : : end if
334 [ + - ]: 65 : if(allocated(Shell1at%neighbours(1,ishell)%sym_in_shell)) then
335 : 65 : ABI_FREE(Shell1at%neighbours(1,ishell)%sym_in_shell)
336 : : end if
337 : : end if
338 : 65 : Shell1at%ncoeff (ishell)=ncoeff
339 : 65 : Shell1at%ncoeff_prev(ishell)=ncoeff_prev
340 : : ncoeff_prev=ncoeff_prev+ncoeff
341 : 65 : write(ab_out,*)' Number of independant coefficients in this shell=',ncoeff
342 : 65 : write(ab_out,*)' Number of interactions in this shell=',Shell1at%neighbours(1,ishell)%n_interactions
343 : 111 : write(ab_out,*) '============================================================================'
344 : : end do
345 : 46 : write(ab_out,*)' >>>>>> Total number of coefficients at the first order=',ncoeff_prev
346 [ + - ]: 46 : if (MPIdata%iam_master) close(16)
347 : 46 : Shell1at%ntotcoeff=ncoeff_prev
348 : :
349 : 46 : end subroutine tdep_init_shell1at
350 : : !!***
351 : :
352 : : !====================================================================================================
353 : :
354 : : !!****f* ABINIT/m_tdep_shell/tdep_init_shell2at
355 : : !! NAME
356 : : !! tdep_init_shell2at
357 : : !!
358 : : !! FUNCTION
359 : : !! Build the shells of second-order interactions: identify the sets of
360 : : !! pairs of atoms that are equivalent by symmetry (within
361 : : !! the cutoff radius Invar%rcut), store for each atom
362 : : !! and each shell the list of interactions and the symmetry operations
363 : : !! relating them to the reference interaction of the shell, and compute
364 : : !! (see tdep_calc_nbcoeff) the number of independent IFC coefficients of
365 : : !! each shell together with the corresponding projectors.
366 : : !! The details are written in the output file and in the
367 : : !! _nbcoeff-phi2.dat file.
368 : : !!
369 : : !! INPUTS
370 : : !! Invar = Input object containing the input variables.
371 : : !! MPIdata = Info on MPI parallelism.
372 : : !! SC = Supercell object, used for the interatomic distances (SC%dij).
373 : : !!
374 : : !! OUTPUT
375 : : !! Shell2at = Shell object for second-order interactions.
376 : : !!
377 : : !! SIDE EFFECTS
378 : : !! Sym = Symmetries object describing all the symmetry operations of the crystal.
379 : : !!
380 : : !! SOURCE
381 : :
382 : 46 : subroutine tdep_init_shell2at(Invar,MPIdata,Shell2at,SC,Sym)
383 : :
384 : : type(Shell_type),intent(inout) :: Shell2at
385 : : type(atdep_dataset_type),intent(in) :: Invar
386 : : type(MPI_enreg_type), intent(in) :: MPIdata
387 : : type(tdep_Supercell_type),intent(in) :: SC
388 : : type(Symmetries_type),intent(inout) :: Sym
389 : :
390 : : integer :: ishell,iatcell,iatom,jatom,eatom,fatom,iatref,jatref
391 : : integer :: natom,natom_unitcell,counter,ncoeff,ncoeff_prev
392 : : integer :: norder,order,nshell_max,nshell,jj
393 : 46 : integer, allocatable :: ref2at(:,:,:),Isym2at(:,:,:),iperm(:)
394 : 46 : real(dp), allocatable :: dist(:)
395 : :
396 : 46 : natom = SC%natom
397 : 46 : natom_unitcell = SC%natom_unitcell
398 : 46 : nshell_max = Invar%nshell_max
399 : 46 : order = 2
400 : 46 : norder = 9
401 : 46 : Shell2at%order = order
402 : 46 : Shell2at%norder = norder
403 : 46 : Shell2at%natom = natom
404 : :
405 : 46 : write(ab_out,*) ' '
406 : 46 : write(ab_out,*) '#############################################################################'
407 : 46 : write(ab_out,*) '###### SECOND ORDER : find the number of coefficients #######################'
408 : 46 : write(ab_out,*) '#############################################################################'
409 : :
410 : : ! - Identify the shells
411 : : ! - Store the index of the atoms included in each shell
412 : : ! - Store the reference atoms for each shell
413 : : ! - Compute the symetry operation between the reference atom and another one
414 : 46 : write(ab_out,*) ' Build the ref2at and Isym2at tables...'
415 [ - + + - : 4676294 : ABI_MALLOC(ref2at ,(natom,natom,3)) ; ref2at (:,:,:)=zero
- + + - +
- + - + -
- + - + +
+ + + +
+ ]
416 [ - + + - : 3117606 : ABI_MALLOC(Isym2at,(natom,natom,2)) ; Isym2at(:,:,:)=zero
- + + - +
- + - + -
- + - + +
+ + + +
+ ]
417 [ - + + - : 248 : ABI_MALLOC(Shell2at%ishell_self,(natom_unitcell)) ; Shell2at%ishell_self(:)=zero
+ - + - -
+ - + - +
+ + ]
418 [ - + + - : 138 : ABI_MALLOC(dist,(natom))
+ - + - -
+ - + ]
419 [ - + + - : 138 : ABI_MALLOC(iperm,(natom))
+ - + - -
+ - + ]
420 : 46 : ishell=0
421 [ + + ]: 156 : do iatcell=1,natom_unitcell
422 : : ! Loop over the atoms by increasing distance, so that the shells are numbered by increasing radius.
423 [ + + ]: 15788 : dist(:)=SC%dij(iatcell,:,1)
424 [ + - + + : 31466 : iperm(:)=(/(jj,jj=1,natom)/)
+ + ]
425 : 110 : call sort_dp(natom,dist,iperm,tol6)
426 [ + + ]: 15834 : do jj=1,natom
427 : 15678 : jatom=iperm(jj)
428 : : ! Interactions are only computed until Rcut in order to have complete shell of neighbours.
429 : : ! Otherwise the symetries are broken.
430 [ + + + + ]: 15678 : if ((ref2at(iatcell,jatom,1).ne.0).or.(SC%dij(iatcell,jatom,1).gt.(Invar%rcut*0.99))) cycle
431 : 522 : ishell=ishell+1
432 [ + + ]: 522 : if (iatcell.eq.jatom) Shell2at%ishell_self(iatcell)=ishell
433 [ + + ]: 83562 : do eatom=1,natom
434 [ + + ]: 18860316 : do fatom=1,natom
435 [ + + + + ]: 18761708 : if ((ref2at(eatom,fatom,1).eq.0).and.&
436 : 82930 : & (abs(SC%dij(iatcell,jatom,1)-SC%dij(eatom,fatom,1)).lt.tol6)) then
437 : 330346 : call tdep_SearchS_2at(SC,iatcell,jatom,eatom,fatom,Isym2at,Sym)
438 [ + + ]: 330346 : if (Isym2at(eatom,fatom,2)==1) then
439 : : ! GA: This polutes the output. Write this in a file if needed.
440 : : !if (Invar%debug) then
441 : : ! write(ab_out,'(a,1x,4(i4,1x),a,i4)') 'For:',iatcell,jatom,eatom,fatom, &
442 : : ! & ' direct transformation with isym=',Isym2at(eatom,fatom,1)
443 : : !end if
444 : 273876 : ref2at(eatom,fatom,1)=iatcell
445 : 273876 : ref2at(eatom,fatom,2)=jatom
446 : 273876 : ref2at(eatom,fatom,3)=ishell
447 : : ! The Phi2 has to be symetric (transposition symetries)
448 : : !if (Invar%debug) then
449 : : ! write(ab_out,'(a,1x,4(i4,1x),a,i4)') 'For:',iatcell,jatom,eatom,fatom, &
450 : : ! & ' transformation+permutation with isym=',Isym2at(eatom,fatom,1)
451 : : !end if
452 : 273876 : ref2at(fatom,eatom,1)=iatcell
453 : 273876 : ref2at(fatom,eatom,2)=jatom
454 : 273876 : ref2at(fatom,eatom,3)=ishell
455 : 273876 : Isym2at(fatom,eatom,1)=Isym2at(eatom,fatom,1)
456 : 273876 : Isym2at(fatom,eatom,2)=2
457 : : else
458 : : !if (Invar%debug) then
459 : : ! write(ab_out,'(a,4(1x,i4))') &
460 : : ! & 'NO SYMETRY OPERATION BETWEEN (iatom,jatom) and (eatom,fatom)=',iatcell,jatom,eatom,fatom
461 : : !end if
462 : : end if
463 : : end if !already treated
464 : : end do !fatom
465 : : end do !eatom
466 : : end do !jj
467 : : end do !iatcell
468 [ - + ]: 46 : ABI_FREE(dist)
469 [ - + ]: 46 : ABI_FREE(iperm)
470 : 46 : Shell2at%nshell = ishell
471 : 46 : nshell = Shell2at%nshell
472 [ - + ]: 46 : if (nshell.gt.nshell_max) then
473 : 0 : write(ab_out,*) ' STOP : The maximum number of shells allowed by the code is:',nshell_max
474 : 0 : write(ab_out,*) ' In the present calculation, the number of shells is:',nshell
475 : 0 : write(ab_out,*) ' Action: increase nshell_max'
476 : 0 : ABI_ERROR('The maximum number of shells allowed by the code is reached')
477 : : end if
478 : :
479 : :
480 : : ! Store all the previous quantities in a better way than in ref2at (without using too memory).
481 : 46 : write(ab_out,*) ' Build the Shell2at datatype...'
482 [ - + + - : 83636 : ABI_MALLOC(Shell2at%neighbours,(natom,nshell))
- + + - +
- + - - +
- + - + +
+ + + +
- ]
483 [ - + + - : 660 : ABI_CALLOC(Shell2at%iatref, (nshell))
+ - + - -
+ - + - +
+ + ]
484 [ - + + - : 614 : ABI_CALLOC(Shell2at%jatref, (nshell))
- + - + -
+ + + ]
485 [ + + ]: 568 : do ishell=1,nshell
486 [ + + ]: 83498 : do iatom=1,natom
487 : 82930 : counter=0
488 [ + + ]: 18844638 : do jatom=1,natom
489 [ + + ]: 18844638 : if (ref2at(iatom,jatom,3).eq.ishell) counter=counter+1
490 : : end do
491 : 82930 : Shell2at%neighbours(iatom,ishell)%n_interactions=counter
492 [ + + ]: 82930 : if (counter.eq.0) cycle
493 [ - + + - : 214230 : ABI_MALLOC(Shell2at%neighbours(iatom,ishell)%atomj_in_shell,(counter))
+ - - + -
+ - + ]
494 [ - + + - : 142820 : ABI_MALLOC(Shell2at%neighbours(iatom,ishell)%sym_in_shell,(counter))
- + - + -
+ ]
495 [ - + + - : 142820 : ABI_MALLOC(Shell2at%neighbours(iatom,ishell)%transpose_in_shell,(counter))
- + - + -
+ ]
496 [ + + ]: 611884 : Shell2at%neighbours(iatom,ishell)%atomj_in_shell(:)=zero
497 [ + + ]: 611884 : Shell2at%neighbours(iatom,ishell)%sym_in_shell(:)=zero
498 [ + + ]: 611884 : Shell2at%neighbours(iatom,ishell)%transpose_in_shell(:)=zero
499 : 71410 : counter=0
500 [ + + ]: 16870824 : do jatom=1,natom
501 [ + + ]: 16881822 : if (ref2at(iatom,jatom,3).eq.ishell) then
502 : 540474 : counter=counter+1
503 : 540474 : Shell2at%neighbours(iatom,ishell)%atomj_in_shell(counter)=jatom
504 : 540474 : Shell2at%iatref (ishell)=ref2at(iatom,jatom,1)
505 : 540474 : Shell2at%jatref (ishell)=ref2at(iatom,jatom,2)
506 : 540474 : Shell2at%neighbours(iatom,ishell)%sym_in_shell (counter)=Isym2at(iatom,jatom,1)
507 : 540474 : Shell2at%neighbours(iatom,ishell)%transpose_in_shell(counter)=Isym2at(iatom,jatom,2)
508 : : end if
509 : : end do
510 : : end do
511 : : end do
512 [ - + ]: 46 : ABI_FREE(ref2at)
513 [ - + ]: 46 : ABI_FREE(Isym2at)
514 : :
515 : : ! Find the number of coefficients of the (3x3) Phi2 for a given shell
516 [ - + + - : 47640 : ABI_CALLOC(Shell2at%proj, (norder,norder,nshell))
+ - + - -
+ - + - +
+ + + + +
+ ]
517 [ - + + - : 614 : ABI_CALLOC(Shell2at%ncoeff, (nshell))
- + - + -
+ + + ]
518 [ - + + - : 614 : ABI_CALLOC(Shell2at%ncoeff_prev, (nshell))
- + - + -
+ + + ]
519 : 46 : write(ab_out,*) ' Number of shells=',nshell
520 : 46 : write(ab_out,*) '============================================================================'
521 [ + - ]: 46 : if (MPIdata%iam_master) open(unit=16,file=trim(Invar%output_prefix)//'_nbcoeff-phi2.dat')
522 : 46 : ncoeff_prev=0
523 [ + + ]: 568 : do ishell=1,nshell
524 : 522 : ncoeff=0
525 : 522 : iatref=Shell2at%iatref(ishell)
526 : 522 : jatref=Shell2at%jatref(ishell)
527 : 522 : write(ab_out,*) 'Shell number:',ishell
528 : 522 : write(ab_out,'(a,i5,a,i5,a,f16.10)') ' Between atom',iatref,' and ',jatref,' the distance is=',SC%dij(iatref,jatref,1)
529 : 522 : call tdep_calc_nbcoeff(SC,MPIdata,iatref,ishell,jatref,1,1,ncoeff,norder,nshell,order,Shell2at%proj,Sym)
530 : 522 : Shell2at%ncoeff (ishell)=ncoeff
531 : 522 : Shell2at%ncoeff_prev(ishell)=ncoeff_prev
532 : 522 : ncoeff_prev=ncoeff_prev+ncoeff
533 : 522 : write(ab_out,*)' Number of independant coefficients in this shell=',ncoeff
534 : 522 : write(ab_out,*)' Number of interactions in this shell=',Shell2at%neighbours(iatref,ishell)%n_interactions
535 : 568 : write(ab_out,*) '============================================================================'
536 : : end do
537 : 46 : write(ab_out,*)' >>>>>> Total number of coefficients at the second order=',ncoeff_prev
538 [ + - ]: 46 : if (MPIdata%iam_master) close(16)
539 : 46 : Shell2at%ntotcoeff=ncoeff_prev
540 : :
541 : : !! Write debugging info
542 : : ! open(unit=91,file='Shell2at.dat')
543 : : ! write(91,*) Shell2at%nshell
544 : : ! do ishell=1,Shell2at%nshell
545 : : ! write(91,*) Shell2at%ncoeff(ishell)
546 : : ! write(91,*) Shell2at%ncoeff_prev(ishell)
547 : : ! write(91,*) Shell2at%iatref(ishell)
548 : : ! write(91,*) Shell2at%jatref(ishell)
549 : : ! do iatom=1,Invar%natom
550 : : ! write(91,*) Shell2at%neighbours(iatom,ishell)%n_interactions
551 : : ! do ii=1,Shell2at%neighbours(iatom,ishell)%n_interactions
552 : : ! write(91,*) Shell2at%neighbours(iatom,ishell)%sym_in_shell(ii)
553 : : ! write(91,*) Shell2at%neighbours(iatom,ishell)%transpose_in_shell(ii)
554 : : ! write(91,*) Shell2at%neighbours(iatom,ishell)%atomj_in_shell(ii)
555 : : ! end do
556 : : ! end do
557 : : ! end do
558 : : ! close(91)
559 : :
560 : 46 : end subroutine tdep_init_shell2at
561 : : !!***
562 : :
563 : : !====================================================================================================
564 : :
565 : : !!****f* ABINIT/m_tdep_shell/tdep_init_shell3at
566 : : !! NAME
567 : : !! tdep_init_shell3at
568 : : !!
569 : : !! FUNCTION
570 : : !! Build the shells of third-order interactions: identify the sets of
571 : : !! triplets of atoms that are equivalent by symmetry (within
572 : : !! the cutoff radius Invar%rcut3), store for each atom
573 : : !! and each shell the list of interactions and the symmetry operations
574 : : !! relating them to the reference interaction of the shell, and compute
575 : : !! (see tdep_calc_nbcoeff) the number of independent IFC coefficients of
576 : : !! each shell together with the corresponding projectors.
577 : : !! The details are written in the output file and in the
578 : : !! _nbcoeff-phi3.dat file.
579 : : !!
580 : : !! INPUTS
581 : : !! Invar = Input object containing the input variables.
582 : : !! MPIdata = Info on MPI parallelism.
583 : : !! SC = Supercell object, used for the interatomic distances (SC%dij).
584 : : !!
585 : : !! OUTPUT
586 : : !! Shell3at = Shell object for third-order interactions.
587 : : !!
588 : : !! SIDE EFFECTS
589 : : !! Sym = Symmetries object describing all the symmetry operations of the crystal.
590 : : !!
591 : : !! SOURCE
592 : :
593 : 16 : subroutine tdep_init_shell3at(Invar,MPIdata,Shell3at,SC,Sym)
594 : :
595 : : type(Shell_type),intent(inout) :: Shell3at
596 : : type(atdep_dataset_type),intent(in) :: Invar
597 : : type(MPI_enreg_type), intent(in) :: MPIdata
598 : : type(tdep_Supercell_type),intent(in) :: SC
599 : : type(Symmetries_type),intent(inout) :: Sym
600 : :
601 : : integer :: ii,ishell,iatom,jatom,katom,eatom,fatom,gatom,iatref,jatref,katref
602 : : integer :: natom,natom_unitcell,watom,xatom,yatom,ninteractions,ncoeff,ncoeff_prev,nshell_tmp
603 : : integer :: find_equivalent,ninter,iat_ref,jat_ref,kat_ref,tmpinter
604 : : integer :: norder,order,nshell_max,nshell
605 : : double precision :: norma,normb,normc
606 : : integer :: Isym3at(2)
607 : 16 : integer, allocatable :: atref(:,:),interactions(:,:)
608 : :
609 : 16 : natom = Invar%natom
610 : 16 : natom_unitcell = Invar%natom_unitcell
611 : 16 : nshell_max = Invar%nshell_max
612 : 16 : order = 3
613 : 16 : norder = 27
614 : 16 : Shell3at%order = order
615 : 16 : Shell3at%norder = norder
616 : 16 : Shell3at%natom = natom
617 : :
618 : 16 : write(ab_out,*) ' '
619 : 16 : write(ab_out,*) '#############################################################################'
620 : 16 : write(ab_out,*) '###### THIRD ORDER : find the number of coefficients ########################'
621 : 16 : write(ab_out,*) '#############################################################################'
622 : :
623 : : ! 1/ Identify the shells
624 [ - + + - : 24064 : ABI_CALLOC(interactions,(natom_unitcell,nshell_max))
- + + - +
- + - - +
- + + + +
+ ]
625 [ - + + - : 24096 : ABI_CALLOC(atref,(nshell_max,3))
+ - + - -
+ - + + +
+ + ]
626 [ - + + - : 80 : ABI_CALLOC(Shell3at%ishell_self,(natom_unitcell))
+ - + - -
+ - + - +
+ + ]
627 : 16 : nshell_tmp=0
628 [ + + ]: 48 : do iatom=1,natom_unitcell
629 [ + + ]: 5728 : do jatom=1,natom
630 [ + + ]: 1064336 : do katom=1,natom
631 : : !FB write(std_out,*) 'NEW COORD1 : iatom,jatom,katom=',iatom,jatom,katom
632 : : ! WARNING: distance(j,k).ne.|djk| due to the inbox procedure when computing distance(j,k).
633 : : ! So, compute |djk| using vec(ij) and vec(ik).
634 : : norma=dsqrt((SC%dij(iatom,katom,2)-SC%dij(iatom,jatom,2))**2+&
635 : : & (SC%dij(iatom,katom,3)-SC%dij(iatom,jatom,3))**2+&
636 : 1058624 : & (SC%dij(iatom,katom,4)-SC%dij(iatom,jatom,4))**2)
637 : : ! Interactions are only computed until Rcut3 in order to have complete shell of neighbours.
638 : : ! Otherwise the symetries are broken.
639 : : if ((SC%dij(iatom,jatom,1).gt.(Invar%rcut3*0.99)).or.&
640 [ + + + + : 1058624 : & (norma .gt.(Invar%rcut3*0.99)).or.&
+ + ]
641 : : & (SC%dij(iatom,katom,1).gt.(Invar%rcut3*0.99))) cycle
642 [ + + ]: 7020 : if (nshell_tmp.eq.0) then
643 [ + + ]: 64 : atref(1,:)=1
644 : 16 : nshell_tmp=nshell_tmp+1
645 [ + - ]: 16 : if ((iatom.eq.jatom).and.(jatom.eq.katom)) Shell3at%ishell_self(iatom)=nshell_tmp
646 : 16 : interactions(iatom,nshell_tmp)=interactions(iatom,nshell_tmp)+1
647 : 16 : cycle
648 : : else
649 : 1324 : find_equivalent=0
650 [ + + ]: 5488 : do ishell=1,nshell_tmp
651 : 5447 : iat_ref=atref(ishell,1) ; jat_ref=atref(ishell,2) ; kat_ref=atref(ishell,3)
652 : : normb=dsqrt((SC%dij(iat_ref,kat_ref,2)-SC%dij(iat_ref,jat_ref,2))**2+&
653 : : & (SC%dij(iat_ref,kat_ref,3)-SC%dij(iat_ref,jat_ref,3))**2+&
654 : 5447 : & (SC%dij(iat_ref,kat_ref,4)-SC%dij(iat_ref,jat_ref,4))**2)
655 [ + + ]: 32479 : do ii=1,6
656 [ + + ]: 28315 : if (ii.eq.1) then ; eatom=iatom ; fatom=jatom ; gatom=katom ; end if
657 [ + + ]: 28315 : if (ii.eq.2) then ; eatom=iatom ; fatom=katom ; gatom=jatom ; end if
658 [ + + ]: 28315 : if (ii.eq.3) then ; eatom=jatom ; fatom=iatom ; gatom=katom ; end if
659 [ + + ]: 28315 : if (ii.eq.4) then ; eatom=jatom ; fatom=katom ; gatom=iatom ; end if
660 [ + + ]: 28315 : if (ii.eq.5) then ; eatom=katom ; fatom=iatom ; gatom=jatom ; end if
661 [ + + ]: 28315 : if (ii.eq.6) then ; eatom=katom ; fatom=jatom ; gatom=iatom ; end if
662 : : normc=dsqrt((SC%dij(eatom,gatom,2)-SC%dij(eatom,fatom,2))**2+&
663 : : & (SC%dij(eatom,gatom,3)-SC%dij(eatom,fatom,3))**2+&
664 : 28315 : & (SC%dij(eatom,gatom,4)-SC%dij(eatom,fatom,4))**2)
665 : : !FB if ((abs(SC%dij(iatom,jatom,1)-SC%dij(eatom,fatom,1)).lt.1.d-3).and.&
666 : : !FB& (abs(norma -normb ).lt.1.d-3).and.&
667 : : !FB& (abs(SC%dij(iatom,katom,1)-SC%dij(eatom,gatom,1)).lt.1.d-3)) then
668 : : if ((abs(SC%dij(iat_ref,jat_ref,1)-SC%dij(eatom,fatom,1)).lt.1.d-6).and.&
669 [ + + + + : 28315 : & (abs(normb -normc ).lt.1.d-6).and.&
+ - ]
670 : 4164 : & (abs(SC%dij(iat_ref,kat_ref,1)-SC%dij(eatom,gatom,1)).lt.1.d-6)) then
671 : 2075 : Isym3at(:)=0
672 : 2075 : call tdep_SearchS_3at(SC,iat_ref,jat_ref,kat_ref,eatom,fatom,gatom,Isym3at,Sym)
673 [ + + ]: 2075 : if (Isym3at(2).eq.1) find_equivalent=1
674 : : if (find_equivalent.eq.1) then
675 : 1283 : interactions(iatom,ishell)=interactions(iatom,ishell)+1
676 : : !FB write(std_out,*) 'The number of interactions in this shell is=',ishell,interactions(iatom,ishell)
677 : 1283 : exit
678 : : end if
679 : : end if
680 : : end do !ii
681 [ + + ]: 5488 : if (find_equivalent.eq.1) exit
682 : : end do !ishell
683 [ + + ]: 1324 : if (find_equivalent.eq.0) then
684 : 41 : nshell_tmp=nshell_tmp+1
685 [ - + ]: 41 : if (nshell_tmp.gt.nshell_max) then
686 : 0 : ABI_ERROR('The shell number index is greater than the shell number max defined in the code')
687 : : end if
688 [ + + ]: 41 : if ((iatom.eq.jatom).and.(jatom.eq.katom)) Shell3at%ishell_self(iatom)=nshell_tmp
689 : 41 : interactions(iatom,nshell_tmp)=interactions(iatom,nshell_tmp)+1
690 : 41 : atref(nshell_tmp,1)=iatom
691 : 41 : atref(nshell_tmp,2)=jatom
692 : 41 : atref(nshell_tmp,3)=katom
693 : : !FB write(std_out,'(a,1x,4(i5,1x))') 'NEW SHELL1 : nshell_tmp,iatom,jatom,katom=',nshell_tmp,iatom,jatom,katom
694 : : end if
695 : : end if
696 : : end do !katom
697 : : end do !jatom
698 : : end do !iatom
699 [ - + ]: 16 : ABI_FREE(atref)
700 : :
701 : : ! 2/ Allocate the datatype Shell3at%...
702 : 16 : Shell3at%nshell=nshell_tmp
703 : 16 : nshell = Shell3at%nshell
704 [ - + + - : 9961 : ABI_MALLOC(Shell3at%neighbours,(natom,nshell))
- + + - +
- + - - +
- + - + +
+ + + +
- ]
705 [ - + + - : 105 : ABI_MALLOC(Shell3at%iatref,(nshell)); Shell3at%iatref(:)=zero
+ - + - -
+ - + - +
+ + ]
706 [ - + + - : 89 : ABI_MALLOC(Shell3at%jatref,(nshell)); Shell3at%jatref(:)=zero
- + - + -
+ + + ]
707 [ - + + - : 89 : ABI_MALLOC(Shell3at%katref,(nshell)); Shell3at%katref(:)=zero
- + - + -
+ + + ]
708 [ + + ]: 73 : do ishell=1,nshell
709 [ + + ]: 9913 : do iatom=1,natom
710 : 9840 : ninteractions=interactions(mod(iatom-1,natom_unitcell)+1,ishell)
711 : 9840 : Shell3at%neighbours(iatom,ishell)%n_interactions=ninteractions
712 [ + + ]: 9840 : if (ninteractions.eq.0) cycle
713 [ - + + - : 25680 : ABI_MALLOC(Shell3at%neighbours(iatom,ishell)%atomj_in_shell,(ninteractions))
+ - + - -
+ - + -
+ ]
714 [ - + + - : 17120 : ABI_MALLOC(Shell3at%neighbours(iatom,ishell)%atomk_in_shell,(ninteractions))
- + - + -
+ ]
715 [ - + + - : 17120 : ABI_MALLOC(Shell3at%neighbours(iatom,ishell)%sym_in_shell,(ninteractions))
- + - + -
+ ]
716 [ - + + - : 17120 : ABI_MALLOC(Shell3at%neighbours(iatom,ishell)%transpose_in_shell,(ninteractions))
- + - + -
+ ]
717 [ + + ]: 106056 : Shell3at%neighbours(iatom,ishell)%atomj_in_shell(:)=zero
718 [ + + ]: 106056 : Shell3at%neighbours(iatom,ishell)%atomk_in_shell(:)=zero
719 [ + + ]: 106056 : Shell3at%neighbours(iatom,ishell)%sym_in_shell(:)=zero
720 [ + + ]: 107393 : Shell3at%neighbours(iatom,ishell)%transpose_in_shell(:)=zero
721 : : end do
722 : : end do
723 [ - + ]: 16 : ABI_FREE(interactions)
724 : :
725 : : ! 3/ Store the index of the (couple of) atoms included in each shell
726 : : ! 4/ Store the reference (couple of) atoms for each shell
727 : : ! 5/ Compute the symetry operation between the reference (couple of) atoms and another one
728 [ - + + - : 9961 : ABI_MALLOC(interactions,(natom,nshell)) ; interactions(:,:)=0
- + + - +
- + - - +
- + + + +
+ ]
729 : 16 : nshell_tmp=0
730 [ + + ]: 2952 : do iatom=1,natom
731 [ + + ]: 575272 : do jatom=1,natom
732 [ + + ]: 118320536 : do katom=1,natom
733 : : !FB write(std_out,*) 'NEW COORD2 : iatom,jatom,katom=',iatom,jatom,katom
734 [ + + ]: 117745280 : if (nshell_tmp.eq.0) then
735 : 16 : nshell_tmp=nshell_tmp+1
736 : 16 : interactions(iatom,nshell_tmp)=1
737 : 16 : Shell3at%iatref(nshell_tmp)=iatom
738 : 16 : Shell3at%jatref(nshell_tmp)=jatom
739 : 16 : Shell3at%katref(nshell_tmp)=katom
740 : 16 : Shell3at%neighbours(iatom,nshell_tmp)%atomj_in_shell(interactions(iatom,nshell_tmp))=jatom
741 : 16 : Shell3at%neighbours(iatom,nshell_tmp)%atomk_in_shell(interactions(iatom,nshell_tmp))=katom
742 : 16 : Shell3at%neighbours(iatom,nshell_tmp)%sym_in_shell(interactions(iatom,nshell_tmp))=1
743 : 16 : Shell3at%neighbours(iatom,nshell_tmp)%transpose_in_shell(interactions(iatom,nshell_tmp))=1
744 : 16 : cycle
745 : : end if
746 : : ! WARNING: SC%dij(j,k).ne.|djk| due to the inbox procedure when computing SC%dij(j,k).
747 : : ! So, compute |djk| using vec(ij) and vec(ik).
748 : : norma=dsqrt((SC%dij(iatom,katom,2)-SC%dij(iatom,jatom,2))**2+&
749 : : & (SC%dij(iatom,katom,3)-SC%dij(iatom,jatom,3))**2+&
750 : 117745264 : & (SC%dij(iatom,katom,4)-SC%dij(iatom,jatom,4))**2)
751 : : ! Interactions are only computed until Rcut3<acell/2 in order to have complete shell of neighbours.
752 : : ! Otherwise the symetries are broken.
753 : : if ((SC%dij(iatom,jatom,1).gt.(Invar%rcut3*0.99)).or.&
754 [ + + + + : 117745264 : & (norma .gt.(Invar%rcut3*0.99)).or.&
+ + ]
755 : : & (SC%dij(iatom,katom,1).gt.(Invar%rcut3*0.99))) cycle
756 : : ! Search if the triplet has already been classified
757 : 97480 : find_equivalent=0
758 [ + + ]: 711119 : do ishell=1,nshell_tmp
759 [ + + ]: 6538231 : do ninter=1,interactions(iatom,ishell)
760 [ + + + + ]: 5827112 : if ((Shell3at%neighbours(iatom,ishell)%atomj_in_shell(ninter).eq.jatom).and.&
761 : 6440751 : & (Shell3at%neighbours(iatom,ishell)%atomk_in_shell(ninter).eq.katom)) find_equivalent=1
762 : : end do
763 : : end do
764 [ + + ]: 97480 : if (find_equivalent.eq.1) cycle
765 : : ! Search if the triplet belongs to a shell already found
766 [ + + ]: 94288 : do ishell=1,nshell_tmp
767 : 94247 : iat_ref=Shell3at%iatref(ishell) ; jat_ref=Shell3at%jatref(ishell) ; kat_ref=Shell3at%katref(ishell)
768 : : normb=dsqrt((SC%dij(iat_ref,kat_ref,2)-SC%dij(iat_ref,jat_ref,2))**2+&
769 : : & (SC%dij(iat_ref,kat_ref,3)-SC%dij(iat_ref,jat_ref,3))**2+&
770 : 94247 : & (SC%dij(iat_ref,kat_ref,4)-SC%dij(iat_ref,jat_ref,4))**2)
771 [ + + ]: 540071 : do ii=1,6
772 [ + + ]: 471967 : if (ii.eq.1) then ; eatom=iatom ; fatom=jatom ; gatom=katom ; end if
773 [ + + ]: 471967 : if (ii.eq.2) then ; eatom=iatom ; fatom=katom ; gatom=jatom ; end if
774 [ + + ]: 471967 : if (ii.eq.3) then ; eatom=jatom ; fatom=iatom ; gatom=katom ; end if
775 [ + + ]: 471967 : if (ii.eq.4) then ; eatom=jatom ; fatom=katom ; gatom=iatom ; end if
776 [ + + ]: 471967 : if (ii.eq.5) then ; eatom=katom ; fatom=iatom ; gatom=jatom ; end if
777 [ + + ]: 471967 : if (ii.eq.6) then ; eatom=katom ; fatom=jatom ; gatom=iatom ; end if
778 : : normc=dsqrt((SC%dij(eatom,gatom,2)-SC%dij(eatom,fatom,2))**2+&
779 : : & (SC%dij(eatom,gatom,3)-SC%dij(eatom,fatom,3))**2+&
780 : 471967 : & (SC%dij(eatom,gatom,4)-SC%dij(eatom,fatom,4))**2)
781 : : !FB if ((abs(SC%dij(iatom,jatom,1)-SC%dij(eatom,fatom,1)).lt.1.d-3).and.&
782 : : !FB& (abs(norma -normb ).lt.1.d-3).and.&
783 : : !FB& (abs(SC%dij(iatom,katom,1)-SC%dij(eatom,gatom,1)).lt.1.d-3)) then
784 : : if ((abs(SC%dij(iat_ref,jat_ref,1)-SC%dij(eatom,fatom,1)).lt.1.d-6).and.&
785 [ + + + + : 471967 : & (abs(normb -normc ).lt.1.d-6).and.&
+ - ]
786 : 68104 : & (abs(SC%dij(iat_ref,kat_ref,1)-SC%dij(eatom,gatom,1)).lt.1.d-6)) then
787 : 41103 : Isym3at(:)=0
788 : 41103 : call tdep_SearchS_3at(SC,iat_ref,jat_ref,kat_ref,eatom,fatom,gatom,Isym3at,Sym)
789 [ + + ]: 41103 : if (Isym3at(2).eq.1) then
790 : : find_equivalent=1
791 : : exit
792 : : end if
793 : : end if
794 : : end do !ii
795 [ + + ]: 94288 : if (find_equivalent.eq.1) exit
796 : : end do !ishell
797 : : ! The triplet belongs to a new shell
798 [ + + ]: 26184 : if (find_equivalent.eq.0) then
799 : 41 : nshell_tmp=nshell_tmp+1
800 : : ! Check that the new shell is allowed
801 [ - + ]: 41 : if (nshell_tmp.gt.nshell) then
802 : 0 : ABI_ERROR('The shell number index is greater than the shell number max computed previously')
803 : : end if
804 : 41 : Shell3at%iatref(nshell_tmp)=iatom
805 : 41 : Shell3at%jatref(nshell_tmp)=jatom
806 : 41 : Shell3at%katref(nshell_tmp)=katom
807 : 41 : eatom=iatom ; fatom=jatom ; gatom=katom
808 [ + + ]: 123 : Isym3at(:)=1
809 : 41 : ishell=nshell_tmp
810 : : !FB write(std_out,'(a,1x,4(i5,1x))') 'NEW SHELL2 : nshell_tmp,iatom,jatom,katom=',nshell_tmp,iatom,jatom,katom
811 : : end if
812 : : ! Classify the informations of the triplet in Shell3at
813 [ + + ]: 755608 : do ii=1,6
814 : : ! The Phi3 has to be symetric (transposition symetries)
815 [ + + ]: 157104 : if (ii==1) then ; watom=eatom ; xatom=fatom ; yatom=gatom ; endif !\Phi3_ijk
816 [ + + ]: 157104 : if (ii==2) then ; watom=eatom ; xatom=gatom ; yatom=fatom ; endif !\Phi3_ikj
817 [ + + ]: 157104 : if (ii==3) then ; watom=fatom ; xatom=eatom ; yatom=gatom ; endif !\Phi3_jik
818 [ + + ]: 157104 : if (ii==4) then ; watom=fatom ; xatom=gatom ; yatom=eatom ; endif !\Phi3_jki
819 [ + + ]: 157104 : if (ii==5) then ; watom=gatom ; xatom=eatom ; yatom=fatom ; endif !\Phi3_kij
820 [ + + ]: 157104 : if (ii==6) then ; watom=gatom ; xatom=fatom ; yatom=eatom ; endif !\Phi3_kji
821 : : ! Do not overwrite the Phi3_iik, Phi3_iji, Phi3_ijj or Phi3_iii IFCs
822 : : ! and avoid double counting of triplet interactions
823 [ + + + + ]: 157104 : if ((eatom.eq.fatom).and.((ii.eq.3).or.(ii.eq.4).or.(ii.eq.6))) cycle
824 [ + + + + ]: 111000 : if ((eatom.eq.gatom).and.((ii.gt.3))) cycle
825 [ + + + + ]: 108080 : if ((fatom.eq.gatom).and.((ii.eq.2).or.(ii.eq.5).or.(ii.eq.6))) cycle
826 [ + + + + : 97480 : if ((eatom.eq.fatom).and.(fatom.eq.gatom).and.(ii.gt.1)) cycle
- + ]
827 : 97480 : interactions(watom,ishell)=interactions(watom,ishell)+1
828 : : !FB write(std_out,*) 'For ishell and eatom=',ishell,watom
829 : : !FB write(std_out,*) ' --> the number of interactions in the shell is=',interactions(watom,ishell)
830 [ - + ]: 97480 : if (interactions(watom,ishell).gt.Shell3at%neighbours(watom,ishell)%n_interactions) then
831 : 0 : write(std_out,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>'
832 : 0 : write(std_out,*) ' >>>>>> Verify that the Rcut used in the input file is lower '
833 : 0 : write(std_out,*) ' >>>>>> than half of the smallest lattice parameter'
834 : 0 : write(std_out,*) ' >>>>>> Solution : Reduce the Rcut parameter'
835 : 0 : write(std_out,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>'
836 : 0 : ABI_ERROR('The interaction number index is greater than the interaction number max computed previously (3rd order)')
837 : : end if
838 : 97480 : Shell3at%neighbours(watom,ishell)%atomj_in_shell(interactions(watom,ishell))=xatom
839 : 97480 : Shell3at%neighbours(watom,ishell)%atomk_in_shell(interactions(watom,ishell))=yatom
840 : 97480 : Shell3at%neighbours(watom,ishell)%sym_in_shell(interactions(watom,ishell))=Isym3at(1)
841 : 117891768 : Shell3at%neighbours(watom,ishell)%transpose_in_shell(interactions(watom,ishell))=ii
842 : : end do !ii
843 : : end do !katom
844 : : end do !jatom
845 : : end do !iatom
846 : : ! Check that each interaction has different symmetry per shell
847 [ + + ]: 73 : do ishell=1,nshell
848 [ + + ]: 9913 : do iatom=1,natom
849 [ + + ]: 9840 : if (Shell3at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
850 [ + + ]: 97553 : do ninter=1,Shell3at%neighbours(iatom,ishell)%n_interactions-1
851 [ + + ]: 1170728 : do tmpinter=ninter+1,Shell3at%neighbours(iatom,ishell)%n_interactions
852 [ + + ]: 1071952 : if (Shell3at%neighbours(iatom,ishell)%sym_in_shell( ninter).eq.&
853 : 88936 : & Shell3at%neighbours(iatom,ishell)%sym_in_shell(tmpinter)) then
854 [ - + ]: 61184 : if (Shell3at%neighbours(iatom,ishell)%transpose_in_shell( ninter).ne.&
855 : : & Shell3at%neighbours(iatom,ishell)%transpose_in_shell(tmpinter)) cycle
856 : 0 : write(std_out,'(a,2(1x,i5))') 'For ishell and iatom =',ishell,iatom
857 : 0 : write(std_out,'(a,i5,a,i5,a,i5)') ' the interactions ',ninter,&
858 : 0 : & ' and ',tmpinter,' have both the same symmetry isym=',Shell3at%neighbours(iatom,ishell)%sym_in_shell( ninter)
859 : 0 : ABI_ERROR('Some interactions are equals due to the symmetry')
860 : : end if
861 : : end do
862 : : end do
863 : : end do
864 : : end do
865 : : ! Check that each equivalent shell has the same set of interactions
866 [ + + ]: 73 : do ishell=1,nshell
867 [ + + ]: 9913 : do iatom=1,natom
868 [ - + ]: 9840 : if (Shell3at%neighbours(mod(iatom-1,natom_unitcell)+1,ishell)%n_interactions.ne.&
869 : 57 : & Shell3at%neighbours( iatom,ishell)%n_interactions) then
870 : 0 : ABI_ERROR('The interaction number index is not equal to the interaction number max computed previously (2)')
871 : : end if
872 : : !BEGIN DEBUG
873 : : ! iatref=Shell3at%iatref(ishell)
874 : : ! jatref=Shell3at%jatref(ishell)
875 : : ! katref=Shell3at%katref(ishell)
876 : : ! if (Shell3at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
877 : : ! do ninter=1,Shell3at%neighbours(iatom,ishell)%n_interactions
878 : : ! jatom=Shell3at%neighbours(iatom,ishell)%atomj_in_shell(ninter)
879 : : ! katom=Shell3at%neighbours(iatom,ishell)%atomk_in_shell(ninter)
880 : : ! isym =Shell3at%neighbours(iatom,ishell)%sym_in_shell(ninter)
881 : : ! itrans=Shell3at%neighbours(iatom,ishell)%transpose_in_shell(ninter)
882 : : ! vectj(:)=zero ; vectk(:)=zero ; vect1(:)=zero ; vect2(:)=zero
883 : : ! do ii=1,3
884 : : ! do jj=1,3
885 : : ! vectj(ii)=vectj(ii)+Sym%S_ref(ii,jj,isym,1)*SC%dij(iatref,jatref,jj+1)
886 : : ! vectk(ii)=vectk(ii)+Sym%S_ref(ii,jj,isym,1)*SC%dij(iatref,katref,jj+1)
887 : : ! end do
888 : : ! end do
889 : : ! if (itrans==1) then ; vect1(:)= vectj(:) ; vect2(:)= vectk(:) ; endif !\Phi3_ijk
890 : : ! if (itrans==2) then ; vect1(:)= vectk(:) ; vect2(:)= vectj(:) ; endif !\Phi3_ikj
891 : : ! if (itrans==3) then ; vect1(:)=-vectj(:) ; vect2(:)= vectk(:)-vectj(:) ; endif !\Phi3_jik
892 : : ! if (itrans==4) then ; vect1(:)= vectk(:)-vectj(:) ; vect2(:)=-vectj(:) ; endif !\Phi3_jki
893 : : ! if (itrans==5) then ; vect1(:)=-vectk(:) ; vect2(:)= vectj(:)-vectk(:) ; endif !\Phi3_kij
894 : : ! if (itrans==6) then ; vect1(:)= vectj(:)-vectk(:) ; vect2(:)=-vectk(:) ; endif !\Phi3_kji
895 : : ! do ii=1,3
896 : : ! if ((abs(SC%dij(iatom,jatom,ii+1)-vect1(ii)).gt.tol8).or.&
897 : : !& (abs(SC%dij(iatom,katom,ii+1)-vect2(ii)).gt.tol8)) then
898 : : ! write(std_out,'(a,4(x,i5))') 'For ishell, iatom, jatom, katom =',ishell,iatom,jatom,katom
899 : : ! write(std_out,'(a,5(x,i5))') ' with isym, itrans, iatref, jatref, katref = ',isym,itrans,iatref,jatref,katref
900 : : ! ABI_ERROR('We do not recover the triplet with the symmetry found')
901 : : ! end if
902 : : ! end do !ii
903 : : ! end do !ninter
904 : : !END DEBUG
905 : : end do !natom
906 : : end do !nshell
907 [ - + ]: 16 : ABI_FREE(interactions)
908 : :
909 : : ! Find the number of coefficients of the (3x3x3) Phi3 for a given shell
910 [ - + + - : 43197 : ABI_CALLOC(Shell3at%proj, (norder,norder,nshell))
+ - + - -
+ - + - +
+ + + + +
+ ]
911 [ - + + - : 89 : ABI_MALLOC(Shell3at%ncoeff ,(nshell)); Shell3at%ncoeff(:)=zero
- + - + -
+ + + ]
912 [ - + + - : 89 : ABI_MALLOC(Shell3at%ncoeff_prev,(nshell)); Shell3at%ncoeff_prev(:)=zero
- + - + -
+ + + ]
913 : 16 : write(ab_out,*) 'Number of shells=',nshell
914 : 16 : write(ab_out,*) '============================================================================'
915 [ + - ]: 16 : if (MPIdata%iam_master) open(unit=16,file=trim(Invar%output_prefix)//'_nbcoeff-phi3.dat')
916 : 16 : ncoeff_prev=0
917 [ + + ]: 73 : do ishell=1,nshell
918 : 57 : ncoeff=0
919 : 57 : iatref=Shell3at%iatref(ishell)
920 : 57 : jatref=Shell3at%jatref(ishell)
921 : 57 : katref=Shell3at%katref(ishell)
922 : 57 : write(ab_out,*) 'Shell number:',ishell
923 : 57 : write(ab_out,'(a,i5,a,i5,a,f16.10)') ' Between atom',iatref,' and ',jatref,' the distance is=',SC%dij(iatref,jatref,1)
924 : 57 : write(ab_out,'(a,i5,a,i5,a,f16.10)') ' Between atom',jatref,' and ',katref,' the distance is=',SC%dij(jatref,katref,1)
925 : 57 : write(ab_out,'(a,i5,a,i5,a,f16.10)') ' Between atom',katref,' and ',iatref,' the distance is=',SC%dij(katref,iatref,1)
926 : 57 : call tdep_calc_nbcoeff(SC,MPIdata,iatref,ishell,jatref,katref,1,ncoeff,norder,nshell,order,Shell3at%proj,Sym)
927 : 57 : Shell3at%ncoeff (ishell)=ncoeff
928 : 57 : Shell3at%ncoeff_prev(ishell)=ncoeff_prev
929 : 57 : ncoeff_prev=ncoeff_prev+ncoeff
930 : 57 : write(ab_out,*)' Number of independant coefficients in this shell=',ncoeff
931 : 57 : write(ab_out,*)' Number of interactions in this shell=',Shell3at%neighbours(iatref,ishell)%n_interactions
932 : : !FB write(ab_out,*)' The ratio is=',dfloat(Shell3at%neighbours(iatref,ishell)%n_interactions)/dfloat(ncoeff)
933 : 73 : write(ab_out,*) '============================================================================'
934 : : end do
935 : 16 : write(ab_out,*)' >>>>>> Total number of coefficients at the third order=',ncoeff_prev
936 [ + - ]: 16 : if (MPIdata%iam_master) close(16)
937 : 16 : Shell3at%ntotcoeff=ncoeff_prev
938 : : !BeginFB
939 : : !FB open(unit=91,file='Shell3at.dat')
940 : : !FB write(91,*) Shell3at%nshell
941 : : !FB do ishell=1,Shell3at%nshell
942 : : !FB write(91,*) Shell3at%ncoeff(ishell)
943 : : !FB write(91,*) Shell3at%ncoeff_prev(ishell)
944 : : !FB write(91,*) Shell3at%iatref(ishell)
945 : : !FB write(91,*) Shell3at%jatref(ishell)
946 : : !FB write(91,*) Shell3at%katref(ishell)
947 : : !FB do iatom=1,Invar%natom
948 : : !FB write(91,*) Shell3at%neighbours(iatom,ishell)%n_interactions
949 : : !FB do ii=1,Shell3at%neighbours(iatom,ishell)%n_interactions
950 : : !FB write(91,*) Shell3at%neighbours(iatom,ishell)%sym_in_shell(ii)
951 : : !FB write(91,*) Shell3at%neighbours(iatom,ishell)%transpose_in_shell(ii)
952 : : !FB write(91,*) Shell3at%neighbours(iatom,ishell)%atomj_in_shell(ii)
953 : : !FB write(91,*) Shell3at%neighbours(iatom,ishell)%atomk_in_shell(ii)
954 : : !FB end do
955 : : !FB end do
956 : : !FB end do
957 : : !FB close(91)
958 : : !EndFB
959 : :
960 : 16 : end subroutine tdep_init_shell3at
961 : : !!***
962 : :
963 : : !====================================================================================================
964 : :
965 : : !!****f* ABINIT/m_tdep_shell/tdep_init_shell4at
966 : : !! NAME
967 : : !! tdep_init_shell4at
968 : : !!
969 : : !! FUNCTION
970 : : !! Build the shells of fourth-order interactions: identify the sets of
971 : : !! quadruplets of atoms that are equivalent by symmetry (within
972 : : !! the cutoff radius Invar%rcut4), store for each atom
973 : : !! and each shell the list of interactions and the symmetry operations
974 : : !! relating them to the reference interaction of the shell, and compute
975 : : !! (see tdep_calc_nbcoeff) the number of independent IFC coefficients of
976 : : !! each shell together with the corresponding projectors.
977 : : !! The details are written in the output file and in the
978 : : !! _nbcoeff-phi4.dat file.
979 : : !!
980 : : !! INPUTS
981 : : !! Invar = Input object containing the input variables.
982 : : !! MPIdata = Info on MPI parallelism.
983 : : !! SC = Supercell object, used for the interatomic distances (SC%dij).
984 : : !!
985 : : !! OUTPUT
986 : : !! Shell4at = Shell object for fourth-order interactions.
987 : : !!
988 : : !! SIDE EFFECTS
989 : : !! Sym = Symmetries object describing all the symmetry operations of the crystal.
990 : : !!
991 : : !! SOURCE
992 : :
993 : 9 : subroutine tdep_init_shell4at(Invar,MPIdata,Shell4at,SC,Sym)
994 : :
995 : : type(Shell_type),intent(inout) :: Shell4at
996 : : type(atdep_dataset_type),intent(in) :: Invar
997 : : type(MPI_enreg_type), intent(in) :: MPIdata
998 : : type(tdep_Supercell_type),intent(in) :: SC
999 : : type(Symmetries_type),intent(inout) :: Sym
1000 : :
1001 : : integer :: ii,ishell,iatom,jatom,katom,latom,eatom,fatom,gatom,hatom,iatref,jatref,katref,latref
1002 : : integer :: natom,natom_unitcell,watom,xatom,yatom,zatom,ninteractions,ncoeff,ncoeff_prev,nshell_tmp
1003 : : integer :: already_found,find_equivalent,ninter,iat_ref,jat_ref,kat_ref,lat_ref,tmpinter
1004 : : integer :: norder,order,nshell_max,nshell
1005 : : double precision :: norma1,norma2,norma3
1006 : : double precision :: normb1,normb2,normb3
1007 : : double precision :: normc1,normc2,normc3
1008 : : integer :: Isym4at(2)
1009 : 9 : integer, allocatable :: atref(:,:),interactions(:,:)
1010 : :
1011 : 9 : natom = Invar%natom
1012 : 9 : natom_unitcell = Invar%natom_unitcell
1013 : 9 : nshell_max = Invar%nshell_max
1014 : 9 : order = 4
1015 : 9 : norder = 81
1016 : 9 : Shell4at%order = order
1017 : 9 : Shell4at%norder = norder
1018 : 9 : Shell4at%natom = natom
1019 : :
1020 : 9 : write(ab_out,*) ' '
1021 : 9 : write(ab_out,*) '#############################################################################'
1022 : 9 : write(ab_out,*) '###### FOURTH ORDER : find the number of coefficients ########################'
1023 : 9 : write(ab_out,*) '#############################################################################'
1024 : :
1025 : : ! 1/ Identify the shells
1026 [ - + + - : 13536 : ABI_CALLOC(interactions,(natom_unitcell,nshell_max))
- + + - +
- + - - +
- + + + +
+ ]
1027 [ - + + - : 18063 : ABI_CALLOC(atref,(nshell_max,4))
+ - + - -
+ - + + +
+ + ]
1028 [ - + + - : 45 : ABI_CALLOC(Shell4at%ishell_self,(natom_unitcell))
+ - + - -
+ - + - +
+ + ]
1029 : 9 : nshell_tmp=0
1030 [ + + ]: 27 : do iatom=1,natom_unitcell
1031 [ + + ]: 3283 : do jatom=1,natom
1032 : : ! Interactions are only computed until Rcut4 in order to have complete shell of neighbours.
1033 : : ! Otherwise the symetries are broken.
1034 [ + + ]: 3256 : if (SC%dij(iatom,jatom,1).gt.(Invar%rcut4*0.99)) cycle
1035 [ + + ]: 23792 : do katom=1,natom
1036 [ + + ]: 23640 : if (SC%dij(iatom,katom,1).gt.(Invar%rcut4*0.99)) cycle
1037 : : ! WARNING: distance(j,k).ne.|djk| due to the inbox procedure when computing distance(j,k).
1038 : : ! So, compute |djk| using vec(ij) and vec(ik).
1039 : : norma1=dsqrt((SC%dij(iatom,katom,2)-SC%dij(iatom,jatom,2))**2+&
1040 : : & (SC%dij(iatom,katom,3)-SC%dij(iatom,jatom,3))**2+&
1041 : 1490 : & (SC%dij(iatom,katom,4)-SC%dij(iatom,jatom,4))**2)
1042 [ + + ]: 1490 : if (norma1 .gt.(Invar%rcut4*0.99)) cycle
1043 [ + + ]: 131774 : do latom=1,natom
1044 : : !FB write(std_out,*) 'NEW COORD1 : iatom,jatom,katom=',iatom,jatom,katom,latom
1045 [ + + ]: 127768 : if (SC%dij(iatom,latom,1).gt.(Invar%rcut4*0.99)) cycle
1046 : : norma2=dsqrt((SC%dij(iatom,latom,2)-SC%dij(iatom,jatom,2))**2+&
1047 : : & (SC%dij(iatom,latom,3)-SC%dij(iatom,jatom,3))**2+&
1048 : 9098 : & (SC%dij(iatom,latom,4)-SC%dij(iatom,jatom,4))**2)
1049 [ + + ]: 9098 : if (norma2 .gt.(Invar%rcut4*0.99)) cycle
1050 : : norma3=dsqrt((SC%dij(iatom,latom,2)-SC%dij(iatom,katom,2))**2+&
1051 : : & (SC%dij(iatom,latom,3)-SC%dij(iatom,katom,3))**2+&
1052 : 5026 : & (SC%dij(iatom,latom,4)-SC%dij(iatom,katom,4))**2)
1053 [ + + ]: 5026 : if (norma3 .gt.(Invar%rcut4*0.99)) cycle
1054 : :
1055 [ + + ]: 3518 : if (nshell_tmp.eq.0) then
1056 [ + + ]: 45 : atref(1,:)=1
1057 : 9 : nshell_tmp=nshell_tmp+1
1058 [ + - + - ]: 9 : if ((iatom.eq.jatom).and.(jatom.eq.katom).and.(katom.eq.latom)) Shell4at%ishell_self(iatom)=nshell_tmp
1059 : 9 : interactions(iatom,nshell_tmp)=interactions(iatom,nshell_tmp)+1
1060 : 9 : cycle
1061 : : end if
1062 : 3509 : find_equivalent=0
1063 [ + + ]: 23130 : do ishell=1,nshell_tmp
1064 : 23090 : iat_ref=atref(ishell,1) ; jat_ref=atref(ishell,2) ; kat_ref=atref(ishell,3) ; lat_ref=atref(ishell,4)
1065 : : normb1=dsqrt((SC%dij(iat_ref,kat_ref,2)-SC%dij(iat_ref,jat_ref,2))**2+&
1066 : : & (SC%dij(iat_ref,kat_ref,3)-SC%dij(iat_ref,jat_ref,3))**2+&
1067 : 23090 : & (SC%dij(iat_ref,kat_ref,4)-SC%dij(iat_ref,jat_ref,4))**2)
1068 : : normb2=dsqrt((SC%dij(iat_ref,lat_ref,2)-SC%dij(iat_ref,jat_ref,2))**2+&
1069 : : & (SC%dij(iat_ref,lat_ref,3)-SC%dij(iat_ref,jat_ref,3))**2+&
1070 : 23090 : & (SC%dij(iat_ref,lat_ref,4)-SC%dij(iat_ref,jat_ref,4))**2)
1071 : : normb3=dsqrt((SC%dij(iat_ref,lat_ref,2)-SC%dij(iat_ref,kat_ref,2))**2+&
1072 : : & (SC%dij(iat_ref,lat_ref,3)-SC%dij(iat_ref,kat_ref,3))**2+&
1073 : 23090 : & (SC%dij(iat_ref,lat_ref,4)-SC%dij(iat_ref,kat_ref,4))**2)
1074 [ + + ]: 519714 : do ii=1,24
1075 [ + + ]: 500093 : if (ii.eq.1 ) then ; eatom=iatom ; fatom=jatom ; gatom=katom ; hatom=latom ; end if !ijkl
1076 [ + + ]: 500093 : if (ii.eq.2 ) then ; eatom=iatom ; fatom=katom ; gatom=jatom ; hatom=latom ; end if !ikjl
1077 [ + + ]: 500093 : if (ii.eq.3 ) then ; eatom=jatom ; fatom=iatom ; gatom=katom ; hatom=latom ; end if !jikl
1078 [ + + ]: 500093 : if (ii.eq.4 ) then ; eatom=jatom ; fatom=katom ; gatom=iatom ; hatom=latom ; end if !jkil
1079 [ + + ]: 500093 : if (ii.eq.5 ) then ; eatom=katom ; fatom=iatom ; gatom=jatom ; hatom=latom ; end if !kijl
1080 [ + + ]: 500093 : if (ii.eq.6 ) then ; eatom=katom ; fatom=jatom ; gatom=iatom ; hatom=latom ; end if !kjil
1081 : :
1082 [ + + ]: 500093 : if (ii.eq.7 ) then ; eatom=iatom ; fatom=jatom ; gatom=latom ; hatom=katom ; end if !ijlk
1083 [ + + ]: 500093 : if (ii.eq.8 ) then ; eatom=iatom ; fatom=katom ; gatom=latom ; hatom=jatom ; end if !iklj
1084 [ + + ]: 500093 : if (ii.eq.9 ) then ; eatom=jatom ; fatom=iatom ; gatom=latom ; hatom=katom ; end if !jilk
1085 [ + + ]: 500093 : if (ii.eq.10) then ; eatom=jatom ; fatom=katom ; gatom=latom ; hatom=iatom ; end if !jkli
1086 [ + + ]: 500093 : if (ii.eq.11) then ; eatom=katom ; fatom=iatom ; gatom=latom ; hatom=jatom ; end if !kilj
1087 [ + + ]: 500093 : if (ii.eq.12) then ; eatom=katom ; fatom=jatom ; gatom=latom ; hatom=iatom ; end if !kjli
1088 : :
1089 [ + + ]: 500093 : if (ii.eq.13) then ; eatom=iatom ; fatom=latom ; gatom=jatom ; hatom=katom ; end if !iljk
1090 [ + + ]: 500093 : if (ii.eq.14) then ; eatom=iatom ; fatom=latom ; gatom=katom ; hatom=jatom ; end if !ilkj
1091 [ + + ]: 500093 : if (ii.eq.15) then ; eatom=jatom ; fatom=latom ; gatom=iatom ; hatom=katom ; end if !jlik
1092 [ + + ]: 500093 : if (ii.eq.16) then ; eatom=jatom ; fatom=latom ; gatom=katom ; hatom=iatom ; end if !jlki
1093 [ + + ]: 500093 : if (ii.eq.17) then ; eatom=katom ; fatom=latom ; gatom=iatom ; hatom=jatom ; end if !klij
1094 [ + + ]: 500093 : if (ii.eq.18) then ; eatom=katom ; fatom=latom ; gatom=jatom ; hatom=iatom ; end if !klji
1095 : :
1096 [ + + ]: 500093 : if (ii.eq.19) then ; eatom=latom ; fatom=iatom ; gatom=jatom ; hatom=katom ; end if !lijk
1097 [ + + ]: 500093 : if (ii.eq.20) then ; eatom=latom ; fatom=iatom ; gatom=katom ; hatom=jatom ; end if !likj
1098 [ + + ]: 500093 : if (ii.eq.21) then ; eatom=latom ; fatom=jatom ; gatom=iatom ; hatom=katom ; end if !ljik
1099 [ + + ]: 500093 : if (ii.eq.22) then ; eatom=latom ; fatom=jatom ; gatom=katom ; hatom=iatom ; end if !ljki
1100 [ + + ]: 500093 : if (ii.eq.23) then ; eatom=latom ; fatom=katom ; gatom=iatom ; hatom=jatom ; end if !lkij
1101 [ + + ]: 500093 : if (ii.eq.24) then ; eatom=latom ; fatom=katom ; gatom=jatom ; hatom=iatom ; end if !lkji
1102 : :
1103 : : normc1=dsqrt((SC%dij(eatom,gatom,2)-SC%dij(eatom,fatom,2))**2+&
1104 : : & (SC%dij(eatom,gatom,3)-SC%dij(eatom,fatom,3))**2+&
1105 : 500093 : & (SC%dij(eatom,gatom,4)-SC%dij(eatom,fatom,4))**2)
1106 : : normc2=dsqrt((SC%dij(eatom,hatom,2)-SC%dij(eatom,fatom,2))**2+&
1107 : : & (SC%dij(eatom,hatom,3)-SC%dij(eatom,fatom,3))**2+&
1108 : 500093 : & (SC%dij(eatom,hatom,4)-SC%dij(eatom,fatom,4))**2)
1109 : : normc3=dsqrt((SC%dij(eatom,hatom,2)-SC%dij(eatom,gatom,2))**2+&
1110 : : & (SC%dij(eatom,hatom,3)-SC%dij(eatom,gatom,3))**2+&
1111 : 500093 : & (SC%dij(eatom,hatom,4)-SC%dij(eatom,gatom,4))**2)
1112 : : if ((abs(SC%dij(iat_ref,jat_ref,1)-SC%dij(eatom,fatom,1)).lt.1.d-6).and.&
1113 : : & (abs(normb1 -normc1 ).lt.1.d-6).and.&
1114 : : & (abs(normb2 -normc2 ).lt.1.d-6).and.&
1115 : : & (abs(normb3 -normc3 ).lt.1.d-6).and.&
1116 [ + + + + : 500093 : & (abs(SC%dij(iat_ref,kat_ref,1)-SC%dij(eatom,gatom,1)).lt.1.d-6).and.&
+ + + + +
- + - ]
1117 : 19621 : & (abs(SC%dij(iat_ref,lat_ref,1)-SC%dij(eatom,hatom,1)).lt.1.d-6)) then
1118 : 6469 : Isym4at(:)=0
1119 : 6469 : call tdep_SearchS_4at(SC,iat_ref,jat_ref,kat_ref,lat_ref,eatom,fatom,gatom,hatom,Isym4at,Sym)
1120 [ + + ]: 6469 : if (Isym4at(2).eq.1) find_equivalent=1
1121 : : if (find_equivalent.eq.1) then
1122 : 3469 : interactions(iatom,ishell)=interactions(iatom,ishell)+1
1123 : : !FB write(std_out,*) 'The number of interactions in this shell is=',ishell,interactions(iatom,ishell)
1124 : 3469 : exit
1125 : : end if
1126 : : end if
1127 : : end do !ii
1128 [ + + ]: 23130 : if (find_equivalent.eq.1) exit
1129 : : end do !ishell
1130 [ + + ]: 27149 : if (find_equivalent.eq.0) then
1131 : 40 : nshell_tmp=nshell_tmp+1
1132 [ - + ]: 40 : if (nshell_tmp.gt.nshell_max) then
1133 : 0 : ABI_ERROR('The shell number index is greater than the shell number max defined in the code')
1134 : : end if
1135 [ + + + + ]: 40 : if ((iatom.eq.jatom).and.(jatom.eq.katom).and.(katom.eq.latom)) Shell4at%ishell_self(iatom)=nshell_tmp
1136 : 40 : interactions(iatom,nshell_tmp)=interactions(iatom,nshell_tmp)+1
1137 : 40 : atref(nshell_tmp,1)=iatom
1138 : 40 : atref(nshell_tmp,2)=jatom
1139 : 40 : atref(nshell_tmp,3)=katom
1140 : 40 : atref(nshell_tmp,4)=latom
1141 : : !FB write(std_out,'(a,1x,5(i5,1x))') 'NEW SHELL1 : nshell_tmp,iatom,jatom,katom,latom=',nshell_tmp,iatom,jatom,katom,latom
1142 : : end if
1143 : : end do !latom
1144 : : end do !katom
1145 : : end do !jatom
1146 : : end do !iatom
1147 [ - + ]: 9 : ABI_FREE(atref)
1148 : :
1149 : : ! 2/ Allocate the datatype Shell4at%...
1150 : 9 : Shell4at%nshell = nshell_tmp
1151 : 9 : nshell = Shell4at%nshell
1152 [ - + + - : 8565 : ABI_MALLOC(Shell4at%neighbours,(natom,nshell))
- + + - +
- + - - +
- + - + +
+ + + +
- ]
1153 [ - + + - : 76 : ABI_CALLOC(Shell4at%iatref,(nshell))
+ - + - -
+ - + - +
+ + ]
1154 [ - + + - : 67 : ABI_CALLOC(Shell4at%jatref,(nshell))
- + - + -
+ + + ]
1155 [ - + + - : 67 : ABI_CALLOC(Shell4at%katref,(nshell))
- + - + -
+ + + ]
1156 [ - + + - : 67 : ABI_CALLOC(Shell4at%latref,(nshell))
- + - + -
+ + + ]
1157 [ + + ]: 58 : do ishell=1,nshell
1158 [ + + ]: 8538 : do iatom=1,natom
1159 : 8480 : ninteractions=interactions(mod(iatom-1,natom_unitcell)+1,ishell)
1160 : 8480 : Shell4at%neighbours(iatom,ishell)%n_interactions=ninteractions
1161 [ + + ]: 8480 : if (ninteractions.eq.0) cycle
1162 [ - + + - : 22368 : ABI_MALLOC(Shell4at%neighbours(iatom,ishell)%atomj_in_shell,(ninteractions))
+ - + - -
+ - + -
+ ]
1163 [ - + + - : 14912 : ABI_MALLOC(Shell4at%neighbours(iatom,ishell)%atomk_in_shell,(ninteractions))
- + - + -
+ ]
1164 [ - + + - : 14912 : ABI_MALLOC(Shell4at%neighbours(iatom,ishell)%atoml_in_shell,(ninteractions))
- + - + -
+ ]
1165 [ - + + - : 14912 : ABI_MALLOC(Shell4at%neighbours(iatom,ishell)%sym_in_shell,(ninteractions))
- + - + -
+ ]
1166 [ - + + - : 14912 : ABI_MALLOC(Shell4at%neighbours(iatom,ishell)%transpose_in_shell,(ninteractions))
- + - + -
+ ]
1167 [ + + ]: 257452 : Shell4at%neighbours(iatom,ishell)%atomj_in_shell(:)=zero
1168 [ + + ]: 257452 : Shell4at%neighbours(iatom,ishell)%atomk_in_shell(:)=zero
1169 [ + + ]: 257452 : Shell4at%neighbours(iatom,ishell)%atoml_in_shell(:)=zero
1170 [ + + ]: 257452 : Shell4at%neighbours(iatom,ishell)%sym_in_shell(:)=zero
1171 [ + + ]: 258525 : Shell4at%neighbours(iatom,ishell)%transpose_in_shell(:)=zero
1172 : : end do
1173 : : end do
1174 [ - + ]: 9 : ABI_FREE(interactions)
1175 : :
1176 : : ! 3/ Store the index of the (couple of) atoms included in each shell
1177 : : ! 4/ Store the reference (couple of) atoms for each shell
1178 : : ! 5/ Compute the symetry operation between the reference (couple of) atoms and another one
1179 [ - + + - : 8565 : ABI_MALLOC(interactions,(natom,nshell)) ; interactions(:,:)=0
- + + - +
- + - - +
- + + + +
+ ]
1180 : 9 : nshell_tmp=0
1181 [ + + ]: 1701 : do iatom=1,natom
1182 [ + + ]: 339061 : do jatom=1,natom
1183 [ + + ]: 337360 : if (SC%dij(iatom,jatom,1).gt.(Invar%rcut4*0.99)) cycle
1184 [ + + ]: 1842072 : do katom=1,natom
1185 [ + + ]: 1830288 : if (SC%dij(iatom,katom,1).gt.(Invar%rcut4*0.99)) cycle
1186 : : ! WARNING: distance(j,k).ne.|djk| due to the inbox procedure when computing distance(j,k).
1187 : : ! So, compute |djk| using vec(ij) and vec(ik).
1188 : : norma1=dsqrt((SC%dij(iatom,katom,2)-SC%dij(iatom,jatom,2))**2+&
1189 : : & (SC%dij(iatom,katom,3)-SC%dij(iatom,jatom,3))**2+&
1190 : 106748 : & (SC%dij(iatom,katom,4)-SC%dij(iatom,jatom,4))**2)
1191 [ + + ]: 106748 : if (norma1 .gt.(Invar%rcut4*0.99)) cycle
1192 [ + + ]: 9711788 : do latom=1,natom
1193 : : !FB write(std_out,*) 'NEW COORD2 : iatom,jatom,katom=',iatom,jatom,katom,latom
1194 [ + + ]: 9320464 : if (SC%dij(iatom,latom,1).gt.(Invar%rcut4*0.99)) cycle
1195 : : norma2=dsqrt((SC%dij(iatom,latom,2)-SC%dij(iatom,jatom,2))**2+&
1196 : : & (SC%dij(iatom,latom,3)-SC%dij(iatom,jatom,3))**2+&
1197 : 645852 : & (SC%dij(iatom,latom,4)-SC%dij(iatom,jatom,4))**2)
1198 [ + + ]: 645852 : if (norma2 .gt.(Invar%rcut4*0.99)) cycle
1199 : : norma3=dsqrt((SC%dij(iatom,latom,2)-SC%dij(iatom,katom,2))**2+&
1200 : : & (SC%dij(iatom,latom,3)-SC%dij(iatom,katom,3))**2+&
1201 : 356924 : & (SC%dij(iatom,latom,4)-SC%dij(iatom,katom,4))**2)
1202 [ + + ]: 356924 : if (norma3 .gt.(Invar%rcut4*0.99)) cycle
1203 [ + + ]: 249996 : if (nshell_tmp.eq.0) then
1204 : 9 : nshell_tmp=nshell_tmp+1
1205 : 9 : interactions(iatom,nshell_tmp)=1
1206 : 9 : Shell4at%iatref(nshell_tmp)=iatom
1207 : 9 : Shell4at%jatref(nshell_tmp)=jatom
1208 : 9 : Shell4at%katref(nshell_tmp)=katom
1209 : 9 : Shell4at%latref(nshell_tmp)=latom
1210 : 9 : Shell4at%neighbours(iatom,nshell_tmp)%atomj_in_shell(interactions(iatom,nshell_tmp))=jatom
1211 : 9 : Shell4at%neighbours(iatom,nshell_tmp)%atomk_in_shell(interactions(iatom,nshell_tmp))=katom
1212 : 9 : Shell4at%neighbours(iatom,nshell_tmp)%atoml_in_shell(interactions(iatom,nshell_tmp))=latom
1213 : 9 : Shell4at%neighbours(iatom,nshell_tmp)%sym_in_shell(interactions(iatom,nshell_tmp))=1
1214 : 9 : Shell4at%neighbours(iatom,nshell_tmp)%transpose_in_shell(interactions(iatom,nshell_tmp))=1
1215 : 9 : cycle
1216 : : end if
1217 : :
1218 : : ! Search if the quadruplet has already been classified
1219 : 249987 : find_equivalent=0
1220 [ + + ]: 2959140 : do ishell=1,nshell_tmp
1221 [ + + ]: 90919726 : do ninter=1,interactions(iatom,ishell)
1222 : : if ((Shell4at%neighbours(iatom,ishell)%atomj_in_shell(ninter).eq.jatom).and.&
1223 [ + + + + : 87960586 : & (Shell4at%neighbours(iatom,ishell)%atomk_in_shell(ninter).eq.katom).and.&
+ + ]
1224 : 90669739 : & (Shell4at%neighbours(iatom,ishell)%atoml_in_shell(ninter).eq.latom)) find_equivalent=1
1225 : : end do
1226 : : end do
1227 [ + + ]: 249987 : if (find_equivalent.eq.1) cycle
1228 : : ! Search if the quadruplet belongs to a shell already found
1229 [ + + ]: 163795 : do ishell=1,nshell_tmp
1230 : 163755 : iat_ref=Shell4at%iatref(ishell) ; jat_ref=Shell4at%jatref(ishell)
1231 : 163755 : kat_ref=Shell4at%katref(ishell) ; lat_ref=Shell4at%latref(ishell)
1232 : : normb1=dsqrt((SC%dij(iat_ref,kat_ref,2)-SC%dij(iat_ref,jat_ref,2))**2+&
1233 : : & (SC%dij(iat_ref,kat_ref,3)-SC%dij(iat_ref,jat_ref,3))**2+&
1234 : 163755 : & (SC%dij(iat_ref,kat_ref,4)-SC%dij(iat_ref,jat_ref,4))**2)
1235 : : normb2=dsqrt((SC%dij(iat_ref,lat_ref,2)-SC%dij(iat_ref,jat_ref,2))**2+&
1236 : : & (SC%dij(iat_ref,lat_ref,3)-SC%dij(iat_ref,jat_ref,3))**2+&
1237 : 163755 : & (SC%dij(iat_ref,lat_ref,4)-SC%dij(iat_ref,jat_ref,4))**2)
1238 : : normb3=dsqrt((SC%dij(iat_ref,lat_ref,2)-SC%dij(iat_ref,kat_ref,2))**2+&
1239 : : & (SC%dij(iat_ref,lat_ref,3)-SC%dij(iat_ref,kat_ref,3))**2+&
1240 : 163755 : & (SC%dij(iat_ref,lat_ref,4)-SC%dij(iat_ref,kat_ref,4))**2)
1241 [ + + ]: 3563337 : do ii=1,24
1242 [ + + ]: 3428489 : if (ii.eq.1 ) then ; eatom=iatom ; fatom=jatom ; gatom=katom ; hatom=latom ; end if !ijkl
1243 [ + + ]: 3428489 : if (ii.eq.2 ) then ; eatom=iatom ; fatom=katom ; gatom=jatom ; hatom=latom ; end if !ikjl
1244 [ + + ]: 3428489 : if (ii.eq.3 ) then ; eatom=jatom ; fatom=iatom ; gatom=katom ; hatom=latom ; end if !jikl
1245 [ + + ]: 3428489 : if (ii.eq.4 ) then ; eatom=jatom ; fatom=katom ; gatom=iatom ; hatom=latom ; end if !jkil
1246 [ + + ]: 3428489 : if (ii.eq.5 ) then ; eatom=katom ; fatom=iatom ; gatom=jatom ; hatom=latom ; end if !kijl
1247 [ + + ]: 3428489 : if (ii.eq.6 ) then ; eatom=katom ; fatom=jatom ; gatom=iatom ; hatom=latom ; end if !kjil
1248 : :
1249 [ + + ]: 3428489 : if (ii.eq.7 ) then ; eatom=iatom ; fatom=jatom ; gatom=latom ; hatom=katom ; end if !ijlk
1250 [ + + ]: 3428489 : if (ii.eq.8 ) then ; eatom=iatom ; fatom=katom ; gatom=latom ; hatom=jatom ; end if !iklj
1251 [ + + ]: 3428489 : if (ii.eq.9 ) then ; eatom=jatom ; fatom=iatom ; gatom=latom ; hatom=katom ; end if !jilk
1252 [ + + ]: 3428489 : if (ii.eq.10) then ; eatom=jatom ; fatom=katom ; gatom=latom ; hatom=iatom ; end if !jkli
1253 [ + + ]: 3428489 : if (ii.eq.11) then ; eatom=katom ; fatom=iatom ; gatom=latom ; hatom=jatom ; end if !kilj
1254 [ + + ]: 3428489 : if (ii.eq.12) then ; eatom=katom ; fatom=jatom ; gatom=latom ; hatom=iatom ; end if !kjli
1255 : :
1256 [ + + ]: 3428489 : if (ii.eq.13) then ; eatom=iatom ; fatom=latom ; gatom=jatom ; hatom=katom ; end if !iljk
1257 [ + + ]: 3428489 : if (ii.eq.14) then ; eatom=iatom ; fatom=latom ; gatom=katom ; hatom=jatom ; end if !ilkj
1258 [ + + ]: 3428489 : if (ii.eq.15) then ; eatom=jatom ; fatom=latom ; gatom=iatom ; hatom=katom ; end if !jlik
1259 [ + + ]: 3428489 : if (ii.eq.16) then ; eatom=jatom ; fatom=latom ; gatom=katom ; hatom=iatom ; end if !jlki
1260 [ + + ]: 3428489 : if (ii.eq.17) then ; eatom=katom ; fatom=latom ; gatom=iatom ; hatom=jatom ; end if !klij
1261 [ + + ]: 3428489 : if (ii.eq.18) then ; eatom=katom ; fatom=latom ; gatom=jatom ; hatom=iatom ; end if !klji
1262 : :
1263 [ + + ]: 3428489 : if (ii.eq.19) then ; eatom=latom ; fatom=iatom ; gatom=jatom ; hatom=katom ; end if !lijk
1264 [ + + ]: 3428489 : if (ii.eq.20) then ; eatom=latom ; fatom=iatom ; gatom=katom ; hatom=jatom ; end if !likj
1265 [ + + ]: 3428489 : if (ii.eq.21) then ; eatom=latom ; fatom=jatom ; gatom=iatom ; hatom=katom ; end if !ljik
1266 [ + + ]: 3428489 : if (ii.eq.22) then ; eatom=latom ; fatom=jatom ; gatom=katom ; hatom=iatom ; end if !ljki
1267 [ + + ]: 3428489 : if (ii.eq.23) then ; eatom=latom ; fatom=katom ; gatom=iatom ; hatom=jatom ; end if !lkij
1268 [ + + ]: 3428489 : if (ii.eq.24) then ; eatom=latom ; fatom=katom ; gatom=jatom ; hatom=iatom ; end if !lkji
1269 : :
1270 : : normc1=dsqrt((SC%dij(eatom,gatom,2)-SC%dij(eatom,fatom,2))**2+&
1271 : : & (SC%dij(eatom,gatom,3)-SC%dij(eatom,fatom,3))**2+&
1272 : 3428489 : & (SC%dij(eatom,gatom,4)-SC%dij(eatom,fatom,4))**2)
1273 : : normc2=dsqrt((SC%dij(eatom,hatom,2)-SC%dij(eatom,fatom,2))**2+&
1274 : : & (SC%dij(eatom,hatom,3)-SC%dij(eatom,fatom,3))**2+&
1275 : 3428489 : & (SC%dij(eatom,hatom,4)-SC%dij(eatom,fatom,4))**2)
1276 : : normc3=dsqrt((SC%dij(eatom,hatom,2)-SC%dij(eatom,gatom,2))**2+&
1277 : : & (SC%dij(eatom,hatom,3)-SC%dij(eatom,gatom,3))**2+&
1278 : 3428489 : & (SC%dij(eatom,hatom,4)-SC%dij(eatom,gatom,4))**2)
1279 : : if ((abs(SC%dij(iat_ref,jat_ref,1)-SC%dij(eatom,fatom,1)).lt.1.d-6).and.&
1280 : : & (abs(normb1 -normc1 ).lt.1.d-6).and.&
1281 : : & (abs(normb2 -normc2 ).lt.1.d-6).and.&
1282 : : & (abs(normb3 -normc3 ).lt.1.d-6).and.&
1283 [ + + + + : 3428489 : & (abs(SC%dij(iat_ref,kat_ref,1)-SC%dij(eatom,gatom,1)).lt.1.d-6).and.&
+ + + + +
- + - ]
1284 : 134848 : & (abs(SC%dij(iat_ref,lat_ref,1)-SC%dij(eatom,hatom,1)).lt.1.d-6)) then
1285 : 61919 : Isym4at(:)=0
1286 : 61919 : call tdep_SearchS_4at(SC,iat_ref,jat_ref,kat_ref,lat_ref,eatom,fatom,gatom,hatom,Isym4at,Sym)
1287 [ + + ]: 61919 : if (Isym4at(2).eq.1) then
1288 : : find_equivalent=1
1289 : : exit
1290 : : end if
1291 : : end if
1292 : : end do !ii
1293 [ + + ]: 163795 : if (find_equivalent.eq.1) exit
1294 : : end do !ishell
1295 : : ! The quadruplet belongs to a new shell
1296 [ + + ]: 28947 : if (find_equivalent.eq.0) then
1297 : 40 : nshell_tmp=nshell_tmp+1
1298 : : ! Check that the new shell is allowed
1299 [ - + ]: 40 : if (nshell_tmp.gt.nshell) then
1300 : 0 : ABI_ERROR('The shell number index is greater than the shell number max computed previously')
1301 : : end if
1302 : 40 : Shell4at%iatref(nshell_tmp)=iatom
1303 : 40 : Shell4at%jatref(nshell_tmp)=jatom
1304 : 40 : Shell4at%katref(nshell_tmp)=katom
1305 : 40 : Shell4at%latref(nshell_tmp)=latom
1306 : 40 : eatom=iatom ; fatom=jatom ; gatom=katom ; hatom=latom
1307 [ + + ]: 120 : Isym4at(:)=1
1308 : 40 : ishell=nshell_tmp
1309 : : !FB write(std_out,'(a,1x,5(i5,1x))') 'NEW SHELL2 : nshell_tmp,iatom,jatom,katom=',nshell_tmp,iatom,jatom,katom,latom
1310 : : end if
1311 : : ! Classify the informations of the quadruplet in Shell4at
1312 [ + + ]: 2553963 : do ii=1,24
1313 : : ! The Phi4 has to be symetric (transposition symetries)
1314 [ + + ]: 694728 : if (ii.eq.1 ) then ; watom=eatom ; xatom=fatom ; yatom=gatom ; zatom=hatom ; end if !ijkl
1315 [ + + ]: 694728 : if (ii.eq.2 ) then ; watom=eatom ; xatom=gatom ; yatom=fatom ; zatom=hatom ; end if !ikjl
1316 [ + + ]: 694728 : if (ii.eq.3 ) then ; watom=fatom ; xatom=eatom ; yatom=gatom ; zatom=hatom ; end if !jikl
1317 [ + + ]: 694728 : if (ii.eq.4 ) then ; watom=fatom ; xatom=gatom ; yatom=eatom ; zatom=hatom ; end if !jkil
1318 [ + + ]: 694728 : if (ii.eq.5 ) then ; watom=gatom ; xatom=eatom ; yatom=fatom ; zatom=hatom ; end if !kijl
1319 [ + + ]: 694728 : if (ii.eq.6 ) then ; watom=gatom ; xatom=fatom ; yatom=eatom ; zatom=hatom ; end if !kjil
1320 : :
1321 [ + + ]: 694728 : if (ii.eq.7 ) then ; watom=eatom ; xatom=fatom ; yatom=hatom ; zatom=gatom ; end if !ijlk
1322 [ + + ]: 694728 : if (ii.eq.8 ) then ; watom=eatom ; xatom=gatom ; yatom=hatom ; zatom=fatom ; end if !iklj
1323 [ + + ]: 694728 : if (ii.eq.9 ) then ; watom=fatom ; xatom=eatom ; yatom=hatom ; zatom=gatom ; end if !jilk
1324 [ + + ]: 694728 : if (ii.eq.10) then ; watom=fatom ; xatom=gatom ; yatom=hatom ; zatom=eatom ; end if !jkli
1325 [ + + ]: 694728 : if (ii.eq.11) then ; watom=gatom ; xatom=eatom ; yatom=hatom ; zatom=fatom ; end if !kilj
1326 [ + + ]: 694728 : if (ii.eq.12) then ; watom=gatom ; xatom=fatom ; yatom=hatom ; zatom=eatom ; end if !kjli
1327 : :
1328 [ + + ]: 694728 : if (ii.eq.13) then ; watom=eatom ; xatom=hatom ; yatom=fatom ; zatom=gatom ; end if !iljk
1329 [ + + ]: 694728 : if (ii.eq.14) then ; watom=eatom ; xatom=hatom ; yatom=gatom ; zatom=fatom ; end if !ilkj
1330 [ + + ]: 694728 : if (ii.eq.15) then ; watom=fatom ; xatom=hatom ; yatom=eatom ; zatom=gatom ; end if !jlik
1331 [ + + ]: 694728 : if (ii.eq.16) then ; watom=fatom ; xatom=hatom ; yatom=gatom ; zatom=eatom ; end if !jlki
1332 [ + + ]: 694728 : if (ii.eq.17) then ; watom=gatom ; xatom=hatom ; yatom=eatom ; zatom=fatom ; end if !klij
1333 [ + + ]: 694728 : if (ii.eq.18) then ; watom=gatom ; xatom=hatom ; yatom=fatom ; zatom=eatom ; end if !klji
1334 : :
1335 [ + + ]: 694728 : if (ii.eq.19) then ; watom=hatom ; xatom=eatom ; yatom=fatom ; zatom=gatom ; end if !lijk
1336 [ + + ]: 694728 : if (ii.eq.20) then ; watom=hatom ; xatom=eatom ; yatom=gatom ; zatom=fatom ; end if !likj
1337 [ + + ]: 694728 : if (ii.eq.21) then ; watom=hatom ; xatom=fatom ; yatom=eatom ; zatom=gatom ; end if !ljik
1338 [ + + ]: 694728 : if (ii.eq.22) then ; watom=hatom ; xatom=fatom ; yatom=gatom ; zatom=eatom ; end if !ljki
1339 [ + + ]: 694728 : if (ii.eq.23) then ; watom=hatom ; xatom=gatom ; yatom=eatom ; zatom=fatom ; end if !lkij
1340 [ + + ]: 694728 : if (ii.eq.24) then ; watom=hatom ; xatom=gatom ; yatom=fatom ; zatom=eatom ; end if !lkji
1341 : : ! Do not overwrite the Phi4_iikl, Phi4_ijil, Phi4_ijjl, Phi4_iiil... IFCs
1342 : : ! and avoid double counting of quadruplet interactions
1343 : 694728 : already_found=0
1344 [ + + ]: 22620552 : do ninter=1,interactions(watom,ishell)
1345 : : if ((Shell4at%neighbours(watom,ishell)%atomj_in_shell(ninter).eq.xatom).and.&
1346 [ + + + + : 22370565 : & (Shell4at%neighbours(watom,ishell)%atomk_in_shell(ninter).eq.yatom).and.&
+ + ]
1347 : 249987 : & (Shell4at%neighbours(watom,ishell)%atoml_in_shell(ninter).eq.zatom)) then
1348 : : !FB write(*,'(a,4(1x,i5))') 'FOR efgh =',eatom,fatom,gatom,hatom
1349 : : !FB write(*,'(a,3(1x,i5))') ' --> ishell,ninter,ninter_tot =',ishell,ninter,interactions(watom,ishell)
1350 : : !FB write(*,'(a,4(1x,i5))') ' --> ALREADY FOUND =',watom,xatom,yatom,zatom
1351 : : already_found=1
1352 : : exit
1353 : : end if
1354 : : end do
1355 [ + + ]: 694728 : if (already_found==1) cycle
1356 : 249987 : interactions(watom,ishell)=interactions(watom,ishell)+1
1357 : : !FB write(std_out,*) 'For ishell and eatom=',ishell,watom
1358 : : !FB write(std_out,*) ' --> the number of interactions in the shell is=',interactions(watom,ishell)
1359 [ - + ]: 249987 : if (interactions(watom,ishell).gt.Shell4at%neighbours(watom,ishell)%n_interactions) then
1360 : 0 : write(std_out,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>'
1361 : 0 : write(std_out,*) ' >>>>>> Verify that the Rcut used in the input file is lower '
1362 : 0 : write(std_out,*) ' >>>>>> than half of the smallest lattice parameter'
1363 : 0 : write(std_out,*) ' >>>>>> Solution : Reduce the Rcut parameter'
1364 : 0 : write(std_out,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>'
1365 : 0 : ABI_ERROR('The interaction number index is greater than the interaction number max computed previously (4th order)')
1366 : : end if
1367 : 249987 : Shell4at%neighbours(watom,ishell)%atomj_in_shell(interactions(watom,ishell))=xatom
1368 : 249987 : Shell4at%neighbours(watom,ishell)%atomk_in_shell(interactions(watom,ishell))=yatom
1369 : 249987 : Shell4at%neighbours(watom,ishell)%atoml_in_shell(interactions(watom,ishell))=zatom
1370 : 249987 : Shell4at%neighbours(watom,ishell)%sym_in_shell(interactions(watom,ishell))=Isym4at(1)
1371 : 10015183 : Shell4at%neighbours(watom,ishell)%transpose_in_shell(interactions(watom,ishell))=ii
1372 : : !DEBUG write(std_out,'(a,9(i5,x))') 'ishell,iatref,jatref,katref,iatom,atomj_in_shell,atomk_in_shell,isym,itrans=',&
1373 : : !DEBUG& ishell,Shell4at%iatref(ishell),Shell4at%jatref(ishell),Shell4at%katref(ishell),watom,xatom,yatom,Isym4at(1),ii
1374 : : end do !ii
1375 : : end do !latom
1376 : : end do !katom
1377 : : end do !jatom
1378 : : end do !iatom
1379 : : ! Check that each interaction has different symmetry per shell
1380 [ + + ]: 58 : do ishell=1,nshell
1381 [ + + ]: 8538 : do iatom=1,natom
1382 [ + + ]: 8480 : if (Shell4at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
1383 [ + + ]: 250045 : do ninter=1,Shell4at%neighbours(iatom,ishell)%n_interactions-1
1384 [ + + ]: 11103716 : do tmpinter=ninter+1,Shell4at%neighbours(iatom,ishell)%n_interactions
1385 [ + + ]: 10852696 : if (Shell4at%neighbours(iatom,ishell)%sym_in_shell( ninter).eq.&
1386 : 242540 : & Shell4at%neighbours(iatom,ishell)%sym_in_shell(tmpinter)) then
1387 [ - + ]: 594612 : if (Shell4at%neighbours(iatom,ishell)%transpose_in_shell( ninter).ne.&
1388 : : & Shell4at%neighbours(iatom,ishell)%transpose_in_shell(tmpinter)) cycle
1389 : 0 : write(std_out,'(a,2(1x,i5))') 'For ishell and iatom =',ishell,iatom
1390 : 0 : write(std_out,'(a,i5,a,i5,a,i5)') ' the interactions ',ninter,&
1391 : 0 : & ' and ',tmpinter,' have both the same symmetry isym=',Shell4at%neighbours(iatom,ishell)%sym_in_shell( ninter)
1392 : 0 : ABI_ERROR('Some interactions are equals due to the symmetry')
1393 : : end if
1394 : : end do
1395 : : end do
1396 : : end do
1397 : : end do
1398 : : ! Check that each equivalent shell has the same set of interactions
1399 [ + + ]: 58 : do ishell=1,nshell
1400 [ + + ]: 8538 : do iatom=1,natom
1401 [ - + ]: 8480 : if (Shell4at%neighbours(mod(iatom-1,natom_unitcell)+1,ishell)%n_interactions.ne.&
1402 : 49 : & Shell4at%neighbours( iatom,ishell)%n_interactions) then
1403 : 0 : ABI_ERROR('The interaction number index is not equal to the interaction number max computed previously (2)')
1404 : : end if
1405 : : !DEBUG iatref=Shell4at%iatref(ishell)
1406 : : !DEBUG jatref=Shell4at%jatref(ishell)
1407 : : !DEBUG katref=Shell4at%katref(ishell)
1408 : : !DEBUG if (Shell4at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
1409 : : !DEBUG do ninter=1,Shell4at%neighbours(iatom,ishell)%n_interactions
1410 : : !DEBUG jatom=Shell4at%neighbours(iatom,ishell)%atomj_in_shell(ninter)
1411 : : !DEBUG katom=Shell4at%neighbours(iatom,ishell)%atomk_in_shell(ninter)
1412 : : !DEBUG isym =Shell4at%neighbours(iatom,ishell)%sym_in_shell(ninter)
1413 : : !DEBUG itrans=Shell4at%neighbours(iatom,ishell)%transpose_in_shell(ninter)
1414 : : !DEBUG vectj(:)=zero ; vectk(:)=zero ; vect1(:)=zero ; vect2(:)=zero
1415 : : !DEBUG do ii=1,3
1416 : : !DEBUG do jj=1,3
1417 : : !DEBUG vectj(ii)=vectj(ii)+Sym%S_ref(ii,jj,isym,1)*SC%dij(iatref,jatref,jj+1)
1418 : : !DEBUG vectk(ii)=vectk(ii)+Sym%S_ref(ii,jj,isym,1)*SC%dij(iatref,katref,jj+1)
1419 : : !DEBUG end do
1420 : : !DEBUG end do
1421 : : !DEBUG if (itrans==1) then ; vect1(:)= vectj(:) ; vect2(:)= vectk(:) ; endif !\Phi3_ijk
1422 : : !DEBUG if (itrans==2) then ; vect1(:)= vectk(:) ; vect2(:)= vectj(:) ; endif !\Phi3_ikj
1423 : : !DEBUG if (itrans==3) then ; vect1(:)=-vectj(:) ; vect2(:)= vectk(:)-vectj(:) ; endif !\Phi3_jik
1424 : : !DEBUG if (itrans==4) then ; vect1(:)= vectk(:)-vectj(:) ; vect2(:)=-vectj(:) ; endif !\Phi3_jki
1425 : : !DEBUG if (itrans==5) then ; vect1(:)=-vectk(:) ; vect2(:)= vectj(:)-vectk(:) ; endif !\Phi3_kij
1426 : : !DEBUG if (itrans==6) then ; vect1(:)= vectj(:)-vectk(:) ; vect2(:)=-vectk(:) ; endif !\Phi3_kji
1427 : : !DEBUG do ii=1,3
1428 : : !DEBUG if ((abs(SC%dij(iatom,jatom,ii+1)-vect1(ii)).gt.tol8).or.&
1429 : : !DEBUG& (abs(SC%dij(iatom,katom,ii+1)-vect2(ii)).gt.tol8)) then
1430 : : !DEBUG write(std_out,'(a,4(x,i5))') 'For ishell, iatom, jatom, katom =',ishell,iatom,jatom,katom
1431 : : !DEBUG write(std_out,'(a,5(x,i5))') ' with isym, itrans, iatref, jatref, katref = ',isym,itrans,iatref,jatref,katref
1432 : : !DEBUG ABI_ERROR('We do not recover the quadruplet with the symmetry found')
1433 : : !DEBUG end if
1434 : : !DEBUG end do !ii
1435 : : !DEBUG end do !ninter
1436 : : end do !natom
1437 : : end do !nshell
1438 [ - + ]: 9 : ABI_FREE(interactions)
1439 : :
1440 : : ! Find the number of coefficients of the (3x3x3x3) Phi4 for a given shell
1441 [ - + + - : 325534 : ABI_CALLOC(Shell4at%proj, (norder,norder,nshell))
+ - + - -
+ - + - +
+ + + + +
+ ]
1442 [ - + + - : 67 : ABI_CALLOC(Shell4at%ncoeff ,(nshell))
- + - + -
+ + + ]
1443 [ - + + - : 67 : ABI_CALLOC(Shell4at%ncoeff_prev,(nshell))
- + - + -
+ + + ]
1444 : 9 : write(ab_out,*) 'Number of shells=',nshell
1445 : 9 : write(ab_out,*) '============================================================================'
1446 [ + - ]: 9 : if (MPIdata%iam_master) open(unit=16,file=trim(Invar%output_prefix)//'_nbcoeff-phi4.dat')
1447 : 9 : ncoeff_prev=0
1448 [ + + ]: 58 : do ishell=1,nshell
1449 : 49 : ncoeff=0
1450 : 49 : iatref=Shell4at%iatref(ishell)
1451 : 49 : jatref=Shell4at%jatref(ishell)
1452 : 49 : katref=Shell4at%katref(ishell)
1453 : 49 : latref=Shell4at%latref(ishell)
1454 : 49 : write(ab_out,*) 'Shell number:',ishell
1455 : 49 : write(ab_out,'(a,i5,a,i5,a,f16.10)') ' Between atom',iatref,' and ',jatref,' the distance is=',SC%dij(iatref,jatref,1)
1456 : 49 : write(ab_out,'(a,i5,a,i5,a,f16.10)') ' Between atom',jatref,' and ',katref,' the distance is=',SC%dij(jatref,katref,1)
1457 : 49 : write(ab_out,'(a,i5,a,i5,a,f16.10)') ' Between atom',katref,' and ',latref,' the distance is=',SC%dij(katref,latref,1)
1458 : 49 : write(ab_out,'(a,i5,a,i5,a,f16.10)') ' Between atom',latref,' and ',iatref,' the distance is=',SC%dij(latref,iatref,1)
1459 : 49 : call tdep_calc_nbcoeff(SC,MPIdata,iatref,ishell,jatref,katref,latref,ncoeff,norder,nshell,order,Shell4at%proj,Sym)
1460 : 49 : Shell4at%ncoeff (ishell)=ncoeff
1461 : 49 : Shell4at%ncoeff_prev(ishell)=ncoeff_prev
1462 : 49 : ncoeff_prev=ncoeff_prev+ncoeff
1463 : 49 : write(ab_out,*)' Number of independant coefficients in this shell=',ncoeff
1464 : 49 : write(ab_out,*)' Number of interactions in this shell=',Shell4at%neighbours(iatref,ishell)%n_interactions
1465 : : !FB write(ab_out,*)' The ratio is=',dfloat(Shell4at%neighbours(iatref,ishell)%n_interactions)/dfloat(ncoeff)
1466 : 58 : write(ab_out,*) '============================================================================'
1467 : : end do
1468 : 9 : write(ab_out,*)' >>>>>> Total number of coefficients at the fourth order=',ncoeff_prev
1469 [ + - ]: 9 : if (MPIdata%iam_master) close(16)
1470 : 9 : Shell4at%ntotcoeff=ncoeff_prev
1471 : : !BeginFB
1472 : : !FB open(unit=91,file='Shell4at.dat')
1473 : : !FB write(91,*) Shell4at%nshell
1474 : : !FB do ishell=1,Shell4at%nshell
1475 : : !FB write(91,*) Shell4at%ncoeff(ishell)
1476 : : !FB write(91,*) Shell4at%ncoeff_prev(ishell)
1477 : : !FB write(91,*) Shell4at%iatref(ishell)
1478 : : !FB write(91,*) Shell4at%jatref(ishell)
1479 : : !FB write(91,*) Shell4at%katref(ishell)
1480 : : !FB do iatom=1,Invar%natom
1481 : : !FB write(91,*) Shell4at%neighbours(iatom,ishell)%n_interactions
1482 : : !FB do ii=1,Shell4at%neighbours(iatom,ishell)%n_interactions
1483 : : !FB write(91,*) Shell4at%neighbours(iatom,ishell)%sym_in_shell(ii)
1484 : : !FB write(91,*) Shell4at%neighbours(iatom,ishell)%transpose_in_shell(ii)
1485 : : !FB write(91,*) Shell4at%neighbours(iatom,ishell)%atomj_in_shell(ii)
1486 : : !FB write(91,*) Shell4at%neighbours(iatom,ishell)%atomk_in_shell(ii)
1487 : : !FB end do
1488 : : !FB end do
1489 : : !FB end do
1490 : : !FB close(91)
1491 : : !EndFB
1492 : :
1493 : 9 : end subroutine tdep_init_shell4at
1494 : : !!***
1495 : :
1496 : : !====================================================================================================
1497 : :
1498 : : !!****f* ABINIT/m_tdep_shell/tdep_destroy_shell
1499 : : !! NAME
1500 : : !! tdep_destroy_shell
1501 : : !!
1502 : : !! FUNCTION
1503 : : !! Free the memory allocated in the Shell object.
1504 : : !!
1505 : : !! SIDE EFFECTS
1506 : : !! Shell = Shell object. All its allocated arrays are deallocated.
1507 : : !!
1508 : : !! SOURCE
1509 : :
1510 : 117 : subroutine tdep_destroy_shell(Shell)
1511 : :
1512 : : class(Shell_type),intent(inout) :: Shell
1513 : :
1514 : : integer :: iatom,ishell,natref
1515 : :
1516 [ - + ]: 117 : ABI_FREE(Shell%ncoeff)
1517 [ - + ]: 117 : ABI_FREE(Shell%ncoeff_prev)
1518 [ - + ]: 117 : ABI_FREE(Shell%iatref)
1519 [ + + ]: 117 : if (Shell%order.gt.1) then
1520 [ - + ]: 71 : ABI_FREE(Shell%jatref)
1521 [ - + ]: 71 : ABI_FREE(Shell%ishell_self)
1522 : : end if
1523 [ + + ]: 117 : if (Shell%order.gt.2) then
1524 [ - + ]: 25 : ABI_FREE(Shell%katref)
1525 : : end if
1526 [ + + ]: 117 : if (Shell%order.gt.3) then
1527 [ - + ]: 9 : ABI_FREE(Shell%latref)
1528 : : end if
1529 [ + + ]: 117 : if (Shell%order.eq.1) then
1530 : : natref=1
1531 : : else
1532 : 71 : natref=Shell%natom
1533 : : end if
1534 [ + + ]: 12069 : do iatom=1,natref
1535 [ + + ]: 113384 : do ishell=1,Shell%nshell
1536 [ + + ]: 113267 : if (Shell%neighbours(iatom,ishell)%n_interactions.ne.0) then
1537 [ - + ]: 87426 : ABI_FREE(Shell%neighbours(iatom,ishell)%atomj_in_shell)
1538 [ - + ]: 87426 : ABI_FREE(Shell%neighbours(iatom,ishell)%sym_in_shell)
1539 [ + - ]: 87426 : if (Shell%order.gt.1) then
1540 [ - + ]: 87426 : ABI_FREE(Shell%neighbours(iatom,ishell)%transpose_in_shell)
1541 : : end if
1542 [ + + ]: 87426 : if (Shell%order.gt.2) then
1543 [ - + ]: 16016 : ABI_FREE(Shell%neighbours(iatom,ishell)%atomk_in_shell)
1544 : : end if
1545 [ + + ]: 87426 : if (Shell%order.gt.3) then
1546 [ - + ]: 7456 : ABI_FREE(Shell%neighbours(iatom,ishell)%atoml_in_shell)
1547 : : end if
1548 : : end if
1549 : : end do
1550 : : end do
1551 [ + - + + : 101432 : ABI_FREE(Shell%neighbours)
- + - + -
+ - + - +
- + ]
1552 [ - + ]: 117 : ABI_FREE(Shell%proj)
1553 : :
1554 : 117 : end subroutine tdep_destroy_shell
1555 : : !!***
1556 : :
1557 : : !====================================================================================================
1558 : :
1559 : : !!****f* ABINIT/m_tdep_shell/tdep_free_shell
1560 : : !! NAME
1561 : : !! tdep_free_shell
1562 : : !!
1563 : : !! FUNCTION
1564 : : !! Free the memory allocated in the Shell object (see tdep_destroy_shell).
1565 : : !! This wrapper is the target of the type-bound procedure Shell%free,
1566 : : !! as a workaround for nvfortran 24.9.
1567 : : !!
1568 : : !! SIDE EFFECTS
1569 : : !! Shell = Shell object. All its allocated arrays are deallocated.
1570 : : !!
1571 : : !! SOURCE
1572 : :
1573 : 117 : subroutine tdep_free_shell(Shell)
1574 : :
1575 : : class(Shell_type),intent(inout) :: Shell
1576 : :
1577 : : ! With nvfortran 24.9 (-fast), the binding free => tdep_destroy_shell is null and calling
1578 : : ! Shell%free() segfaults (jump to address 0), whereas a binding to a wrapper works.
1579 : 117 : call tdep_destroy_shell(Shell)
1580 : :
1581 : 117 : end subroutine tdep_free_shell
1582 : : !!***
1583 : :
1584 : : !====================================================================================================
1585 : :
1586 : : !!****f* ABINIT/m_tdep_shell/tdep_calc_nbcoeff
1587 : : !! NAME
1588 : : !! tdep_calc_nbcoeff
1589 : : !!
1590 : : !! FUNCTION
1591 : : !! Compute the number of independent coefficients of the IFC (at a given
1592 : : !! order) for a given shell, and the projector onto these coefficients.
1593 : : !! For each symmetry operation (combined with each permutation of the atoms)
1594 : : !! that leaves the reference interaction of the shell invariant, the
1595 : : !! corresponding linear constraints on the norder IFC coefficients are
1596 : : !! built. At the third and fourth orders, when at least two atoms of the
1597 : : !! reference interaction coincide, the constraints coming from the
1598 : : !! permutation symmetry of the IFC are added. The set of constraint vectors is then
1599 : : !! orthonormalized (Gram-Schmidt), and the (norder-ncount) vectors
1600 : : !! orthogonal to the ncount independent constraints span the space of the
1601 : : !! independent coefficients.
1602 : : !! The details are written (by the master process) in the file opened on
1603 : : !! unit 16 by the caller (_nbcoeff-phi*.dat).
1604 : : !!
1605 : : !! INPUTS
1606 : : !! SC = Supercell object, used for the interatomic distances (SC%dij).
1607 : : !! MPIdata = Info on MPI parallelism.
1608 : : !! iatcell = Index of the reference atom of the shell, in the unitcell.
1609 : : !! ishell = Index of the shell.
1610 : : !! jatom = Index of the second atom of the reference interaction (order >= 2).
1611 : : !! katom = Index of the third atom of the reference interaction (order >= 3).
1612 : : !! latom = Index of the fourth atom of the reference interaction (order == 4).
1613 : : !! norder = Dimension of the IFC at this order (3**order).
1614 : : !! nshell = Number of shells.
1615 : : !! order = Order of the interactions (1, 2, 3 or 4).
1616 : : !! Sym = Symmetries object describing all the symmetry operations of the crystal.
1617 : : !!
1618 : : !! OUTPUT
1619 : : !! proj(norder,norder,nshell) = proj(:,1:ncoeff,ishell) contains the
1620 : : !! projectors onto the independent
1621 : : !! coefficients of the shell ishell.
1622 : : !!
1623 : : !! SIDE EFFECTS
1624 : : !! ncoeff = Number of independent coefficients of the shell.
1625 : : !! Unchanged if the routine returns early (order 1 with iatcell=1,
1626 : : !! or order 2 with jatom=iatcell).
1627 : : !!
1628 : : !! SOURCE
1629 : :
1630 : 693 : subroutine tdep_calc_nbcoeff(SC,MPIdata,iatcell,ishell,jatom,katom,latom,&
1631 : 693 : & ncoeff,norder,nshell,order,proj,Sym)
1632 : :
1633 : : integer,intent(in) :: iatcell,ishell,jatom,katom,latom,nshell,order,norder
1634 : : integer,intent(inout) :: ncoeff
1635 : : type(tdep_Supercell_type),intent(in) :: SC
1636 : : type(MPI_enreg_type), intent(in) :: MPIdata
1637 : : type(Symmetries_type),intent(in) :: Sym
1638 : : double precision,intent(out) :: proj(norder,norder,nshell)
1639 : :
1640 : : integer :: ii,jj,kk,ll,isym,LWORK,INFO,const_tot,itemp,nconst_perm,nconst_loc
1641 : : integer :: ncount,icoeff,jatcell,katcell,latcell,mu,nu,xi,zeta
1642 : : integer :: inv,watom,xatom,yatom,zatom,isyminv,nsyminv,facorder,iseed
1643 : 693 : integer, allocatable :: iconst(:)
1644 : : double precision :: prod_scal,drandom
1645 : : double precision :: eigvec(3,3)
1646 : : double precision :: vect_trial(3),vect_trial1(3),vect_trial2(3),vect_trial3(3)
1647 : : double precision :: vect_trial4(3),vect_trial5(3),vect_trial6(3)
1648 : : double precision :: WR(3),WI(3),VL(3,3),VR(3,3)
1649 : 693 : double precision, allocatable :: WORK(:)
1650 : : double complex :: eigenvectors(3,3),eigenvalues(3)
1651 : : double complex :: pp(3,3),ppp(3,3,3),pppp(3,3,3,3),lambda
1652 : 693 : double complex, allocatable :: tab_vec(:,:),temp(:,:),alphaij(:,:,:),constraints(:,:,:)
1653 : : logical :: ok
1654 : 693 : logical, allocatable :: unchanged(:)
1655 : : character(len=500) :: message
1656 : :
1657 [ + + + + ]: 693 : if (iatcell==1.and.order==1) return
1658 [ + + + + ]: 647 : if (jatom==iatcell.and.order==2) return
1659 : : !FB if (katom==iatcell.and.jatom==iatcell.and.order==3) return
1660 : :
1661 [ + + ]: 582 : if (order==1) then
1662 : : facorder=1
1663 [ + + ]: 563 : else if (order==2) then
1664 : : facorder=2
1665 [ + + ]: 106 : else if (order==3) then
1666 : : facorder=6
1667 [ + - ]: 49 : else if (order==4) then
1668 : 49 : facorder=24
1669 : : end if
1670 : :
1671 : : ! If we want to remove the constraints coming from the symetries
1672 : : !FB if (order==3) then
1673 : : !FB do ii=1,norder
1674 : : !FB proj(ii,ii,ishell)=1.d0
1675 : : !FB end do
1676 : : !FB ncoeff=norder
1677 : : !FB return
1678 : : !FB end if
1679 : :
1680 : 582 : nconst_loc=0
1681 : 582 : const_tot=0
1682 : 582 : nsyminv=Sym%nsym*facorder
1683 [ - + + - : 332988576 : ABI_MALLOC(alphaij,(nsyminv,norder,norder)); alphaij(:,:,:)=czero
- + + - -
+ + - + -
+ - - + -
+ + + + +
+ + ]
1684 [ - + + - : 88752 : ABI_MALLOC(iconst,(nsyminv)) ; iconst(:)=0
+ - + - -
+ - + +
+ ]
1685 [ - + + - : 88170 : ABI_MALLOC(unchanged,(nsyminv)) ; unchanged(:)=.false.
- + - + +
+ ]
1686 : :
1687 : : ! ================================================================================================
1688 : : ! ================ Big loop over symmetries and invariance (nsym*facorder) =======================
1689 : : ! ================================================================================================
1690 [ + - ]: 582 : if (MPIdata%iam_master) write(16,'(a)') ' '
1691 [ + - ]: 582 : if (MPIdata%iam_master) write(16,'(a,i4)') 'For shell number=',ishell
1692 [ + + ]: 87306 : do isyminv=1,nsyminv
1693 : 86730 : isym=(isyminv-1)/facorder+1
1694 : 86730 : inv=isyminv-(isym-1)*facorder
1695 [ + + ]: 86730 : if (isym==1) cycle
1696 : :
1697 : : ! For the 1st order: Search if the atom is let invariant
1698 [ + + ]: 84279 : if (order==1) then
1699 [ + + ]: 223 : if (Sym%indsym(4,isym,iatcell)==iatcell) then
1700 [ + - ]: 107 : if (MPIdata%iam_master) then
1701 : 107 : write(16,'(a,1x,i3)')'===========The atom is kept invariant for isym=',isym
1702 : : end if
1703 : : else
1704 : : cycle
1705 : : end if
1706 : : end if
1707 : :
1708 : : ! For the 2nd order: Search if the bond is kept invariant or reversed
1709 [ + + ]: 84163 : if (order==2) then
1710 : 22934 : vect_trial(:)=zero
1711 [ + + ]: 22934 : if (inv==1) then ; watom=iatcell ; xatom=jatom ; endif !\Phi_ij
1712 [ + + ]: 22934 : if (inv==2) then ; watom=jatom ; xatom=iatcell ; endif !\Phi_ji
1713 [ + + ]: 91736 : do ii=1,3
1714 [ + + ]: 298142 : do jj=1,3
1715 : 275208 : vect_trial(ii)=vect_trial(ii)+Sym%S_ref(ii,jj,isym,1)*SC%dij(watom,xatom,jj+1)
1716 : : end do
1717 : : end do
1718 : 22934 : jatcell=mod(jatom-1,SC%natom_unitcell)+1
1719 : : if ((sum(abs(vect_trial(:)-SC%dij(iatcell,jatom,2:4))).lt.tol8).and.&
1720 [ + + + + : 91736 : & (Sym%indsym(4,isym,watom)==iatcell).and.&
+ + + - ]
1721 : : & (Sym%indsym(4,isym,xatom)==jatcell)) then
1722 [ + - ]: 1738 : if (MPIdata%iam_master) then
1723 [ + + ]: 1738 : if (inv==1) write(16,'(a,1x,i3)')'===========The bond is kept invariant for isym=',isym
1724 [ + + ]: 1738 : if (inv==2) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j) --> (j,i) for isym=',isym
1725 : : end if
1726 : : else
1727 : : cycle
1728 : : end if
1729 : : end if
1730 : :
1731 : : ! For the 3rd order : 6 permutations at all
1732 [ + + ]: 62967 : if (order==3) then
1733 : 13434 : vect_trial1(:)=zero
1734 : 13434 : vect_trial2(:)=zero
1735 : 13434 : vect_trial3(:)=zero
1736 [ + + ]: 13434 : if (inv==1) then ; watom=iatcell ; xatom=jatom ; yatom=katom ; endif !\Phi3_ijk
1737 [ + + ]: 13434 : if (inv==2) then ; watom=iatcell ; xatom=katom ; yatom=jatom ; endif !\Phi3_ikj
1738 [ + + ]: 13434 : if (inv==3) then ; watom=jatom ; xatom=iatcell ; yatom=katom ; endif !\Phi3_jik
1739 [ + + ]: 13434 : if (inv==4) then ; watom=jatom ; xatom=katom ; yatom=iatcell ; endif !\Phi3_jki
1740 [ + + ]: 13434 : if (inv==5) then ; watom=katom ; xatom=iatcell ; yatom=jatom ; endif !\Phi3_kij
1741 [ + + ]: 13434 : if (inv==6) then ; watom=katom ; xatom=jatom ; yatom=iatcell ; endif !\Phi3_kji
1742 [ + + ]: 53736 : do ii=1,3
1743 [ + + ]: 174642 : do jj=1,3
1744 : 120906 : vect_trial1(ii)=vect_trial1(ii)+Sym%S_ref(ii,jj,isym,1)*SC%dij(watom,xatom,jj+1)
1745 : 120906 : vect_trial2(ii)=vect_trial2(ii)+Sym%S_ref(ii,jj,isym,1)*SC%dij(xatom,yatom,jj+1)
1746 : 161208 : vect_trial3(ii)=vect_trial3(ii)+Sym%S_ref(ii,jj,isym,1)*SC%dij(yatom,watom,jj+1)
1747 : : end do
1748 : : end do
1749 : 13434 : jatcell=mod(jatom-1,SC%natom_unitcell)+1
1750 : 13434 : katcell=mod(katom-1,SC%natom_unitcell)+1
1751 : : if ((sum(abs(vect_trial1(:)-SC%dij(iatcell,jatom ,2:4))).lt.tol8).and.&
1752 : : & (sum(abs(vect_trial2(:)-SC%dij(jatom ,katom ,2:4))).lt.tol8).and.&
1753 : : & (sum(abs(vect_trial3(:)-SC%dij(katom ,iatcell,2:4))).lt.tol8).and.&
1754 : : & (Sym%indsym(4,isym,watom)==iatcell).and.&
1755 [ + + + + : 134340 : & (Sym%indsym(4,isym,xatom)==jatcell).and.&
+ + + + +
+ + - + +
+ - + - ]
1756 : : & (Sym%indsym(4,isym,yatom)==katcell)) then
1757 [ + - ]: 4007 : if (MPIdata%iam_master) then
1758 [ + + ]: 4007 : if (inv==1) write(16,'(a,1x,i3)')'===========The bond is kept invariant for isym=',isym
1759 [ + + ]: 4007 : if (inv==2) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k) --> (i,k,j) for isym=',isym
1760 [ + + ]: 4007 : if (inv==3) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k) --> (j,i,k) for isym=',isym
1761 [ + + ]: 4007 : if (inv==4) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k) --> (j,k,i) for isym=',isym
1762 [ + + ]: 4007 : if (inv==5) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k) --> (k,i,j) for isym=',isym
1763 [ + + ]: 4007 : if (inv==6) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k) --> (k,j,i) for isym=',isym
1764 : : end if
1765 : : else
1766 : : cycle
1767 : : end if
1768 : : end if
1769 : :
1770 : : ! For the 4th order : 24 permutations at all
1771 [ + + ]: 53540 : if (order==4) then
1772 : 47688 : vect_trial1(:)=zero
1773 : 47688 : vect_trial2(:)=zero
1774 : 47688 : vect_trial3(:)=zero
1775 : 47688 : vect_trial4(:)=zero
1776 : 47688 : vect_trial5(:)=zero
1777 : 47688 : vect_trial6(:)=zero
1778 [ + + ]: 47688 : if (inv==1) then ; watom=iatcell ; xatom=jatom ; yatom=katom ; zatom=latom ; endif !\Phi4_ijkl
1779 [ + + ]: 47688 : if (inv==2) then ; watom=iatcell ; xatom=katom ; yatom=jatom ; zatom=latom ; endif !\Phi4_ikjl
1780 [ + + ]: 47688 : if (inv==3) then ; watom=jatom ; xatom=iatcell ; yatom=katom ; zatom=latom ; endif !\Phi4_jikl
1781 [ + + ]: 47688 : if (inv==4) then ; watom=jatom ; xatom=katom ; yatom=iatcell ; zatom=latom ; endif !\Phi4_jkil
1782 [ + + ]: 47688 : if (inv==5) then ; watom=katom ; xatom=iatcell ; yatom=jatom ; zatom=latom ; endif !\Phi4_kijl
1783 [ + + ]: 47688 : if (inv==6) then ; watom=katom ; xatom=jatom ; yatom=iatcell ; zatom=latom ; endif !\Phi4_kjil
1784 : :
1785 [ + + ]: 47688 : if (inv==7 ) then ; watom=iatcell ; xatom=jatom ; yatom=latom ; zatom=katom ; endif !\Phi4_ijlk
1786 [ + + ]: 47688 : if (inv==8 ) then ; watom=iatcell ; xatom=katom ; yatom=latom ; zatom=jatom ; endif !\Phi4_iklj
1787 [ + + ]: 47688 : if (inv==9 ) then ; watom=jatom ; xatom=iatcell ; yatom=latom ; zatom=katom ; endif !\Phi4_jilk
1788 [ + + ]: 47688 : if (inv==10) then ; watom=jatom ; xatom=katom ; yatom=latom ; zatom=iatcell ; endif !\Phi4_jkli
1789 [ + + ]: 47688 : if (inv==11) then ; watom=katom ; xatom=iatcell ; yatom=latom ; zatom=jatom ; endif !\Phi4_kilj
1790 [ + + ]: 47688 : if (inv==12) then ; watom=katom ; xatom=jatom ; yatom=latom ; zatom=iatcell ; endif !\Phi4_kjli
1791 : :
1792 [ + + ]: 47688 : if (inv==13) then ; watom=iatcell ; xatom=latom ; yatom=jatom ; zatom=katom ; endif !\Phi4_iljk
1793 [ + + ]: 47688 : if (inv==14) then ; watom=iatcell ; xatom=latom ; yatom=katom ; zatom=jatom ; endif !\Phi4_ilkj
1794 [ + + ]: 47688 : if (inv==15) then ; watom=jatom ; xatom=latom ; yatom=iatcell ; zatom=katom ; endif !\Phi4_jlik
1795 [ + + ]: 47688 : if (inv==16) then ; watom=jatom ; xatom=latom ; yatom=katom ; zatom=iatcell ; endif !\Phi4_jlki
1796 [ + + ]: 47688 : if (inv==17) then ; watom=katom ; xatom=latom ; yatom=iatcell ; zatom=jatom ; endif !\Phi4_klij
1797 [ + + ]: 47688 : if (inv==18) then ; watom=katom ; xatom=latom ; yatom=jatom ; zatom=iatcell ; endif !\Phi4_klji
1798 : :
1799 [ + + ]: 47688 : if (inv==19) then ; watom=latom ; xatom=iatcell ; yatom=jatom ; zatom=katom ; endif !\Phi4_lijk
1800 [ + + ]: 47688 : if (inv==20) then ; watom=latom ; xatom=iatcell ; yatom=katom ; zatom=jatom ; endif !\Phi4_likj
1801 [ + + ]: 47688 : if (inv==21) then ; watom=latom ; xatom=jatom ; yatom=iatcell ; zatom=katom ; endif !\Phi4_ljik
1802 [ + + ]: 47688 : if (inv==22) then ; watom=latom ; xatom=jatom ; yatom=katom ; zatom=iatcell ; endif !\Phi4_ljki
1803 [ + + ]: 47688 : if (inv==23) then ; watom=latom ; xatom=katom ; yatom=iatcell ; zatom=jatom ; endif !\Phi4_lkij
1804 [ + + ]: 47688 : if (inv==24) then ; watom=latom ; xatom=katom ; yatom=jatom ; zatom=iatcell ; endif !\Phi4_lkji
1805 : :
1806 [ + + ]: 190752 : do ii=1,3
1807 [ + + ]: 619944 : do jj=1,3
1808 : 429192 : vect_trial1(ii)=vect_trial1(ii)+Sym%S_ref(ii,jj,isym,1)*SC%dij(watom,xatom,jj+1)
1809 : 429192 : vect_trial2(ii)=vect_trial2(ii)+Sym%S_ref(ii,jj,isym,1)*SC%dij(watom,yatom,jj+1)
1810 : 429192 : vect_trial3(ii)=vect_trial3(ii)+Sym%S_ref(ii,jj,isym,1)*SC%dij(watom,zatom,jj+1)
1811 : 429192 : vect_trial4(ii)=vect_trial4(ii)+Sym%S_ref(ii,jj,isym,1)*SC%dij(xatom,yatom,jj+1)
1812 : 429192 : vect_trial5(ii)=vect_trial5(ii)+Sym%S_ref(ii,jj,isym,1)*SC%dij(xatom,zatom,jj+1)
1813 : 572256 : vect_trial6(ii)=vect_trial6(ii)+Sym%S_ref(ii,jj,isym,1)*SC%dij(yatom,zatom,jj+1)
1814 : : end do
1815 : : end do
1816 : 47688 : jatcell=mod(jatom-1,SC%natom_unitcell)+1
1817 : 47688 : katcell=mod(katom-1,SC%natom_unitcell)+1
1818 : 47688 : latcell=mod(latom-1,SC%natom_unitcell)+1
1819 : : if ((sum(abs(vect_trial1(:)-SC%dij(iatcell,jatom,2:4))).lt.tol8).and.&
1820 : : & (sum(abs(vect_trial2(:)-SC%dij(iatcell,katom,2:4))).lt.tol8).and.&
1821 : : & (sum(abs(vect_trial3(:)-SC%dij(iatcell,latom,2:4))).lt.tol8).and.&
1822 : : & (sum(abs(vect_trial4(:)-SC%dij(jatom ,katom,2:4))).lt.tol8).and.&
1823 : : & (sum(abs(vect_trial5(:)-SC%dij(jatom ,latom,2:4))).lt.tol8).and.&
1824 : : & (sum(abs(vect_trial6(:)-SC%dij(katom ,latom,2:4))).lt.tol8).and.&
1825 : : & (Sym%indsym(4,isym,watom)==iatcell).and.&
1826 : : & (Sym%indsym(4,isym,xatom)==jatcell).and.&
1827 [ + + + + : 906072 : & (Sym%indsym(4,isym,yatom)==katcell).and.&
+ + + + +
+ + + + +
+ + + + +
- + - + -
+ + + - +
- + - ]
1828 : : & (Sym%indsym(4,isym,zatom)==latcell)) then
1829 [ + - ]: 9580 : if (MPIdata%iam_master) then
1830 [ + + ]: 9580 : if (inv==1 ) write(16,'(a,1x,i3)')'===========The bond is kept invariant for isym=',isym
1831 [ + + ]: 9580 : if (inv==2 ) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k,l) --> (i,k,j,l) for isym=',isym !\Phi4_ikjl
1832 [ + + ]: 9580 : if (inv==3 ) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k,l) --> (j,i,k,l) for isym=',isym !\Phi4_jikl
1833 [ + + ]: 9580 : if (inv==4 ) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k,l) --> (j,k,i,l) for isym=',isym !\Phi4_jkil
1834 [ + + ]: 9580 : if (inv==5 ) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k,l) --> (k,i,j,l) for isym=',isym !\Phi4_kijl
1835 [ + + ]: 9580 : if (inv==6 ) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k,l) --> (k,j,i,l) for isym=',isym !\Phi4_kjil
1836 : :
1837 [ + + ]: 9580 : if (inv==7 ) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k,l) --> (i,j,l,k) for isym=',isym !\Phi4_ijlk
1838 [ + + ]: 9580 : if (inv==8 ) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k,l) --> (i,k,l,j) for isym=',isym !\Phi4_iklj
1839 [ + + ]: 9580 : if (inv==9 ) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k,l) --> (j,i,l,k) for isym=',isym !\Phi4_jilk
1840 [ + + ]: 9580 : if (inv==10) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k,l) --> (j,k,l,i) for isym=',isym !\Phi4_jkli
1841 [ + + ]: 9580 : if (inv==11) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k,l) --> (k,i,l,j) for isym=',isym !\Phi4_kilj
1842 [ + + ]: 9580 : if (inv==12) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k,l) --> (k,j,l,i) for isym=',isym !\Phi4_kjli
1843 : :
1844 [ + + ]: 9580 : if (inv==13) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k,l) --> (i,l,j,k) for isym=',isym !\Phi4_iljk
1845 [ + + ]: 9580 : if (inv==14) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k,l) --> (i,l,k,j) for isym=',isym !\Phi4_ilkj
1846 [ + + ]: 9580 : if (inv==15) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k,l) --> (j,l,i,k) for isym=',isym !\Phi4_jlik
1847 [ + + ]: 9580 : if (inv==16) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k,l) --> (j,l,k,i) for isym=',isym !\Phi4_jlki
1848 [ + + ]: 9580 : if (inv==17) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k,l) --> (k,l,i,j) for isym=',isym !\Phi4_klij
1849 [ + + ]: 9580 : if (inv==18) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k,l) --> (k,l,j,i) for isym=',isym !\Phi4_klji
1850 : :
1851 [ + + ]: 9580 : if (inv==19) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k,l) --> (l,i,j,k) for isym=',isym !\Phi4_lijk
1852 [ + + ]: 9580 : if (inv==20) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k,l) --> (l,i,k,j) for isym=',isym !\Phi4_likj
1853 [ + + ]: 9580 : if (inv==21) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k,l) --> (l,j,i,k) for isym=',isym !\Phi4_ljik
1854 [ + + ]: 9580 : if (inv==22) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k,l) --> (l,j,k,i) for isym=',isym !\Phi4_ljki
1855 [ + + ]: 9580 : if (inv==23) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k,l) --> (l,k,i,j) for isym=',isym !\Phi4_lkij
1856 [ + + ]: 9580 : if (inv==24) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k,l) --> (l,k,j,i) for isym=',isym !\Phi4_lkji
1857 : : end if
1858 : : else
1859 : : cycle
1860 : : end if
1861 : : end if
1862 : :
1863 : : ! Write the S_ref matrix
1864 : : !FB write(16,'(3(f16.12,1x))') Sym%S_ref(1,1,isym,1),Sym%S_ref(1,2,isym,1),Sym%S_ref(1,3,isym,1)
1865 : : !FB write(16,'(3(f16.12,1x))') Sym%S_ref(2,1,isym,1),Sym%S_ref(2,2,isym,1),Sym%S_ref(2,3,isym,1)
1866 : : !FB write(16,'(3(f16.12,1x))') Sym%S_ref(3,1,isym,1),Sym%S_ref(3,2,isym,1),Sym%S_ref(3,3,isym,1)
1867 : :
1868 : : ! Diagonalize the S_ref matrix
1869 [ + + ]: 61728 : do ii=1,3
1870 [ + + ]: 200616 : do jj=1,3
1871 : 185184 : eigvec(ii,jj)=Sym%S_ref(jj,ii,isym,1)
1872 : : end do
1873 : : end do
1874 : 15432 : LWORK=4*3
1875 [ - + + + ]: 200616 : ABI_MALLOC(WORK,(LWORK)); WORK(:)=zero
1876 : : ! This one is real and could be non-symmetric
1877 : 15432 : call dgeev( 'N', 'V', 3, eigvec, 3, WR, WI, VL, 3, VR, 3, WORK, LWORK, INFO)
1878 [ - + ]: 15432 : ABI_FREE(WORK)
1879 : :
1880 : : ! Build the real and imaginary parts of the eigenvectors and eigenvalues
1881 : 15432 : jj=0
1882 [ + + ]: 61728 : do ii=1,3
1883 : 46296 : eigenvalues(ii)=dcmplx(WR(ii),-WI(ii))
1884 [ + + + + ]: 61728 : if (WI(ii).ne.zero.and.jj==0) then
1885 [ + + ]: 31936 : do kk=1,3
1886 : 31936 : eigenvectors(kk,ii)=dcmplx(VR(kk,ii),VR(kk,ii+1))
1887 : : end do
1888 : : jj=jj+1
1889 [ + + - + ]: 38312 : else if (WI(ii).ne.zero.and.jj==1) then
1890 [ + + ]: 31936 : do kk=1,3
1891 : 31936 : eigenvectors(kk,ii)=dcmplx(VR(kk,ii-1),-VR(kk,ii))
1892 : : end do
1893 : : jj=jj+1
1894 : : else
1895 [ + + ]: 121312 : do kk=1,3
1896 : 121312 : eigenvectors(kk,ii)=dcmplx(VR(kk,ii),zero)
1897 : : end do
1898 : : end if
1899 : : end do
1900 : :
1901 : : ! Write the eigenvalues and eigenvectors
1902 : 15432 : ok=.true.
1903 [ + + ]: 61728 : do ii=1,3
1904 [ + + + + : 61728 : if ((aimag(eigenvalues(1)).ne.0).or.(aimag(eigenvalues(2)).ne.0).or.(aimag(eigenvalues(3)).ne.0)) then
- + ]
1905 : 23952 : ok=.false.
1906 : : end if
1907 : : end do
1908 [ + + + - ]: 15432 : if (.not.ok.and.MPIdata%iam_master) write(16,'(a)') ' WARNING: THERE IS COMPLEX EIGENVALUES'
1909 : :
1910 : : ! If the transformation matrix keeps the bond invariant:
1911 : : ! Phi_{\alpha\beta}=\sum_{\mu\nu} S_{\alpha\mu}.S_{\beta\nu}.Phi_{\mu\nu}
1912 : : ! If lambda and p are the eigenvectors and eigenvalues of the S matrix, then:
1913 : : ! \sum_{\alpha\beta} p_{\alpha}^l.p_{\beta}^k Phi_{\alpha\beta}
1914 : : ! = \sum_{\mu\nu,\alpha\beta} p_{\alpha}^l.p_{\beta}^k.S_{\alpha\mu}.S_{\beta\nu}.Phi_{\mu\nu}
1915 : : ! = lambda^{*l}.lambda^{*k} \sum_{\mu\nu} p_{\mu}^l.p_{\nu}^k.Phi_{\mu\nu}
1916 : : ! So, if lambda^{*l}.lambda^{*k} = -1, we must have:
1917 : : ! \sum_{\alpha\beta} p_{\alpha}^l.p_{\beta}^k.Phi_{\alpha\beta}= 0
1918 : : !
1919 : : ! In the case of the reversed bond, one obtains the following constraint:
1920 : : ! \sum_{\alpha\beta} (lambda^{*l}.lambda^{*k}.p_{\alpha}^l.p_{\beta}^k-p_{\beta}^l.p_{\alpha}^k).Phi_{\alpha\beta}= 0
1921 : : ! which applies whether lambda^{*l}.lambda^{*k} = \pm 1
1922 : : !
1923 : : ! We obtain n vectors with norder coefficients (defined in the R^norder space).
1924 : : ! The space of the independent solutions are in the R^(norder-n) space, orthogonal
1925 : : ! to the space spanned by the starting n vectors.
1926 [ + + ]: 15432 : if (order==1) then
1927 [ + + ]: 428 : do ii=1,3
1928 : 321 : lambda=eigenvalues(ii)
1929 [ + + - + ]: 321 : if ((abs(real(lambda)-1.d0).lt.tol6).and.(abs(aimag(lambda)).lt.tol6)) cycle
1930 : 211 : unchanged(isyminv)=.true.
1931 : 211 : iconst(isyminv)=iconst(isyminv)+1
1932 : : !FB const_tot=const_tot+1
1933 : : !FB write(16,*)' The eigenvalue',ii
1934 : : !FB write(16,*)' is equal to ',lambda
1935 [ + + ]: 951 : do mu=1,3
1936 : 954 : alphaij(isyminv,mu,iconst(isyminv))=eigenvectors(mu,ii)
1937 : : end do
1938 : : !FB write(16,*)' Real & imaginary parts of the eigenvectors product:'
1939 : : !FB write(16,'(3(f16.12,1x))') real(alphaij(isyminv,:,iconst(isyminv)))
1940 : : !FB write(16,'(3(f16.12,1x))') aimag(alphaij(isyminv,:,iconst(isyminv)))
1941 : : end do !ii
1942 [ + + ]: 15325 : else if (order==2) then
1943 [ + + ]: 6952 : do ii=1,3
1944 [ + + ]: 22594 : do jj=1,3
1945 : 15642 : lambda=eigenvalues(ii)*eigenvalues(jj)
1946 [ + - + + : 9014 : if (((abs(real(lambda)-1.d0).lt.tol6).and.(abs(aimag(lambda)).lt.tol6).and.(inv==1)).or.&
+ + ]
1947 [ + + ]: 15642 : & ((abs(real(lambda)-1.d0).lt.tol6).and.(abs(aimag(lambda)).lt.tol6).and.(inv==2).and.(ii==jj))) cycle
1948 : 9456 : unchanged(isyminv)=.true.
1949 : 9456 : iconst(isyminv)=iconst(isyminv)+1
1950 : : !FB const_tot=const_tot+1
1951 : : !FB write(16,*)' The product of eigenvalues',ii,jj
1952 : : !FB write(16,*)' is equal to ',lambda
1953 [ + + ]: 37824 : do mu=1,3
1954 [ + + ]: 122928 : do nu=1,3
1955 : 113472 : pp(mu,nu)=eigenvectors(mu,ii)*eigenvectors(nu,jj)
1956 : : end do
1957 : : end do
1958 [ + + ]: 43038 : do mu=1,3
1959 [ + + ]: 129114 : do nu=1,3
1960 [ + + ]: 113472 : if (inv==1) then
1961 : 31500 : alphaij(isyminv,(mu-1)*3+nu,iconst(isyminv))=pp(mu,nu)
1962 [ + - ]: 53604 : else if (inv==2) then
1963 : 53604 : alphaij(isyminv,(mu-1)*3+nu,iconst(isyminv))=lambda*pp(mu,nu)-pp(nu,mu)
1964 : : else
1965 : 0 : ABI_BUG('This symetry is neither Keptinvariant nor Reversed')
1966 : : end if
1967 : : end do
1968 : : end do
1969 : : !FB write(16,*)' Real & imaginary parts of the eigenvectors product:'
1970 : : !FB write(16,'(9(f16.12,1x))') real(alphaij(isyminv,:,iconst(isyminv)))
1971 : : !FB write(16,'(9(f16.12,1x))') aimag(alphaij(isyminv,:,iconst(isyminv)))
1972 : : end do !jj
1973 : : end do !ii
1974 [ + + ]: 13587 : else if (order==3) then
1975 [ + + ]: 16028 : do ii=1,3
1976 [ + + ]: 52091 : do jj=1,3
1977 [ + + ]: 156273 : do kk=1,3
1978 : 108189 : lambda=eigenvalues(ii)*eigenvalues(jj)*eigenvalues(kk)
1979 : : if (((abs(real(lambda)-1.d0).lt.tol6).and.(abs(aimag(lambda)).lt.tol6).and.(inv==1)).or.&
1980 : : & ((abs(real(lambda)-1.d0).lt.tol6).and.(abs(aimag(lambda)).lt.tol6).and.(inv==2).and.(jj==kk)).or.&
1981 : : & ((abs(real(lambda)-1.d0).lt.tol6).and.(abs(aimag(lambda)).lt.tol6).and.(inv==3).and.(ii==jj)).or.&
1982 : : & ((abs(real(lambda)-1.d0).lt.tol6).and.(abs(aimag(lambda)).lt.tol6).and.(inv==4).and.(ii==jj).and.(jj==kk)).or.&
1983 [ + - + + : 108189 : & ((abs(real(lambda)-1.d0).lt.tol6).and.(abs(aimag(lambda)).lt.tol6).and.(inv==5).and.(ii==jj).and.(jj==kk)).or.&
+ + + + +
- + + + -
+ + + + +
+ + - + +
+ + + - +
+ ]
1984 [ + + ]: 108189 : & ((abs(real(lambda)-1.d0).lt.tol6).and.(abs(aimag(lambda)).lt.tol6).and.(inv==6).and.(ii==kk))) cycle
1985 : 93508 : unchanged(isyminv)=.true.
1986 : 93508 : iconst(isyminv)=iconst(isyminv)+1
1987 : : !FB const_tot=const_tot+1
1988 : : !FB write(16,*)' The product of eigenvalues',ii,jj
1989 : : !FB write(16,*)' is equal to ',lambda
1990 [ + + ]: 374032 : do mu=1,3
1991 [ + + ]: 1215604 : do nu=1,3
1992 [ + + ]: 3646812 : do xi=1,3
1993 : 3366288 : ppp(mu,nu,xi)=eigenvectors(mu,ii)*eigenvectors(nu,jj)*eigenvectors(xi,kk)
1994 : : end do !xi
1995 : : end do !nu
1996 : : end do !mu
1997 [ + + ]: 410095 : do mu=1,3
1998 [ + + ]: 1230285 : do nu=1,3
1999 [ + + ]: 3646812 : do xi=1,3
2000 : 841572 : if (inv==1) then
2001 : 349353 : alphaij(isyminv,(mu-1)*9+(nu-1)*3+xi,iconst(isyminv))=ppp(mu,nu,xi)
2002 : : else if (inv==2) then
2003 : 432648 : alphaij(isyminv,(mu-1)*9+(nu-1)*3+xi,iconst(isyminv))=lambda*ppp(mu,nu,xi)-ppp(mu,xi,nu)
2004 : : else if (inv==3) then
2005 : 456327 : alphaij(isyminv,(mu-1)*9+(nu-1)*3+xi,iconst(isyminv))=lambda*ppp(mu,nu,xi)-ppp(nu,mu,xi)
2006 : : else if (inv==4) then
2007 : 436752 : alphaij(isyminv,(mu-1)*9+(nu-1)*3+xi,iconst(isyminv))=lambda*ppp(mu,nu,xi)-ppp(nu,xi,mu)
2008 : : else if (inv==5) then
2009 : 436752 : alphaij(isyminv,(mu-1)*9+(nu-1)*3+xi,iconst(isyminv))=lambda*ppp(mu,nu,xi)-ppp(xi,mu,nu)
2010 : : else if (inv==6) then
2011 : 412884 : alphaij(isyminv,(mu-1)*9+(nu-1)*3+xi,iconst(isyminv))=lambda*ppp(mu,nu,xi)-ppp(xi,nu,mu)
2012 : : else
2013 : 0 : ABI_BUG('This symetry is neither Keptinvariant nor Reversed')
2014 : : end if
2015 : : end do !xi
2016 : : end do !nu
2017 : : end do !mu
2018 : : !FB write(16,*)' Real & imaginary parts of the eigenvectors product:'
2019 : : !FB write(16,'(27(f16.12,1x))') real(alphaij(isyminv,:,iconst(isyminv)))
2020 : : !FB write(16,'(27(f16.12,1x))') aimag(alphaij(isyminv,:,iconst(isyminv)))
2021 : : end do !kk
2022 : : end do !jj
2023 : : end do !ii
2024 [ + - ]: 9580 : else if (order==4) then
2025 [ + + ]: 38320 : do ii=1,3
2026 [ + + ]: 124540 : do jj=1,3
2027 [ + + ]: 373620 : do kk=1,3
2028 [ + + ]: 1120860 : do ll=1,3
2029 : 775980 : lambda=eigenvalues(ii)*eigenvalues(jj)*eigenvalues(kk)*eigenvalues(ll)
2030 [ + + + - ]: 775980 : if ((abs(real(lambda)-1.d0).lt.tol6).and.(abs(aimag(lambda)).lt.tol6)) then
2031 : : if ((inv==1 ) .or.& !\Phi4_ijkl
2032 : : & ((inv==2 ).and.(jj==kk)) .or.& !\Phi4_ikjl
2033 : : & ((inv==3 ).and.(ii==jj)) .or.& !\Phi4_jikl
2034 : : & ((inv==4 ).and.(ii==jj).and.(jj==kk)) .or.& !\Phi4_jkil
2035 : : & ((inv==5 ).and.(ii==jj).and.(jj==kk)) .or.& !\Phi4_kijl
2036 : : & ((inv==6 ).and.(ii==kk)) .or.& !\Phi4_kjil
2037 : :
2038 : : & ((inv==7 ).and.(kk==ll)) .or.& !\Phi4_ijlk
2039 : : & ((inv==8 ).and.(jj==kk).and.(kk==ll)) .or.& !\Phi4_iklj
2040 : : & ((inv==9 ).and.(ii==jj).and.(kk==ll)) .or.& !\Phi4_jilk
2041 : : & ((inv==10).and.(ii==jj).and.(jj==kk).and.(kk==ll)).or.& !\Phi4_jkli
2042 : : & ((inv==11).and.(ii==jj).and.(jj==kk).and.(kk==ll)).or.& !\Phi4_kilj
2043 : : & ((inv==12).and.(ii==kk).and.(kk==ll)) .or.& !\Phi4_kjli
2044 : :
2045 : : & ((inv==13).and.(jj==kk).and.(kk==ll)) .or.& !\Phi4_iljk
2046 : : & ((inv==14).and.(jj==ll)) .or.& !\Phi4_ilkj
2047 : : & ((inv==15).and.(ii==jj).and.(jj==kk).and.(kk==ll)).or.& !\Phi4_jlik
2048 : : & ((inv==16).and.(ii==jj).and.(jj==ll)) .or.& !\Phi4_jlki
2049 : : & ((inv==17).and.(ii==kk).and.(jj==ll)) .or.& !\Phi4_klij
2050 : : & ((inv==18).and.(ii==jj).and.(jj==kk).and.(kk==ll)).or.& !\Phi4_klji
2051 : :
2052 : : & ((inv==19).and.(ii==jj).and.(jj==kk).and.(kk==ll)).or.& !\Phi4_lijk
2053 : : & ((inv==20).and.(ii==jj).and.(jj==ll)) .or.& !\Phi4_likj
2054 : : & ((inv==21).and.(ii==kk).and.(kk==ll)) .or.& !\Phi4_ljik
2055 : : & ((inv==22).and.(ii==ll)) .or.& !\Phi4_ljki
2056 [ + + + + : 310696 : & ((inv==23).and.(ii==jj).and.(jj==kk).and.(kk==ll)).or.& !\Phi4_lkij
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + ]
2057 : : & ((inv==24).and.(ii==ll).and.(jj==kk))) cycle !\Phi4_lkji
2058 : : end if
2059 : 706594 : unchanged(isyminv)=.true.
2060 : 706594 : iconst(isyminv)=iconst(isyminv)+1
2061 : : !FB const_tot=const_tot+1
2062 : : !FB write(16,*)' The product of eigenvalues',ii,jj
2063 : : !FB write(16,*)' is equal to ',lambda
2064 [ + + ]: 2826376 : do mu=1,3
2065 [ + + ]: 9185722 : do nu=1,3
2066 [ + + ]: 27557166 : do xi=1,3
2067 [ + + ]: 82671498 : do zeta=1,3
2068 : 76312152 : pppp(mu,nu,xi,zeta)=eigenvectors(mu,ii)*eigenvectors(nu,jj)*eigenvectors(xi,kk)*eigenvectors(zeta,ll)
2069 : : end do !zeta
2070 : : end do !xi
2071 : : end do !nu
2072 : : end do !mu
2073 [ + + ]: 3085036 : do mu=1,3
2074 [ + + ]: 9255108 : do nu=1,3
2075 [ + + ]: 27557166 : do xi=1,3
2076 [ + + ]: 82671498 : do zeta=1,3
2077 [ + + + + : 57234114 : itemp=(mu-1)*27+(nu-1)*9+(xi-1)*3+zeta
+ + + + +
+ + + + +
+ + + + +
+ + + + +
- ]
2078 : 20938446 : if (inv==1) then ; alphaij(isyminv,itemp,iconst(isyminv))=pppp(mu,nu,xi,zeta)
2079 : 2433240 : else if (inv==2 ) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(mu,xi,nu,zeta)
2080 : 2586492 : else if (inv==3 ) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(nu,mu,xi,zeta)
2081 : 2571588 : else if (inv==4 ) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(nu,xi,mu,zeta)
2082 : 2571588 : else if (inv==5 ) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(xi,mu,nu,zeta)
2083 : 2328750 : else if (inv==6 ) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(xi,nu,mu,zeta)
2084 : :
2085 : 2470662 : else if (inv==7 ) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(mu,nu,zeta,xi)
2086 : 2410722 : else if (inv==8 ) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(mu,xi,zeta,nu)
2087 : 2562516 : else if (inv==9 ) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(nu,mu,zeta,xi)
2088 : 2358072 : else if (inv==10) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(nu,xi,zeta,mu)
2089 : 2358072 : else if (inv==11) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(xi,mu,zeta,nu)
2090 : 2311092 : else if (inv==12) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(xi,nu,zeta,mu)
2091 : :
2092 : 2410722 : else if (inv==13) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(mu,zeta,nu,xi)
2093 : 2186190 : else if (inv==14) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(mu,zeta,xi,nu)
2094 : 2358072 : else if (inv==15) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(nu,zeta,mu,xi)
2095 : 2311092 : else if (inv==16) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(nu,zeta,xi,mu)
2096 : 2454624 : else if (inv==17) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(xi,zeta,mu,nu)
2097 : 2573208 : else if (inv==18) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(xi,zeta,nu,mu)
2098 : :
2099 : 2358072 : else if (inv==19) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(zeta,mu,nu,xi)
2100 : 2311092 : else if (inv==20) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(zeta,mu,xi,nu)
2101 : 2311092 : else if (inv==21) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(zeta,nu,mu,xi)
2102 : 2103084 : else if (inv==22) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(zeta,nu,xi,mu)
2103 : 2573208 : else if (inv==23) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(zeta,xi,mu,nu)
2104 : 2460456 : else if (inv==24) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(zeta,xi,nu,mu)
2105 : 0 : else ; ABI_BUG('This symetry is neither Keptinvariant nor Reversed')
2106 : : end if
2107 : : end do !zeta
2108 : : end do !xi
2109 : : end do !nu
2110 : : end do !mu
2111 : : !FB write(16,*)' Real & imaginary parts of the eigenvectors product:'
2112 : : !FB write(16,'(81(f16.12,1x))') real(alphaij(isyminv,:,iconst(isyminv)))
2113 : : !FB write(16,'(81(f16.12,1x))') aimag(alphaij(isyminv,:,iconst(isyminv)))
2114 : : end do !ll
2115 : : end do !kk
2116 : : end do !jj
2117 : : end do !ii
2118 : : else
2119 : 0 : ABI_BUG('Only the first, second, third and fourth order are allowed')
2120 : : end if
2121 : :
2122 : :
2123 : : !FB=================================================================
2124 : : !FB======== TO CLEAN ===============================================
2125 : : !FB=================================================================
2126 : 15432 : nconst_loc=const_tot+iconst(isyminv)
2127 : 15432 : ii=0
2128 [ - + + - : 123884634 : ABI_MALLOC(tab_vec,(norder,nconst_loc)); tab_vec(:,:)=czero
+ + + - +
- + + - +
- + + + +
+ ]
2129 [ + + ]: 6045393 : do itemp=1,isyminv
2130 [ + + ]: 6045393 : if (unchanged(itemp)) then
2131 [ + + ]: 6447159 : do jj=1,iconst(itemp)
2132 : 1655358 : ii=ii+1
2133 [ + + ]: 128614707 : tab_vec(:,ii)=alphaij(itemp,:,jj)
2134 : : end do
2135 : : end if
2136 : : end do
2137 : :
2138 [ + + ]: 1655364 : do kk=2,nconst_loc
2139 [ + + ]: 114878467 : do jj=1,kk-1
2140 [ + + ]: 9017348536 : prod_scal=sum( real(tab_vec(:,jj))* real(tab_vec(:,jj))+aimag(tab_vec(:,jj))*aimag(tab_vec(:,jj)))
2141 [ + + ]: 114863035 : if (abs(prod_scal).gt.tol8) then
2142 [ + + + + ]: 13645954863 : tab_vec(:,kk)=tab_vec(:,kk)-sum(tab_vec(:,kk)*conjg(tab_vec(:,jj)))/dcmplx(prod_scal,zero)*tab_vec(:,jj)
2143 [ + + ]: 6866115828 : do ii=1,norder
2144 [ + + ]: 6779839035 : if (abs( real(tab_vec(ii,kk))).lt.tol8) tab_vec(ii,kk)=dcmplx(zero,aimag(tab_vec(ii,kk)))
2145 [ + + ]: 6866115828 : if (abs(aimag(tab_vec(ii,kk))).lt.tol8) tab_vec(ii,kk)=dcmplx( real(tab_vec(ii,kk)),zero)
2146 : : end do
2147 : : end if
2148 : : end do
2149 : : end do
2150 : :
2151 : : ! On stocke les vecteurs non-nuls
2152 [ - + + - : 123869202 : ABI_MALLOC(temp ,(norder,nconst_loc)); temp(:,:) =czero
+ + + - -
+ - + + +
+ + ]
2153 : 15432 : ii=0
2154 [ + + ]: 1670790 : do kk=1,nconst_loc
2155 [ + + ]: 123822906 : prod_scal=sum( real(tab_vec(:,kk))* real(tab_vec(:,kk))+aimag(tab_vec(:,kk))*aimag(tab_vec(:,kk)))
2156 [ + + ]: 1670790 : if (abs(prod_scal).gt.tol8) then
2157 : 852331 : ii=ii+1
2158 [ + + ]: 63503200 : temp(:,ii)=tab_vec(:,kk)/dsqrt(prod_scal)
2159 : : end if
2160 : : end do
2161 [ - + ]: 15432 : ABI_FREE(tab_vec)
2162 : 15432 : iconst(isyminv)=ii-const_tot
2163 : 15432 : const_tot=const_tot+iconst(isyminv)
2164 : :
2165 : 15432 : ii=0
2166 [ + + + + : 70953568626 : alphaij(:,:,:)=czero
+ + ]
2167 [ + + ]: 6045393 : do itemp=1,isyminv
2168 [ + + ]: 6045393 : if (unchanged(itemp)) then
2169 [ + + ]: 5644132 : do jj=1,iconst(itemp)
2170 : 852331 : ii=ii+1
2171 [ + + ]: 68295001 : alphaij(itemp,:,jj)=temp(:,ii)
2172 : : end do
2173 : : end if
2174 : : end do
2175 [ - + ]: 15432 : ABI_FREE(temp)
2176 : : !FB=================================================================
2177 : : !FB======== TO CLEAN ===============================================
2178 : : !FB=================================================================
2179 : :
2180 : :
2181 : :
2182 : :
2183 : :
2184 : :
2185 : : ! WARNING: There are some minimum and maximum of constraints
2186 [ + + + + ]: 15432 : if (order==1.and.(iconst(isyminv).eq.3)) then
2187 : 6 : ncoeff=0
2188 [ + + + + ]: 78 : proj(:,:,ishell)=zero
2189 [ - + ]: 6 : ABI_FREE(unchanged)
2190 [ - + ]: 6 : ABI_FREE(alphaij)
2191 : 6 : ABI_FREE(iconst)
2192 : 6 : return
2193 [ + + - + ]: 15426 : else if (order==1.and.(iconst(isyminv).gt.3)) then
2194 : 0 : ABI_BUG(' First order : There are more than 3 constraints')
2195 : : end if
2196 [ + + - + ]: 15426 : if (order==2.and.(iconst(isyminv).gt.8)) then
2197 : 0 : ABI_BUG(' Second order : There are more than 8 constraints')
2198 : : end if
2199 [ + + - + ]: 15426 : if (order==3.and.(iconst(isyminv).gt.27)) then
2200 : 0 : ABI_BUG(' Third order : There are more than 27 constraints')
2201 : : end if
2202 [ + + + - ]: 16002 : if (order==4.and.(iconst(isyminv).gt.81)) then
2203 : 0 : ABI_BUG(' Fourth order : There are more than 81 constraints')
2204 : : end if
2205 : : end do !isyminv
2206 : : ! ================================================================================================
2207 : : ! =========== End big loop over symetries and facorder ===========================================
2208 : : ! ================================================================================================
2209 : 576 : nconst_perm=0
2210 : : ! The (iik, iji, ijj and iii) third order IFCs are symmetric with respect to some permutations.
2211 : : ! Some constraints have to be added :
2212 [ + + ]: 576 : if (order.eq.3) then
2213 [ + + + - : 57 : if ((iatcell.eq.jatom).or.(iatcell.eq.katom).or.(jatom.eq.katom)) then
+ + ]
2214 : 44 : nconst_perm=5
2215 [ + - ]: 44 : if (MPIdata%iam_master) write(16,'(a)')'=========== The IFCs are symmetric'
2216 : 44 : const_tot=const_tot+nconst_perm*norder
2217 [ - + + - : 193820 : ABI_MALLOC(constraints,(nconst_perm,norder,norder)) ; constraints(:,:,:)=czero
- + + - +
- + - - +
- + + + +
+ + + ]
2218 : 44 : ii=0
2219 [ + + ]: 176 : do mu=1,3
2220 [ + + ]: 585 : do nu=1,3
2221 [ + + ]: 1716 : do xi=1,3
2222 : 1188 : ii=ii+1
2223 [ + + ]: 1188 : if (iatcell.eq.jatom) then
2224 [ + + ]: 1053 : if (mu.eq.nu) cycle
2225 : 702 : constraints(1,(mu-1)*9+(nu-1)*3+xi,ii)= cone
2226 : 702 : constraints(1,(nu-1)*9+(mu-1)*3+xi,ii)=-cone
2227 : : end if
2228 [ + + ]: 837 : if (iatcell.eq.katom) then
2229 [ + + ]: 378 : if (mu.eq.xi) cycle
2230 : 252 : constraints(2,(mu-1)*9+(nu-1)*3+xi,ii)= cone
2231 : 252 : constraints(2,(xi-1)*9+(nu-1)*3+mu,ii)=-cone
2232 : : end if
2233 [ + + ]: 711 : if (jatom.eq.katom) then
2234 [ + + ]: 387 : if (nu.eq.xi) cycle
2235 : 216 : constraints(3,(mu-1)*9+(nu-1)*3+xi,ii)= cone
2236 : 216 : constraints(3,(mu-1)*9+(xi-1)*3+nu,ii)=-cone
2237 : : end if
2238 [ + + + + ]: 540 : if ((iatcell.eq.jatom).and.(jatom.eq.katom)) then
2239 [ + - ]: 126 : if ((nu.eq.xi).and.(nu.eq.mu)) cycle
2240 : 126 : constraints(4,(mu-1)*9+(nu-1)*3+xi,ii)= cone
2241 : 126 : constraints(4,(xi-1)*9+(mu-1)*3+nu,ii)=-cone
2242 : : end if
2243 [ + + + + ]: 936 : if ((iatcell.eq.jatom).and.(jatom.eq.katom)) then
2244 [ + - ]: 126 : if ((nu.eq.xi).and.(nu.eq.mu)) cycle
2245 : 126 : constraints(5,(mu-1)*9+(nu-1)*3+xi,ii)= cone
2246 : 126 : constraints(5,(nu-1)*9+(xi-1)*3+mu,ii)=-cone
2247 : : end if
2248 : : end do
2249 : : end do
2250 : : end do
2251 : : end if
2252 : : end if
2253 : :
2254 : : ! The (iikl, ijil, ijki, ijjl, ijkj, ijkk, iiil, iiki, ijii, ijjj, iiii)
2255 : : ! fourth order IFCs are symmetric with respect to some permutations.
2256 : : ! Some constraints have to be added :
2257 [ + + ]: 576 : if (order.eq.4) then
2258 : : if ((iatcell.eq.jatom).or.(iatcell.eq.katom).or.(iatcell.eq.latom)&
2259 [ + + + - : 49 : & .or.(jatom.eq.katom).or.(jatom.eq.latom).or.(katom.eq.latom)) then
+ - + + +
- - + ]
2260 : 47 : nconst_perm=17
2261 [ + - ]: 47 : if (MPIdata%iam_master) write(16,'(a)')'=========== The IFCs are symmetric'
2262 : 47 : const_tot=const_tot+nconst_perm*norder
2263 [ - + + - : 5554601 : ABI_MALLOC(constraints,(nconst_perm,norder,norder)) ; constraints(:,:,:)=czero
- + + - +
- + - - +
- + - + +
+ + + +
+ ]
2264 : 47 : ii=0
2265 [ + + ]: 188 : do mu=1,3
2266 [ + + ]: 613 : do nu=1,3
2267 [ + + ]: 1833 : do xi=1,3
2268 [ + + ]: 5499 : do zeta=1,3
2269 : 3807 : ii=ii+1
2270 [ + + ]: 3807 : if (iatcell.eq.jatom) then
2271 [ + + ]: 3240 : if (mu.eq.nu) cycle
2272 : 2160 : constraints(1,(mu-1)*27+(nu-1)*9+(xi-1)*3+zeta,ii)= cone
2273 : 2160 : constraints(1,(nu-1)*27+(mu-1)*9+(xi-1)*3+zeta,ii)=-cone
2274 : : end if
2275 [ + + ]: 2727 : if (iatcell.eq.katom) then
2276 [ + + ]: 1188 : if (mu.eq.xi) cycle
2277 : 792 : constraints(2,(mu-1)*27+(nu-1)*9+(xi-1)*3+zeta,ii)= cone
2278 : 792 : constraints(2,(xi-1)*27+(nu-1)*9+(mu-1)*3+zeta,ii)=-cone
2279 : : end if
2280 [ + + ]: 2331 : if (iatcell.eq.latom) then
2281 [ + + ]: 432 : if (mu.eq.zeta) cycle
2282 : 288 : constraints(3,(mu -1)*27+(nu-1)*9+(xi-1)*3+zeta,ii)= cone
2283 : 288 : constraints(3,(zeta-1)*27+(nu-1)*9+(xi-1)*3+mu ,ii)=-cone
2284 : : end if
2285 [ + + ]: 2187 : if (jatom.eq.katom) then
2286 [ + + ]: 1215 : if (nu.eq.xi) cycle
2287 : 702 : constraints(4,(mu-1)*27+(nu-1)*9+(xi-1)*3+zeta,ii)= cone
2288 : 702 : constraints(4,(mu-1)*27+(xi-1)*9+(nu-1)*3+zeta,ii)=-cone
2289 : : end if
2290 [ + + ]: 1674 : if (jatom.eq.latom) then
2291 [ + + ]: 306 : if (nu.eq.zeta) cycle
2292 : 180 : constraints(5,(mu-1)*27+(nu -1)*9+(xi-1)*3+zeta,ii)= cone
2293 : 180 : constraints(5,(mu-1)*27+(zeta-1)*9+(xi-1)*3+nu ,ii)=-cone
2294 : : end if
2295 [ + + ]: 1548 : if (katom.eq.latom) then
2296 [ + + ]: 720 : if (xi.eq.zeta) cycle
2297 : 414 : constraints(6,(mu-1)*27+(nu-1)*9+(xi -1)*3+zeta,ii)= cone
2298 : 414 : constraints(6,(mu-1)*27+(nu-1)*9+(zeta-1)*3+xi ,ii)=-cone
2299 : : end if
2300 : :
2301 [ + + + + ]: 1242 : if ((iatcell.eq.jatom).and.(jatom.eq.katom)) then
2302 [ + - ]: 180 : if ((mu.eq.nu).and.(nu.eq.xi)) cycle
2303 : 180 : constraints(7,(mu-1)*27+(nu-1)*9+(xi-1)*3+zeta,ii)= cone
2304 : 180 : constraints(7,(xi-1)*27+(mu-1)*9+(nu-1)*3+zeta,ii)=-cone
2305 : : end if
2306 [ + + + + ]: 1242 : if ((iatcell.eq.jatom).and.(jatom.eq.katom)) then
2307 [ + - ]: 180 : if ((mu.eq.nu).and.(nu.eq.xi)) cycle
2308 : 180 : constraints(8,(mu-1)*27+(nu-1)*9+(xi-1)*3+zeta,ii)= cone
2309 : 180 : constraints(8,(nu-1)*27+(xi-1)*9+(mu-1)*3+zeta,ii)=-cone
2310 : : end if
2311 : :
2312 [ + + - + ]: 1242 : if ((iatcell.eq.jatom).and.(jatom.eq.latom)) then
2313 [ # # ]: 0 : if ((mu.eq.nu).and.(nu.eq.zeta)) cycle
2314 : 0 : constraints(9,(mu-1)*27+(nu -1)*9+(xi-1)*3+zeta,ii)= cone
2315 : 0 : constraints(9,(nu-1)*27+(zeta-1)*9+(xi-1)*3+mu ,ii)=-cone
2316 : : end if
2317 [ + + - + ]: 1242 : if ((iatcell.eq.jatom).and.(jatom.eq.latom)) then
2318 [ # # ]: 0 : if ((mu.eq.nu).and.(nu.eq.zeta)) cycle
2319 : 0 : constraints(10,(mu -1)*27+(nu-1)*9+(xi-1)*3+zeta,ii)= cone
2320 : 0 : constraints(10,(zeta-1)*27+(mu-1)*9+(xi-1)*3+nu ,ii)=-cone
2321 : : end if
2322 : :
2323 [ + + - + ]: 1242 : if ((iatcell.eq.katom).and.(katom.eq.latom)) then
2324 [ # # ]: 0 : if ((mu.eq.xi).and.(xi.eq.zeta)) cycle
2325 : 0 : constraints(11,(mu-1)*27+(nu-1)*9+(xi -1)*3+zeta,ii)= cone
2326 : 0 : constraints(11,(xi-1)*27+(nu-1)*9+(zeta-1)*3+mu ,ii)=-cone
2327 : : end if
2328 [ + + - + ]: 1242 : if ((iatcell.eq.katom).and.(katom.eq.latom)) then
2329 [ # # ]: 0 : if ((mu.eq.xi).and.(xi.eq.zeta)) cycle
2330 : 0 : constraints(12,(mu -1)*27+(nu-1)*9+(xi-1)*3+zeta,ii)= cone
2331 : 0 : constraints(12,(zeta-1)*27+(nu-1)*9+(mu-1)*3+xi ,ii)=-cone
2332 : : end if
2333 : :
2334 [ + + + + ]: 1242 : if ((jatom.eq.katom).and.(katom.eq.latom)) then
2335 [ + - ]: 54 : if ((nu.eq.xi).and.(xi.eq.zeta)) cycle
2336 : 54 : constraints(13,(mu-1)*27+(nu-1)*9+(xi -1)*3+zeta,ii)= cone
2337 : 54 : constraints(13,(mu-1)*27+(xi-1)*9+(zeta-1)*3+nu ,ii)=-cone
2338 : : end if
2339 [ + + + + ]: 1242 : if ((jatom.eq.katom).and.(katom.eq.latom)) then
2340 [ + - ]: 54 : if ((nu.eq.xi).and.(xi.eq.zeta)) cycle
2341 : 54 : constraints(14,(mu-1)*27+(nu -1)*9+(xi-1)*3+zeta,ii)= cone
2342 : 54 : constraints(14,(mu-1)*27+(zeta-1)*9+(nu-1)*3+xi ,ii)=-cone
2343 : : end if
2344 : :
2345 [ + + + + : 1242 : if ((iatcell.eq.jatom).and.(jatom.eq.katom).and.(katom.eq.latom)) then
- + ]
2346 [ # # # # ]: 0 : if ((mu.eq.nu).and.(nu.eq.xi).and.(xi.eq.zeta)) cycle
2347 : 0 : constraints(15,(mu-1)*27+(nu-1)*9+(xi -1)*3+zeta,ii)= cone
2348 : 0 : constraints(15,(nu-1)*27+(xi-1)*9+(zeta-1)*3+mu ,ii)=-cone
2349 : : end if
2350 [ + + + + : 1242 : if ((iatcell.eq.jatom).and.(jatom.eq.katom).and.(katom.eq.latom)) then
- + ]
2351 [ # # # # ]: 0 : if ((mu.eq.nu).and.(nu.eq.xi).and.(xi.eq.zeta)) cycle
2352 : 0 : constraints(16,(mu-1)*27+(nu -1)*9+(xi-1)*3+zeta,ii)= cone
2353 : 0 : constraints(16,(xi-1)*27+(zeta-1)*9+(mu-1)*3+nu ,ii)=-cone
2354 : : end if
2355 [ + + + + : 2511 : if ((iatcell.eq.jatom).and.(jatom.eq.katom).and.(katom.eq.latom)) then
- + ]
2356 [ # # # # ]: 0 : if ((mu.eq.nu).and.(nu.eq.xi).and.(xi.eq.zeta)) cycle
2357 : 0 : constraints(17,(mu -1)*27+(nu-1)*9+(xi-1)*3+zeta,ii)= cone
2358 : 0 : constraints(17,(zeta-1)*27+(mu-1)*9+(nu-1)*3+xi ,ii)=-cone
2359 : : end if
2360 : : end do
2361 : : end do
2362 : : end do
2363 : : end do
2364 : : end if
2365 : : end if
2366 : :
2367 : : ! In the case where the matrix has norder**2 inequivalent and non-zero elements
2368 [ + + ]: 576 : if (const_tot==0) then
2369 : 49 : write(message,'(a,1x,i3,1x,a)') 'For shell number=',ishell,'there is no symetry operation reducing the number of coefficients'
2370 : 49 : ABI_WARNING(message)
2371 [ + + + + ]: 4459 : proj(:,:,ishell)=0.d0
2372 [ + + ]: 490 : do ii=1,norder
2373 : 490 : proj(ii,ii,ishell)=1.d0
2374 : : end do
2375 : 49 : ncoeff=norder
2376 [ - + ]: 49 : ABI_FREE(unchanged)
2377 [ - + ]: 49 : ABI_FREE(alphaij)
2378 [ - + ]: 49 : ABI_FREE(iconst)
2379 : 49 : return
2380 : : end if
2381 : :
2382 : : ! When some constraints have been found
2383 : 527 : ncount=const_tot
2384 [ + - ]: 527 : if (MPIdata%iam_master) then
2385 : 527 : write(16,'(a,1x,i7,1x,a)') 'There is a total of ',ncount,' non-independant constraints for this shell'
2386 : : end if
2387 : 527 : ii=0
2388 [ - + + - : 5809944 : ABI_MALLOC(tab_vec,(norder,ncount)); tab_vec(:,:)=czero
- + + - +
- + - - +
- + + + +
+ ]
2389 [ - + + - : 5809417 : ABI_MALLOC(temp ,(norder,ncount)); temp(:,:) =czero
- + + - -
+ - + + +
+ + ]
2390 [ + + ]: 86733 : do isyminv=1,nsyminv
2391 [ + + ]: 86733 : if (unchanged(isyminv)) then
2392 [ + + ]: 22144 : do jj=1,iconst(isyminv)
2393 : 6724 : ii=ii+1
2394 [ + + ]: 349978 : tab_vec(:,ii)=alphaij(isyminv,:,jj)
2395 : : end do
2396 : : end if
2397 : : end do
2398 [ - + ]: 527 : ABI_FREE(unchanged)
2399 [ - + ]: 527 : ABI_FREE(alphaij)
2400 [ - + ]: 527 : ABI_FREE(iconst)
2401 : : ! Add the constraints coming from the symmetry of the IFCs (at the 3rd order)
2402 [ + + ]: 527 : if (nconst_perm.gt.0) then
2403 [ + + ]: 5086 : do jj=1,norder
2404 [ + + ]: 75745 : do kk=1,nconst_perm
2405 : 70659 : ii=ii+1
2406 [ + + ]: 5478273 : tab_vec(:,ii)=constraints(kk,:,jj)
2407 : : end do
2408 : : end do
2409 [ - + ]: 91 : ABI_FREE(constraints)
2410 : : end if
2411 [ - + ]: 527 : if (ii.ne.ncount) then
2412 : 0 : write(message,'(i7,1x,a,1x,i7)') ii,' non equal to ',ncount
2413 : 0 : ABI_BUG(message)
2414 : : end if
2415 [ + + ]: 9746 : do ii=1,norder
2416 [ + + ]: 5740199 : do jj=1,ncount
2417 [ + + ]: 5730453 : if (abs( real(tab_vec(ii,jj))).lt.tol8) tab_vec(ii,jj)=dcmplx(zero,aimag(tab_vec(ii,jj)))
2418 [ + + ]: 5739672 : if (abs(aimag(tab_vec(ii,jj))).lt.tol8) tab_vec(ii,jj)=dcmplx( real(tab_vec(ii,jj)),zero)
2419 : : end do
2420 : : end do
2421 : :
2422 : : ! On stocke les vecteurs non-nuls
2423 : 527 : ii=0
2424 [ + + ]: 77910 : do kk=1,ncount
2425 [ + + ]: 5807836 : prod_scal=sum( real(tab_vec(:,kk))* real(tab_vec(:,kk))+aimag(tab_vec(:,kk))*aimag(tab_vec(:,kk)))
2426 [ + + ]: 77910 : if (abs(prod_scal).gt.tol8) then
2427 : 13150 : ii=ii+1
2428 [ + + ]: 784702 : temp(:,ii)=tab_vec(:,kk)/dsqrt(prod_scal)
2429 : : end if
2430 : : end do
2431 : 527 : ncount=ii
2432 [ - + ]: 527 : ABI_FREE(tab_vec)
2433 [ - + + - : 786810 : ABI_MALLOC(tab_vec,(norder,ncount)); tab_vec(:,1:ncount)=temp(:,1:ncount)
- + + - +
- + - - +
- + + + +
+ ]
2434 [ - + ]: 527 : ABI_FREE(temp)
2435 [ - + + - : 786283 : ABI_MALLOC(temp ,(norder,ncount)); temp(:,:) =czero
- + + - -
+ - + + +
+ + ]
2436 : :
2437 : : ! L'ensemble des vecteurs reduisants l'espace de R^norder a R^n ne forment pas une base
2438 : : ! independante. Il faut donc trouver les vecteurs independants.
2439 : : ! --> Orthogonalisation de Gram-Schmidt
2440 [ + + ]: 13150 : do kk=2,ncount
2441 [ + + ]: 870340 : do jj=1,kk-1
2442 [ + + ]: 65795250 : prod_scal=sum( real(tab_vec(:,jj))* real(tab_vec(:,jj))+aimag(tab_vec(:,jj))*aimag(tab_vec(:,jj)))
2443 [ + + ]: 869813 : if (abs(prod_scal).gt.tol8) then
2444 [ + + + + ]: 81343776 : tab_vec(:,kk)=tab_vec(:,kk)-sum(tab_vec(:,kk)*conjg(tab_vec(:,jj)))/dcmplx(prod_scal,zero)*tab_vec(:,jj)
2445 [ + + ]: 40939968 : do ii=1,norder
2446 [ + + ]: 40403808 : if (abs( real(tab_vec(ii,kk))).lt.tol8) tab_vec(ii,kk)=dcmplx(zero,aimag(tab_vec(ii,kk)))
2447 [ + + ]: 40939968 : if (abs(aimag(tab_vec(ii,kk))).lt.tol8) tab_vec(ii,kk)=dcmplx( real(tab_vec(ii,kk)),zero)
2448 : : end do
2449 : : !FB else
2450 : : !FB write(ab_out,*)'One prod_scal equals zero'
2451 : : end if
2452 : : end do
2453 : : end do
2454 : :
2455 : : ! On stocke les vecteurs non-nuls
2456 : 527 : ii=0
2457 [ + + ]: 13677 : do kk=1,ncount
2458 [ + + ]: 784702 : prod_scal=sum( real(tab_vec(:,kk))* real(tab_vec(:,kk))+aimag(tab_vec(:,kk))*aimag(tab_vec(:,kk)))
2459 [ + + ]: 13677 : if (abs(prod_scal).gt.tol8) then
2460 : 6751 : ii=ii+1
2461 [ + + ]: 336772 : temp(:,ii)=tab_vec(:,kk)/dsqrt(prod_scal)
2462 : : end if
2463 : : end do
2464 : 527 : ncount=ii
2465 [ - + ]: 527 : ABI_FREE(tab_vec)
2466 : :
2467 : : ! On ecrit les vecteurs non-nuls
2468 : : !FB write(16,*) ' '
2469 : : !FB write(16,*) ' ========The final set of vectors is:'
2470 : : !FB do kk=1,ncount
2471 : : !FB write(16,'(81(f16.12,1x))') real(temp(:,kk))
2472 : : !FB write(16,'(81(f16.12,1x))') aimag(temp(:,kk))
2473 : : !FB end do
2474 [ + - ]: 527 : if (MPIdata%iam_master) then
2475 : 527 : write(16,'(a,1x,i7,1x,a)') ' ======= Finally, there are ',ncount,' independent vectors'
2476 : : end if
2477 [ + + - + ]: 527 : if (ncount.gt.8.and.order==2) then
2478 : 0 : ABI_ERROR(' Order 2 : There are too many independent vectors')
2479 : : end if
2480 [ + + - + ]: 527 : if (ncount.gt.27.and.order==3) then
2481 : 0 : ABI_ERROR(' Order 3 : There are too many independent vectors')
2482 : : end if
2483 [ - + - - ]: 527 : if (ncount.gt.81.and.order==4) then
2484 : 0 : ABI_ERROR(' Order 4 : There are too many independent vectors')
2485 : : end if
2486 : :
2487 : : ! On cherche les (norder-ncount) vecteurs orthogonaux aux vecteurs non-nuls
2488 : : ! --> Orthogonalisation de Gram-Schmidt
2489 [ - + + - : 407534 : ABI_MALLOC(tab_vec,(norder,norder)); tab_vec(:,:)=czero
- + + - +
- + - - +
- + + + +
+ ]
2490 : 527 : iseed=-5
2491 [ + + ]: 9746 : do kk=1,norder
2492 [ + + ]: 9746 : if (kk.le.ncount) then
2493 [ + + ]: 336772 : tab_vec(:,kk)=temp(:,kk)
2494 : : else
2495 [ + + ]: 68654 : do jj=1,norder
2496 : 66186 : drandom=uniformrandom(iseed)
2497 : 68654 : tab_vec(jj,kk)=dcmplx(drandom,zero)
2498 : : end do
2499 [ + + ]: 57116 : do jj=1,kk-1
2500 [ + + ]: 3453084 : prod_scal=sum( real(tab_vec(:,jj))* real(tab_vec(:,jj))+aimag(tab_vec(:,jj))*aimag(tab_vec(:,jj)))
2501 [ + - ]: 54648 : if (abs(prod_scal).gt.tol8) then
2502 [ + + + + ]: 6851520 : tab_vec(:,kk)=tab_vec(:,kk)-sum(tab_vec(:,kk)*conjg(tab_vec(:,jj)))/prod_scal*tab_vec(:,jj)
2503 [ + + ]: 3453084 : do ii=1,norder
2504 [ + + ]: 3398436 : if (abs( real(tab_vec(ii,kk))).lt.tol8) tab_vec(ii,kk)=dcmplx(zero,aimag(tab_vec(ii,kk)))
2505 [ + + ]: 3453084 : if (abs(aimag(tab_vec(ii,kk))).lt.tol8) tab_vec(ii,kk)=dcmplx( real(tab_vec(ii,kk)),zero)
2506 : : end do
2507 : : end if
2508 [ + + ]: 3453084 : prod_scal=sum( real(tab_vec(:,kk))* real(tab_vec(:,kk))+aimag(tab_vec(:,kk))*aimag(tab_vec(:,kk)))
2509 [ + + ]: 3455552 : tab_vec(:,kk)=tab_vec(:,kk)/dsqrt(prod_scal)
2510 : : end do
2511 : : end if
2512 : : end do
2513 [ - + ]: 527 : ABI_FREE(temp)
2514 : :
2515 : : ! On ecrit les vecteurs non-nuls
2516 : : !FB write(16,*) ' '
2517 : : !FB write(16,*) ' ========The orthogonal set of vectors is:'
2518 [ + + ]: 2995 : do kk=ncount+1,norder
2519 : : !FB write(16,'(81(f16.12,1x))') real(tab_vec(:,kk))
2520 : : !FB write(16,'(81(f16.12,1x))') aimag(tab_vec(:,kk))
2521 : : if ((abs(aimag(tab_vec(1,kk))).gt.tol6).or.&
2522 : 2468 : & (abs(aimag(tab_vec(1,kk))).gt.tol6).or.&
2523 [ - + ]: 2995 : & (abs(aimag(tab_vec(1,kk))).gt.tol6)) then
2524 : 0 : ABI_ERROR('the constraint has an imaginary part')
2525 : : end if
2526 : : end do
2527 : 527 : ncoeff=norder-ncount
2528 [ + - ]: 527 : if (MPIdata%iam_master) then
2529 : 527 : write(16,'(a,1x,i7,1x,a)') ' ======= Finally, there are ',ncoeff,' coefficients'
2530 : : end if
2531 : :
2532 : : ! On copie tab_vec dans proj
2533 [ + + ]: 2995 : do icoeff=1,ncoeff
2534 [ + + ]: 69181 : proj(:,icoeff,ishell)=tab_vec(:,ncount+icoeff)
2535 : : end do
2536 [ - + ]: 527 : ABI_FREE(tab_vec)
2537 : :
2538 [ - + ]: 693 : end subroutine tdep_calc_nbcoeff
2539 : : !!***
2540 : :
2541 : : !====================================================================================================
2542 : :
2543 [ # # # # : 0 : end module m_tdep_shell
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # #
# ]
2544 : : !!***
|
