Branch data 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
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # #
# ]
|
