Line data Source code
1 :
2 : #if defined HAVE_CONFIG_H
3 : #include "config.h"
4 : #endif
5 :
6 : #include "abi_common.h"
7 :
8 : module m_tdep_shell
9 :
10 : use defs_basis
11 : use m_errors
12 : use m_abicore
13 : use m_xmpi
14 : use m_io_tools
15 : use m_numeric_tools, only : uniformrandom
16 : use m_tdep_dataset, only : atdep_dataset_type, MPI_enreg_type
17 : use m_tdep_sym, only : Symmetries_type, tdep_SearchS_2at, tdep_SearchS_3at, tdep_SearchS_4at
18 : use m_tdep_sampling, only : tdep_Sampling_type
19 :
20 : type List_of_neighbours
21 : integer :: n_interactions
22 : integer, allocatable :: atomj_in_shell(:)
23 : integer, allocatable :: atomk_in_shell(:)
24 : integer, allocatable :: atoml_in_shell(:)
25 : integer, allocatable :: sym_in_shell(:)
26 : integer, allocatable :: transpose_in_shell(:)
27 : end type List_of_neighbours
28 :
29 : type Shell_type
30 :
31 : integer :: order
32 : ! Order of interaction (1, 2, 3, or 4)
33 :
34 : integer :: norder
35 : ! Dimension of IFC matrix at that order (3**order)
36 :
37 : integer :: natom
38 : ! Number of atoms
39 :
40 : integer :: ntotcoeff
41 : ! Total number of coefficients
42 :
43 : integer :: nshell
44 : ! Number of shells
45 :
46 : integer, allocatable :: ncoeff(:)
47 : ! ncoeff(nshell)
48 : ! Number of coefficients in each shell
49 :
50 : integer, allocatable :: ncoeff_prev(:)
51 : ! ncoeff_prev(nshell)
52 :
53 : integer, allocatable :: ishell_self(:)
54 : ! ncoeff_prev(natom_unitcell)
55 :
56 : integer, allocatable :: iatref(:)
57 : ! iatref(nshell)
58 :
59 : integer, allocatable :: jatref(:)
60 : ! jatref(nshell)
61 :
62 : integer, allocatable :: katref(:)
63 : ! katref(nshell)
64 :
65 : integer, allocatable :: latref(:)
66 : ! latref(nshell)
67 :
68 : double precision, allocatable :: proj(:,:,:)
69 : ! proj(norder, norder, nshell)
70 : ! Projector onto the subset of non-zero coefficients
71 :
72 : type(List_of_neighbours),allocatable :: neighbours(:,:)
73 : ! neighbours(natom, nshell)
74 :
75 : end type Shell_type
76 :
77 : public :: tdep_init_shell1at
78 : public :: tdep_init_shell2at
79 : public :: tdep_init_shell3at
80 : public :: tdep_init_shell4at
81 : public :: tdep_destroy_shell
82 : public :: tdep_calc_nbcoeff
83 :
84 : contains
85 :
86 : !====================================================================================================
87 44 : subroutine tdep_init_shell1at(Shell1at,Invar,MD,Sym,MPIdata)
88 :
89 : type(Shell_type),intent(out) :: Shell1at
90 : type(atdep_dataset_type),intent(in) :: Invar
91 : type(tdep_Sampling_type),intent(in) :: MD
92 : type(Symmetries_type),intent(inout) :: Sym
93 : type(MPI_enreg_type), intent(in) :: MPIdata
94 :
95 : integer :: ishell,iatcell,iatom,eatom,iatref,isym
96 : integer :: natom,natom_unitcell,counter,ncoeff,ncoeff_prev
97 : integer :: norder,order,nshell_max,nshell
98 44 : integer, allocatable :: ref1at(:,:),Isym1at(:,:)
99 :
100 44 : natom = Invar%natom
101 44 : natom_unitcell = Invar%natom_unitcell
102 44 : nshell_max = Invar%nshell_max
103 44 : order = 1
104 44 : norder = 3
105 44 : Shell1at%order = order
106 44 : Shell1at%norder = norder
107 44 : Shell1at%natom = natom
108 :
109 44 : write(Invar%stdout,*) ' '
110 44 : write(Invar%stdout,*) '#############################################################################'
111 44 : write(Invar%stdout,*) '####### FIRST ORDER : find the number of coefficients #######################'
112 44 : write(Invar%stdout,*) '#############################################################################'
113 :
114 :
115 : ! - Identify the shells
116 : ! - Store the index of the atoms included in each shell
117 : ! - Store the reference atoms for each shell
118 : ! - Compute the symetry operation between the reference atom and another one
119 44 : write(Invar%stdout,*) ' Build the ref1at and Isym1at tables...'
120 13984 : ABI_MALLOC(ref1at ,(natom,2)) ; ref1at (:,:)=zero
121 7058 : ABI_MALLOC(Isym1at,(natom,1)) ; Isym1at(:,:)=zero
122 44 : ishell=0
123 150 : do iatcell=1,natom_unitcell
124 106 : if (ref1at(iatcell,1).ne.0) cycle
125 63 : ishell=ishell+1
126 9405 : do eatom=1,natom
127 9404 : if (ref1at(eatom,1).eq.0) then
128 52172 : do isym=1,Sym%nsym
129 : !FB write(Invar%stdlog,'(4(i5,x))') Sym%indsym(4,isym,eatom),eatom,iatcell,isym
130 52172 : if (Sym%indsym(4,isym,eatom).eq.iatcell) then
131 6882 : Isym1at(eatom,1)=isym
132 6882 : ref1at(eatom,1)=iatcell
133 6882 : ref1at(eatom,2)=ishell
134 6882 : if (Invar%debug) write(Invar%stdout,'(a,1x,2(i4,1x),a,i4)') &
135 224 : & 'For:',iatcell,eatom,' direct transformation with isym=',Isym1at(eatom,1)
136 : exit
137 : end if
138 : end do !isym
139 : end if !already treated
140 : end do !eatom
141 : end do !iatcell
142 44 : Shell1at%nshell = ishell
143 44 : nshell = Shell1at%nshell
144 44 : if (nshell.gt.nshell_max) then
145 0 : write(Invar%stdout,*) ' STOP : The maximum number of shells allowed by the code is:',nshell_max
146 0 : write(Invar%stdout,*) ' In the present calculation, the number of shells is:',nshell
147 0 : write(Invar%stdout,*) ' Action: increase nshell_max'
148 0 : ABI_ERROR('The maximum number of shells allowed by the code is reached')
149 : end if
150 :
151 :
152 : ! Store all the previous quantities in a better way than in ref1at (without using too memory).
153 44 : write(Invar%stdout,*) ' Build the Shell1at datatype...'
154 258 : ABI_MALLOC(Shell1at%neighbours,(1,nshell))
155 195 : ABI_CALLOC(Shell1at%iatref, (nshell))
156 107 : do ishell=1,Shell1at%nshell
157 : counter=0
158 9361 : do iatom=1,natom
159 9361 : if (ref1at(iatom,2).eq.ishell) counter=counter+1
160 : end do
161 63 : Shell1at%neighbours(1,ishell)%n_interactions=counter
162 63 : if (counter.eq.0) then
163 : cycle
164 : end if
165 7071 : ABI_CALLOC(Shell1at%neighbours(1,ishell)%atomj_in_shell,(counter))
166 7008 : ABI_CALLOC(Shell1at%neighbours(1,ishell)%sym_in_shell,(counter))
167 : counter=0
168 9405 : do iatom=1,natom
169 9361 : if (ref1at(iatom,2).eq.ishell) then
170 6882 : counter=counter+1
171 6882 : Shell1at%neighbours(1,ishell)%atomj_in_shell(counter)=iatom
172 6882 : Shell1at%iatref(ishell)=ref1at(iatom,1)
173 6882 : Shell1at%neighbours(1,ishell)%sym_in_shell(counter)=Isym1at(iatom,1)
174 : end if
175 : end do
176 : end do
177 44 : ABI_FREE(ref1at)
178 44 : ABI_FREE(Isym1at)
179 :
180 : ! Find the number of coefficients of the (3x3) Phi2 for a given shell
181 951 : ABI_CALLOC(Shell1at%proj, (norder,norder,nshell))
182 151 : ABI_CALLOC(Shell1at%ncoeff, (nshell))
183 151 : ABI_CALLOC(Shell1at%ncoeff_prev, (nshell))
184 44 : write(Invar%stdout,*) ' Number of shells=',nshell
185 44 : write(Invar%stdout,*) '============================================================================'
186 44 : if (MPIdata%iam_master) open(unit=16,file=trim(Invar%output_prefix)//'_nbcoeff-phi1.dat')
187 44 : ncoeff_prev=0
188 107 : do ishell=1,nshell
189 63 : ncoeff=0
190 63 : iatref=Shell1at%iatref(ishell)
191 63 : write(Invar%stdout,*) 'Shell number:',ishell
192 63 : write(Invar%stdout,'(a,i5,a)') ' For atom',iatref,':'
193 63 : call tdep_calc_nbcoeff(MD%distance,iatref,Invar,ishell,1,1,1,MPIdata,ncoeff,norder,nshell,order,Shell1at%proj,Sym)
194 63 : ncoeff=0
195 : if (ncoeff.eq.0) then
196 63 : Shell1at%neighbours(1,ishell)%n_interactions=0
197 63 : if(allocated(Shell1at%neighbours(1,ishell)%atomj_in_shell)) then
198 63 : ABI_FREE(Shell1at%neighbours(1,ishell)%atomj_in_shell)
199 : end if
200 63 : if(allocated(Shell1at%neighbours(1,ishell)%sym_in_shell)) then
201 63 : ABI_FREE(Shell1at%neighbours(1,ishell)%sym_in_shell)
202 : end if
203 : end if
204 63 : Shell1at%ncoeff (ishell)=ncoeff
205 63 : Shell1at%ncoeff_prev(ishell)=ncoeff_prev
206 : ncoeff_prev=ncoeff_prev+ncoeff
207 63 : write(Invar%stdout,*)' Number of independant coefficients in this shell=',ncoeff
208 63 : write(Invar%stdout,*)' Number of interactions in this shell=',Shell1at%neighbours(1,ishell)%n_interactions
209 : !FB write(Invar%stdout,*)' The ratio is=',dfloat(Shell1at%neighbours(iatref,ishell)%n_interactions)/dfloat(ncoeff)
210 107 : write(Invar%stdout,*) '============================================================================'
211 : end do
212 44 : write(Invar%stdout,*)' >>>>>> Total number of coefficients at the first order=',ncoeff_prev
213 44 : if (MPIdata%iam_master) close(16)
214 44 : Shell1at%ntotcoeff=ncoeff_prev
215 :
216 44 : end subroutine tdep_init_shell1at
217 :
218 : !====================================================================================================
219 44 : subroutine tdep_init_shell2at(Shell2at,Invar,MD,Sym,MPIdata)
220 :
221 : type(Shell_type),intent(out) :: Shell2at
222 : type(atdep_dataset_type),intent(in) :: Invar
223 : type(tdep_Sampling_type),intent(in) :: MD
224 : type(Symmetries_type),intent(inout) :: Sym
225 : type(MPI_enreg_type), intent(in) :: MPIdata
226 :
227 : integer :: ishell,iatcell,iatom,jatom,eatom,fatom,iatref,jatref
228 : integer :: natom,natom_unitcell,counter,ncoeff,ncoeff_prev
229 : integer :: norder,order,nshell_max,nshell
230 44 : integer, allocatable :: ref2at(:,:,:),Isym2at(:,:,:)
231 :
232 44 : natom = Invar%natom
233 44 : natom_unitcell = Invar%natom_unitcell
234 44 : nshell_max = Invar%nshell_max
235 44 : order = 2
236 44 : norder = 9
237 44 : Shell2at%order = order
238 44 : Shell2at%norder = norder
239 44 : Shell2at%natom = natom
240 :
241 44 : write(Invar%stdout,*) ' '
242 44 : write(Invar%stdout,*) '#############################################################################'
243 44 : write(Invar%stdout,*) '###### SECOND ORDER : find the number of coefficients #######################'
244 44 : write(Invar%stdout,*) '#############################################################################'
245 :
246 : ! - Identify the shells
247 : ! - Store the index of the atoms included in each shell
248 : ! - Store the reference atoms for each shell
249 : ! - Compute the symetry operation between the reference atom and another one
250 44 : write(Invar%stdout,*) ' Build the ref2at and Isym2at tables...'
251 4437922 : ABI_MALLOC(ref2at ,(natom,natom,3)) ; ref2at (:,:,:)=zero
252 2958688 : ABI_MALLOC(Isym2at,(natom,natom,2)) ; Isym2at(:,:,:)=zero
253 238 : ABI_MALLOC(Shell2at%ishell_self,(natom_unitcell)) ; Shell2at%ishell_self(:)=zero
254 44 : ishell=0
255 150 : do iatcell=1,natom_unitcell
256 15036 : do jatom=1,natom
257 : ! Interactions are only computed until Rcut in order to have complete shell of neighbours.
258 : ! Otherwise the symetries are broken.
259 14886 : if ((ref2at(iatcell,jatom,1).ne.0).or.(MD%distance(iatcell,jatom,1).gt.(Invar%rcut*0.99))) cycle
260 500 : ishell=ishell+1
261 500 : if (iatcell.eq.jatom) Shell2at%ishell_self(iatcell)=ishell
262 79180 : do eatom=1,natom
263 17985552 : do fatom=1,natom
264 17892092 : if ((ref2at(eatom,fatom,1).eq.0).and.&
265 : !FB& (abs(MD%distance(iatcell,jatom,1)-MD%distance(eatom,fatom,1)).lt.1.d-3)) then
266 78574 : & (abs(MD%distance(iatcell,jatom,1)-MD%distance(eatom,fatom,1)).lt.tol6)) then
267 306801 : call tdep_SearchS_2at(Invar,iatcell,jatom,eatom,fatom,Isym2at,Sym,MD%xred_ideal)
268 306801 : if (Isym2at(eatom,fatom,2)==1) then
269 253104 : if (Invar%debug) write(Invar%stdout,'(a,1x,4(i4,1x),a,i4)') &
270 5440 : & 'For:',iatcell,jatom,eatom,fatom,' direct transformation with isym=',Isym2at(eatom,fatom,1)
271 253104 : ref2at(eatom,fatom,1)=iatcell
272 253104 : ref2at(eatom,fatom,2)=jatom
273 253104 : ref2at(eatom,fatom,3)=ishell
274 : ! The Phi2 has to be symetric (transposition symetries)
275 253104 : if (Invar%debug) write(Invar%stdout,'(a,1x,4(i4,1x),a,i4)') &
276 5440 : & 'For:',iatcell,jatom,eatom,fatom,' transformation+permutation with isym=',Isym2at(eatom,fatom,1)
277 253104 : ref2at(fatom,eatom,1)=iatcell
278 253104 : ref2at(fatom,eatom,2)=jatom
279 253104 : ref2at(fatom,eatom,3)=ishell
280 253104 : Isym2at(fatom,eatom,1)=Isym2at(eatom,fatom,1)
281 253104 : Isym2at(fatom,eatom,2)=2
282 : else
283 53697 : if (Invar%debug) write(Invar%stdout,'(a,4(1x,i4))') &
284 3640 : & 'NO SYMETRY OPERATION BETWEEN (iatom,jatom) and (eatom,fatom)=',iatcell,jatom,eatom,fatom
285 : end if
286 : end if !already treated
287 : end do !fatom
288 : end do !eatom
289 : end do !jatom
290 : end do !iatcell
291 44 : Shell2at%nshell = ishell
292 44 : nshell = Shell2at%nshell
293 44 : if (nshell.gt.nshell_max) then
294 0 : write(Invar%stdout,*) ' STOP : The maximum number of shells allowed by the code is:',nshell_max
295 0 : write(Invar%stdout,*) ' In the present calculation, the number of shells is:',nshell
296 0 : write(Invar%stdout,*) ' Action: increase nshell_max'
297 0 : ABI_ERROR('The maximum number of shells allowed by the code is reached')
298 : end if
299 :
300 :
301 : ! Store all the previous quantities in a better way than in ref2at (without using too memory).
302 44 : write(Invar%stdout,*) ' Build the Shell2at datatype...'
303 79250 : ABI_MALLOC(Shell2at%neighbours,(natom,nshell))
304 632 : ABI_CALLOC(Shell2at%iatref, (nshell))
305 588 : ABI_CALLOC(Shell2at%jatref, (nshell))
306 544 : do ishell=1,nshell
307 79118 : do iatom=1,natom
308 78574 : counter=0
309 17970666 : do jatom=1,natom
310 17970666 : if (ref2at(iatom,jatom,3).eq.ishell) counter=counter+1
311 : end do
312 78574 : Shell2at%neighbours(iatom,ishell)%n_interactions=counter
313 78574 : if (counter.eq.0) cycle
314 201162 : ABI_MALLOC(Shell2at%neighbours(iatom,ishell)%atomj_in_shell,(counter))
315 134108 : ABI_MALLOC(Shell2at%neighbours(iatom,ishell)%sym_in_shell,(counter))
316 134108 : ABI_MALLOC(Shell2at%neighbours(iatom,ishell)%transpose_in_shell,(counter))
317 566380 : Shell2at%neighbours(iatom,ishell)%atomj_in_shell(:)=zero
318 566380 : Shell2at%neighbours(iatom,ishell)%sym_in_shell(:)=zero
319 566380 : Shell2at%neighbours(iatom,ishell)%transpose_in_shell(:)=zero
320 67054 : counter=0
321 15996830 : do jatom=1,natom
322 16007850 : if (ref2at(iatom,jatom,3).eq.ishell) then
323 499326 : counter=counter+1
324 499326 : Shell2at%neighbours(iatom,ishell)%atomj_in_shell(counter)=jatom
325 499326 : Shell2at%iatref (ishell)=ref2at(iatom,jatom,1)
326 499326 : Shell2at%jatref (ishell)=ref2at(iatom,jatom,2)
327 499326 : Shell2at%neighbours(iatom,ishell)%sym_in_shell (counter)=Isym2at(iatom,jatom,1)
328 499326 : Shell2at%neighbours(iatom,ishell)%transpose_in_shell(counter)=Isym2at(iatom,jatom,2)
329 : end if
330 : end do
331 : end do
332 : end do
333 44 : ABI_FREE(ref2at)
334 44 : ABI_FREE(Isym2at)
335 :
336 : ! Find the number of coefficients of the (3x3) Phi2 for a given shell
337 45632 : ABI_CALLOC(Shell2at%proj, (norder,norder,nshell))
338 588 : ABI_CALLOC(Shell2at%ncoeff, (nshell))
339 588 : ABI_CALLOC(Shell2at%ncoeff_prev, (nshell))
340 44 : write(Invar%stdout,*) ' Number of shells=',nshell
341 44 : write(Invar%stdout,*) '============================================================================'
342 44 : if (MPIdata%iam_master) open(unit=16,file=trim(Invar%output_prefix)//'_nbcoeff-phi2.dat')
343 44 : ncoeff_prev=0
344 544 : do ishell=1,nshell
345 500 : ncoeff=0
346 500 : iatref=Shell2at%iatref(ishell)
347 500 : jatref=Shell2at%jatref(ishell)
348 500 : write(Invar%stdout,*) 'Shell number:',ishell
349 500 : write(Invar%stdout,'(a,i5,a,i5,a,f16.10)') ' Between atom',iatref,' and ',jatref,' the distance is=',MD%distance(iatref,jatref,1)
350 500 : call tdep_calc_nbcoeff(MD%distance,iatref,Invar,ishell,jatref,1,1,MPIdata,ncoeff,norder,nshell,order,Shell2at%proj,Sym)
351 500 : Shell2at%ncoeff (ishell)=ncoeff
352 500 : Shell2at%ncoeff_prev(ishell)=ncoeff_prev
353 500 : ncoeff_prev=ncoeff_prev+ncoeff
354 500 : write(Invar%stdout,*)' Number of independant coefficients in this shell=',ncoeff
355 500 : write(Invar%stdout,*)' Number of interactions in this shell=',Shell2at%neighbours(iatref,ishell)%n_interactions
356 : !FB write(Invar%stdout,*)' The ratio is=',dfloat(Shell2at%neighbours(iatref,ishell)%n_interactions)/dfloat(ncoeff)
357 544 : write(Invar%stdout,*) '============================================================================'
358 : end do
359 44 : write(Invar%stdout,*)' >>>>>> Total number of coefficients at the second order=',ncoeff_prev
360 44 : if (MPIdata%iam_master) close(16)
361 44 : Shell2at%ntotcoeff=ncoeff_prev
362 : !BeginFB
363 : !FB open(unit=91,file='Shell2at.dat')
364 : !FB write(91,*) Shell2at%nshell
365 : !FB do ishell=1,Shell2at%nshell
366 : !FB write(91,*) Shell2at%ncoeff(ishell)
367 : !FB write(91,*) Shell2at%ncoeff_prev(ishell)
368 : !FB write(91,*) Shell2at%iatref(ishell)
369 : !FB write(91,*) Shell2at%jatref(ishell)
370 : !FB do iatom=1,Invar%natom
371 : !FB write(91,*) Shell2at%neighbours(iatom,ishell)%n_interactions
372 : !FB do ii=1,Shell2at%neighbours(iatom,ishell)%n_interactions
373 : !FB write(91,*) Shell2at%neighbours(iatom,ishell)%sym_in_shell(ii)
374 : !FB write(91,*) Shell2at%neighbours(iatom,ishell)%transpose_in_shell(ii)
375 : !FB write(91,*) Shell2at%neighbours(iatom,ishell)%atomj_in_shell(ii)
376 : !FB end do
377 : !FB end do
378 : !FB end do
379 : !FB close(91)
380 : !EndFB
381 :
382 44 : end subroutine tdep_init_shell2at
383 :
384 : !====================================================================================================
385 14 : subroutine tdep_init_shell3at(Shell3at,Invar,MD,Sym,MPIdata)
386 :
387 : type(Shell_type),intent(out) :: Shell3at
388 : type(atdep_dataset_type),intent(in) :: Invar
389 : type(tdep_Sampling_type),intent(in) :: MD
390 : type(Symmetries_type),intent(inout) :: Sym
391 : type(MPI_enreg_type), intent(in) :: MPIdata
392 :
393 : integer :: ii,ishell,iatom,jatom,katom,eatom,fatom,gatom,iatref,jatref,katref
394 : integer :: natom,natom_unitcell,watom,xatom,yatom,ninteractions,ncoeff,ncoeff_prev,nshell_tmp
395 : integer :: find_equivalent,ninter,iat_ref,jat_ref,kat_ref,tmpinter
396 : integer :: norder,order,nshell_max,nshell
397 : double precision :: norma,normb,normc
398 : integer :: Isym3at(2)
399 14 : integer, allocatable :: atref(:,:),interactions(:,:)
400 :
401 14 : natom = Invar%natom
402 14 : natom_unitcell = Invar%natom_unitcell
403 14 : nshell_max = Invar%nshell_max
404 14 : order = 3
405 14 : norder = 27
406 14 : Shell3at%order = order
407 14 : Shell3at%norder = norder
408 14 : Shell3at%natom = natom
409 :
410 14 : write(Invar%stdout,*) ' '
411 14 : write(Invar%stdout,*) '#############################################################################'
412 14 : write(Invar%stdout,*) '###### THIRD ORDER : find the number of coefficients ########################'
413 14 : write(Invar%stdout,*) '#############################################################################'
414 :
415 : ! 1/ Identify the shells
416 21056 : ABI_CALLOC(interactions,(natom_unitcell,nshell_max))
417 21084 : ABI_CALLOC(atref,(nshell_max,3))
418 70 : ABI_CALLOC(Shell3at%ishell_self,(natom_unitcell))
419 14 : nshell_tmp=0
420 42 : do iatom=1,natom_unitcell
421 4930 : do jatom=1,natom
422 905428 : do katom=1,natom
423 : !FB write(Invar%stdlog,*) 'NEW COORD1 : iatom,jatom,katom=',iatom,jatom,katom
424 : ! WARNING: distance(j,k).ne.|djk| due to the inbox procedure when computing distance(j,k).
425 : ! So, compute |djk| using vec(ij) and vec(ik).
426 : norma=dsqrt((MD%distance(iatom,katom,2)-MD%distance(iatom,jatom,2))**2+&
427 : & (MD%distance(iatom,katom,3)-MD%distance(iatom,jatom,3))**2+&
428 900512 : & (MD%distance(iatom,katom,4)-MD%distance(iatom,jatom,4))**2)
429 : ! Interactions are only computed until Rcut3 in order to have complete shell of neighbours.
430 : ! Otherwise the symetries are broken.
431 : if ((MD%distance(iatom,jatom,1).gt.(Invar%rcut3*0.99)).or.&
432 900512 : & (norma .gt.(Invar%rcut3*0.99)).or.&
433 : & (MD%distance(iatom,katom,1).gt.(Invar%rcut3*0.99))) cycle
434 6032 : if (nshell_tmp.eq.0) then
435 56 : atref(1,:)=1
436 14 : nshell_tmp=nshell_tmp+1
437 14 : if ((iatom.eq.jatom).and.(jatom.eq.katom)) Shell3at%ishell_self(iatom)=nshell_tmp
438 14 : interactions(iatom,nshell_tmp)=interactions(iatom,nshell_tmp)+1
439 14 : cycle
440 : else
441 1130 : find_equivalent=0
442 4838 : do ishell=1,nshell_tmp
443 4803 : iat_ref=atref(ishell,1) ; jat_ref=atref(ishell,2) ; kat_ref=atref(ishell,3)
444 : normb=dsqrt((MD%distance(iat_ref,kat_ref,2)-MD%distance(iat_ref,jat_ref,2))**2+&
445 : & (MD%distance(iat_ref,kat_ref,3)-MD%distance(iat_ref,jat_ref,3))**2+&
446 4803 : & (MD%distance(iat_ref,kat_ref,4)-MD%distance(iat_ref,jat_ref,4))**2)
447 28755 : do ii=1,6
448 25047 : if (ii.eq.1) then ; eatom=iatom ; fatom=jatom ; gatom=katom ; end if
449 25047 : if (ii.eq.2) then ; eatom=iatom ; fatom=katom ; gatom=jatom ; end if
450 25047 : if (ii.eq.3) then ; eatom=jatom ; fatom=iatom ; gatom=katom ; end if
451 25047 : if (ii.eq.4) then ; eatom=jatom ; fatom=katom ; gatom=iatom ; end if
452 25047 : if (ii.eq.5) then ; eatom=katom ; fatom=iatom ; gatom=jatom ; end if
453 25047 : if (ii.eq.6) then ; eatom=katom ; fatom=jatom ; gatom=iatom ; end if
454 : normc=dsqrt((MD%distance(eatom,gatom,2)-MD%distance(eatom,fatom,2))**2+&
455 : & (MD%distance(eatom,gatom,3)-MD%distance(eatom,fatom,3))**2+&
456 25047 : & (MD%distance(eatom,gatom,4)-MD%distance(eatom,fatom,4))**2)
457 : !FB if ((abs(MD%distance(iatom,jatom,1)-MD%distance(eatom,fatom,1)).lt.1.d-3).and.&
458 : !FB& (abs(norma -normb ).lt.1.d-3).and.&
459 : !FB& (abs(MD%distance(iatom,katom,1)-MD%distance(eatom,gatom,1)).lt.1.d-3)) then
460 : if ((abs(MD%distance(iat_ref,jat_ref,1)-MD%distance(eatom,fatom,1)).lt.1.d-6).and.&
461 25047 : & (abs(normb -normc ).lt.1.d-6).and.&
462 3708 : & (abs(MD%distance(iat_ref,kat_ref,1)-MD%distance(eatom,gatom,1)).lt.1.d-6)) then
463 1815 : Isym3at(:)=0
464 1815 : call tdep_SearchS_3at(Invar,iat_ref,jat_ref,kat_ref,eatom,fatom,gatom,Isym3at,Sym,MD%xred_ideal)
465 1815 : if (Isym3at(2).eq.1) find_equivalent=1
466 : if (find_equivalent.eq.1) then
467 1095 : interactions(iatom,ishell)=interactions(iatom,ishell)+1
468 : !FB write(Invar%stdlog,*) 'The number of interactions in this shell is=',ishell,interactions(iatom,ishell)
469 1095 : exit
470 : end if
471 : end if
472 : end do !ii
473 4838 : if (find_equivalent.eq.1) exit
474 : end do !ishell
475 1130 : if (find_equivalent.eq.0) then
476 35 : nshell_tmp=nshell_tmp+1
477 35 : if (nshell_tmp.gt.nshell_max) then
478 0 : ABI_ERROR('The shell number index is greater than the shell number max defined in the code')
479 : end if
480 35 : if ((iatom.eq.jatom).and.(jatom.eq.katom)) Shell3at%ishell_self(iatom)=nshell_tmp
481 35 : interactions(iatom,nshell_tmp)=interactions(iatom,nshell_tmp)+1
482 35 : atref(nshell_tmp,1)=iatom
483 35 : atref(nshell_tmp,2)=jatom
484 35 : atref(nshell_tmp,3)=katom
485 : !FB write(Invar%stdlog,'(a,1x,4(i5,1x))') 'NEW SHELL1 : nshell_tmp,iatom,jatom,katom=',nshell_tmp,iatom,jatom,katom
486 : end if
487 : end if
488 : end do !katom
489 : end do !jatom
490 : end do !iatom
491 14 : ABI_FREE(atref)
492 :
493 : ! 2/ Allocate the datatype Shell3at%...
494 14 : Shell3at%nshell=nshell_tmp
495 14 : nshell = Shell3at%nshell
496 8433 : ABI_MALLOC(Shell3at%neighbours,(natom,nshell))
497 91 : ABI_MALLOC(Shell3at%iatref,(nshell)); Shell3at%iatref(:)=zero
498 77 : ABI_MALLOC(Shell3at%jatref,(nshell)); Shell3at%jatref(:)=zero
499 77 : ABI_MALLOC(Shell3at%katref,(nshell)); Shell3at%katref(:)=zero
500 63 : do ishell=1,nshell
501 8391 : do iatom=1,natom
502 8328 : ninteractions=interactions(mod(iatom-1,natom_unitcell)+1,ishell)
503 8328 : Shell3at%neighbours(iatom,ishell)%n_interactions=ninteractions
504 8328 : if (ninteractions.eq.0) cycle
505 21144 : ABI_MALLOC(Shell3at%neighbours(iatom,ishell)%atomj_in_shell,(ninteractions))
506 14096 : ABI_MALLOC(Shell3at%neighbours(iatom,ishell)%atomk_in_shell,(ninteractions))
507 14096 : ABI_MALLOC(Shell3at%neighbours(iatom,ishell)%sym_in_shell,(ninteractions))
508 14096 : ABI_MALLOC(Shell3at%neighbours(iatom,ishell)%transpose_in_shell,(ninteractions))
509 86436 : Shell3at%neighbours(iatom,ishell)%atomj_in_shell(:)=zero
510 86436 : Shell3at%neighbours(iatom,ishell)%atomk_in_shell(:)=zero
511 86436 : Shell3at%neighbours(iatom,ishell)%sym_in_shell(:)=zero
512 87765 : Shell3at%neighbours(iatom,ishell)%transpose_in_shell(:)=zero
513 : end do
514 : end do
515 14 : ABI_FREE(interactions)
516 :
517 : ! 3/ Store the index of the (couple of) atoms included in each shell
518 : ! 4/ Store the reference (couple of) atoms for each shell
519 : ! 5/ Compute the symetry operation between the reference (couple of) atoms and another one
520 8433 : ABI_MALLOC(interactions,(natom,nshell)) ; interactions(:,:)=0
521 14 : nshell_tmp=0
522 2554 : do iatom=1,natom
523 495818 : do jatom=1,natom
524 102331388 : do katom=1,natom
525 : !FB write(Invar%stdlog,*) 'NEW COORD2 : iatom,jatom,katom=',iatom,jatom,katom
526 101835584 : if (nshell_tmp.eq.0) then
527 14 : nshell_tmp=nshell_tmp+1
528 14 : interactions(iatom,nshell_tmp)=1
529 14 : Shell3at%iatref(nshell_tmp)=iatom
530 14 : Shell3at%jatref(nshell_tmp)=jatom
531 14 : Shell3at%katref(nshell_tmp)=katom
532 14 : Shell3at%neighbours(iatom,nshell_tmp)%atomj_in_shell(interactions(iatom,nshell_tmp))=jatom
533 14 : Shell3at%neighbours(iatom,nshell_tmp)%atomk_in_shell(interactions(iatom,nshell_tmp))=katom
534 14 : Shell3at%neighbours(iatom,nshell_tmp)%sym_in_shell(interactions(iatom,nshell_tmp))=1
535 14 : Shell3at%neighbours(iatom,nshell_tmp)%transpose_in_shell(interactions(iatom,nshell_tmp))=1
536 14 : cycle
537 : end if
538 : ! WARNING: MD%distance(j,k).ne.|djk| due to the inbox procedure when computing MD%distance(j,k).
539 : ! So, compute |djk| using vec(ij) and vec(ik).
540 : norma=dsqrt((MD%distance(iatom,katom,2)-MD%distance(iatom,jatom,2))**2+&
541 : & (MD%distance(iatom,katom,3)-MD%distance(iatom,jatom,3))**2+&
542 101835570 : & (MD%distance(iatom,katom,4)-MD%distance(iatom,jatom,4))**2)
543 : ! Interactions are only computed until Rcut3<acell/2 in order to have complete shell of neighbours.
544 : ! Otherwise the symetries are broken.
545 : if ((MD%distance(iatom,jatom,1).gt.(Invar%rcut3*0.99)).or.&
546 101835570 : & (norma .gt.(Invar%rcut3*0.99)).or.&
547 : & (MD%distance(iatom,katom,1).gt.(Invar%rcut3*0.99))) cycle
548 : ! Search if the triplet has already been classified
549 79374 : find_equivalent=0
550 595656 : do ishell=1,nshell_tmp
551 5404412 : do ninter=1,interactions(iatom,ishell)
552 4808756 : if ((Shell3at%neighbours(iatom,ishell)%atomj_in_shell(ninter).eq.jatom).and.&
553 5325038 : & (Shell3at%neighbours(iatom,ishell)%atomk_in_shell(ninter).eq.katom)) find_equivalent=1
554 : end do
555 : end do
556 79374 : if (find_equivalent.eq.1) cycle
557 : ! Search if the triplet belongs to a shell already found
558 80466 : do ishell=1,nshell_tmp
559 80431 : iat_ref=Shell3at%iatref(ishell) ; jat_ref=Shell3at%jatref(ishell) ; kat_ref=Shell3at%katref(ishell)
560 : normb=dsqrt((MD%distance(iat_ref,kat_ref,2)-MD%distance(iat_ref,jat_ref,2))**2+&
561 : & (MD%distance(iat_ref,kat_ref,3)-MD%distance(iat_ref,jat_ref,3))**2+&
562 80431 : & (MD%distance(iat_ref,kat_ref,4)-MD%distance(iat_ref,jat_ref,4))**2)
563 464491 : do ii=1,6
564 405351 : if (ii.eq.1) then ; eatom=iatom ; fatom=jatom ; gatom=katom ; end if
565 405351 : if (ii.eq.2) then ; eatom=iatom ; fatom=katom ; gatom=jatom ; end if
566 405351 : if (ii.eq.3) then ; eatom=jatom ; fatom=iatom ; gatom=katom ; end if
567 405351 : if (ii.eq.4) then ; eatom=jatom ; fatom=katom ; gatom=iatom ; end if
568 405351 : if (ii.eq.5) then ; eatom=katom ; fatom=iatom ; gatom=jatom ; end if
569 405351 : if (ii.eq.6) then ; eatom=katom ; fatom=jatom ; gatom=iatom ; end if
570 : normc=dsqrt((MD%distance(eatom,gatom,2)-MD%distance(eatom,fatom,2))**2+&
571 : & (MD%distance(eatom,gatom,3)-MD%distance(eatom,fatom,3))**2+&
572 405351 : & (MD%distance(eatom,gatom,4)-MD%distance(eatom,fatom,4))**2)
573 : !FB if ((abs(MD%distance(iatom,jatom,1)-MD%distance(eatom,fatom,1)).lt.1.d-3).and.&
574 : !FB& (abs(norma -normb ).lt.1.d-3).and.&
575 : !FB& (abs(MD%distance(iatom,katom,1)-MD%distance(eatom,gatom,1)).lt.1.d-3)) then
576 : if ((abs(MD%distance(iat_ref,jat_ref,1)-MD%distance(eatom,fatom,1)).lt.1.d-6).and.&
577 405351 : & (abs(normb -normc ).lt.1.d-6).and.&
578 59140 : & (abs(MD%distance(iat_ref,kat_ref,1)-MD%distance(eatom,gatom,1)).lt.1.d-6)) then
579 35171 : Isym3at(:)=0
580 35171 : call tdep_SearchS_3at(Invar,iat_ref,jat_ref,kat_ref,eatom,fatom,gatom,Isym3at,Sym,MD%xred_ideal)
581 35171 : if (Isym3at(2).eq.1) then
582 : find_equivalent=1
583 : exit
584 : end if
585 : end if
586 : end do !ii
587 80466 : if (find_equivalent.eq.1) exit
588 : end do !ishell
589 : ! The triplet belongs to a new shell
590 21326 : if (find_equivalent.eq.0) then
591 35 : nshell_tmp=nshell_tmp+1
592 : ! Check that the new shell is allowed
593 35 : if (nshell_tmp.gt.nshell) then
594 0 : ABI_ERROR('The shell number index is greater than the shell number max computed previously')
595 : end if
596 35 : Shell3at%iatref(nshell_tmp)=iatom
597 35 : Shell3at%jatref(nshell_tmp)=jatom
598 35 : Shell3at%katref(nshell_tmp)=katom
599 35 : eatom=iatom ; fatom=jatom ; gatom=katom
600 105 : Isym3at(:)=1
601 35 : ishell=nshell_tmp
602 : !FB write(Invar%stdlog,'(a,1x,4(i5,1x))') 'NEW SHELL2 : nshell_tmp,iatom,jatom,katom=',nshell_tmp,iatom,jatom,katom
603 : end if
604 : ! Classify the informations of the triplet in Shell3at
605 642546 : do ii=1,6
606 : ! The Phi3 has to be symetric (transposition symetries)
607 127956 : if (ii==1) then ; watom=eatom ; xatom=fatom ; yatom=gatom ; endif !\Phi3_ijk
608 127956 : if (ii==2) then ; watom=eatom ; xatom=gatom ; yatom=fatom ; endif !\Phi3_ikj
609 127956 : if (ii==3) then ; watom=fatom ; xatom=eatom ; yatom=gatom ; endif !\Phi3_jik
610 127956 : if (ii==4) then ; watom=fatom ; xatom=gatom ; yatom=eatom ; endif !\Phi3_jki
611 127956 : if (ii==5) then ; watom=gatom ; xatom=eatom ; yatom=fatom ; endif !\Phi3_kij
612 127956 : if (ii==6) then ; watom=gatom ; xatom=fatom ; yatom=eatom ; endif !\Phi3_kji
613 : ! Do not overwrite the Phi3_iik, Phi3_iji, Phi3_ijj or Phi3_iii IFCs
614 : ! and avoid double counting of triplet interactions
615 127956 : if ((eatom.eq.fatom).and.((ii.eq.3).or.(ii.eq.4).or.(ii.eq.6))) cycle
616 92106 : if ((eatom.eq.gatom).and.((ii.gt.3))) cycle
617 89580 : if ((fatom.eq.gatom).and.((ii.eq.2).or.(ii.eq.5).or.(ii.eq.6))) cycle
618 79374 : if ((eatom.eq.fatom).and.(fatom.eq.gatom).and.(ii.gt.1)) cycle
619 79374 : interactions(watom,ishell)=interactions(watom,ishell)+1
620 : !FB write(Invar%stdlog,*) 'For ishell and eatom=',ishell,watom
621 : !FB write(Invar%stdlog,*) ' --> the number of interactions in the shell is=',interactions(watom,ishell)
622 79374 : if (interactions(watom,ishell).gt.Shell3at%neighbours(watom,ishell)%n_interactions) then
623 0 : write(Invar%stdlog,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>'
624 0 : write(Invar%stdlog,*) ' >>>>>> Verify that the Rcut used in the input file is lower '
625 0 : write(Invar%stdlog,*) ' >>>>>> than half of the smallest lattice parameter'
626 0 : write(Invar%stdlog,*) ' >>>>>> Solution : Reduce the Rcut parameter'
627 0 : write(Invar%stdlog,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>'
628 0 : ABI_ERROR('The interaction number index is greater than the interaction number max computed previously (3rd order)')
629 : end if
630 79374 : Shell3at%neighbours(watom,ishell)%atomj_in_shell(interactions(watom,ishell))=xatom
631 79374 : Shell3at%neighbours(watom,ishell)%atomk_in_shell(interactions(watom,ishell))=yatom
632 79374 : Shell3at%neighbours(watom,ishell)%sym_in_shell(interactions(watom,ishell))=Isym3at(1)
633 101953320 : Shell3at%neighbours(watom,ishell)%transpose_in_shell(interactions(watom,ishell))=ii
634 : !DEBUG write(Invar%stdlog,'(a,9(i5,x))') 'ishell,iatref,jatref,katref,iatom,atomj_in_shell,atomk_in_shell,isym,itrans=',&
635 : !DEBUG& ishell,Shell3at%iatref(ishell),Shell3at%jatref(ishell),Shell3at%katref(ishell),watom,xatom,yatom,Isym3at(1),ii
636 : end do !ii
637 : end do !katom
638 : end do !jatom
639 : end do !iatom
640 : ! Check that each interaction has different symmetry per shell
641 63 : do ishell=1,nshell
642 8391 : do iatom=1,natom
643 8328 : if (Shell3at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
644 79437 : do ninter=1,Shell3at%neighbours(iatom,ishell)%n_interactions-1
645 964484 : do tmpinter=ninter+1,Shell3at%neighbours(iatom,ishell)%n_interactions
646 883816 : if (Shell3at%neighbours(iatom,ishell)%sym_in_shell( ninter).eq.&
647 72340 : & Shell3at%neighbours(iatom,ishell)%sym_in_shell(tmpinter)) then
648 48896 : if (Shell3at%neighbours(iatom,ishell)%transpose_in_shell( ninter).ne.&
649 : & Shell3at%neighbours(iatom,ishell)%transpose_in_shell(tmpinter)) cycle
650 0 : write(std_out,'(a,2(1x,i5))') 'For ishell and iatom =',ishell,iatom
651 0 : write(std_out,'(a,i5,a,i5,a,i5)') ' the interactions ',ninter,&
652 0 : & ' and ',tmpinter,' have both the same symmetry isym=',Shell3at%neighbours(iatom,ishell)%sym_in_shell( ninter)
653 0 : ABI_ERROR('Some interactions are equals due to the symmetry')
654 : end if
655 : end do
656 : end do
657 : end do
658 : end do
659 : ! Check that each equivalent shell has the same set of interactions
660 63 : do ishell=1,nshell
661 8391 : do iatom=1,natom
662 8328 : if (Shell3at%neighbours(mod(iatom-1,natom_unitcell)+1,ishell)%n_interactions.ne.&
663 49 : & Shell3at%neighbours( iatom,ishell)%n_interactions) then
664 0 : ABI_ERROR('The interaction number index is not equal to the interaction number max computed previously (2)')
665 : end if
666 : !DEBUG iatref=Shell3at%iatref(ishell)
667 : !DEBUG jatref=Shell3at%jatref(ishell)
668 : !DEBUG katref=Shell3at%katref(ishell)
669 : !DEBUG if (Shell3at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
670 : !DEBUG do ninter=1,Shell3at%neighbours(iatom,ishell)%n_interactions
671 : !DEBUG jatom=Shell3at%neighbours(iatom,ishell)%atomj_in_shell(ninter)
672 : !DEBUG katom=Shell3at%neighbours(iatom,ishell)%atomk_in_shell(ninter)
673 : !DEBUG isym =Shell3at%neighbours(iatom,ishell)%sym_in_shell(ninter)
674 : !DEBUG itrans=Shell3at%neighbours(iatom,ishell)%transpose_in_shell(ninter)
675 : !DEBUG vectj(:)=zero ; vectk(:)=zero ; vect1(:)=zero ; vect2(:)=zero
676 : !DEBUG do ii=1,3
677 : !DEBUG do jj=1,3
678 : !DEBUG vectj(ii)=vectj(ii)+Sym%S_ref(ii,jj,isym,1)*MD%distance(iatref,jatref,jj+1)
679 : !DEBUG vectk(ii)=vectk(ii)+Sym%S_ref(ii,jj,isym,1)*MD%distance(iatref,katref,jj+1)
680 : !DEBUG end do
681 : !DEBUG end do
682 : !DEBUG if (itrans==1) then ; vect1(:)= vectj(:) ; vect2(:)= vectk(:) ; endif !\Phi3_ijk
683 : !DEBUG if (itrans==2) then ; vect1(:)= vectk(:) ; vect2(:)= vectj(:) ; endif !\Phi3_ikj
684 : !DEBUG if (itrans==3) then ; vect1(:)=-vectj(:) ; vect2(:)= vectk(:)-vectj(:) ; endif !\Phi3_jik
685 : !DEBUG if (itrans==4) then ; vect1(:)= vectk(:)-vectj(:) ; vect2(:)=-vectj(:) ; endif !\Phi3_jki
686 : !DEBUG if (itrans==5) then ; vect1(:)=-vectk(:) ; vect2(:)= vectj(:)-vectk(:) ; endif !\Phi3_kij
687 : !DEBUG if (itrans==6) then ; vect1(:)= vectj(:)-vectk(:) ; vect2(:)=-vectk(:) ; endif !\Phi3_kji
688 : !DEBUG do ii=1,3
689 : !DEBUG if ((abs(MD%distance(iatom,jatom,ii+1)-vect1(ii)).gt.tol8).or.&
690 : !DEBUG& (abs(MD%distance(iatom,katom,ii+1)-vect2(ii)).gt.tol8)) then
691 : !DEBUG write(std_out,'(a,4(x,i5))') 'For ishell, iatom, jatom, katom =',ishell,iatom,jatom,katom
692 : !DEBUG write(std_out,'(a,5(x,i5))') ' with isym, itrans, iatref, jatref, katref = ',isym,itrans,iatref,jatref,katref
693 : !DEBUG ABI_ERROR('We do not recover the triplet with the symmetry found')
694 : !DEBUG end if
695 : !DEBUG end do !ii
696 : !DEBUG end do !ninter
697 : end do !natom
698 : end do !nshell
699 14 : ABI_FREE(interactions)
700 :
701 : ! Find the number of coefficients of the (3x3x3) Phi3 for a given shell
702 37135 : ABI_CALLOC(Shell3at%proj, (norder,norder,nshell))
703 77 : ABI_MALLOC(Shell3at%ncoeff ,(nshell)); Shell3at%ncoeff(:)=zero
704 77 : ABI_MALLOC(Shell3at%ncoeff_prev,(nshell)); Shell3at%ncoeff_prev(:)=zero
705 14 : write(Invar%stdout,*) 'Number of shells=',nshell
706 14 : write(Invar%stdout,*) '============================================================================'
707 14 : if (MPIdata%iam_master) open(unit=16,file=trim(Invar%output_prefix)//'_nbcoeff-phi3.dat')
708 14 : ncoeff_prev=0
709 63 : do ishell=1,nshell
710 49 : ncoeff=0
711 49 : iatref=Shell3at%iatref(ishell)
712 49 : jatref=Shell3at%jatref(ishell)
713 49 : katref=Shell3at%katref(ishell)
714 49 : write(Invar%stdout,*) 'Shell number:',ishell
715 49 : write(Invar%stdout,'(a,i5,a,i5,a,f16.10)') ' Between atom',iatref,' and ',jatref,' the distance is=',MD%distance(iatref,jatref,1)
716 49 : write(Invar%stdout,'(a,i5,a,i5,a,f16.10)') ' Between atom',jatref,' and ',katref,' the distance is=',MD%distance(jatref,katref,1)
717 49 : write(Invar%stdout,'(a,i5,a,i5,a,f16.10)') ' Between atom',katref,' and ',iatref,' the distance is=',MD%distance(katref,iatref,1)
718 49 : call tdep_calc_nbcoeff(MD%distance,iatref,Invar,ishell,jatref,katref,1,MPIdata,ncoeff,norder,nshell,order,Shell3at%proj,Sym)
719 49 : Shell3at%ncoeff (ishell)=ncoeff
720 49 : Shell3at%ncoeff_prev(ishell)=ncoeff_prev
721 49 : ncoeff_prev=ncoeff_prev+ncoeff
722 49 : write(Invar%stdout,*)' Number of independant coefficients in this shell=',ncoeff
723 49 : write(Invar%stdout,*)' Number of interactions in this shell=',Shell3at%neighbours(iatref,ishell)%n_interactions
724 : !FB write(Invar%stdout,*)' The ratio is=',dfloat(Shell3at%neighbours(iatref,ishell)%n_interactions)/dfloat(ncoeff)
725 63 : write(Invar%stdout,*) '============================================================================'
726 : end do
727 14 : write(Invar%stdout,*)' >>>>>> Total number of coefficients at the third order=',ncoeff_prev
728 14 : if (MPIdata%iam_master) close(16)
729 14 : Shell3at%ntotcoeff=ncoeff_prev
730 : !BeginFB
731 : !FB open(unit=91,file='Shell3at.dat')
732 : !FB write(91,*) Shell3at%nshell
733 : !FB do ishell=1,Shell3at%nshell
734 : !FB write(91,*) Shell3at%ncoeff(ishell)
735 : !FB write(91,*) Shell3at%ncoeff_prev(ishell)
736 : !FB write(91,*) Shell3at%iatref(ishell)
737 : !FB write(91,*) Shell3at%jatref(ishell)
738 : !FB write(91,*) Shell3at%katref(ishell)
739 : !FB do iatom=1,Invar%natom
740 : !FB write(91,*) Shell3at%neighbours(iatom,ishell)%n_interactions
741 : !FB do ii=1,Shell3at%neighbours(iatom,ishell)%n_interactions
742 : !FB write(91,*) Shell3at%neighbours(iatom,ishell)%sym_in_shell(ii)
743 : !FB write(91,*) Shell3at%neighbours(iatom,ishell)%transpose_in_shell(ii)
744 : !FB write(91,*) Shell3at%neighbours(iatom,ishell)%atomj_in_shell(ii)
745 : !FB write(91,*) Shell3at%neighbours(iatom,ishell)%atomk_in_shell(ii)
746 : !FB end do
747 : !FB end do
748 : !FB end do
749 : !FB close(91)
750 : !EndFB
751 :
752 14 : end subroutine tdep_init_shell3at
753 :
754 : !====================================================================================================
755 :
756 8 : subroutine tdep_init_shell4at(Shell4at,Invar,MD,Sym,MPIdata)
757 :
758 : type(Shell_type),intent(out) :: Shell4at
759 : type(atdep_dataset_type),intent(in) :: Invar
760 : type(tdep_Sampling_type),intent(in) :: MD
761 : type(Symmetries_type),intent(inout) :: Sym
762 : type(MPI_enreg_type), intent(in) :: MPIdata
763 :
764 : integer :: ii,ishell,iatom,jatom,katom,latom,eatom,fatom,gatom,hatom,iatref,jatref,katref,latref
765 : integer :: natom,natom_unitcell,watom,xatom,yatom,zatom,ninteractions,ncoeff,ncoeff_prev,nshell_tmp
766 : integer :: already_found,find_equivalent,ninter,iat_ref,jat_ref,kat_ref,lat_ref,tmpinter
767 : integer :: norder,order,nshell_max,nshell
768 : double precision :: norma1,norma2,norma3
769 : double precision :: normb1,normb2,normb3
770 : double precision :: normc1,normc2,normc3
771 : integer :: Isym4at(2)
772 8 : integer, allocatable :: atref(:,:),interactions(:,:)
773 :
774 8 : natom = Invar%natom
775 8 : natom_unitcell = Invar%natom_unitcell
776 8 : nshell_max = Invar%nshell_max
777 8 : order = 4
778 8 : norder = 81
779 8 : Shell4at%order = order
780 8 : Shell4at%norder = norder
781 8 : Shell4at%natom = natom
782 :
783 8 : write(Invar%stdout,*) ' '
784 8 : write(Invar%stdout,*) '#############################################################################'
785 8 : write(Invar%stdout,*) '###### FOURTH ORDER : find the number of coefficients ########################'
786 8 : write(Invar%stdout,*) '#############################################################################'
787 :
788 : ! 1/ Identify the shells
789 12032 : ABI_CALLOC(interactions,(natom_unitcell,nshell_max))
790 16056 : ABI_CALLOC(atref,(nshell_max,4))
791 40 : ABI_CALLOC(Shell4at%ishell_self,(natom_unitcell))
792 8 : nshell_tmp=0
793 24 : do iatom=1,natom_unitcell
794 2848 : do jatom=1,natom
795 : ! Interactions are only computed until Rcut4 in order to have complete shell of neighbours.
796 : ! Otherwise the symetries are broken.
797 2824 : if (MD%distance(iatom,jatom,1).gt.(Invar%rcut4*0.99)) cycle
798 21620 : do katom=1,natom
799 21480 : if (MD%distance(iatom,katom,1).gt.(Invar%rcut4*0.99)) cycle
800 : ! WARNING: distance(j,k).ne.|djk| due to the inbox procedure when computing distance(j,k).
801 : ! So, compute |djk| using vec(ij) and vec(ik).
802 : norma1=dsqrt((MD%distance(iatom,katom,2)-MD%distance(iatom,jatom,2))**2+&
803 : & (MD%distance(iatom,katom,3)-MD%distance(iatom,jatom,3))**2+&
804 1440 : & (MD%distance(iatom,katom,4)-MD%distance(iatom,jatom,4))**2)
805 1440 : if (norma1 .gt.(Invar%rcut4*0.99)) cycle
806 125700 : do latom=1,natom
807 : !FB write(Invar%stdlog,*) 'NEW COORD1 : iatom,jatom,katom=',iatom,jatom,katom,latom
808 122152 : if (MD%distance(iatom,latom,1).gt.(Invar%rcut4*0.99)) cycle
809 : norma2=dsqrt((MD%distance(iatom,latom,2)-MD%distance(iatom,jatom,2))**2+&
810 : & (MD%distance(iatom,latom,3)-MD%distance(iatom,jatom,3))**2+&
811 8968 : & (MD%distance(iatom,latom,4)-MD%distance(iatom,jatom,4))**2)
812 8968 : if (norma2 .gt.(Invar%rcut4*0.99)) cycle
813 : norma3=dsqrt((MD%distance(iatom,latom,2)-MD%distance(iatom,katom,2))**2+&
814 : & (MD%distance(iatom,latom,3)-MD%distance(iatom,katom,3))**2+&
815 4944 : & (MD%distance(iatom,latom,4)-MD%distance(iatom,katom,4))**2)
816 4944 : if (norma3 .gt.(Invar%rcut4*0.99)) cycle
817 :
818 3460 : if (nshell_tmp.eq.0) then
819 40 : atref(1,:)=1
820 8 : nshell_tmp=nshell_tmp+1
821 8 : if ((iatom.eq.jatom).and.(jatom.eq.katom).and.(katom.eq.latom)) Shell4at%ishell_self(iatom)=nshell_tmp
822 8 : interactions(iatom,nshell_tmp)=interactions(iatom,nshell_tmp)+1
823 8 : cycle
824 : end if
825 3452 : find_equivalent=0
826 22849 : do ishell=1,nshell_tmp
827 22811 : iat_ref=atref(ishell,1) ; jat_ref=atref(ishell,2) ; kat_ref=atref(ishell,3) ; lat_ref=atref(ishell,4)
828 : normb1=dsqrt((MD%distance(iat_ref,kat_ref,2)-MD%distance(iat_ref,jat_ref,2))**2+&
829 : & (MD%distance(iat_ref,kat_ref,3)-MD%distance(iat_ref,jat_ref,3))**2+&
830 22811 : & (MD%distance(iat_ref,kat_ref,4)-MD%distance(iat_ref,jat_ref,4))**2)
831 : normb2=dsqrt((MD%distance(iat_ref,lat_ref,2)-MD%distance(iat_ref,jat_ref,2))**2+&
832 : & (MD%distance(iat_ref,lat_ref,3)-MD%distance(iat_ref,jat_ref,3))**2+&
833 22811 : & (MD%distance(iat_ref,lat_ref,4)-MD%distance(iat_ref,jat_ref,4))**2)
834 : normb3=dsqrt((MD%distance(iat_ref,lat_ref,2)-MD%distance(iat_ref,kat_ref,2))**2+&
835 : & (MD%distance(iat_ref,lat_ref,3)-MD%distance(iat_ref,kat_ref,3))**2+&
836 22811 : & (MD%distance(iat_ref,lat_ref,4)-MD%distance(iat_ref,kat_ref,4))**2)
837 513635 : do ii=1,24
838 494238 : if (ii.eq.1 ) then ; eatom=iatom ; fatom=jatom ; gatom=katom ; hatom=latom ; end if !ijkl
839 494238 : if (ii.eq.2 ) then ; eatom=iatom ; fatom=katom ; gatom=jatom ; hatom=latom ; end if !ikjl
840 494238 : if (ii.eq.3 ) then ; eatom=jatom ; fatom=iatom ; gatom=katom ; hatom=latom ; end if !jikl
841 494238 : if (ii.eq.4 ) then ; eatom=jatom ; fatom=katom ; gatom=iatom ; hatom=latom ; end if !jkil
842 494238 : if (ii.eq.5 ) then ; eatom=katom ; fatom=iatom ; gatom=jatom ; hatom=latom ; end if !kijl
843 494238 : if (ii.eq.6 ) then ; eatom=katom ; fatom=jatom ; gatom=iatom ; hatom=latom ; end if !kjil
844 :
845 494238 : if (ii.eq.7 ) then ; eatom=iatom ; fatom=jatom ; gatom=latom ; hatom=katom ; end if !ijlk
846 494238 : if (ii.eq.8 ) then ; eatom=iatom ; fatom=katom ; gatom=latom ; hatom=jatom ; end if !iklj
847 494238 : if (ii.eq.9 ) then ; eatom=jatom ; fatom=iatom ; gatom=latom ; hatom=katom ; end if !jilk
848 494238 : if (ii.eq.10) then ; eatom=jatom ; fatom=katom ; gatom=latom ; hatom=iatom ; end if !jkli
849 494238 : if (ii.eq.11) then ; eatom=katom ; fatom=iatom ; gatom=latom ; hatom=jatom ; end if !kilj
850 494238 : if (ii.eq.12) then ; eatom=katom ; fatom=jatom ; gatom=latom ; hatom=iatom ; end if !kjli
851 :
852 494238 : if (ii.eq.13) then ; eatom=iatom ; fatom=latom ; gatom=jatom ; hatom=katom ; end if !iljk
853 494238 : if (ii.eq.14) then ; eatom=iatom ; fatom=latom ; gatom=katom ; hatom=jatom ; end if !ilkj
854 494238 : if (ii.eq.15) then ; eatom=jatom ; fatom=latom ; gatom=iatom ; hatom=katom ; end if !jlik
855 494238 : if (ii.eq.16) then ; eatom=jatom ; fatom=latom ; gatom=katom ; hatom=iatom ; end if !jlki
856 494238 : if (ii.eq.17) then ; eatom=katom ; fatom=latom ; gatom=iatom ; hatom=jatom ; end if !klij
857 494238 : if (ii.eq.18) then ; eatom=katom ; fatom=latom ; gatom=jatom ; hatom=iatom ; end if !klji
858 :
859 494238 : if (ii.eq.19) then ; eatom=latom ; fatom=iatom ; gatom=jatom ; hatom=katom ; end if !lijk
860 494238 : if (ii.eq.20) then ; eatom=latom ; fatom=iatom ; gatom=katom ; hatom=jatom ; end if !likj
861 494238 : if (ii.eq.21) then ; eatom=latom ; fatom=jatom ; gatom=iatom ; hatom=katom ; end if !ljik
862 494238 : if (ii.eq.22) then ; eatom=latom ; fatom=jatom ; gatom=katom ; hatom=iatom ; end if !ljki
863 494238 : if (ii.eq.23) then ; eatom=latom ; fatom=katom ; gatom=iatom ; hatom=jatom ; end if !lkij
864 494238 : if (ii.eq.24) then ; eatom=latom ; fatom=katom ; gatom=jatom ; hatom=iatom ; end if !lkji
865 :
866 : normc1=dsqrt((MD%distance(eatom,gatom,2)-MD%distance(eatom,fatom,2))**2+&
867 : & (MD%distance(eatom,gatom,3)-MD%distance(eatom,fatom,3))**2+&
868 494238 : & (MD%distance(eatom,gatom,4)-MD%distance(eatom,fatom,4))**2)
869 : normc2=dsqrt((MD%distance(eatom,hatom,2)-MD%distance(eatom,fatom,2))**2+&
870 : & (MD%distance(eatom,hatom,3)-MD%distance(eatom,fatom,3))**2+&
871 494238 : & (MD%distance(eatom,hatom,4)-MD%distance(eatom,fatom,4))**2)
872 : normc3=dsqrt((MD%distance(eatom,hatom,2)-MD%distance(eatom,gatom,2))**2+&
873 : & (MD%distance(eatom,hatom,3)-MD%distance(eatom,gatom,3))**2+&
874 494238 : & (MD%distance(eatom,hatom,4)-MD%distance(eatom,gatom,4))**2)
875 : if ((abs(MD%distance(iat_ref,jat_ref,1)-MD%distance(eatom,fatom,1)).lt.1.d-6).and.&
876 : & (abs(normb1 -normc1 ).lt.1.d-6).and.&
877 : & (abs(normb2 -normc2 ).lt.1.d-6).and.&
878 : & (abs(normb3 -normc3 ).lt.1.d-6).and.&
879 494238 : & (abs(MD%distance(iat_ref,kat_ref,1)-MD%distance(eatom,gatom,1)).lt.1.d-6).and.&
880 19397 : & (abs(MD%distance(iat_ref,lat_ref,1)-MD%distance(eatom,hatom,1)).lt.1.d-6)) then
881 6414 : Isym4at(:)=0
882 6414 : call tdep_SearchS_4at(Invar,iat_ref,jat_ref,kat_ref,lat_ref,eatom,fatom,gatom,hatom,Isym4at,Sym,MD%xred_ideal)
883 6414 : if (Isym4at(2).eq.1) find_equivalent=1
884 : if (find_equivalent.eq.1) then
885 3414 : interactions(iatom,ishell)=interactions(iatom,ishell)+1
886 : !FB write(Invar%stdlog,*) 'The number of interactions in this shell is=',ishell,interactions(iatom,ishell)
887 3414 : exit
888 : end if
889 : end if
890 : end do !ii
891 22849 : if (find_equivalent.eq.1) exit
892 : end do !ishell
893 24932 : if (find_equivalent.eq.0) then
894 38 : nshell_tmp=nshell_tmp+1
895 38 : if (nshell_tmp.gt.nshell_max) then
896 0 : ABI_ERROR('The shell number index is greater than the shell number max defined in the code')
897 : end if
898 38 : if ((iatom.eq.jatom).and.(jatom.eq.katom).and.(katom.eq.latom)) Shell4at%ishell_self(iatom)=nshell_tmp
899 38 : interactions(iatom,nshell_tmp)=interactions(iatom,nshell_tmp)+1
900 38 : atref(nshell_tmp,1)=iatom
901 38 : atref(nshell_tmp,2)=jatom
902 38 : atref(nshell_tmp,3)=katom
903 38 : atref(nshell_tmp,4)=latom
904 : !FB write(Invar%stdlog,'(a,1x,5(i5,1x))') 'NEW SHELL1 : nshell_tmp,iatom,jatom,katom,latom=',nshell_tmp,iatom,jatom,katom,latom
905 : end if
906 : end do !latom
907 : end do !katom
908 : end do !jatom
909 : end do !iatom
910 8 : ABI_FREE(atref)
911 :
912 : ! 2/ Allocate the datatype Shell4at%...
913 8 : Shell4at%nshell = nshell_tmp
914 8 : nshell = Shell4at%nshell
915 7910 : ABI_MALLOC(Shell4at%neighbours,(natom,nshell))
916 70 : ABI_CALLOC(Shell4at%iatref,(nshell))
917 62 : ABI_CALLOC(Shell4at%jatref,(nshell))
918 62 : ABI_CALLOC(Shell4at%katref,(nshell))
919 62 : ABI_CALLOC(Shell4at%latref,(nshell))
920 54 : do ishell=1,nshell
921 7886 : do iatom=1,natom
922 7832 : ninteractions=interactions(mod(iatom-1,natom_unitcell)+1,ishell)
923 7832 : Shell4at%neighbours(iatom,ishell)%n_interactions=ninteractions
924 7832 : if (ninteractions.eq.0) cycle
925 20424 : ABI_MALLOC(Shell4at%neighbours(iatom,ishell)%atomj_in_shell,(ninteractions))
926 13616 : ABI_MALLOC(Shell4at%neighbours(iatom,ishell)%atomk_in_shell,(ninteractions))
927 13616 : ABI_MALLOC(Shell4at%neighbours(iatom,ishell)%atoml_in_shell,(ninteractions))
928 13616 : ABI_MALLOC(Shell4at%neighbours(iatom,ishell)%sym_in_shell,(ninteractions))
929 13616 : ABI_MALLOC(Shell4at%neighbours(iatom,ishell)%transpose_in_shell,(ninteractions))
930 250540 : Shell4at%neighbours(iatom,ishell)%atomj_in_shell(:)=zero
931 250540 : Shell4at%neighbours(iatom,ishell)%atomk_in_shell(:)=zero
932 250540 : Shell4at%neighbours(iatom,ishell)%atoml_in_shell(:)=zero
933 250540 : Shell4at%neighbours(iatom,ishell)%sym_in_shell(:)=zero
934 251610 : Shell4at%neighbours(iatom,ishell)%transpose_in_shell(:)=zero
935 : end do
936 : end do
937 8 : ABI_FREE(interactions)
938 :
939 : ! 3/ Store the index of the (couple of) atoms included in each shell
940 : ! 4/ Store the reference (couple of) atoms for each shell
941 : ! 5/ Compute the symetry operation between the reference (couple of) atoms and another one
942 7910 : ABI_MALLOC(interactions,(natom,nshell)) ; interactions(:,:)=0
943 8 : nshell_tmp=0
944 1484 : do iatom=1,natom
945 292188 : do jatom=1,natom
946 290704 : if (MD%distance(iatom,jatom,1).gt.(Invar%rcut4*0.99)) cycle
947 1607496 : do katom=1,natom
948 1597008 : if (MD%distance(iatom,katom,1).gt.(Invar%rcut4*0.99)) cycle
949 : ! WARNING: distance(j,k).ne.|djk| due to the inbox procedure when computing distance(j,k).
950 : ! So, compute |djk| using vec(ij) and vec(ik).
951 : norma1=dsqrt((MD%distance(iatom,katom,2)-MD%distance(iatom,jatom,2))**2+&
952 : & (MD%distance(iatom,katom,3)-MD%distance(iatom,jatom,3))**2+&
953 101348 : & (MD%distance(iatom,katom,4)-MD%distance(iatom,jatom,4))**2)
954 101348 : if (norma1 .gt.(Invar%rcut4*0.99)) cycle
955 9055796 : do latom=1,natom
956 : !FB write(Invar%stdlog,*) 'NEW COORD2 : iatom,jatom,katom=',iatom,jatom,katom,latom
957 8713936 : if (MD%distance(iatom,latom,1).gt.(Invar%rcut4*0.99)) cycle
958 : norma2=dsqrt((MD%distance(iatom,latom,2)-MD%distance(iatom,jatom,2))**2+&
959 : & (MD%distance(iatom,latom,3)-MD%distance(iatom,jatom,3))**2+&
960 631812 : & (MD%distance(iatom,latom,4)-MD%distance(iatom,jatom,4))**2)
961 631812 : if (norma2 .gt.(Invar%rcut4*0.99)) cycle
962 : norma3=dsqrt((MD%distance(iatom,latom,2)-MD%distance(iatom,katom,2))**2+&
963 : & (MD%distance(iatom,latom,3)-MD%distance(iatom,katom,3))**2+&
964 348068 : & (MD%distance(iatom,latom,4)-MD%distance(iatom,katom,4))**2)
965 348068 : if (norma3 .gt.(Invar%rcut4*0.99)) cycle
966 243732 : if (nshell_tmp.eq.0) then
967 8 : nshell_tmp=nshell_tmp+1
968 8 : interactions(iatom,nshell_tmp)=1
969 8 : Shell4at%iatref(nshell_tmp)=iatom
970 8 : Shell4at%jatref(nshell_tmp)=jatom
971 8 : Shell4at%katref(nshell_tmp)=katom
972 8 : Shell4at%latref(nshell_tmp)=latom
973 8 : Shell4at%neighbours(iatom,nshell_tmp)%atomj_in_shell(interactions(iatom,nshell_tmp))=jatom
974 8 : Shell4at%neighbours(iatom,nshell_tmp)%atomk_in_shell(interactions(iatom,nshell_tmp))=katom
975 8 : Shell4at%neighbours(iatom,nshell_tmp)%atoml_in_shell(interactions(iatom,nshell_tmp))=latom
976 8 : Shell4at%neighbours(iatom,nshell_tmp)%sym_in_shell(interactions(iatom,nshell_tmp))=1
977 8 : Shell4at%neighbours(iatom,nshell_tmp)%transpose_in_shell(interactions(iatom,nshell_tmp))=1
978 8 : cycle
979 : end if
980 :
981 : ! Search if the quadruplet has already been classified
982 243724 : find_equivalent=0
983 2934077 : do ishell=1,nshell_tmp
984 90581033 : do ninter=1,interactions(iatom,ishell)
985 : if ((Shell4at%neighbours(iatom,ishell)%atomj_in_shell(ninter).eq.jatom).and.&
986 87646956 : & (Shell4at%neighbours(iatom,ishell)%atomk_in_shell(ninter).eq.katom).and.&
987 90337309 : & (Shell4at%neighbours(iatom,ishell)%atoml_in_shell(ninter).eq.latom)) find_equivalent=1
988 : end do
989 : end do
990 243724 : if (find_equivalent.eq.1) cycle
991 : ! Search if the quadruplet belongs to a shell already found
992 159476 : do ishell=1,nshell_tmp
993 159438 : iat_ref=Shell4at%iatref(ishell) ; jat_ref=Shell4at%jatref(ishell)
994 159438 : kat_ref=Shell4at%katref(ishell) ; lat_ref=Shell4at%latref(ishell)
995 : normb1=dsqrt((MD%distance(iat_ref,kat_ref,2)-MD%distance(iat_ref,jat_ref,2))**2+&
996 : & (MD%distance(iat_ref,kat_ref,3)-MD%distance(iat_ref,jat_ref,3))**2+&
997 159438 : & (MD%distance(iat_ref,kat_ref,4)-MD%distance(iat_ref,jat_ref,4))**2)
998 : normb2=dsqrt((MD%distance(iat_ref,lat_ref,2)-MD%distance(iat_ref,jat_ref,2))**2+&
999 : & (MD%distance(iat_ref,lat_ref,3)-MD%distance(iat_ref,jat_ref,3))**2+&
1000 159438 : & (MD%distance(iat_ref,lat_ref,4)-MD%distance(iat_ref,jat_ref,4))**2)
1001 : normb3=dsqrt((MD%distance(iat_ref,lat_ref,2)-MD%distance(iat_ref,kat_ref,2))**2+&
1002 : & (MD%distance(iat_ref,lat_ref,3)-MD%distance(iat_ref,kat_ref,3))**2+&
1003 159438 : & (MD%distance(iat_ref,lat_ref,4)-MD%distance(iat_ref,kat_ref,4))**2)
1004 3484080 : do ii=1,24
1005 3352040 : if (ii.eq.1 ) then ; eatom=iatom ; fatom=jatom ; gatom=katom ; hatom=latom ; end if !ijkl
1006 3352040 : if (ii.eq.2 ) then ; eatom=iatom ; fatom=katom ; gatom=jatom ; hatom=latom ; end if !ikjl
1007 3352040 : if (ii.eq.3 ) then ; eatom=jatom ; fatom=iatom ; gatom=katom ; hatom=latom ; end if !jikl
1008 3352040 : if (ii.eq.4 ) then ; eatom=jatom ; fatom=katom ; gatom=iatom ; hatom=latom ; end if !jkil
1009 3352040 : if (ii.eq.5 ) then ; eatom=katom ; fatom=iatom ; gatom=jatom ; hatom=latom ; end if !kijl
1010 3352040 : if (ii.eq.6 ) then ; eatom=katom ; fatom=jatom ; gatom=iatom ; hatom=latom ; end if !kjil
1011 :
1012 3352040 : if (ii.eq.7 ) then ; eatom=iatom ; fatom=jatom ; gatom=latom ; hatom=katom ; end if !ijlk
1013 3352040 : if (ii.eq.8 ) then ; eatom=iatom ; fatom=katom ; gatom=latom ; hatom=jatom ; end if !iklj
1014 3352040 : if (ii.eq.9 ) then ; eatom=jatom ; fatom=iatom ; gatom=latom ; hatom=katom ; end if !jilk
1015 3352040 : if (ii.eq.10) then ; eatom=jatom ; fatom=katom ; gatom=latom ; hatom=iatom ; end if !jkli
1016 3352040 : if (ii.eq.11) then ; eatom=katom ; fatom=iatom ; gatom=latom ; hatom=jatom ; end if !kilj
1017 3352040 : if (ii.eq.12) then ; eatom=katom ; fatom=jatom ; gatom=latom ; hatom=iatom ; end if !kjli
1018 :
1019 3352040 : if (ii.eq.13) then ; eatom=iatom ; fatom=latom ; gatom=jatom ; hatom=katom ; end if !iljk
1020 3352040 : if (ii.eq.14) then ; eatom=iatom ; fatom=latom ; gatom=katom ; hatom=jatom ; end if !ilkj
1021 3352040 : if (ii.eq.15) then ; eatom=jatom ; fatom=latom ; gatom=iatom ; hatom=katom ; end if !jlik
1022 3352040 : if (ii.eq.16) then ; eatom=jatom ; fatom=latom ; gatom=katom ; hatom=iatom ; end if !jlki
1023 3352040 : if (ii.eq.17) then ; eatom=katom ; fatom=latom ; gatom=iatom ; hatom=jatom ; end if !klij
1024 3352040 : if (ii.eq.18) then ; eatom=katom ; fatom=latom ; gatom=jatom ; hatom=iatom ; end if !klji
1025 :
1026 3352040 : if (ii.eq.19) then ; eatom=latom ; fatom=iatom ; gatom=jatom ; hatom=katom ; end if !lijk
1027 3352040 : if (ii.eq.20) then ; eatom=latom ; fatom=iatom ; gatom=katom ; hatom=jatom ; end if !likj
1028 3352040 : if (ii.eq.21) then ; eatom=latom ; fatom=jatom ; gatom=iatom ; hatom=katom ; end if !ljik
1029 3352040 : if (ii.eq.22) then ; eatom=latom ; fatom=jatom ; gatom=katom ; hatom=iatom ; end if !ljki
1030 3352040 : if (ii.eq.23) then ; eatom=latom ; fatom=katom ; gatom=iatom ; hatom=jatom ; end if !lkij
1031 3352040 : if (ii.eq.24) then ; eatom=latom ; fatom=katom ; gatom=jatom ; hatom=iatom ; end if !lkji
1032 :
1033 : normc1=dsqrt((MD%distance(eatom,gatom,2)-MD%distance(eatom,fatom,2))**2+&
1034 : & (MD%distance(eatom,gatom,3)-MD%distance(eatom,fatom,3))**2+&
1035 3352040 : & (MD%distance(eatom,gatom,4)-MD%distance(eatom,fatom,4))**2)
1036 : normc2=dsqrt((MD%distance(eatom,hatom,2)-MD%distance(eatom,fatom,2))**2+&
1037 : & (MD%distance(eatom,hatom,3)-MD%distance(eatom,fatom,3))**2+&
1038 3352040 : & (MD%distance(eatom,hatom,4)-MD%distance(eatom,fatom,4))**2)
1039 : normc3=dsqrt((MD%distance(eatom,hatom,2)-MD%distance(eatom,gatom,2))**2+&
1040 : & (MD%distance(eatom,hatom,3)-MD%distance(eatom,gatom,3))**2+&
1041 3352040 : & (MD%distance(eatom,hatom,4)-MD%distance(eatom,gatom,4))**2)
1042 : if ((abs(MD%distance(iat_ref,jat_ref,1)-MD%distance(eatom,fatom,1)).lt.1.d-6).and.&
1043 : & (abs(normb1 -normc1 ).lt.1.d-6).and.&
1044 : & (abs(normb2 -normc2 ).lt.1.d-6).and.&
1045 : & (abs(normb3 -normc3 ).lt.1.d-6).and.&
1046 3352040 : & (abs(MD%distance(iat_ref,kat_ref,1)-MD%distance(eatom,gatom,1)).lt.1.d-6).and.&
1047 132040 : & (abs(MD%distance(iat_ref,lat_ref,1)-MD%distance(eatom,hatom,1)).lt.1.d-6)) then
1048 60410 : Isym4at(:)=0
1049 60410 : call tdep_SearchS_4at(Invar,iat_ref,jat_ref,kat_ref,lat_ref,eatom,fatom,gatom,hatom,Isym4at,Sym,MD%xred_ideal)
1050 60410 : if (Isym4at(2).eq.1) then
1051 : find_equivalent=1
1052 : exit
1053 : end if
1054 : end if
1055 : end do !ii
1056 159476 : if (find_equivalent.eq.1) exit
1057 : end do !ishell
1058 : ! The quadruplet belongs to a new shell
1059 27436 : if (find_equivalent.eq.0) then
1060 38 : nshell_tmp=nshell_tmp+1
1061 : ! Check that the new shell is allowed
1062 38 : if (nshell_tmp.gt.nshell) then
1063 0 : ABI_ERROR('The shell number index is greater than the shell number max computed previously')
1064 : end if
1065 38 : Shell4at%iatref(nshell_tmp)=iatom
1066 38 : Shell4at%jatref(nshell_tmp)=jatom
1067 38 : Shell4at%katref(nshell_tmp)=katom
1068 38 : Shell4at%latref(nshell_tmp)=latom
1069 38 : eatom=iatom ; fatom=jatom ; gatom=katom ; hatom=latom
1070 114 : Isym4at(:)=1
1071 38 : ishell=nshell_tmp
1072 : !FB write(Invar%stdlog,'(a,1x,5(i5,1x))') 'NEW SHELL2 : nshell_tmp,iatom,jatom,katom=',nshell_tmp,iatom,jatom,katom,latom
1073 : end if
1074 : ! Classify the informations of the quadruplet in Shell4at
1075 2282908 : do ii=1,24
1076 : ! The Phi4 has to be symetric (transposition symetries)
1077 658464 : if (ii.eq.1 ) then ; watom=eatom ; xatom=fatom ; yatom=gatom ; zatom=hatom ; end if !ijkl
1078 658464 : if (ii.eq.2 ) then ; watom=eatom ; xatom=gatom ; yatom=fatom ; zatom=hatom ; end if !ikjl
1079 658464 : if (ii.eq.3 ) then ; watom=fatom ; xatom=eatom ; yatom=gatom ; zatom=hatom ; end if !jikl
1080 658464 : if (ii.eq.4 ) then ; watom=fatom ; xatom=gatom ; yatom=eatom ; zatom=hatom ; end if !jkil
1081 658464 : if (ii.eq.5 ) then ; watom=gatom ; xatom=eatom ; yatom=fatom ; zatom=hatom ; end if !kijl
1082 658464 : if (ii.eq.6 ) then ; watom=gatom ; xatom=fatom ; yatom=eatom ; zatom=hatom ; end if !kjil
1083 :
1084 658464 : if (ii.eq.7 ) then ; watom=eatom ; xatom=fatom ; yatom=hatom ; zatom=gatom ; end if !ijlk
1085 658464 : if (ii.eq.8 ) then ; watom=eatom ; xatom=gatom ; yatom=hatom ; zatom=fatom ; end if !iklj
1086 658464 : if (ii.eq.9 ) then ; watom=fatom ; xatom=eatom ; yatom=hatom ; zatom=gatom ; end if !jilk
1087 658464 : if (ii.eq.10) then ; watom=fatom ; xatom=gatom ; yatom=hatom ; zatom=eatom ; end if !jkli
1088 658464 : if (ii.eq.11) then ; watom=gatom ; xatom=eatom ; yatom=hatom ; zatom=fatom ; end if !kilj
1089 658464 : if (ii.eq.12) then ; watom=gatom ; xatom=fatom ; yatom=hatom ; zatom=eatom ; end if !kjli
1090 :
1091 658464 : if (ii.eq.13) then ; watom=eatom ; xatom=hatom ; yatom=fatom ; zatom=gatom ; end if !iljk
1092 658464 : if (ii.eq.14) then ; watom=eatom ; xatom=hatom ; yatom=gatom ; zatom=fatom ; end if !ilkj
1093 658464 : if (ii.eq.15) then ; watom=fatom ; xatom=hatom ; yatom=eatom ; zatom=gatom ; end if !jlik
1094 658464 : if (ii.eq.16) then ; watom=fatom ; xatom=hatom ; yatom=gatom ; zatom=eatom ; end if !jlki
1095 658464 : if (ii.eq.17) then ; watom=gatom ; xatom=hatom ; yatom=eatom ; zatom=fatom ; end if !klij
1096 658464 : if (ii.eq.18) then ; watom=gatom ; xatom=hatom ; yatom=fatom ; zatom=eatom ; end if !klji
1097 :
1098 658464 : if (ii.eq.19) then ; watom=hatom ; xatom=eatom ; yatom=fatom ; zatom=gatom ; end if !lijk
1099 658464 : if (ii.eq.20) then ; watom=hatom ; xatom=eatom ; yatom=gatom ; zatom=fatom ; end if !likj
1100 658464 : if (ii.eq.21) then ; watom=hatom ; xatom=fatom ; yatom=eatom ; zatom=gatom ; end if !ljik
1101 658464 : if (ii.eq.22) then ; watom=hatom ; xatom=fatom ; yatom=gatom ; zatom=eatom ; end if !ljki
1102 658464 : if (ii.eq.23) then ; watom=hatom ; xatom=gatom ; yatom=eatom ; zatom=fatom ; end if !lkij
1103 658464 : if (ii.eq.24) then ; watom=hatom ; xatom=gatom ; yatom=fatom ; zatom=eatom ; end if !lkji
1104 : ! Do not overwrite the Phi4_iikl, Phi4_ijil, Phi4_ijjl, Phi4_iiil... IFCs
1105 : ! and avoid double counting of quadruplet interactions
1106 658464 : already_found=0
1107 22371744 : do ninter=1,interactions(watom,ishell)
1108 : if ((Shell4at%neighbours(watom,ishell)%atomj_in_shell(ninter).eq.xatom).and.&
1109 22128020 : & (Shell4at%neighbours(watom,ishell)%atomk_in_shell(ninter).eq.yatom).and.&
1110 243724 : & (Shell4at%neighbours(watom,ishell)%atoml_in_shell(ninter).eq.zatom)) then
1111 : !FB write(*,'(a,4(1x,i5))') 'FOR efgh =',eatom,fatom,gatom,hatom
1112 : !FB write(*,'(a,3(1x,i5))') ' --> ishell,ninter,ninter_tot =',ishell,ninter,interactions(watom,ishell)
1113 : !FB write(*,'(a,4(1x,i5))') ' --> ALREADY FOUND =',watom,xatom,yatom,zatom
1114 : already_found=1
1115 : exit
1116 : end if
1117 : end do
1118 658464 : if (already_found==1) cycle
1119 243724 : interactions(watom,ishell)=interactions(watom,ishell)+1
1120 : !FB write(Invar%stdlog,*) 'For ishell and eatom=',ishell,watom
1121 : !FB write(Invar%stdlog,*) ' --> the number of interactions in the shell is=',interactions(watom,ishell)
1122 243724 : if (interactions(watom,ishell).gt.Shell4at%neighbours(watom,ishell)%n_interactions) then
1123 0 : write(Invar%stdlog,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>'
1124 0 : write(Invar%stdlog,*) ' >>>>>> Verify that the Rcut used in the input file is lower '
1125 0 : write(Invar%stdlog,*) ' >>>>>> than half of the smallest lattice parameter'
1126 0 : write(Invar%stdlog,*) ' >>>>>> Solution : Reduce the Rcut parameter'
1127 0 : write(Invar%stdlog,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>'
1128 0 : ABI_ERROR('The interaction number index is greater than the interaction number max computed previously (4th order)')
1129 : end if
1130 243724 : Shell4at%neighbours(watom,ishell)%atomj_in_shell(interactions(watom,ishell))=xatom
1131 243724 : Shell4at%neighbours(watom,ishell)%atomk_in_shell(interactions(watom,ishell))=yatom
1132 243724 : Shell4at%neighbours(watom,ishell)%atoml_in_shell(interactions(watom,ishell))=zatom
1133 243724 : Shell4at%neighbours(watom,ishell)%sym_in_shell(interactions(watom,ishell))=Isym4at(1)
1134 9372392 : Shell4at%neighbours(watom,ishell)%transpose_in_shell(interactions(watom,ishell))=ii
1135 : !DEBUG write(Invar%stdlog,'(a,9(i5,x))') 'ishell,iatref,jatref,katref,iatom,atomj_in_shell,atomk_in_shell,isym,itrans=',&
1136 : !DEBUG& ishell,Shell4at%iatref(ishell),Shell4at%jatref(ishell),Shell4at%katref(ishell),watom,xatom,yatom,Isym4at(1),ii
1137 : end do !ii
1138 : end do !latom
1139 : end do !katom
1140 : end do !jatom
1141 : end do !iatom
1142 : ! Check that each interaction has different symmetry per shell
1143 54 : do ishell=1,nshell
1144 7886 : do iatom=1,natom
1145 7832 : if (Shell4at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
1146 243778 : do ninter=1,Shell4at%neighbours(iatom,ishell)%n_interactions-1
1147 11057276 : do tmpinter=ninter+1,Shell4at%neighbours(iatom,ishell)%n_interactions
1148 10812520 : if (Shell4at%neighbours(iatom,ishell)%sym_in_shell( ninter).eq.&
1149 236924 : & Shell4at%neighbours(iatom,ishell)%sym_in_shell(tmpinter)) then
1150 599796 : if (Shell4at%neighbours(iatom,ishell)%transpose_in_shell( ninter).ne.&
1151 : & Shell4at%neighbours(iatom,ishell)%transpose_in_shell(tmpinter)) cycle
1152 0 : write(std_out,'(a,2(1x,i5))') 'For ishell and iatom =',ishell,iatom
1153 0 : write(std_out,'(a,i5,a,i5,a,i5)') ' the interactions ',ninter,&
1154 0 : & ' and ',tmpinter,' have both the same symmetry isym=',Shell4at%neighbours(iatom,ishell)%sym_in_shell( ninter)
1155 0 : ABI_ERROR('Some interactions are equals due to the symmetry')
1156 : end if
1157 : end do
1158 : end do
1159 : end do
1160 : end do
1161 : ! Check that each equivalent shell has the same set of interactions
1162 54 : do ishell=1,nshell
1163 7886 : do iatom=1,natom
1164 7832 : if (Shell4at%neighbours(mod(iatom-1,natom_unitcell)+1,ishell)%n_interactions.ne.&
1165 46 : & Shell4at%neighbours( iatom,ishell)%n_interactions) then
1166 0 : ABI_ERROR('The interaction number index is not equal to the interaction number max computed previously (2)')
1167 : end if
1168 : !DEBUG iatref=Shell4at%iatref(ishell)
1169 : !DEBUG jatref=Shell4at%jatref(ishell)
1170 : !DEBUG katref=Shell4at%katref(ishell)
1171 : !DEBUG if (Shell4at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
1172 : !DEBUG do ninter=1,Shell4at%neighbours(iatom,ishell)%n_interactions
1173 : !DEBUG jatom=Shell4at%neighbours(iatom,ishell)%atomj_in_shell(ninter)
1174 : !DEBUG katom=Shell4at%neighbours(iatom,ishell)%atomk_in_shell(ninter)
1175 : !DEBUG isym =Shell4at%neighbours(iatom,ishell)%sym_in_shell(ninter)
1176 : !DEBUG itrans=Shell4at%neighbours(iatom,ishell)%transpose_in_shell(ninter)
1177 : !DEBUG vectj(:)=zero ; vectk(:)=zero ; vect1(:)=zero ; vect2(:)=zero
1178 : !DEBUG do ii=1,3
1179 : !DEBUG do jj=1,3
1180 : !DEBUG vectj(ii)=vectj(ii)+Sym%S_ref(ii,jj,isym,1)*MD%distance(iatref,jatref,jj+1)
1181 : !DEBUG vectk(ii)=vectk(ii)+Sym%S_ref(ii,jj,isym,1)*MD%distance(iatref,katref,jj+1)
1182 : !DEBUG end do
1183 : !DEBUG end do
1184 : !DEBUG if (itrans==1) then ; vect1(:)= vectj(:) ; vect2(:)= vectk(:) ; endif !\Phi3_ijk
1185 : !DEBUG if (itrans==2) then ; vect1(:)= vectk(:) ; vect2(:)= vectj(:) ; endif !\Phi3_ikj
1186 : !DEBUG if (itrans==3) then ; vect1(:)=-vectj(:) ; vect2(:)= vectk(:)-vectj(:) ; endif !\Phi3_jik
1187 : !DEBUG if (itrans==4) then ; vect1(:)= vectk(:)-vectj(:) ; vect2(:)=-vectj(:) ; endif !\Phi3_jki
1188 : !DEBUG if (itrans==5) then ; vect1(:)=-vectk(:) ; vect2(:)= vectj(:)-vectk(:) ; endif !\Phi3_kij
1189 : !DEBUG if (itrans==6) then ; vect1(:)= vectj(:)-vectk(:) ; vect2(:)=-vectk(:) ; endif !\Phi3_kji
1190 : !DEBUG do ii=1,3
1191 : !DEBUG if ((abs(MD%distance(iatom,jatom,ii+1)-vect1(ii)).gt.tol8).or.&
1192 : !DEBUG& (abs(MD%distance(iatom,katom,ii+1)-vect2(ii)).gt.tol8)) then
1193 : !DEBUG write(std_out,'(a,4(x,i5))') 'For ishell, iatom, jatom, katom =',ishell,iatom,jatom,katom
1194 : !DEBUG write(std_out,'(a,5(x,i5))') ' with isym, itrans, iatref, jatref, katref = ',isym,itrans,iatref,jatref,katref
1195 : !DEBUG ABI_ERROR('We do not recover the quadruplet with the symmetry found')
1196 : !DEBUG end if
1197 : !DEBUG end do !ii
1198 : !DEBUG end do !ninter
1199 : end do !natom
1200 : end do !nshell
1201 8 : ABI_FREE(interactions)
1202 :
1203 : ! Find the number of coefficients of the (3x3x3x3) Phi4 for a given shell
1204 305602 : ABI_CALLOC(Shell4at%proj, (norder,norder,nshell))
1205 62 : ABI_CALLOC(Shell4at%ncoeff ,(nshell))
1206 62 : ABI_CALLOC(Shell4at%ncoeff_prev,(nshell))
1207 8 : write(Invar%stdout,*) 'Number of shells=',nshell
1208 8 : write(Invar%stdout,*) '============================================================================'
1209 8 : if (MPIdata%iam_master) open(unit=16,file=trim(Invar%output_prefix)//'_nbcoeff-phi4.dat')
1210 8 : ncoeff_prev=0
1211 54 : do ishell=1,nshell
1212 46 : ncoeff=0
1213 46 : iatref=Shell4at%iatref(ishell)
1214 46 : jatref=Shell4at%jatref(ishell)
1215 46 : katref=Shell4at%katref(ishell)
1216 46 : latref=Shell4at%latref(ishell)
1217 46 : write(Invar%stdout,*) 'Shell number:',ishell
1218 46 : write(Invar%stdout,'(a,i5,a,i5,a,f16.10)') ' Between atom',iatref,' and ',jatref,' the distance is=',MD%distance(iatref,jatref,1)
1219 46 : write(Invar%stdout,'(a,i5,a,i5,a,f16.10)') ' Between atom',jatref,' and ',katref,' the distance is=',MD%distance(jatref,katref,1)
1220 46 : write(Invar%stdout,'(a,i5,a,i5,a,f16.10)') ' Between atom',katref,' and ',latref,' the distance is=',MD%distance(katref,latref,1)
1221 46 : write(Invar%stdout,'(a,i5,a,i5,a,f16.10)') ' Between atom',latref,' and ',iatref,' the distance is=',MD%distance(latref,iatref,1)
1222 46 : call tdep_calc_nbcoeff(MD%distance,iatref,Invar,ishell,jatref,katref,latref,MPIdata,ncoeff,norder,nshell,order,Shell4at%proj,Sym)
1223 46 : Shell4at%ncoeff (ishell)=ncoeff
1224 46 : Shell4at%ncoeff_prev(ishell)=ncoeff_prev
1225 46 : ncoeff_prev=ncoeff_prev+ncoeff
1226 46 : write(Invar%stdout,*)' Number of independant coefficients in this shell=',ncoeff
1227 46 : write(Invar%stdout,*)' Number of interactions in this shell=',Shell4at%neighbours(iatref,ishell)%n_interactions
1228 : !FB write(Invar%stdout,*)' The ratio is=',dfloat(Shell4at%neighbours(iatref,ishell)%n_interactions)/dfloat(ncoeff)
1229 54 : write(Invar%stdout,*) '============================================================================'
1230 : end do
1231 8 : write(Invar%stdout,*)' >>>>>> Total number of coefficients at the fourth order=',ncoeff_prev
1232 8 : if (MPIdata%iam_master) close(16)
1233 8 : Shell4at%ntotcoeff=ncoeff_prev
1234 : !BeginFB
1235 : !FB open(unit=91,file='Shell4at.dat')
1236 : !FB write(91,*) Shell4at%nshell
1237 : !FB do ishell=1,Shell4at%nshell
1238 : !FB write(91,*) Shell4at%ncoeff(ishell)
1239 : !FB write(91,*) Shell4at%ncoeff_prev(ishell)
1240 : !FB write(91,*) Shell4at%iatref(ishell)
1241 : !FB write(91,*) Shell4at%jatref(ishell)
1242 : !FB write(91,*) Shell4at%katref(ishell)
1243 : !FB do iatom=1,Invar%natom
1244 : !FB write(91,*) Shell4at%neighbours(iatom,ishell)%n_interactions
1245 : !FB do ii=1,Shell4at%neighbours(iatom,ishell)%n_interactions
1246 : !FB write(91,*) Shell4at%neighbours(iatom,ishell)%sym_in_shell(ii)
1247 : !FB write(91,*) Shell4at%neighbours(iatom,ishell)%transpose_in_shell(ii)
1248 : !FB write(91,*) Shell4at%neighbours(iatom,ishell)%atomj_in_shell(ii)
1249 : !FB write(91,*) Shell4at%neighbours(iatom,ishell)%atomk_in_shell(ii)
1250 : !FB end do
1251 : !FB end do
1252 : !FB end do
1253 : !FB close(91)
1254 : !EndFB
1255 :
1256 8 : end subroutine tdep_init_shell4at
1257 :
1258 : !====================================================================================================
1259 110 : subroutine tdep_destroy_shell(Shell)
1260 :
1261 : type(Shell_type),intent(inout) :: Shell
1262 :
1263 : integer :: iatom,ishell,natref
1264 :
1265 110 : ABI_FREE(Shell%ncoeff)
1266 110 : ABI_FREE(Shell%ncoeff_prev)
1267 110 : ABI_FREE(Shell%iatref)
1268 110 : if (Shell%order.gt.1) then
1269 66 : ABI_FREE(Shell%jatref)
1270 66 : ABI_FREE(Shell%ishell_self)
1271 : end if
1272 110 : if (Shell%order.gt.2) then
1273 22 : ABI_FREE(Shell%katref)
1274 : end if
1275 110 : if (Shell%order.gt.3) then
1276 8 : ABI_FREE(Shell%latref)
1277 : end if
1278 110 : if (Shell%order.eq.1) then
1279 : natref=1
1280 : else
1281 66 : natref=Shell%natom
1282 : end if
1283 11052 : do iatom=1,natref
1284 105849 : do ishell=1,Shell%nshell
1285 105739 : if (Shell%neighbours(iatom,ishell)%n_interactions.ne.0) then
1286 80910 : ABI_FREE(Shell%neighbours(iatom,ishell)%atomj_in_shell)
1287 80910 : ABI_FREE(Shell%neighbours(iatom,ishell)%sym_in_shell)
1288 80910 : if (Shell%order.gt.1) then
1289 80910 : ABI_FREE(Shell%neighbours(iatom,ishell)%transpose_in_shell)
1290 : end if
1291 80910 : if (Shell%order.gt.2) then
1292 13856 : ABI_FREE(Shell%neighbours(iatom,ishell)%atomk_in_shell)
1293 : end if
1294 80910 : if (Shell%order.gt.3) then
1295 6808 : ABI_FREE(Shell%neighbours(iatom,ishell)%atoml_in_shell)
1296 : end if
1297 : end if
1298 : end do
1299 : end do
1300 94907 : ABI_FREE(Shell%neighbours)
1301 110 : ABI_FREE(Shell%proj)
1302 :
1303 110 : end subroutine tdep_destroy_shell
1304 :
1305 : !====================================================================================================
1306 :
1307 658 : subroutine tdep_calc_nbcoeff(distance,iatcell,Invar,ishell,jatom,katom,latom,MPIdata,&
1308 658 : & ncoeff,norder,nshell,order,proj,Sym)
1309 :
1310 : integer,intent(in) :: iatcell,ishell,jatom,katom,latom,nshell,order,norder
1311 : integer,intent(inout) :: ncoeff
1312 : type(atdep_dataset_type),intent(in) :: Invar
1313 : type(Symmetries_type),intent(in) :: Sym
1314 : type(MPI_enreg_type), intent(in) :: MPIdata
1315 : double precision,intent(in) :: distance(Invar%natom,Invar%natom,4)
1316 : double precision,intent(out) :: proj(norder,norder,nshell)
1317 :
1318 : integer :: ii,jj,kk,ll,isym,LWORK,INFO,const_tot,itemp,nconst_perm,nconst_loc
1319 : integer :: ncount,icoeff,jatcell,katcell,latcell,mu,nu,xi,zeta
1320 : integer :: inv,watom,xatom,yatom,zatom,isyminv,nsyminv,facorder,iseed
1321 658 : integer, allocatable :: iconst(:)
1322 : double precision :: prod_scal,drandom
1323 : double precision :: eigvec(3,3)
1324 : double precision :: vect_trial(3),vect_trial1(3),vect_trial2(3),vect_trial3(3)
1325 : double precision :: vect_trial4(3),vect_trial5(3),vect_trial6(3)
1326 : double precision :: WR(3),WI(3),VL(3,3),VR(3,3)
1327 658 : double precision, allocatable :: WORK(:)
1328 : double complex :: eigenvectors(3,3),eigenvalues(3)
1329 : double complex :: pp(3,3),ppp(3,3,3),pppp(3,3,3,3),lambda
1330 658 : double complex, allocatable :: tab_vec(:,:),temp(:,:),alphaij(:,:,:),constraints(:,:,:)
1331 : logical :: ok
1332 658 : logical, allocatable :: unchanged(:)
1333 : character(len=500) :: message
1334 :
1335 658 : if (iatcell==1.and.order==1) return
1336 614 : if (jatom==iatcell.and.order==2) return
1337 : !FB if (katom==iatcell.and.jatom==iatcell.and.order==3) return
1338 :
1339 551 : if (order==1) then
1340 : facorder=1
1341 532 : else if (order==2) then
1342 : facorder=2
1343 95 : else if (order==3) then
1344 : facorder=6
1345 46 : else if (order==4) then
1346 46 : facorder=24
1347 : end if
1348 :
1349 : ! If we want to remove the constraints coming from the symetries
1350 : !FB if (order==3) then
1351 : !FB do ii=1,norder
1352 : !FB proj(ii,ii,ishell)=1.d0
1353 : !FB end do
1354 : !FB ncoeff=norder
1355 : !FB return
1356 : !FB end if
1357 :
1358 551 : nconst_loc=0
1359 551 : const_tot=0
1360 551 : nsyminv=Sym%nsym*facorder
1361 309119431 : ABI_MALLOC(alphaij,(nsyminv,norder,norder)); alphaij(:,:,:)=czero
1362 82323 : ABI_MALLOC(iconst,(nsyminv)) ; iconst(:)=0
1363 81772 : ABI_MALLOC(unchanged,(nsyminv)) ; unchanged(:)=.false.
1364 :
1365 : ! ================================================================================================
1366 : ! ================ Big loop over symmetries and invariance (nsym*facorder) =======================
1367 : ! ================================================================================================
1368 551 : if (MPIdata%iam_master) write(16,'(a)') ' '
1369 551 : if (MPIdata%iam_master) write(16,'(a,i4)') 'For shell number=',ishell
1370 80939 : do isyminv=1,nsyminv
1371 80394 : isym=(isyminv-1)/facorder+1
1372 80394 : inv=isyminv-(isym-1)*facorder
1373 80394 : if (isym==1) cycle
1374 :
1375 : ! For the 1st order: Search if the atom is let invariant
1376 78103 : if (order==1) then
1377 223 : if (Sym%indsym(4,isym,iatcell)==iatcell) then
1378 107 : if (MPIdata%iam_master) then
1379 107 : write(16,'(a,1x,i3)')'===========The atom is kept invariant for isym=',isym
1380 : end if
1381 : else
1382 : cycle
1383 : end if
1384 : end if
1385 :
1386 : ! For the 2nd order: Search if the bond is kept invariant or reversed
1387 77987 : if (order==2) then
1388 21534 : vect_trial(:)=zero
1389 21534 : if (inv==1) then ; watom=iatcell ; xatom=jatom ; endif !\Phi_ij
1390 21534 : if (inv==2) then ; watom=jatom ; xatom=iatcell ; endif !\Phi_ji
1391 86136 : do ii=1,3
1392 279942 : do jj=1,3
1393 258408 : vect_trial(ii)=vect_trial(ii)+Sym%S_ref(ii,jj,isym,1)*distance(watom,xatom,jj+1)
1394 : end do
1395 : end do
1396 21534 : jatcell=mod(jatom-1,Invar%natom_unitcell)+1
1397 : if ((sum(abs(vect_trial(:)-distance(iatcell,jatom,2:4))).lt.tol8).and.&
1398 86136 : & (Sym%indsym(4,isym,watom)==iatcell).and.&
1399 : & (Sym%indsym(4,isym,xatom)==jatcell)) then
1400 1664 : if (MPIdata%iam_master) then
1401 1664 : if (inv==1) write(16,'(a,1x,i3)')'===========The bond is kept invariant for isym=',isym
1402 1664 : if (inv==2) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j) --> (j,i) for isym=',isym
1403 : end if
1404 : else
1405 : cycle
1406 : end if
1407 : end if
1408 :
1409 : ! For the 3rd order : 6 permutations at all
1410 58117 : if (order==3) then
1411 12042 : vect_trial1(:)=zero
1412 12042 : vect_trial2(:)=zero
1413 12042 : vect_trial3(:)=zero
1414 12042 : if (inv==1) then ; watom=iatcell ; xatom=jatom ; yatom=katom ; endif !\Phi3_ijk
1415 12042 : if (inv==2) then ; watom=iatcell ; xatom=katom ; yatom=jatom ; endif !\Phi3_ikj
1416 12042 : if (inv==3) then ; watom=jatom ; xatom=iatcell ; yatom=katom ; endif !\Phi3_jik
1417 12042 : if (inv==4) then ; watom=jatom ; xatom=katom ; yatom=iatcell ; endif !\Phi3_jki
1418 12042 : if (inv==5) then ; watom=katom ; xatom=iatcell ; yatom=jatom ; endif !\Phi3_kij
1419 12042 : if (inv==6) then ; watom=katom ; xatom=jatom ; yatom=iatcell ; endif !\Phi3_kji
1420 48168 : do ii=1,3
1421 156546 : do jj=1,3
1422 108378 : vect_trial1(ii)=vect_trial1(ii)+Sym%S_ref(ii,jj,isym,1)*distance(watom,xatom,jj+1)
1423 108378 : vect_trial2(ii)=vect_trial2(ii)+Sym%S_ref(ii,jj,isym,1)*distance(xatom,yatom,jj+1)
1424 144504 : vect_trial3(ii)=vect_trial3(ii)+Sym%S_ref(ii,jj,isym,1)*distance(yatom,watom,jj+1)
1425 : end do
1426 : end do
1427 12042 : jatcell=mod(jatom-1,Invar%natom_unitcell)+1
1428 12042 : katcell=mod(katom-1,Invar%natom_unitcell)+1
1429 : if ((sum(abs(vect_trial1(:)-distance(iatcell,jatom ,2:4))).lt.tol8).and.&
1430 : & (sum(abs(vect_trial2(:)-distance(jatom ,katom ,2:4))).lt.tol8).and.&
1431 : & (sum(abs(vect_trial3(:)-distance(katom ,iatcell,2:4))).lt.tol8).and.&
1432 : & (Sym%indsym(4,isym,watom)==iatcell).and.&
1433 120420 : & (Sym%indsym(4,isym,xatom)==jatcell).and.&
1434 : & (Sym%indsym(4,isym,yatom)==katcell)) then
1435 3766 : if (MPIdata%iam_master) then
1436 3766 : if (inv==1) write(16,'(a,1x,i3)')'===========The bond is kept invariant for isym=',isym
1437 3766 : if (inv==2) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k) --> (i,k,j) for isym=',isym
1438 3766 : if (inv==3) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k) --> (j,i,k) for isym=',isym
1439 3766 : if (inv==4) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k) --> (j,k,i) for isym=',isym
1440 3766 : if (inv==5) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k) --> (k,i,j) for isym=',isym
1441 3766 : if (inv==6) write(16,'(a,1x,i3)')'===========The bond is reversed with (i,j,k) --> (k,j,i) for isym=',isym
1442 : end if
1443 : else
1444 : cycle
1445 : end if
1446 : end if
1447 :
1448 : ! For the 4th order : 24 permutations at all
1449 49841 : if (order==4) then
1450 44304 : vect_trial1(:)=zero
1451 44304 : vect_trial2(:)=zero
1452 44304 : vect_trial3(:)=zero
1453 44304 : vect_trial4(:)=zero
1454 44304 : vect_trial5(:)=zero
1455 44304 : vect_trial6(:)=zero
1456 44304 : if (inv==1) then ; watom=iatcell ; xatom=jatom ; yatom=katom ; zatom=latom ; endif !\Phi4_ijkl
1457 44304 : if (inv==2) then ; watom=iatcell ; xatom=katom ; yatom=jatom ; zatom=latom ; endif !\Phi4_ikjl
1458 44304 : if (inv==3) then ; watom=jatom ; xatom=iatcell ; yatom=katom ; zatom=latom ; endif !\Phi4_jikl
1459 44304 : if (inv==4) then ; watom=jatom ; xatom=katom ; yatom=iatcell ; zatom=latom ; endif !\Phi4_jkil
1460 44304 : if (inv==5) then ; watom=katom ; xatom=iatcell ; yatom=jatom ; zatom=latom ; endif !\Phi4_kijl
1461 44304 : if (inv==6) then ; watom=katom ; xatom=jatom ; yatom=iatcell ; zatom=latom ; endif !\Phi4_kjil
1462 :
1463 44304 : if (inv==7 ) then ; watom=iatcell ; xatom=jatom ; yatom=latom ; zatom=katom ; endif !\Phi4_ijlk
1464 44304 : if (inv==8 ) then ; watom=iatcell ; xatom=katom ; yatom=latom ; zatom=jatom ; endif !\Phi4_iklj
1465 44304 : if (inv==9 ) then ; watom=jatom ; xatom=iatcell ; yatom=latom ; zatom=katom ; endif !\Phi4_jilk
1466 44304 : if (inv==10) then ; watom=jatom ; xatom=katom ; yatom=latom ; zatom=iatcell ; endif !\Phi4_jkli
1467 44304 : if (inv==11) then ; watom=katom ; xatom=iatcell ; yatom=latom ; zatom=jatom ; endif !\Phi4_kilj
1468 44304 : if (inv==12) then ; watom=katom ; xatom=jatom ; yatom=latom ; zatom=iatcell ; endif !\Phi4_kjli
1469 :
1470 44304 : if (inv==13) then ; watom=iatcell ; xatom=latom ; yatom=jatom ; zatom=katom ; endif !\Phi4_iljk
1471 44304 : if (inv==14) then ; watom=iatcell ; xatom=latom ; yatom=katom ; zatom=jatom ; endif !\Phi4_ilkj
1472 44304 : if (inv==15) then ; watom=jatom ; xatom=latom ; yatom=iatcell ; zatom=katom ; endif !\Phi4_jlik
1473 44304 : if (inv==16) then ; watom=jatom ; xatom=latom ; yatom=katom ; zatom=iatcell ; endif !\Phi4_jlki
1474 44304 : if (inv==17) then ; watom=katom ; xatom=latom ; yatom=iatcell ; zatom=jatom ; endif !\Phi4_klij
1475 44304 : if (inv==18) then ; watom=katom ; xatom=latom ; yatom=jatom ; zatom=iatcell ; endif !\Phi4_klji
1476 :
1477 44304 : if (inv==19) then ; watom=latom ; xatom=iatcell ; yatom=jatom ; zatom=katom ; endif !\Phi4_lijk
1478 44304 : if (inv==20) then ; watom=latom ; xatom=iatcell ; yatom=katom ; zatom=jatom ; endif !\Phi4_likj
1479 44304 : if (inv==21) then ; watom=latom ; xatom=jatom ; yatom=iatcell ; zatom=katom ; endif !\Phi4_ljik
1480 44304 : if (inv==22) then ; watom=latom ; xatom=jatom ; yatom=katom ; zatom=iatcell ; endif !\Phi4_ljki
1481 44304 : if (inv==23) then ; watom=latom ; xatom=katom ; yatom=iatcell ; zatom=jatom ; endif !\Phi4_lkij
1482 44304 : if (inv==24) then ; watom=latom ; xatom=katom ; yatom=jatom ; zatom=iatcell ; endif !\Phi4_lkji
1483 :
1484 177216 : do ii=1,3
1485 575952 : do jj=1,3
1486 398736 : vect_trial1(ii)=vect_trial1(ii)+Sym%S_ref(ii,jj,isym,1)*distance(watom,xatom,jj+1)
1487 398736 : vect_trial2(ii)=vect_trial2(ii)+Sym%S_ref(ii,jj,isym,1)*distance(watom,yatom,jj+1)
1488 398736 : vect_trial3(ii)=vect_trial3(ii)+Sym%S_ref(ii,jj,isym,1)*distance(watom,zatom,jj+1)
1489 398736 : vect_trial4(ii)=vect_trial4(ii)+Sym%S_ref(ii,jj,isym,1)*distance(xatom,yatom,jj+1)
1490 398736 : vect_trial5(ii)=vect_trial5(ii)+Sym%S_ref(ii,jj,isym,1)*distance(xatom,zatom,jj+1)
1491 531648 : vect_trial6(ii)=vect_trial6(ii)+Sym%S_ref(ii,jj,isym,1)*distance(yatom,zatom,jj+1)
1492 : end do
1493 : end do
1494 44304 : jatcell=mod(jatom-1,Invar%natom_unitcell)+1
1495 44304 : katcell=mod(katom-1,Invar%natom_unitcell)+1
1496 44304 : latcell=mod(latom-1,Invar%natom_unitcell)+1
1497 : if ((sum(abs(vect_trial1(:)-distance(iatcell,jatom,2:4))).lt.tol8).and.&
1498 : & (sum(abs(vect_trial2(:)-distance(iatcell,katom,2:4))).lt.tol8).and.&
1499 : & (sum(abs(vect_trial3(:)-distance(iatcell,latom,2:4))).lt.tol8).and.&
1500 : & (sum(abs(vect_trial4(:)-distance(jatom ,katom,2:4))).lt.tol8).and.&
1501 : & (sum(abs(vect_trial5(:)-distance(jatom ,latom,2:4))).lt.tol8).and.&
1502 : & (sum(abs(vect_trial6(:)-distance(katom ,latom,2:4))).lt.tol8).and.&
1503 : & (Sym%indsym(4,isym,watom)==iatcell).and.&
1504 : & (Sym%indsym(4,isym,xatom)==jatcell).and.&
1505 841776 : & (Sym%indsym(4,isym,yatom)==katcell).and.&
1506 : & (Sym%indsym(4,isym,zatom)==latcell)) then
1507 8954 : if (MPIdata%iam_master) then
1508 8954 : if (inv==1 ) write(16,'(a,1x,i3)')'===========The bond is kept invariant for isym=',isym
1509 8954 : 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
1510 8954 : 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
1511 8954 : 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
1512 8954 : 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
1513 8954 : 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
1514 :
1515 8954 : 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
1516 8954 : 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
1517 8954 : 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
1518 8954 : 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
1519 8954 : 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
1520 8954 : 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
1521 :
1522 8954 : 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
1523 8954 : 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
1524 8954 : 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
1525 8954 : 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
1526 8954 : 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
1527 8954 : 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
1528 :
1529 8954 : 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
1530 8954 : 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
1531 8954 : 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
1532 8954 : 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
1533 8954 : 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
1534 8954 : 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
1535 : end if
1536 : else
1537 : cycle
1538 : end if
1539 : end if
1540 :
1541 : ! Write the S_ref matrix
1542 : !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)
1543 : !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)
1544 : !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)
1545 :
1546 : ! Diagonalize the S_ref matrix
1547 57964 : do ii=1,3
1548 188383 : do jj=1,3
1549 173892 : eigvec(ii,jj)=Sym%S_ref(jj,ii,isym,1)
1550 : end do
1551 : end do
1552 14491 : LWORK=4*3
1553 188383 : ABI_MALLOC(WORK,(LWORK)); WORK(:)=zero
1554 : ! This one is real and could be non-symmetric
1555 14491 : call dgeev( 'N', 'V', 3, eigvec, 3, WR, WI, VL, 3, VR, 3, WORK, LWORK, INFO)
1556 14491 : ABI_FREE(WORK)
1557 :
1558 : ! Build the real and imaginary parts of the eigenvectors and eigenvalues
1559 14491 : jj=0
1560 57964 : do ii=1,3
1561 43473 : eigenvalues(ii)=dcmplx(WR(ii),-WI(ii))
1562 57964 : if (WI(ii).ne.zero.and.jj==0) then
1563 29944 : do kk=1,3
1564 29944 : eigenvectors(kk,ii)=dcmplx(VR(kk,ii),VR(kk,ii+1))
1565 : end do
1566 : jj=jj+1
1567 35987 : else if (WI(ii).ne.zero.and.jj==1) then
1568 29944 : do kk=1,3
1569 29944 : eigenvectors(kk,ii)=dcmplx(VR(kk,ii-1),-VR(kk,ii))
1570 : end do
1571 : jj=jj+1
1572 : else
1573 114004 : do kk=1,3
1574 114004 : eigenvectors(kk,ii)=dcmplx(VR(kk,ii),zero)
1575 : end do
1576 : end if
1577 : end do
1578 :
1579 : ! Write the eigenvalues and eigenvectors
1580 14491 : ok=.true.
1581 57964 : do ii=1,3
1582 57964 : if ((aimag(eigenvalues(1)).ne.0).or.(aimag(eigenvalues(2)).ne.0).or.(aimag(eigenvalues(3)).ne.0)) then
1583 22458 : ok=.false.
1584 : end if
1585 : end do
1586 14491 : if (.not.ok.and.MPIdata%iam_master) write(16,'(a)') ' WARNING: THERE IS COMPLEX EIGENVALUES'
1587 :
1588 : ! If the transformation matrix keeps the bond invariant:
1589 : ! Phi_{\alpha\beta}=\sum_{\mu\nu} S_{\alpha\mu}.S_{\beta\nu}.Phi_{\mu\nu}
1590 : ! If lambda and p are the eigenvectors and eigenvalues of the S matrix, then:
1591 : ! \sum_{\alpha\beta} p_{\alpha}^l.p_{\beta}^k Phi_{\alpha\beta}
1592 : ! = \sum_{\mu\nu,\alpha\beta} p_{\alpha}^l.p_{\beta}^k.S_{\alpha\mu}.S_{\beta\nu}.Phi_{\mu\nu}
1593 : ! = lambda^{*l}.lambda^{*k} \sum_{\mu\nu} p_{\mu}^l.p_{\nu}^k.Phi_{\mu\nu}
1594 : ! So, if lambda^{*l}.lambda^{*k} = -1, we must have:
1595 : ! \sum_{\alpha\beta} p_{\alpha}^l.p_{\beta}^k.Phi_{\alpha\beta}= 0
1596 : !
1597 : ! In the case of the reversed bond, one obtains the following constraint:
1598 : ! \sum_{\alpha\beta} (lambda^{*l}.lambda^{*k}.p_{\alpha}^l.p_{\beta}^k-p_{\beta}^l.p_{\alpha}^k).Phi_{\alpha\beta}= 0
1599 : ! which applies whether lambda^{*l}.lambda^{*k} = \pm 1
1600 : !
1601 : ! We obtain n vectors with norder coefficients (defined in the R^norder space).
1602 : ! The space of the independent solutions are in the R^(norder-n) space, orthogonal
1603 : ! to the space spanned by the starting n vectors.
1604 14491 : if (order==1) then
1605 428 : do ii=1,3
1606 321 : lambda=eigenvalues(ii)
1607 321 : if ((abs(real(lambda)-1.d0).lt.tol6).and.(abs(aimag(lambda)).lt.tol6)) cycle
1608 211 : unchanged(isyminv)=.true.
1609 211 : iconst(isyminv)=iconst(isyminv)+1
1610 : !FB const_tot=const_tot+1
1611 : !FB write(16,*)' The eigenvalue',ii
1612 : !FB write(16,*)' is equal to ',lambda
1613 951 : do mu=1,3
1614 954 : alphaij(isyminv,mu,iconst(isyminv))=eigenvectors(mu,ii)
1615 : end do
1616 : !FB write(16,*)' Real & imaginary parts of the eigenvectors product:'
1617 : !FB write(16,'(3(f16.12,1x))') real(alphaij(isyminv,:,iconst(isyminv)))
1618 : !FB write(16,'(3(f16.12,1x))') aimag(alphaij(isyminv,:,iconst(isyminv)))
1619 : end do !ii
1620 14384 : else if (order==2) then
1621 6656 : do ii=1,3
1622 21632 : do jj=1,3
1623 14976 : lambda=eigenvalues(ii)*eigenvalues(jj)
1624 8632 : if (((abs(real(lambda)-1.d0).lt.tol6).and.(abs(aimag(lambda)).lt.tol6).and.(inv==1)).or.&
1625 14976 : & ((abs(real(lambda)-1.d0).lt.tol6).and.(abs(aimag(lambda)).lt.tol6).and.(inv==2).and.(ii==jj))) cycle
1626 9044 : unchanged(isyminv)=.true.
1627 9044 : iconst(isyminv)=iconst(isyminv)+1
1628 : !FB const_tot=const_tot+1
1629 : !FB write(16,*)' The product of eigenvalues',ii,jj
1630 : !FB write(16,*)' is equal to ',lambda
1631 36176 : do mu=1,3
1632 117572 : do nu=1,3
1633 108528 : pp(mu,nu)=eigenvectors(mu,ii)*eigenvectors(nu,jj)
1634 : end do
1635 : end do
1636 41168 : do mu=1,3
1637 123504 : do nu=1,3
1638 108528 : if (inv==1) then
1639 30312 : alphaij(isyminv,(mu-1)*3+nu,iconst(isyminv))=pp(mu,nu)
1640 51084 : else if (inv==2) then
1641 51084 : alphaij(isyminv,(mu-1)*3+nu,iconst(isyminv))=lambda*pp(mu,nu)-pp(nu,mu)
1642 : else
1643 0 : ABI_BUG('This symetry is neither Keptinvariant nor Reversed')
1644 : end if
1645 : end do
1646 : end do
1647 : !FB write(16,*)' Real & imaginary parts of the eigenvectors product:'
1648 : !FB write(16,'(9(f16.12,1x))') real(alphaij(isyminv,:,iconst(isyminv)))
1649 : !FB write(16,'(9(f16.12,1x))') aimag(alphaij(isyminv,:,iconst(isyminv)))
1650 : end do !jj
1651 : end do !ii
1652 12720 : else if (order==3) then
1653 15064 : do ii=1,3
1654 48958 : do jj=1,3
1655 146874 : do kk=1,3
1656 101682 : lambda=eigenvalues(ii)*eigenvalues(jj)*eigenvalues(kk)
1657 : if (((abs(real(lambda)-1.d0).lt.tol6).and.(abs(aimag(lambda)).lt.tol6).and.(inv==1)).or.&
1658 : & ((abs(real(lambda)-1.d0).lt.tol6).and.(abs(aimag(lambda)).lt.tol6).and.(inv==2).and.(jj==kk)).or.&
1659 : & ((abs(real(lambda)-1.d0).lt.tol6).and.(abs(aimag(lambda)).lt.tol6).and.(inv==3).and.(ii==jj)).or.&
1660 : & ((abs(real(lambda)-1.d0).lt.tol6).and.(abs(aimag(lambda)).lt.tol6).and.(inv==4).and.(ii==jj).and.(jj==kk)).or.&
1661 101682 : & ((abs(real(lambda)-1.d0).lt.tol6).and.(abs(aimag(lambda)).lt.tol6).and.(inv==5).and.(ii==jj).and.(jj==kk)).or.&
1662 101682 : & ((abs(real(lambda)-1.d0).lt.tol6).and.(abs(aimag(lambda)).lt.tol6).and.(inv==6).and.(ii==kk))) cycle
1663 88073 : unchanged(isyminv)=.true.
1664 88073 : iconst(isyminv)=iconst(isyminv)+1
1665 : !FB const_tot=const_tot+1
1666 : !FB write(16,*)' The product of eigenvalues',ii,jj
1667 : !FB write(16,*)' is equal to ',lambda
1668 352292 : do mu=1,3
1669 1144949 : do nu=1,3
1670 3434847 : do xi=1,3
1671 3170628 : ppp(mu,nu,xi)=eigenvectors(mu,ii)*eigenvectors(nu,jj)*eigenvectors(xi,kk)
1672 : end do !xi
1673 : end do !nu
1674 : end do !mu
1675 386186 : do mu=1,3
1676 1158558 : do nu=1,3
1677 3434847 : do xi=1,3
1678 792657 : if (inv==1) then
1679 330534 : alphaij(isyminv,(mu-1)*9+(nu-1)*3+xi,iconst(isyminv))=ppp(mu,nu,xi)
1680 : else if (inv==2) then
1681 407646 : alphaij(isyminv,(mu-1)*9+(nu-1)*3+xi,iconst(isyminv))=lambda*ppp(mu,nu,xi)-ppp(mu,xi,nu)
1682 : else if (inv==3) then
1683 428328 : alphaij(isyminv,(mu-1)*9+(nu-1)*3+xi,iconst(isyminv))=lambda*ppp(mu,nu,xi)-ppp(nu,mu,xi)
1684 : else if (inv==4) then
1685 410940 : alphaij(isyminv,(mu-1)*9+(nu-1)*3+xi,iconst(isyminv))=lambda*ppp(mu,nu,xi)-ppp(nu,xi,mu)
1686 : else if (inv==5) then
1687 410940 : alphaij(isyminv,(mu-1)*9+(nu-1)*3+xi,iconst(isyminv))=lambda*ppp(mu,nu,xi)-ppp(xi,mu,nu)
1688 : else if (inv==6) then
1689 389583 : alphaij(isyminv,(mu-1)*9+(nu-1)*3+xi,iconst(isyminv))=lambda*ppp(mu,nu,xi)-ppp(xi,nu,mu)
1690 : else
1691 0 : ABI_BUG('This symetry is neither Keptinvariant nor Reversed')
1692 : end if
1693 : end do !xi
1694 : end do !nu
1695 : end do !mu
1696 : !FB write(16,*)' Real & imaginary parts of the eigenvectors product:'
1697 : !FB write(16,'(27(f16.12,1x))') real(alphaij(isyminv,:,iconst(isyminv)))
1698 : !FB write(16,'(27(f16.12,1x))') aimag(alphaij(isyminv,:,iconst(isyminv)))
1699 : end do !kk
1700 : end do !jj
1701 : end do !ii
1702 8954 : else if (order==4) then
1703 35816 : do ii=1,3
1704 116402 : do jj=1,3
1705 349206 : do kk=1,3
1706 1047618 : do ll=1,3
1707 725274 : lambda=eigenvalues(ii)*eigenvalues(jj)*eigenvalues(kk)*eigenvalues(ll)
1708 725274 : if ((abs(real(lambda)-1.d0).lt.tol6).and.(abs(aimag(lambda)).lt.tol6)) then
1709 : if ((inv==1 ) .or.& !\Phi4_ijkl
1710 : & ((inv==2 ).and.(jj==kk)) .or.& !\Phi4_ikjl
1711 : & ((inv==3 ).and.(ii==jj)) .or.& !\Phi4_jikl
1712 : & ((inv==4 ).and.(ii==jj).and.(jj==kk)) .or.& !\Phi4_jkil
1713 : & ((inv==5 ).and.(ii==jj).and.(jj==kk)) .or.& !\Phi4_kijl
1714 : & ((inv==6 ).and.(ii==kk)) .or.& !\Phi4_kjil
1715 :
1716 : & ((inv==7 ).and.(kk==ll)) .or.& !\Phi4_ijlk
1717 : & ((inv==8 ).and.(jj==kk).and.(kk==ll)) .or.& !\Phi4_iklj
1718 : & ((inv==9 ).and.(ii==jj).and.(kk==ll)) .or.& !\Phi4_jilk
1719 : & ((inv==10).and.(ii==jj).and.(jj==kk).and.(kk==ll)).or.& !\Phi4_jkli
1720 : & ((inv==11).and.(ii==jj).and.(jj==kk).and.(kk==ll)).or.& !\Phi4_kilj
1721 : & ((inv==12).and.(ii==kk).and.(kk==ll)) .or.& !\Phi4_kjli
1722 :
1723 : & ((inv==13).and.(jj==kk).and.(kk==ll)) .or.& !\Phi4_iljk
1724 : & ((inv==14).and.(jj==ll)) .or.& !\Phi4_ilkj
1725 : & ((inv==15).and.(ii==jj).and.(jj==kk).and.(kk==ll)).or.& !\Phi4_jlik
1726 : & ((inv==16).and.(ii==jj).and.(jj==ll)) .or.& !\Phi4_jlki
1727 : & ((inv==17).and.(ii==kk).and.(jj==ll)) .or.& !\Phi4_klij
1728 : & ((inv==18).and.(ii==jj).and.(jj==kk).and.(kk==ll)).or.& !\Phi4_klji
1729 :
1730 : & ((inv==19).and.(ii==jj).and.(jj==kk).and.(kk==ll)).or.& !\Phi4_lijk
1731 : & ((inv==20).and.(ii==jj).and.(jj==ll)) .or.& !\Phi4_likj
1732 : & ((inv==21).and.(ii==kk).and.(kk==ll)) .or.& !\Phi4_ljik
1733 : & ((inv==22).and.(ii==ll)) .or.& !\Phi4_ljki
1734 290830 : & ((inv==23).and.(ii==jj).and.(jj==kk).and.(kk==ll)).or.& !\Phi4_lkij
1735 : & ((inv==24).and.(ii==ll).and.(jj==kk))) cycle !\Phi4_lkji
1736 : end if
1737 660368 : unchanged(isyminv)=.true.
1738 660368 : iconst(isyminv)=iconst(isyminv)+1
1739 : !FB const_tot=const_tot+1
1740 : !FB write(16,*)' The product of eigenvalues',ii,jj
1741 : !FB write(16,*)' is equal to ',lambda
1742 2641472 : do mu=1,3
1743 8584784 : do nu=1,3
1744 25754352 : do xi=1,3
1745 77263056 : do zeta=1,3
1746 71319744 : pppp(mu,nu,xi,zeta)=eigenvectors(mu,ii)*eigenvectors(nu,jj)*eigenvectors(xi,kk)*eigenvectors(zeta,ll)
1747 : end do !zeta
1748 : end do !xi
1749 : end do !nu
1750 : end do !mu
1751 2883230 : do mu=1,3
1752 8649690 : do nu=1,3
1753 25754352 : do xi=1,3
1754 77263056 : do zeta=1,3
1755 53489808 : itemp=(mu-1)*27+(nu-1)*9+(xi-1)*3+zeta
1756 19560096 : if (inv==1) then ; alphaij(isyminv,itemp,iconst(isyminv))=pppp(mu,nu,xi,zeta)
1757 2269458 : else if (inv==2 ) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(mu,xi,nu,zeta)
1758 2400354 : else if (inv==3 ) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(nu,mu,xi,zeta)
1759 2396628 : else if (inv==4 ) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(nu,xi,mu,zeta)
1760 2396628 : else if (inv==5 ) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(xi,mu,nu,zeta)
1761 2175660 : else if (inv==6 ) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(xi,nu,mu,zeta)
1762 :
1763 2312226 : else if (inv==7 ) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(mu,nu,zeta,xi)
1764 2260548 : else if (inv==8 ) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(mu,xi,zeta,nu)
1765 2392416 : else if (inv==9 ) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(nu,mu,zeta,xi)
1766 2211462 : else if (inv==10) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(nu,xi,zeta,mu)
1767 2211462 : else if (inv==11) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(xi,mu,zeta,nu)
1768 2173554 : else if (inv==12) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(xi,nu,zeta,mu)
1769 :
1770 2260548 : else if (inv==13) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(mu,zeta,nu,xi)
1771 2050110 : else if (inv==14) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(mu,zeta,xi,nu)
1772 2211462 : else if (inv==15) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(nu,zeta,mu,xi)
1773 2167236 : else if (inv==16) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(nu,zeta,xi,mu)
1774 2278692 : else if (inv==17) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(xi,zeta,mu,nu)
1775 2388366 : else if (inv==18) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(xi,zeta,nu,mu)
1776 :
1777 2211462 : else if (inv==19) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(zeta,mu,nu,xi)
1778 2167236 : else if (inv==20) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(zeta,mu,xi,nu)
1779 2173554 : else if (inv==21) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(zeta,nu,mu,xi)
1780 1977696 : else if (inv==22) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(zeta,nu,xi,mu)
1781 2388366 : else if (inv==23) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(zeta,xi,mu,nu)
1782 2284524 : else if (inv==24) then ; alphaij(isyminv,itemp,iconst(isyminv))=lambda*pppp(mu,nu,xi,zeta)-pppp(zeta,xi,nu,mu)
1783 0 : else ; ABI_BUG('This symetry is neither Keptinvariant nor Reversed')
1784 : end if
1785 : end do !zeta
1786 : end do !xi
1787 : end do !nu
1788 : end do !mu
1789 : !FB write(16,*)' Real & imaginary parts of the eigenvectors product:'
1790 : !FB write(16,'(81(f16.12,1x))') real(alphaij(isyminv,:,iconst(isyminv)))
1791 : !FB write(16,'(81(f16.12,1x))') aimag(alphaij(isyminv,:,iconst(isyminv)))
1792 : end do !ll
1793 : end do !kk
1794 : end do !jj
1795 : end do !ii
1796 : else
1797 0 : ABI_BUG('Only the first, second, third and fourth order are allowed')
1798 : end if
1799 :
1800 :
1801 : !FB=================================================================
1802 : !FB======== TO CLEAN ===============================================
1803 : !FB=================================================================
1804 14491 : nconst_loc=const_tot+iconst(isyminv)
1805 14491 : ii=0
1806 115825604 : ABI_MALLOC(tab_vec,(norder,nconst_loc)); tab_vec(:,:)=czero
1807 5620087 : do itemp=1,isyminv
1808 5620087 : if (unchanged(itemp)) then
1809 6174477 : do jj=1,iconst(itemp)
1810 1549045 : ii=ii+1
1811 120393072 : tab_vec(:,ii)=alphaij(itemp,:,jj)
1812 : end do
1813 : end if
1814 : end do
1815 :
1816 1549051 : do kk=2,nconst_loc
1817 107409811 : do jj=1,kk-1
1818 8428220064 : prod_scal=sum( real(tab_vec(:,jj))* real(tab_vec(:,jj))+aimag(tab_vec(:,jj))*aimag(tab_vec(:,jj)))
1819 107395320 : if (abs(prod_scal).gt.tol8) then
1820 12754247719 : tab_vec(:,kk)=tab_vec(:,kk)-sum(tab_vec(:,kk)*conjg(tab_vec(:,jj)))/dcmplx(prod_scal,zero)*tab_vec(:,jj)
1821 6417457066 : do ii=1,norder
1822 6336790653 : if (abs( real(tab_vec(ii,kk))).lt.tol8) tab_vec(ii,kk)=dcmplx(zero,aimag(tab_vec(ii,kk)))
1823 6417457066 : if (abs(aimag(tab_vec(ii,kk))).lt.tol8) tab_vec(ii,kk)=dcmplx( real(tab_vec(ii,kk)),zero)
1824 : end do
1825 : end if
1826 : end do
1827 : end do
1828 :
1829 : ! On stocke les vecteurs non-nuls
1830 115811113 : ABI_MALLOC(temp ,(norder,nconst_loc)); temp(:,:) =czero
1831 14491 : ii=0
1832 1563536 : do kk=1,nconst_loc
1833 115767640 : prod_scal=sum( real(tab_vec(:,kk))* real(tab_vec(:,kk))+aimag(tab_vec(:,kk))*aimag(tab_vec(:,kk)))
1834 1563536 : if (abs(prod_scal).gt.tol8) then
1835 797602 : ii=ii+1
1836 59370346 : temp(:,ii)=tab_vec(:,kk)/dsqrt(prod_scal)
1837 : end if
1838 : end do
1839 14491 : ABI_FREE(tab_vec)
1840 14491 : iconst(isyminv)=ii-const_tot
1841 14491 : const_tot=const_tot+iconst(isyminv)
1842 :
1843 14491 : ii=0
1844 66176433313 : alphaij(:,:,:)=czero
1845 5620087 : do itemp=1,isyminv
1846 5620087 : if (unchanged(itemp)) then
1847 5423034 : do jj=1,iconst(itemp)
1848 797602 : ii=ii+1
1849 63995778 : alphaij(itemp,:,jj)=temp(:,ii)
1850 : end do
1851 : end if
1852 : end do
1853 14491 : ABI_FREE(temp)
1854 : !FB=================================================================
1855 : !FB======== TO CLEAN ===============================================
1856 : !FB=================================================================
1857 :
1858 :
1859 :
1860 :
1861 :
1862 :
1863 : ! WARNING: There are some minimum and maximum of constraints
1864 14491 : if (order==1.and.(iconst(isyminv).eq.3)) then
1865 6 : ncoeff=0
1866 78 : proj(:,:,ishell)=zero
1867 6 : ABI_FREE(unchanged)
1868 6 : ABI_FREE(alphaij)
1869 6 : ABI_FREE(iconst)
1870 6 : return
1871 14485 : else if (order==1.and.(iconst(isyminv).gt.3)) then
1872 0 : ABI_BUG(' First order : There are more than 3 constraints')
1873 : end if
1874 14485 : if (order==2.and.(iconst(isyminv).gt.8)) then
1875 0 : ABI_BUG(' Second order : There are more than 8 constraints')
1876 : end if
1877 14485 : if (order==3.and.(iconst(isyminv).gt.27)) then
1878 0 : ABI_BUG(' Third order : There are more than 27 constraints')
1879 : end if
1880 15030 : if (order==4.and.(iconst(isyminv).gt.81)) then
1881 0 : ABI_BUG(' Fourth order : There are more than 81 constraints')
1882 : end if
1883 : end do !isyminv
1884 : ! ================================================================================================
1885 : ! =========== End big loop over symetries and facorder ===========================================
1886 : ! ================================================================================================
1887 545 : nconst_perm=0
1888 : ! The (iik, iji, ijj and iii) third order IFCs are symmetric with respect to some permutations.
1889 : ! Some constraints have to be added :
1890 545 : if (order.eq.3) then
1891 49 : if ((iatcell.eq.jatom).or.(iatcell.eq.katom).or.(jatom.eq.katom)) then
1892 39 : nconst_perm=5
1893 39 : if (MPIdata%iam_master) write(16,'(a)')'=========== The IFCs are symmetric'
1894 39 : const_tot=const_tot+nconst_perm*norder
1895 171795 : ABI_MALLOC(constraints,(nconst_perm,norder,norder)) ; constraints(:,:,:)=czero
1896 39 : ii=0
1897 156 : do mu=1,3
1898 517 : do nu=1,3
1899 1521 : do xi=1,3
1900 1053 : ii=ii+1
1901 1053 : if (iatcell.eq.jatom) then
1902 918 : if (mu.eq.nu) cycle
1903 612 : constraints(1,(mu-1)*9+(nu-1)*3+xi,ii)= cone
1904 612 : constraints(1,(nu-1)*9+(mu-1)*3+xi,ii)=-cone
1905 : end if
1906 747 : if (iatcell.eq.katom) then
1907 342 : if (mu.eq.xi) cycle
1908 228 : constraints(2,(mu-1)*9+(nu-1)*3+xi,ii)= cone
1909 228 : constraints(2,(xi-1)*9+(nu-1)*3+mu,ii)=-cone
1910 : end if
1911 633 : if (jatom.eq.katom) then
1912 363 : if (nu.eq.xi) cycle
1913 204 : constraints(3,(mu-1)*9+(nu-1)*3+xi,ii)= cone
1914 204 : constraints(3,(mu-1)*9+(xi-1)*3+nu,ii)=-cone
1915 : end if
1916 474 : if ((iatcell.eq.jatom).and.(jatom.eq.katom)) then
1917 114 : if ((nu.eq.xi).and.(nu.eq.mu)) cycle
1918 114 : constraints(4,(mu-1)*9+(nu-1)*3+xi,ii)= cone
1919 114 : constraints(4,(xi-1)*9+(mu-1)*3+nu,ii)=-cone
1920 : end if
1921 825 : if ((iatcell.eq.jatom).and.(jatom.eq.katom)) then
1922 114 : if ((nu.eq.xi).and.(nu.eq.mu)) cycle
1923 114 : constraints(5,(mu-1)*9+(nu-1)*3+xi,ii)= cone
1924 114 : constraints(5,(nu-1)*9+(xi-1)*3+mu,ii)=-cone
1925 : end if
1926 : end do
1927 : end do
1928 : end do
1929 : end if
1930 : end if
1931 :
1932 : ! The (iikl, ijil, ijki, ijjl, ijkj, ijkk, iiil, iiki, ijii, ijjj, iiii)
1933 : ! fourth order IFCs are symmetric with respect to some permutations.
1934 : ! Some constraints have to be added :
1935 545 : if (order.eq.4) then
1936 : if ((iatcell.eq.jatom).or.(iatcell.eq.katom).or.(iatcell.eq.latom)&
1937 46 : & .or.(jatom.eq.katom).or.(jatom.eq.latom).or.(katom.eq.latom)) then
1938 44 : nconst_perm=17
1939 44 : if (MPIdata%iam_master) write(16,'(a)')'=========== The IFCs are symmetric'
1940 44 : const_tot=const_tot+nconst_perm*norder
1941 5200052 : ABI_MALLOC(constraints,(nconst_perm,norder,norder)) ; constraints(:,:,:)=czero
1942 44 : ii=0
1943 176 : do mu=1,3
1944 574 : do nu=1,3
1945 1716 : do xi=1,3
1946 5148 : do zeta=1,3
1947 3564 : ii=ii+1
1948 3564 : if (iatcell.eq.jatom) then
1949 2997 : if (mu.eq.nu) cycle
1950 1998 : constraints(1,(mu-1)*27+(nu-1)*9+(xi-1)*3+zeta,ii)= cone
1951 1998 : constraints(1,(nu-1)*27+(mu-1)*9+(xi-1)*3+zeta,ii)=-cone
1952 : end if
1953 2565 : if (iatcell.eq.katom) then
1954 1080 : if (mu.eq.xi) cycle
1955 720 : constraints(2,(mu-1)*27+(nu-1)*9+(xi-1)*3+zeta,ii)= cone
1956 720 : constraints(2,(xi-1)*27+(nu-1)*9+(mu-1)*3+zeta,ii)=-cone
1957 : end if
1958 2205 : if (iatcell.eq.latom) then
1959 396 : if (mu.eq.zeta) cycle
1960 264 : constraints(3,(mu -1)*27+(nu-1)*9+(xi-1)*3+zeta,ii)= cone
1961 264 : constraints(3,(zeta-1)*27+(nu-1)*9+(xi-1)*3+mu ,ii)=-cone
1962 : end if
1963 2073 : if (jatom.eq.katom) then
1964 1155 : if (nu.eq.xi) cycle
1965 672 : constraints(4,(mu-1)*27+(nu-1)*9+(xi-1)*3+zeta,ii)= cone
1966 672 : constraints(4,(mu-1)*27+(xi-1)*9+(nu-1)*3+zeta,ii)=-cone
1967 : end if
1968 1590 : if (jatom.eq.latom) then
1969 294 : if (nu.eq.zeta) cycle
1970 174 : constraints(5,(mu-1)*27+(nu -1)*9+(xi-1)*3+zeta,ii)= cone
1971 174 : constraints(5,(mu-1)*27+(zeta-1)*9+(xi-1)*3+nu ,ii)=-cone
1972 : end if
1973 1470 : if (katom.eq.latom) then
1974 660 : if (xi.eq.zeta) cycle
1975 378 : constraints(6,(mu-1)*27+(nu-1)*9+(xi -1)*3+zeta,ii)= cone
1976 378 : constraints(6,(mu-1)*27+(nu-1)*9+(zeta-1)*3+xi ,ii)=-cone
1977 : end if
1978 :
1979 1188 : if ((iatcell.eq.jatom).and.(jatom.eq.katom)) then
1980 162 : if ((mu.eq.nu).and.(nu.eq.xi)) cycle
1981 162 : constraints(7,(mu-1)*27+(nu-1)*9+(xi-1)*3+zeta,ii)= cone
1982 162 : constraints(7,(xi-1)*27+(mu-1)*9+(nu-1)*3+zeta,ii)=-cone
1983 : end if
1984 1188 : if ((iatcell.eq.jatom).and.(jatom.eq.katom)) then
1985 162 : if ((mu.eq.nu).and.(nu.eq.xi)) cycle
1986 162 : constraints(8,(mu-1)*27+(nu-1)*9+(xi-1)*3+zeta,ii)= cone
1987 162 : constraints(8,(nu-1)*27+(xi-1)*9+(mu-1)*3+zeta,ii)=-cone
1988 : end if
1989 :
1990 1188 : if ((iatcell.eq.jatom).and.(jatom.eq.latom)) then
1991 0 : if ((mu.eq.nu).and.(nu.eq.zeta)) cycle
1992 0 : constraints(9,(mu-1)*27+(nu -1)*9+(xi-1)*3+zeta,ii)= cone
1993 0 : constraints(9,(nu-1)*27+(zeta-1)*9+(xi-1)*3+mu ,ii)=-cone
1994 : end if
1995 1188 : if ((iatcell.eq.jatom).and.(jatom.eq.latom)) then
1996 0 : if ((mu.eq.nu).and.(nu.eq.zeta)) cycle
1997 0 : constraints(10,(mu -1)*27+(nu-1)*9+(xi-1)*3+zeta,ii)= cone
1998 0 : constraints(10,(zeta-1)*27+(mu-1)*9+(xi-1)*3+nu ,ii)=-cone
1999 : end if
2000 :
2001 1188 : if ((iatcell.eq.katom).and.(katom.eq.latom)) then
2002 0 : if ((mu.eq.xi).and.(xi.eq.zeta)) cycle
2003 0 : constraints(11,(mu-1)*27+(nu-1)*9+(xi -1)*3+zeta,ii)= cone
2004 0 : constraints(11,(xi-1)*27+(nu-1)*9+(zeta-1)*3+mu ,ii)=-cone
2005 : end if
2006 1188 : if ((iatcell.eq.katom).and.(katom.eq.latom)) then
2007 0 : if ((mu.eq.xi).and.(xi.eq.zeta)) cycle
2008 0 : constraints(12,(mu -1)*27+(nu-1)*9+(xi-1)*3+zeta,ii)= cone
2009 0 : constraints(12,(zeta-1)*27+(nu-1)*9+(mu-1)*3+xi ,ii)=-cone
2010 : end if
2011 :
2012 1188 : if ((jatom.eq.katom).and.(katom.eq.latom)) then
2013 54 : if ((nu.eq.xi).and.(xi.eq.zeta)) cycle
2014 54 : constraints(13,(mu-1)*27+(nu-1)*9+(xi -1)*3+zeta,ii)= cone
2015 54 : constraints(13,(mu-1)*27+(xi-1)*9+(zeta-1)*3+nu ,ii)=-cone
2016 : end if
2017 1188 : if ((jatom.eq.katom).and.(katom.eq.latom)) then
2018 54 : if ((nu.eq.xi).and.(xi.eq.zeta)) cycle
2019 54 : constraints(14,(mu-1)*27+(nu -1)*9+(xi-1)*3+zeta,ii)= cone
2020 54 : constraints(14,(mu-1)*27+(zeta-1)*9+(nu-1)*3+xi ,ii)=-cone
2021 : end if
2022 :
2023 1188 : if ((iatcell.eq.jatom).and.(jatom.eq.katom).and.(katom.eq.latom)) then
2024 0 : if ((mu.eq.nu).and.(nu.eq.xi).and.(xi.eq.zeta)) cycle
2025 0 : constraints(15,(mu-1)*27+(nu-1)*9+(xi -1)*3+zeta,ii)= cone
2026 0 : constraints(15,(nu-1)*27+(xi-1)*9+(zeta-1)*3+mu ,ii)=-cone
2027 : end if
2028 1188 : if ((iatcell.eq.jatom).and.(jatom.eq.katom).and.(katom.eq.latom)) then
2029 0 : if ((mu.eq.nu).and.(nu.eq.xi).and.(xi.eq.zeta)) cycle
2030 0 : constraints(16,(mu-1)*27+(nu -1)*9+(xi-1)*3+zeta,ii)= cone
2031 0 : constraints(16,(xi-1)*27+(zeta-1)*9+(mu-1)*3+nu ,ii)=-cone
2032 : end if
2033 2376 : if ((iatcell.eq.jatom).and.(jatom.eq.katom).and.(katom.eq.latom)) then
2034 0 : if ((mu.eq.nu).and.(nu.eq.xi).and.(xi.eq.zeta)) cycle
2035 0 : constraints(17,(mu -1)*27+(nu-1)*9+(xi-1)*3+zeta,ii)= cone
2036 0 : constraints(17,(zeta-1)*27+(mu-1)*9+(nu-1)*3+xi ,ii)=-cone
2037 : end if
2038 : end do
2039 : end do
2040 : end do
2041 : end do
2042 : end if
2043 : end if
2044 :
2045 : ! In the case where the matrix has norder**2 inequivalent and non-zero elements
2046 545 : if (const_tot==0) then
2047 47 : write(message,'(a,1x,i3,1x,a)') 'For shell number=',ishell,'there is no symetry operation reducing the number of coefficients'
2048 47 : ABI_WARNING(message)
2049 4277 : proj(:,:,ishell)=0.d0
2050 470 : do ii=1,norder
2051 470 : proj(ii,ii,ishell)=1.d0
2052 : end do
2053 47 : ncoeff=norder
2054 47 : ABI_FREE(unchanged)
2055 47 : ABI_FREE(alphaij)
2056 47 : ABI_FREE(iconst)
2057 47 : return
2058 : end if
2059 :
2060 : ! When some constraints have been found
2061 498 : ncount=const_tot
2062 498 : if (MPIdata%iam_master) then
2063 498 : write(16,'(a,1x,i7,1x,a)') 'There is a total of ',ncount,' non-independant constraints for this shell'
2064 : end if
2065 498 : ii=0
2066 5427766 : ABI_MALLOC(tab_vec,(norder,ncount)); tab_vec(:,:)=czero
2067 5427268 : ABI_MALLOC(temp ,(norder,ncount)); temp(:,:) =czero
2068 80376 : do isyminv=1,nsyminv
2069 80376 : if (unchanged(isyminv)) then
2070 20714 : do jj=1,iconst(isyminv)
2071 6235 : ii=ii+1
2072 324617 : tab_vec(:,ii)=alphaij(isyminv,:,jj)
2073 : end do
2074 : end if
2075 : end do
2076 498 : ABI_FREE(unchanged)
2077 498 : ABI_FREE(alphaij)
2078 498 : ABI_FREE(iconst)
2079 : ! Add the constraints coming from the symmetry of the IFCs (at the 3rd order)
2080 498 : if (nconst_perm.gt.0) then
2081 4700 : do jj=1,norder
2082 70553 : do kk=1,nconst_perm
2083 65853 : ii=ii+1
2084 5120253 : tab_vec(:,ii)=constraints(kk,:,jj)
2085 : end do
2086 : end do
2087 83 : ABI_FREE(constraints)
2088 : end if
2089 498 : if (ii.ne.ncount) then
2090 0 : write(message,'(i7,1x,a,1x,i7)') ii,' non equal to ',ncount
2091 0 : ABI_BUG(message)
2092 : end if
2093 9096 : do ii=1,norder
2094 5362782 : do jj=1,ncount
2095 5353686 : if (abs( real(tab_vec(ii,jj))).lt.tol8) tab_vec(ii,jj)=dcmplx(zero,aimag(tab_vec(ii,jj)))
2096 5362284 : if (abs(aimag(tab_vec(ii,jj))).lt.tol8) tab_vec(ii,jj)=dcmplx( real(tab_vec(ii,jj)),zero)
2097 : end do
2098 : end do
2099 :
2100 : ! On stocke les vecteurs non-nuls
2101 498 : ii=0
2102 72586 : do kk=1,ncount
2103 5425774 : prod_scal=sum( real(tab_vec(:,kk))* real(tab_vec(:,kk))+aimag(tab_vec(:,kk))*aimag(tab_vec(:,kk)))
2104 72586 : if (abs(prod_scal).gt.tol8) then
2105 12145 : ii=ii+1
2106 726070 : temp(:,ii)=tab_vec(:,kk)/dsqrt(prod_scal)
2107 : end if
2108 : end do
2109 498 : ncount=ii
2110 498 : ABI_FREE(tab_vec)
2111 728062 : ABI_MALLOC(tab_vec,(norder,ncount)); tab_vec(:,1:ncount)=temp(:,1:ncount)
2112 498 : ABI_FREE(temp)
2113 727564 : ABI_MALLOC(temp ,(norder,ncount)); temp(:,:) =czero
2114 :
2115 : ! L'ensemble des vecteurs reduisants l'espace de R^norder a R^n ne forment pas une base
2116 : ! independante. Il faut donc trouver les vecteurs independants.
2117 : ! --> Orthogonalisation de Gram-Schmidt
2118 12145 : do kk=2,ncount
2119 801548 : do jj=1,kk-1
2120 60684808 : prod_scal=sum( real(tab_vec(:,jj))* real(tab_vec(:,jj))+aimag(tab_vec(:,jj))*aimag(tab_vec(:,jj)))
2121 801050 : if (abs(prod_scal).gt.tol8) then
2122 75080365 : tab_vec(:,kk)=tab_vec(:,kk)-sum(tab_vec(:,kk)*conjg(tab_vec(:,jj)))/dcmplx(prod_scal,zero)*tab_vec(:,jj)
2123 37787200 : do ii=1,norder
2124 37293165 : if (abs( real(tab_vec(ii,kk))).lt.tol8) tab_vec(ii,kk)=dcmplx(zero,aimag(tab_vec(ii,kk)))
2125 37787200 : if (abs(aimag(tab_vec(ii,kk))).lt.tol8) tab_vec(ii,kk)=dcmplx( real(tab_vec(ii,kk)),zero)
2126 : end do
2127 : !FB else
2128 : !FB write(Invar%stdout,*)'One prod_scal equals zero'
2129 : end if
2130 : end do
2131 : end do
2132 :
2133 : ! On stocke les vecteurs non-nuls
2134 498 : ii=0
2135 12643 : do kk=1,ncount
2136 726070 : prod_scal=sum( real(tab_vec(:,kk))* real(tab_vec(:,kk))+aimag(tab_vec(:,kk))*aimag(tab_vec(:,kk)))
2137 12643 : if (abs(prod_scal).gt.tol8) then
2138 6262 : ii=ii+1
2139 312352 : temp(:,ii)=tab_vec(:,kk)/dsqrt(prod_scal)
2140 : end if
2141 : end do
2142 498 : ncount=ii
2143 498 : ABI_FREE(tab_vec)
2144 :
2145 : ! On ecrit les vecteurs non-nuls
2146 : !FB write(16,*) ' '
2147 : !FB write(16,*) ' ========The final set of vectors is:'
2148 : !FB do kk=1,ncount
2149 : !FB write(16,'(81(f16.12,1x))') real(temp(:,kk))
2150 : !FB write(16,'(81(f16.12,1x))') aimag(temp(:,kk))
2151 : !FB end do
2152 498 : if (MPIdata%iam_master) then
2153 498 : write(16,'(a,1x,i7,1x,a)') ' ======= Finally, there are ',ncount,' independent vectors'
2154 : end if
2155 498 : if (ncount.gt.8.and.order==2) then
2156 0 : ABI_ERROR(' Order 2 : There are too many independent vectors')
2157 : end if
2158 498 : if (ncount.gt.27.and.order==3) then
2159 0 : ABI_ERROR(' Order 3 : There are too many independent vectors')
2160 : end if
2161 498 : if (ncount.gt.81.and.order==4) then
2162 0 : ABI_ERROR(' Order 4 : There are too many independent vectors')
2163 : end if
2164 :
2165 : ! On cherche les (norder-ncount) vecteurs orthogonaux aux vecteurs non-nuls
2166 : ! --> Orthogonalisation de Gram-Schmidt
2167 379824 : ABI_MALLOC(tab_vec,(norder,norder)); tab_vec(:,:)=czero
2168 498 : iseed=-5
2169 9096 : do kk=1,norder
2170 9096 : if (kk.le.ncount) then
2171 312352 : tab_vec(:,kk)=temp(:,kk)
2172 : else
2173 65480 : do jj=1,norder
2174 63144 : drandom=uniformrandom(iseed)
2175 65480 : tab_vec(jj,kk)=dcmplx(drandom,zero)
2176 : end do
2177 54433 : do jj=1,kk-1
2178 3330896 : prod_scal=sum( real(tab_vec(:,jj))* real(tab_vec(:,jj))+aimag(tab_vec(:,jj))*aimag(tab_vec(:,jj)))
2179 52097 : if (abs(prod_scal).gt.tol8) then
2180 6609695 : tab_vec(:,kk)=tab_vec(:,kk)-sum(tab_vec(:,kk)*conjg(tab_vec(:,jj)))/prod_scal*tab_vec(:,jj)
2181 3330896 : do ii=1,norder
2182 3278799 : if (abs( real(tab_vec(ii,kk))).lt.tol8) tab_vec(ii,kk)=dcmplx(zero,aimag(tab_vec(ii,kk)))
2183 3330896 : if (abs(aimag(tab_vec(ii,kk))).lt.tol8) tab_vec(ii,kk)=dcmplx( real(tab_vec(ii,kk)),zero)
2184 : end do
2185 : end if
2186 3330896 : prod_scal=sum( real(tab_vec(:,kk))* real(tab_vec(:,kk))+aimag(tab_vec(:,kk))*aimag(tab_vec(:,kk)))
2187 3333232 : tab_vec(:,kk)=tab_vec(:,kk)/dsqrt(prod_scal)
2188 : end do
2189 : end if
2190 : end do
2191 498 : ABI_FREE(temp)
2192 :
2193 : ! On ecrit les vecteurs non-nuls
2194 : !FB write(16,*) ' '
2195 : !FB write(16,*) ' ========The orthogonal set of vectors is:'
2196 2834 : do kk=ncount+1,norder
2197 : !FB write(16,'(81(f16.12,1x))') real(tab_vec(:,kk))
2198 : !FB write(16,'(81(f16.12,1x))') aimag(tab_vec(:,kk))
2199 : if ((abs(aimag(tab_vec(1,kk))).gt.tol6).or.&
2200 2336 : & (abs(aimag(tab_vec(1,kk))).gt.tol6).or.&
2201 2834 : & (abs(aimag(tab_vec(1,kk))).gt.tol6)) then
2202 0 : ABI_ERROR('the constraint has an imaginary part')
2203 : end if
2204 : end do
2205 498 : ncoeff=norder-ncount
2206 498 : if (MPIdata%iam_master) then
2207 498 : write(16,'(a,1x,i7,1x,a)') ' ======= Finally, there are ',ncoeff,' coefficients'
2208 : end if
2209 :
2210 : ! On copie tab_vec dans proj
2211 2834 : do icoeff=1,ncoeff
2212 65978 : proj(:,icoeff,ishell)=tab_vec(:,ncount+icoeff)
2213 : end do
2214 498 : ABI_FREE(tab_vec)
2215 :
2216 658 : end subroutine tdep_calc_nbcoeff
2217 :
2218 : !====================================================================================================
2219 :
2220 0 : end module m_tdep_shell
|
