Branch data Line data Source code
1 : :
2 : : !****f* ABINIT/m_polynomal_coeff
3 : : !!
4 : : !! NAME
5 : : !! m_polynomial_coeff
6 : : !!
7 : : !! FUNCTION
8 : : !! Module with the datatype polynomial coefficients
9 : : !!
10 : : !! COPYRIGHT
11 : : !! Copyright (C) 2010-2026 ABINIT group (AM)
12 : : !! This file is distributed under the terms of the
13 : : !! GNU General Public Licence, see ~abinit/COPYING
14 : : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
16 : : !!
17 : : !! SOURCE
18 : :
19 : :
20 : : #if defined HAVE_CONFIG_H
21 : : #include "config.h"
22 : : #endif
23 : :
24 : : #include "abi_common.h"
25 : : !#include "trace.h"
26 : :
27 : :
28 : : module m_polynomial_coeff
29 : :
30 : : use defs_basis
31 : : use m_errors
32 : : use m_abicore
33 : : use m_polynomial_term
34 : : USE_MPI
35 : : use m_xmpi
36 : :
37 : : use m_sort, only : sort_dp
38 : : use m_io_tools, only : open_file, get_unit
39 : : use m_symtk, only : symchk, symatm
40 : : use m_crystal, only : crystal_t,symbols_crystal
41 : : use m_supercell, only : getPBCIndexes_supercell,distance_supercell,findBound_supercell
42 : : !use m_geometry, only : xcart2xred,metric
43 : : use m_dtfil, only : isfile
44 : : use m_hashtable_strval, only: hash_table_t
45 : : use m_dynamic_array, only: int2d_array_type
46 : :
47 : :
48 : : implicit none
49 : : private
50 : :
51 : : ! symmetric pairs and information about it.
52 : : type,public :: SymPairs_t
53 : : integer :: ncoeff_sym,nstr_sym
54 : : integer,allocatable :: list_symcoeff(:,:,:),list_symstr(:,:,:)
55 : : integer :: natom,nrpt, nsym
56 : : real(dp) :: cutoff
57 : : !arrays
58 : : type(crystal_t), pointer :: crystal
59 : : integer, allocatable :: cell(:,:)
60 : : real(dp), allocatable:: dist(:, :, :, :)
61 : : character(len=5),allocatable :: symbols(:)
62 : : integer :: sc_size(3)
63 : : real(dp):: range_ifc(3)
64 : : integer :: fit_iatom=-1
65 : : contains
66 : : procedure :: init => SymPairs_t_init
67 : : procedure :: free => SymPairs_t_free
68 : : procedure :: generateTerms => SymPairs_t_generateTerms
69 : : final :: SymPairs_t_final
70 : : end type SymPairs_T
71 : :
72 : : type, private :: symlist_t
73 : : integer(dp) :: nsym=0, power=0
74 : : integer(dp) :: max=0
75 : : integer(dp) :: counter=0
76 : : integer, allocatable :: list(:)
77 : : contains
78 : : procedure :: init=> symlist_init
79 : : procedure :: next => symlist_next
80 : : procedure :: free => symlist_free
81 : : end type symlist_t
82 : :
83 : : type, private :: IrreducibleCombinations_t
84 : : ! use the hash table to see if the term is already there. And the values of the table is -1 if the term is not bounding term.
85 : : type(hash_table_t) :: table
86 : : type(int2d_array_type) :: array
87 : : contains
88 : : procedure :: init => IrreducibleCombinations_init
89 : : procedure :: free => IrreducibleCombinations_free
90 : : procedure :: reset_array => IrreducibleCombinations_reset_array
91 : : procedure :: add_irr => IrreducibleCombinations_add_irr
92 : : end type IrreducibleCombinations_T
93 : :
94 : : type, private :: polyform_t
95 : : integer :: n=0
96 : : integer, allocatable :: ind(:) ! ind(n)
97 : : integer, allocatable :: order(:) ! orders(n)
98 : : integer, allocatable :: expansion(:) ! orders(power), where power=sum(order)
99 : : contains
100 : : procedure:: init => polyform_init
101 : : procedure:: free => polyform_free
102 : : procedure:: expand => polyform_expand
103 : : procedure:: from_expansion=> polyform_from_expansion
104 : : end type polyform_t
105 : :
106 : : public :: polynomial_coeff_broadcast
107 : : public :: polynomial_coeff_evaluate
108 : : public :: polynomial_coeff_free
109 : : public :: polynomial_coeff_init
110 : : public :: prepare_for_getList
111 : : public :: get_crystal_cutoff
112 : : public :: polynomial_coeff_getList
113 : : public :: polynomial_coeff_getName
114 : : public :: polynomial_coeff_getNorder
115 : : public :: polynomial_coeff_getOrder1
116 : : public :: polynomial_coeff_MPIrecv
117 : : public :: polynomial_coeff_MPIsend
118 : : public :: polynomial_coeff_setName
119 : : public :: polynomial_coeff_setCoefficient
120 : : public :: polynomial_coeff_writeXML
121 : : public :: polynomial_coeff_getEvenAnhaStrain
122 : : public :: coeffs_list_copy
123 : : public :: coeffs_list_conc
124 : : public :: coeffs_list_conc_onsite
125 : : public :: coeffs_list_append
126 : : public :: coeffs_list_truncate
127 : : public :: polynomial_coeff_list_free
128 : : public :: generateTermsFromList
129 : : public :: find_irpt
130 : : public :: coeffs_compare
131 : : public :: coeffs_list_reduce_duplicate
132 : : public :: operator(==)
133 : : public :: assignment(=)
134 : : public :: operator(+)
135 : :
136 : :
137 : : private :: computeNorder
138 : : private :: computeCombinationFromList
139 : : private :: computeSymmetricCombinations
140 : : !private :: computeSymmetricCombinations_old
141 : : private :: getCoeffFromList
142 : : private :: reduce_zero_combinations
143 : : !private :: check_irreducibility
144 : : private :: sort_combination_list
145 : : private :: sort_combination
146 : :
147 : :
148 : : !!***
149 : :
150 : : !!****t* m_polynomial_coeff/polynomial_coeff_type
151 : : !! NAME
152 : : !! polynomial_coeff_type
153 : : !!
154 : : !! FUNCTION
155 : : !! structure for a polynomial coefficient
156 : : !! contains the value of the coefficient and a
157 : : !! list of terms (displacements and/or strain) relating to the coefficient by symmetry
158 : : !!
159 : : !! SOURCE
160 : :
161 : : type, public :: polynomial_coeff_type
162 : :
163 : : character(len=200) :: name = ""
164 : : ! Name of the polynomial_coeff (Sr_y-O1_y)^3) for example
165 : :
166 : : integer :: nterm = 0
167 : : ! Number of terms (short range interaction) for this polynomial_coeff
168 : :
169 : : real(dp) :: coefficient = zero
170 : : ! coefficient = value of the coefficient of this term
171 : : ! \frac{\partial E^{k}}{\partial \tau^{k}}
172 : :
173 : : type(polynomial_term_type),dimension(:),allocatable :: terms
174 : : ! polynomial_term(nterm)<type(polynomial_term)>
175 : : ! contains all the displacements for this coefficient
176 : :
177 : : character(len=500) :: debug_str="uninitialized"
178 : :
179 : : contains
180 : : !final :: polynomial_coeff_final
181 : : end type polynomial_coeff_type
182 : : !!***
183 : :
184 : :
185 : : interface operator (==)
186 : : module procedure coeffs_compare
187 : : end interface operator (==)
188 : :
189 : : interface operator (+)
190 : : module procedure coeffs_list_conc
191 : : end interface operator (+)
192 : :
193 : : interface assignment (=)
194 : : module procedure coeffs_list_copy
195 : : end interface assignment (=)
196 : :
197 : :
198 : :
199 : : CONTAINS !===========================================================================================
200 : :
201 : :
202 : : !!****f* m_polynomial_coeff/polynomial_coeff_init
203 : : !!
204 : : !! NAME
205 : : !! polynomial_coeff_init
206 : : !!
207 : : !! FUNCTION
208 : : !! Initialize a polynomial_coeff datatype
209 : : !!
210 : : !! INPUTS
211 : : !! name = Name of the polynomial_coeff (Sr_y-O1_y)^3) for example
212 : : !! nterm = Number of terms (short range interaction) for this polynomial_coeff
213 : : !! coefficient = Value of the coefficient of this term
214 : : !! terms(nterm)<type(polynomial_term)> = array of polynomial_term_type
215 : : !! check = TRUE if this list of terms has to be check. We remove the symetric of equivalent terms
216 : : !! for example: ((Sr_y-O1_y)^1 and -1*(Sr_y-O1_y)^1 => zero
217 : : !! FALSE, defaut, do nothing
218 : : !! OUTPUT
219 : : !! polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype to be initialized
220 : : !!
221 : : !! SOURCE
222 : :
223 [ + + + - : 65 : subroutine polynomial_coeff_init(coefficient,nterm,polynomial_coeff,terms,name, check, debug_str)
+ - + - +
- + - + -
+ - ]
224 : :
225 : : !Arguments ------------------------------------
226 : : !scalars
227 : : integer, intent(in) :: nterm
228 : : real(dp),intent(in) :: coefficient
229 : : logical,optional,intent(in) :: check
230 : : !arrays
231 : : character(len=200),optional,intent(in) :: name
232 : : type(polynomial_term_type),intent(in) :: terms(nterm)
233 : : type(polynomial_coeff_type), intent(out) :: polynomial_coeff
234 : : !Local variables-------------------------------
235 : : !scalar
236 : : integer :: iterm1,iterm2
237 : : integer :: ii,nterm_tmp
238 : : real(dp):: coefficient_tmp
239 : : logical :: check_in
240 : : !arrays
241 : 53310 : real(dp) :: weights(nterm)
242 : : character(len=200) :: name_tmp
243 : : character(*), optional :: debug_str
244 : : ! *************************************************************************
245 : : !First free before initilisation
246 : : !call polynomial_coeff_free(polynomial_coeff)
247 : 26655 : check_in = .false.
248 [ + + ]: 26655 : if(present(check)) check_in = check
249 : :
250 [ + + ]: 26655 : if (present(debug_str)) then
251 [ + - ]: 1770 : polynomial_coeff%debug_str = debug_str
252 : : else
253 : 24885 : polynomial_coeff%debug_str = "initialized"
254 : : end if
255 : :
256 [ + + ]: 26655 : if(check_in)then
257 : : ! Check if the list of term is available or contains identical terms
258 : : ! in this case, remove all the not needed terms
259 : 800823 : nterm_tmp = 0
260 [ + + ]: 800823 : weights(:) = one
261 [ + + ]: 800823 : do iterm1=1,nterm
262 [ + + ]: 781301 : if(abs(weights(iterm1)) < tol16)cycle ! FIXME: do nothing?
263 : 663057 : weights(iterm1) = terms(iterm1)%weight
264 [ + + ]: 15005735 : do iterm2=iterm1+1,nterm
265 [ + + ]: 14342678 : if(abs(weights(iterm2)) < tol16)cycle
266 : : ! if the terms are identical we check the weight
267 [ + + ]: 14857002 : if(terms(iterm1)==terms(iterm2))then
268 : 118244 : weights(iterm1) = weights(iterm1) + terms(iterm2)%weight
269 : 118244 : weights(iterm2) = 0
270 : : end if
271 : : end do
272 [ + + ]: 682579 : if(abs(weights(iterm1)) > tol16) then
273 : 619796 : weights(iterm1)= anint(weights(iterm1)/weights(iterm1))
274 : : end if
275 : : end do
276 : :
277 : :
278 : : ! Count the number of terms
279 [ + + + + ]: 800823 : nterm_tmp=count(abs(weights) > tol16)
280 [ + + ]: 19522 : if (nterm_tmp ==0)then
281 : : coefficient_tmp = 0.0
282 : : else
283 : 17327 : coefficient_tmp = coefficient
284 : : end if
285 : : else
286 : 7133 : nterm_tmp = nterm
287 : 7133 : coefficient_tmp = coefficient
288 [ + + ]: 249545 : weights(:) = terms(:)%weight
289 : : end if!end Check
290 : :
291 [ + + ]: 26655 : if(present(name))then
292 : 17613 : name_tmp = name
293 : : else
294 : 9042 : name_tmp = ""
295 : : end if
296 : :
297 : : !Initilisation
298 : 26655 : polynomial_coeff%name = name_tmp
299 : 26655 : polynomial_coeff%nterm = nterm_tmp
300 : 26655 : polynomial_coeff%coefficient = coefficient_tmp
301 [ + + + - : 942173 : ABI_MALLOC(polynomial_coeff%terms,(polynomial_coeff%nterm))
+ - + + -
+ - + - +
+ + + - -
+ - + - +
- + - + -
+ - + ]
302 : 26655 : iterm1 = 0
303 [ + + ]: 1050368 : do ii = 1,nterm
304 [ + + ]: 1050368 : if(abs(weights(ii)) > tol16)then
305 : 862208 : iterm1 = iterm1 + 1
306 : 862208 : call polynomial_term_copy(terms(ii), polynomial_coeff%terms(iterm1))
307 : 862208 : if(sum(polynomial_coeff%terms(iterm1)%power_strain)>1) then
308 : : endif
309 : : end if
310 : : end do
311 [ + + ]: 26655 : end subroutine polynomial_coeff_init
312 : : !!***
313 : :
314 : : !!****f* m_polynomial_coeff/polynomial_coeff_free
315 : : !!
316 : : !! NAME
317 : : !! polynomial_coeff_free
318 : : !!
319 : : !! FUNCTION
320 : : !! Free polynomial_coeff datatype
321 : : !!
322 : : !! INPUTS
323 : : !! polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype
324 : : !!
325 : : !! OUTPUT
326 : : !! polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype
327 : : !!
328 : : !! SOURCE
329 : :
330 : 29362 : subroutine polynomial_coeff_free(polynomial_coeff)
331 : :
332 : : !Arguments ------------------------------------
333 : : !scalars
334 : : !arrays
335 : : type(polynomial_coeff_type), intent(inout) :: polynomial_coeff
336 : : !Local variables-------------------------------
337 : : !scalar
338 : : integer :: ii
339 : : !arrays
340 : :
341 : : ! *************************************************************************
342 : :
343 [ + + ]: 29362 : if(allocated(polynomial_coeff%terms))then
344 [ + + ]: 895935 : do ii = 1,polynomial_coeff%nterm
345 : 895935 : call polynomial_term_free(polynomial_coeff%terms(ii))
346 : : end do
347 : : end if
348 [ + + + + : 898458 : ABI_SFREE(polynomial_coeff%terms)
- + - + -
+ - + - +
- + - + -
+ ]
349 : 29362 : polynomial_coeff%name = ""
350 : 29362 : polynomial_coeff%nterm = 0
351 : 29362 : polynomial_coeff%coefficient = zero
352 : :
353 : : !if(trim(polynomial_coeff%debug_str) == "uninitialized") then
354 : : ! ABI_WARNING("Polynomial coeff: is freed before initialization")
355 : : !else
356 : : !if (trim(polynomial_coeff%debug_str) == "freed") then
357 : : ! print *, "- Polynomial coeff: ", trim(polynomial_coeff%debug_str), " ->freed"
358 : : !ABI_ERROR("Polynomial coeff: is freed twice")
359 : : !endif
360 : : !else
361 : : ! print *, "- Polynomial coeff: ", trim(polynomial_coeff%debug_str), " ->freed"
362 : : ! end if
363 : : ! if debug_str starts with freed
364 : : !if(trim(polynomial_coeff%debug_str(1:5)) == "freed") then
365 : : ! print *, "- Polynomial coeff: ", trim(polynomial_coeff%debug_str), " ->freed"
366 : : !end if
367 : :
368 : : !polynomial_coeff%debug_str = "freed "//polynomial_coeff%debug_str
369 : :
370 : 29362 : end subroutine polynomial_coeff_free
371 : : !!***
372 : :
373 : : !!****f* m_polynomial_coeff/polynomial_coeff_list_free
374 : : !!
375 : : !! NAME
376 : : !! polynomial_coeff_list_free
377 : : !!
378 : : !! FUNCTION
379 : : !! Free polynomial_coeff datatype
380 : : !!
381 : : !! INPUTS
382 : : !! polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype
383 : : !!
384 : : !! OUTPUT
385 : : !! polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype
386 : : !!
387 : : !! SOURCE
388 : :
389 : 155 : subroutine polynomial_coeff_list_free(polynomial_coeff_list)
390 : :
391 : : !Arguments ------------------------------------
392 : : !scalars
393 : : !arrays
394 : : type(polynomial_coeff_type),allocatable, intent(inout) :: polynomial_coeff_list(:)
395 : : !Local variables-------------------------------
396 : : !scalar
397 : : integer :: i,ncoeff
398 : : !arrays
399 : :
400 : : ! *************************************************************************
401 : :
402 : : !Free output
403 [ + + ]: 155 : if(allocated(polynomial_coeff_list))then
404 : 117 : ncoeff = size(polynomial_coeff_list)
405 [ + + ]: 12796 : do i=1,ncoeff
406 : 12796 : call polynomial_coeff_free(polynomial_coeff_list(i))
407 : : enddo
408 : : endif
409 [ + + + + : 12834 : ABI_SFREE(polynomial_coeff_list)
- + - - -
- - - - -
- - - - -
- - - -
+ ]
410 : :
411 : 155 : end subroutine polynomial_coeff_list_free
412 : : !!***
413 : :
414 : :
415 : :
416 : : !!****f* m_polynomial_coeff/polynomial_coeff_setCoefficient
417 : : !!
418 : : !! NAME
419 : : !! polynomial_coeff_setCoefficient
420 : : !!
421 : : !! FUNCTION
422 : : !! set the coefficient for of polynomial_coeff
423 : : !!
424 : : !! INPUTS
425 : : !! coefficient = coefficient of this coefficient
426 : : !!
427 : : !! OUTPUT
428 : : !! polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype
429 : : !!
430 : : !! SOURCE
431 : :
432 : 32 : subroutine polynomial_coeff_setCoefficient(coefficient,polynomial_coeff)
433 : :
434 : : !Arguments ------------------------------------
435 : : !scalars
436 : : real(dp),intent(in) :: coefficient
437 : : !arrays
438 : : type(polynomial_coeff_type), intent(inout) :: polynomial_coeff
439 : : !Local variables-------------------------------
440 : : !scalar
441 : : !arrays
442 : : ! *************************************************************************
443 : :
444 : 32 : polynomial_coeff%coefficient = coefficient
445 : :
446 : 32 : end subroutine polynomial_coeff_setCoefficient
447 : : !!***
448 : :
449 : : !!****f* m_polynomial_coeff/polynomial_coeff_setName
450 : : !!
451 : : !! NAME
452 : : !! polynomial_coeff_setName
453 : : !!
454 : : !! FUNCTION
455 : : !! set the name of a polynomial_coeff type
456 : : !!
457 : : !! INPUTS
458 : : !! name = name of the coeff
459 : : !!
460 : : !! OUTPUT
461 : : !! polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype
462 : : !!
463 : : !! SOURCE
464 : :
465 : 1863 : subroutine polynomial_coeff_setName(name,polynomial_coeff)
466 : :
467 : : !Arguments ------------------------------------
468 : : !scalars
469 : : !arrays
470 : : character(len=200),intent(in) :: name
471 : : type(polynomial_coeff_type), intent(inout) :: polynomial_coeff
472 : : !Local variables-------------------------------
473 : : !scalar
474 : : !arrays
475 : : ! *************************************************************************
476 : :
477 : 1863 : polynomial_coeff%name = name
478 : :
479 : 1863 : end subroutine polynomial_coeff_setName
480 : : !!***
481 : :
482 : :
483 : : !!****f* m_polynomial_coeff/polynomial_coeff_getName
484 : : !!
485 : : !! NAME
486 : : !! polynomial_coeff_getName
487 : : !!
488 : : !! FUNCTION
489 : : !! get the name of a polynomial coefficient
490 : : !!
491 : : !! INPUTS
492 : : !! natom = number of atoms
493 : : !! polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype
494 : : !! symbols(natom) = array with the atomic symbol:["Sr","Ru","O1","O2","O3"]
495 : : !! recompute = (optional) flag to set if the name has to be recomputed
496 : : !! iterm = (optional) number of the term used for the name
497 : : !!
498 : : !! OUTPUT
499 : : !! name = name xof the coefficients
500 : : !!
501 : : !! SOURCE
502 : :
503 [ - + ]: 6841 : subroutine polynomial_coeff_getName(name,polynomial_coeff,symbols,recompute,iterm)
504 : :
505 : : !Arguments ------------------------------------
506 : : !scalars
507 : : integer,optional,intent(in) :: iterm
508 : : !arrays
509 : : character(len=5),intent(in) :: symbols(:)
510 : : character(len=200),intent(out):: name
511 : : type(polynomial_coeff_type),optional, intent(in) :: polynomial_coeff
512 : : logical,optional,intent(in) :: recompute
513 : : !Local variables-------------------------------
514 : : !scalar
515 : : integer :: ii,idisp,iterm_in
516 : : logical :: need_recompute
517 : : !arrays
518 : : integer :: cell_atm1(3),cell_atm2(3)
519 : : character(len=1) :: mutodir(9) = (/"x","y","z","1","2","3","4","5","6"/)
520 : : character(len=1) :: dir
521 : : character(len=2) :: power_disp,power_dispchar
522 : : character(len=20) :: atm1,atm2
523 : : character(len=100):: atm1_tmp,atm2_tmp
524 : : character(len=200):: text
525 : : character(len=500):: message
526 : : ! *************************************************************************
527 : :
528 : : !Reset output
529 : 6841 : name=""
530 : 6841 : iterm_in = 1
531 : :
532 : : !Set the optional arguments
533 : 6841 : need_recompute = .FALSE.
534 [ + - ]: 6841 : if(present(recompute)) need_recompute = recompute
535 [ - + ]: 6841 : if(present(iterm)) then
536 : 0 : iterm_in = iterm
537 : : else
538 [ + - ]: 6841 : if(need_recompute)then
539 : 6841 : iterm_in = -1
540 [ + + ]: 157125 : do ii=1,polynomial_coeff%nterm
541 : : ! Find the index of the ref
542 [ + + ]: 157125 : if(iterm_in==-1) then !Need to find the reference term
543 [ + + ]: 208611 : do idisp=1,polynomial_coeff%terms(ii)%ndisp
544 [ + - ]: 208611 : if(polynomial_coeff%terms(ii)%direction(idisp) > 0) then
545 : 205088 : iterm_in = ii
546 [ + + + + : 1112463 : if(any(polynomial_coeff%terms(ii)%cell(:,1,idisp) /= 0).or.&
+ + + + +
+ ]
547 : : & any(polynomial_coeff%terms(ii)%cell(:,2,idisp) /= 0)) then
548 : : iterm_in = -1
549 : : exit
550 : : end if
551 : : end if
552 : : end do!end do disp
553 : : else
554 : : exit
555 : : end if
556 : : end do!end do term
557 : : ! If not find, we set to the first element
558 [ + + ]: 6841 : if(iterm_in==-1) iterm_in = 1
559 : : else
560 : : iterm_in = 1
561 : : end if
562 : : end if
563 : : !Do check
564 [ + - - + ]: 6841 : if(iterm_in > polynomial_coeff%nterm.or.iterm_in < 0) then
565 : : write(message, '(5a)')&
566 : 0 : & ' The number of the requested term for the generation of',ch10,&
567 : 0 : & 'the name of the coefficient is not possible.',ch10,&
568 : 0 : & 'Action: Contact Abinit group.'
569 : 0 : ABI_BUG(message)
570 : : end if
571 : :
572 [ - + - - ]: 6841 : if(polynomial_coeff%name /= "".and..not.need_recompute)then
573 : 0 : name = polynomial_coeff%name
574 : : else
575 : : ! Nedd to recompute
576 [ + + ]: 22904 : do idisp=1,polynomial_coeff%terms(iterm_in)%ndisp
577 : 16063 : text = ""
578 : : !Fill variables for this displacement
579 : 16063 : write(power_dispchar,'(I0)') polynomial_coeff%terms(iterm_in)%power_disp(idisp)
580 : 16063 : power_disp=trim(power_dispchar)
581 : :
582 : 16063 : atm1=symbols(polynomial_coeff%terms(iterm_in)%atindx(1,idisp))
583 : 16063 : atm2=symbols(polynomial_coeff%terms(iterm_in)%atindx(2,idisp))
584 : 16063 : dir=mutodir(polynomial_coeff%terms(iterm_in)%direction(idisp))
585 [ + + ]: 64252 : cell_atm1=polynomial_coeff%terms(iterm_in)%cell(:,1,idisp)
586 [ + + ]: 64252 : cell_atm2=polynomial_coeff%terms(iterm_in)%cell(:,2,idisp)
587 : : ! Construct ATM1
588 [ + + + - : 64252 : if (any(cell_atm1(:) /= 0) )then
- + ]
589 : 0 : write(atm1_tmp,'(4a,I0,a,I0,a,I0,a)') trim(atm1),"_",dir,"[",cell_atm1(1)," ",&
590 : 0 : & cell_atm1(2)," ",cell_atm1(3),"]"
591 : : else
592 [ + - ]: 16063 : atm1_tmp = trim(atm1)//"_"//dir
593 : : end if
594 : : ! Construct ATM2
595 [ + + + + : 52604 : if(any(cell_atm2(:) /= 0))then
+ + ]
596 : 4876 : write(atm2_tmp,'(4a,I0,a,I0,a,I0,a)') trim(atm2),"_",dir,"[",cell_atm2(1)," ",&
597 : 9752 : & cell_atm2(2)," ",cell_atm2(3),"]"
598 : : else
599 [ + - ]: 11187 : atm2_tmp = trim(atm2)//"_"//dir
600 : : end if
601 [ + - ]: 16063 : text="("//trim(atm1_tmp)//"-"//trim(atm2_tmp)//")^"//power_disp
602 : 22904 : name = trim(name)//trim(text)
603 : : end do
604 : : !Strain case
605 [ + + ]: 9594 : do idisp=1,polynomial_coeff%terms(iterm_in)%nstrain
606 : 2753 : write(power_dispchar,'(I0)') polynomial_coeff%terms(iterm_in)%power_strain(idisp)
607 : 2753 : power_disp=trim(power_dispchar)
608 : 2753 : dir=mutodir(3+polynomial_coeff%terms(iterm_in)%strain(idisp))
609 [ + - ]: 2753 : text="("//"eta_"//trim(dir)//")^"//power_disp
610 : 9594 : name = trim(name)//trim(text)
611 : : end do
612 : : end if
613 : :
614 : 6841 : end subroutine polynomial_coeff_getName
615 : : !!***
616 : :
617 : : !!****f* m_polynomial_coeff/polynomial_coeff_broadcast
618 : : !! NAME
619 : : !! polynomial_coeff_broadcast
620 : : !!
621 : : !! FUNCTION
622 : : !! MPI broadcast polynomial_coefficent datatype
623 : : !!
624 : : !! INPUTS
625 : : !! source = rank of source
626 : : !! comm = MPI communicator
627 : : !!
628 : : !! SIDE EFFECTS
629 : : !! coefficients<type(polynomial_coefficent_type)>= Input if node is source,
630 : : !! other nodes returns with a completely initialized instance.
631 : : !!
632 : : !! SOURCE
633 : :
634 : 2500 : subroutine polynomial_coeff_broadcast(coefficients, source, comm)
635 : :
636 : : !Arguments ------------------------------------
637 : : !array
638 : : type(polynomial_coeff_type),intent(inout) :: coefficients
639 : : integer, intent(in) :: source,comm
640 : :
641 : : !Local variables-------------------------------
642 : : !scalars
643 : : integer :: ierr,ii
644 : : !arrays
645 : :
646 : : ! *************************************************************************
647 : :
648 : :
649 [ + + ]: 2068 : if (xmpi_comm_size(comm) == 1) return
650 : :
651 : : ! Free the output
652 [ + + ]: 144 : if (xmpi_comm_rank(comm) /= source) then
653 : 96 : call polynomial_coeff_free(coefficients)
654 : : end if
655 : :
656 : : ! Transmit variables
657 : 144 : call xmpi_bcast(coefficients%name, source, comm, ierr)
658 : 144 : call xmpi_bcast(coefficients%nterm, source, comm, ierr)
659 : 144 : call xmpi_bcast(coefficients%coefficient, source, comm, ierr)
660 : :
661 : : !Allocate arrays on the other nodes.
662 [ + + ]: 144 : if (xmpi_comm_rank(comm) /= source) then
663 [ - + + - : 3504 : ABI_MALLOC(coefficients%terms,(coefficients%nterm))
+ - + - -
+ - + - +
+ + + - ]
664 [ + + ]: 3312 : do ii=1,coefficients%nterm
665 : 3312 : call polynomial_term_free(coefficients%terms(ii))
666 : : end do
667 : : end if
668 : : ! Set the number of term on each node (needed for allocations of array)
669 [ + + ]: 4968 : do ii = 1,coefficients%nterm
670 : 4824 : call xmpi_bcast(coefficients%terms(ii)%ndisp, source, comm, ierr)
671 : 4824 : call xmpi_bcast(coefficients%terms(ii)%nstrain, source, comm, ierr)
672 : 4968 : call xmpi_bcast(coefficients%terms(ii)%nindex, source, comm, ierr)
673 : : end do
674 : :
675 : : ! Allocate arrays on the other nodes
676 [ + + ]: 144 : if (xmpi_comm_rank(comm) /= source) then
677 [ + + ]: 3312 : do ii = 1,coefficients%nterm
678 [ - + + - : 9648 : ABI_MALLOC(coefficients%terms(ii)%atindx,(2,coefficients%terms(ii)%ndisp))
+ - + - -
+ - + -
+ ]
679 [ + + + + ]: 27984 : coefficients%terms(ii)%atindx = 0
680 [ - + + - : 9648 : ABI_MALLOC(coefficients%terms(ii)%direction,(coefficients%terms(ii)%ndisp))
+ - + - -
+ - + -
+ ]
681 [ - + + - : 9648 : ABI_MALLOC(coefficients%terms(ii)%cell,(3,2,coefficients%terms(ii)%ndisp))
+ - + - -
+ - + -
+ ]
682 [ - + + - : 9648 : ABI_MALLOC(coefficients%terms(ii)%power_disp,(coefficients%terms(ii)%ndisp))
+ - + - -
+ - + -
+ ]
683 [ + + + - : 9648 : ABI_MALLOC(coefficients%terms(ii)%power_strain,(coefficients%terms(ii)%nstrain))
+ - + + -
+ - + -
+ ]
684 [ + + + - : 9648 : ABI_MALLOC(coefficients%terms(ii)%strain,(coefficients%terms(ii)%nstrain))
+ - + + -
+ - + -
+ ]
685 [ + + ]: 3312 : if(coefficients%terms(ii)%nindex>-1) then
686 [ - + + - : 5472 : ABI_MALLOC(coefficients%terms(ii)%index_coeff,(coefficients%terms(ii)%nindex))
+ - + - -
+ - + -
+ ]
687 : : end if
688 : : end do
689 : : end if
690 : :
691 : : ! Transfert value
692 [ + + ]: 4968 : do ii = 1,coefficients%nterm
693 : 4824 : call xmpi_bcast(coefficients%terms(ii)%weight, source, comm, ierr)
694 : 4824 : call xmpi_bcast(coefficients%terms(ii)%atindx, source, comm, ierr)
695 : 4824 : call xmpi_bcast(coefficients%terms(ii)%direction, source, comm, ierr)
696 : 4824 : call xmpi_bcast(coefficients%terms(ii)%cell, source, comm, ierr)
697 : 4824 : call xmpi_bcast(coefficients%terms(ii)%power_disp, source, comm, ierr)
698 : 4824 : call xmpi_bcast(coefficients%terms(ii)%power_strain, source, comm, ierr)
699 : 4824 : call xmpi_bcast(coefficients%terms(ii)%strain, source, comm, ierr)
700 [ + + ]: 4968 : if(coefficients%terms(ii)%nindex>-1) then
701 : 2736 : call xmpi_bcast(coefficients%terms(ii)%index_coeff, source, comm, ierr)
702 : : end if
703 : : end do
704 : : end subroutine polynomial_coeff_broadcast
705 : : !!***
706 : :
707 : : !!****f* m_polynomial_coeff/polynomial_coeff_MPIsend
708 : : !! NAME
709 : : !! polynomial_coeff_MPIsend
710 : : !!
711 : : !! FUNCTION
712 : : !! MPI send the polynomial_coefficent datatype
713 : : !!
714 : : !! INPUTS
715 : : !! tag = tag of the message to send
716 : : !! dest= rank of Dest
717 : : !! comm= MPI communicator
718 : : !!
719 : : !! SIDE EFFECTS
720 : : !! polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype
721 : : !!
722 : : !! SOURCE
723 : :
724 : 88 : subroutine polynomial_coeff_MPIsend(coefficients, tag, dest, comm)
725 : :
726 : : !Arguments ------------------------------------
727 : : !array
728 : : type(polynomial_coeff_type),intent(inout) :: coefficients
729 : : integer, intent(in) :: dest,comm,tag
730 : :
731 : : !Local variables-------------------------------
732 : : !scalars
733 : : integer :: ierr,ii
734 : : integer :: my_rank
735 : : !arrays
736 : :
737 : : ! *************************************************************************
738 : :
739 [ - + ]: 88 : if (xmpi_comm_size(comm) == 1) return
740 : :
741 : 88 : my_rank = xmpi_comm_rank(comm)
742 : : ! Transmit variables
743 : 88 : call xmpi_send(coefficients%name, dest, 9*tag+0, comm, ierr)
744 : 88 : call xmpi_send(coefficients%nterm, dest, 9*tag+1, comm, ierr)
745 : 88 : call xmpi_send(coefficients%coefficient, dest, 9*tag+2, comm, ierr)
746 : :
747 : : ! Set the number of term on each node (needed for allocations of array)
748 [ + + ]: 3760 : do ii = 1,coefficients%nterm
749 : 3672 : call xmpi_send(coefficients%terms(ii)%ndisp, dest, 9*tag+3, comm, ierr)
750 : 3672 : call xmpi_send(coefficients%terms(ii)%nstrain, dest, 9*tag+4, comm, ierr)
751 : 3760 : call xmpi_send(coefficients%terms(ii)%nindex, dest, 9*tag+13, comm, ierr)
752 : : end do
753 : :
754 : : ! Transfert value
755 [ + + ]: 3760 : do ii = 1,coefficients%nterm
756 : 3672 : call xmpi_send(coefficients%terms(ii)%weight, dest, 9*tag+5, comm, ierr)
757 : 3672 : call xmpi_send(coefficients%terms(ii)%atindx, dest, 9*tag+6, comm, ierr)
758 : 3672 : call xmpi_send(coefficients%terms(ii)%direction, dest, 9*tag+7, comm, ierr)
759 : 3672 : call xmpi_send(coefficients%terms(ii)%cell, dest, 9*tag+8, comm, ierr)
760 : 3672 : call xmpi_send(coefficients%terms(ii)%power_disp, dest, 9*tag+9, comm, ierr)
761 : 3672 : call xmpi_send(coefficients%terms(ii)%power_strain, dest, 9*tag+10, comm, ierr)
762 : 3672 : call xmpi_send(coefficients%terms(ii)%strain, dest, 9*tag+11, comm, ierr)
763 [ + - ]: 3760 : if (coefficients%terms(ii)%nindex>-1) then
764 : 3672 : call xmpi_send(coefficients%terms(ii)%index_coeff, dest, 9*tag+14, comm, ierr)
765 : : end if
766 : : end do
767 : : end subroutine polynomial_coeff_MPIsend
768 : : !!***
769 : :
770 : : !!****f* m_polynomial_coeff/polynomial_coeff_MPIrecv
771 : : !! NAME
772 : : !! polynomial_coeff_MPIrecv
773 : : !!
774 : : !! FUNCTION
775 : : !! MPI receive the polynomial_coefficent datatype
776 : : !!
777 : : !! INPUTS
778 : : !! tag = tag of the message to receive
779 : : !! source = rank of Source
780 : : !! comm = MPI communicator
781 : : !!
782 : : !! SIDE EFFECTS
783 : : !! coefficients<type(polynomial_coefficent_type)>= polynomial_coeff datatype
784 : : !!
785 : : !! SOURCE
786 : :
787 : 88 : subroutine polynomial_coeff_MPIrecv(coefficients, tag, source, comm)
788 : :
789 : : !Arguments ------------------------------------
790 : : !array
791 : : type(polynomial_coeff_type),intent(inout) :: coefficients
792 : : integer, intent(in) :: source,comm,tag
793 : :
794 : : !Local variables-------------------------------
795 : : !scalars
796 : : integer :: ierr,ii
797 : : !arrays
798 : :
799 : : ! *************************************************************************
800 : :
801 [ - + ]: 88 : if (xmpi_comm_size(comm) == 1) return
802 : :
803 : :
804 : : ! Free the output
805 : 88 : call polynomial_coeff_free(coefficients)
806 : :
807 : : ! Transmit variables
808 : 88 : call xmpi_recv(coefficients%name, source, 9*tag+0, comm, ierr)
809 : 88 : call xmpi_recv(coefficients%nterm, source, 9*tag+1, comm, ierr)
810 : 88 : call xmpi_recv(coefficients%coefficient, source, 9*tag+2, comm, ierr)
811 : :
812 : : !Allocate arrays on the other nodes.
813 [ - + + - : 3936 : ABI_MALLOC(coefficients%terms,(coefficients%nterm))
+ - + - -
+ - + - +
+ + + - ]
814 [ + + ]: 3760 : do ii=1,coefficients%nterm
815 : 3760 : call polynomial_term_free(coefficients%terms(ii))
816 : : end do
817 : :
818 : : ! Set the number of term on each node (needed for allocations of array)
819 [ + + ]: 3760 : do ii = 1,coefficients%nterm
820 : 3672 : call xmpi_recv(coefficients%terms(ii)%ndisp, source, 9*tag+3, comm, ierr)
821 : 3672 : call xmpi_recv(coefficients%terms(ii)%nstrain, source, 9*tag+4, comm, ierr)
822 : 3760 : call xmpi_recv(coefficients%terms(ii)%nindex, source, 9*tag+13, comm, ierr)
823 : : end do
824 : :
825 : : ! Allocate arrays on the other nodes
826 [ + + ]: 3760 : do ii = 1,coefficients%nterm
827 [ - + + - : 11016 : ABI_MALLOC(coefficients%terms(ii)%atindx,(2,coefficients%terms(ii)%ndisp))
+ - + - -
+ - + -
+ ]
828 [ + + + + ]: 32256 : coefficients%terms(ii)%atindx = 0
829 [ - + + - : 11016 : ABI_MALLOC(coefficients%terms(ii)%direction,(coefficients%terms(ii)%ndisp))
+ - + - -
+ - + -
+ ]
830 [ - + + - : 11016 : ABI_MALLOC(coefficients%terms(ii)%cell,(3,2,coefficients%terms(ii)%ndisp))
+ - + - -
+ - + -
+ ]
831 [ - + + - : 11016 : ABI_MALLOC(coefficients%terms(ii)%power_disp,(coefficients%terms(ii)%ndisp))
+ - + - -
+ - + -
+ ]
832 [ + + + - : 11016 : ABI_MALLOC(coefficients%terms(ii)%power_strain,(coefficients%terms(ii)%nstrain))
+ - + + -
+ - + -
+ ]
833 [ + + + - : 11016 : ABI_MALLOC(coefficients%terms(ii)%strain,(coefficients%terms(ii)%nstrain))
+ - + + -
+ - + -
+ ]
834 [ + - ]: 3760 : if (coefficients%terms(ii)%nindex>-1) then
835 [ - + + - : 11016 : ABI_MALLOC(coefficients%terms(ii)%index_coeff,(coefficients%terms(ii)%nindex))
+ - + - -
+ - + -
+ ]
836 : : end if
837 : : end do
838 : :
839 : : ! Transfert value
840 [ + + ]: 3760 : do ii = 1,coefficients%nterm
841 : 3672 : call xmpi_recv(coefficients%terms(ii)%weight, source, 9*tag+5, comm, ierr)
842 : 3672 : call xmpi_recv(coefficients%terms(ii)%atindx, source, 9*tag+6, comm, ierr)
843 : 3672 : call xmpi_recv(coefficients%terms(ii)%direction, source, 9*tag+7, comm, ierr)
844 : 3672 : call xmpi_recv(coefficients%terms(ii)%cell, source, 9*tag+8, comm, ierr)
845 : 3672 : call xmpi_recv(coefficients%terms(ii)%power_disp, source, 9*tag+9, comm, ierr)
846 : 3672 : call xmpi_recv(coefficients%terms(ii)%power_strain, source, 9*tag+10, comm, ierr)
847 : 3672 : call xmpi_recv(coefficients%terms(ii)%strain, source, 9*tag+11, comm, ierr)
848 [ + - ]: 3760 : if (coefficients%terms(ii)%nindex>-1) then
849 : 3672 : call xmpi_recv(coefficients%terms(ii)%index_coeff, source , 9*tag+14, comm, ierr)
850 : : end if
851 : : end do
852 : :
853 : : end subroutine polynomial_coeff_MPIrecv
854 : : !!***
855 : :
856 : :
857 : :
858 : :
859 : : !!****f*m_polynomial_coeff/polynomial_coeff_writeXML
860 : : !! NAME
861 : : !! polynomial_coeff_writeXML
862 : : !!
863 : : !! FUNCTION
864 : : !! This routine print the coefficents into XML format
865 : : !!
866 : : !! COPYRIGHT
867 : : !! Copyright (C) 2000-2026 ABINIT group (AM)
868 : : !! This file is distributed under the terms of the
869 : : !! GNU General Public License, see ~abinit/COPYING
870 : : !! or http://www.gnu.org/copyleft/gpl.txt .
871 : : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
872 : : !!
873 : : !! INPUTS
874 : : !! coeffs(ncoeffs)<type(polynomial_coeff)> = array of polynomial_coeff datatype
875 : : !! ncoeff = number of coeffs to print
876 : : !! filename = optional,the name of output file
877 : : !! default is coefficients.xml
878 : : !! unit = optional,unit of the output file
879 : : !! newfile = optional, TRUE the coefficients are print in new XML (print the headers)
880 : : !! FALSE (requieres unit) will not print the headers
881 : : !! replace = optional, TRUE replace filename if filename exists
882 : : !! FALSE, default not replace if filename exists
883 : : !!
884 : : !! OUTPUT
885 : : !!
886 : : !! SOURCE
887 : :
888 : 15 : subroutine polynomial_coeff_writeXML(coeffs,ncoeff,filename,unit,newfile,replace)
889 : :
890 : : !Arguments ------------------------------------
891 : : !scalars
892 : : integer, intent(in) :: ncoeff
893 : : integer,optional,intent(in) :: unit
894 : : logical,optional,intent(in) :: newfile,replace
895 : : !arrays
896 : : type(polynomial_coeff_type), intent(in) :: coeffs(ncoeff)
897 : : character(len=fnlen),optional,intent(in) :: filename
898 : : !Local variables-------------------------------
899 : : !scalar
900 : : integer :: icoeff,idisp,iterm
901 : : integer :: unit_xml
902 : : logical :: need_header,need_to_replace
903 : : character(len=500) :: message
904 : : character(len=fnlen) :: namefile
905 : : character(len=1) :: direction
906 : :
907 : : !arrays
908 : :
909 : : ! *************************************************************************
910 : :
911 : : !fill the default
912 : 15 : unit_xml = get_unit()
913 : :
914 : : !Check the inputs
915 [ + - ]: 15 : if(present(filename))then
916 : 15 : namefile=trim(filename)
917 : : else
918 : 0 : namefile='coefficients.xml'
919 : : end if
920 : :
921 : 15 : need_to_replace = .FALSE.
922 [ - + ]: 15 : if(present(replace))then
923 : 0 : need_to_replace = replace
924 : : end if
925 : :
926 : 15 : need_header = .TRUE.
927 [ - + ]: 15 : if(present(newfile))then
928 [ + + ]: 15 : if (newfile) then
929 : 12 : unit_xml = get_unit()
930 : 12 : need_header = .TRUE.
931 : :
932 [ - + ]: 12 : if(.not. need_to_replace) call isfile(namefile,'new')
933 : : else
934 [ - + ]: 3 : if(.not.present(unit))then
935 : 0 : write(message,'(a,a)')' You need to specified the unit'
936 : 0 : ABI_ERROR(message)
937 : : else
938 : 3 : need_header = .FALSE.
939 : 3 : unit_xml = unit
940 : : end if
941 : : end if
942 : : end if
943 [ - + ]: 15 : if (size(coeffs) /= ncoeff) then
944 : 0 : write(message,'(a,a)')' The number of coeffs does not correspond to ncoeff'
945 : 0 : ABI_ERROR(message)
946 : : end if
947 : :
948 : : !Print the coefficients into XML file
949 [ + - ]: 15 : if(ncoeff>0)then
950 [ + + ]: 15 : if(need_header)then
951 : : ! open new file
952 [ - + ]: 12 : if (open_file(namefile,message,unit=unit_xml,form="formatted",&
953 : : & status="new",action="write") /= 0) then
954 : 0 : ABI_ERROR(message)
955 : : end if
956 : : else
957 : : ! just open the file to append the coefficient
958 : 3 : open(unit=unit_xml,file=namefile,position="append")
959 : : end if
960 : :
961 : : ! Write header
962 [ + + ]: 15 : if (need_header)then
963 : 12 : write(message,'(a,a,a)')ch10,&
964 : 24 : & ' Generation of the xml file for the fitted polynomial in ',trim(namefile)
965 : :
966 : 12 : call wrtout(ab_out,message,'COLL')
967 : 12 : call wrtout(std_out,message,'COLL')
968 : 12 : WRITE(unit_xml,'("<?xml version=""1.0"" ?>")')
969 : : end if
970 : 15 : WRITE(unit_xml,'("<Heff_definition>")')
971 : : ! Close header
972 [ + + ]: 95 : do icoeff = 1, ncoeff
973 : : WRITE(unit_xml,'(" <coefficient number=""",I0,""" value=""",E19.10,""" text=""",a,""">")') &
974 [ + - ]: 80 : icoeff,coeffs(icoeff)%coefficient,trim(coeffs(icoeff)%name)
975 [ + + ]: 2438 : do iterm = 1,coeffs(icoeff)%nterm
976 : : WRITE(unit_xml,'(" <term weight=""",F9.6,""">")') &
977 : 2358 : coeffs(icoeff)%terms(iterm)%weight
978 [ + + ]: 8022 : do idisp=1,coeffs(icoeff)%terms(iterm)%ndisp
979 : : ! Atomic displacement case
980 : 7552 : select case(coeffs(icoeff)%terms(iterm)%direction(idisp))
981 : : case(1)
982 : 1888 : direction ="x"
983 : : case(2)
984 : 1888 : direction ="y"
985 : : case(3)
986 [ + + + - ]: 5664 : direction ="z"
987 : : end select
988 : 5664 : WRITE(unit_xml,'(a,I0,a,I0,3a,I0,a)') " <displacement_diff atom_a=""",&
989 : 5664 : & coeffs(icoeff)%terms(iterm)%atindx(1,idisp)-1,""" atom_b=""",&
990 : 5664 : & coeffs(icoeff)%terms(iterm)%atindx(2,idisp)-1,""" direction=""",&
991 : 5664 : & direction,""" power=""",coeffs(icoeff)%terms(iterm)%power_disp(idisp),&
992 : 11328 : & """>"
993 : 5664 : WRITE(unit_xml,'(" <cell_a>")',advance='no')
994 : : WRITE(unit_xml,'(3(I0,a,I0,a,I0))',advance='no')&
995 : 5664 : & coeffs(icoeff)%terms(iterm)%cell(1,1,idisp)," ",&
996 : 5664 : & coeffs(icoeff)%terms(iterm)%cell(2,1,idisp)," ",&
997 : 11328 : & coeffs(icoeff)%terms(iterm)%cell(3,1,idisp)
998 : 5664 : WRITE(unit_xml,'("</cell_a>")')
999 : 5664 : WRITE(unit_xml,'(" <cell_b>")',advance='no')
1000 : : WRITE(unit_xml,'(3(I0,a,I0,a,I0))',advance='no')&
1001 : 5664 : & coeffs(icoeff)%terms(iterm)%cell(1,2,idisp)," ",&
1002 : 5664 : & coeffs(icoeff)%terms(iterm)%cell(2,2,idisp)," ",&
1003 : 11328 : & coeffs(icoeff)%terms(iterm)%cell(3,2,idisp)
1004 : 5664 : WRITE(unit_xml,'("</cell_b>")')
1005 : 8022 : WRITE(unit_xml,'(" </displacement_diff>")')
1006 : : end do
1007 [ + + ]: 2892 : do idisp=1,coeffs(icoeff)%terms(iterm)%nstrain
1008 : : ! Strain case
1009 : : WRITE(unit_xml,'(" <strain power=""",i2,""" voigt=""",i2,"""/>")')&
1010 : 534 : & coeffs(icoeff)%terms(iterm)%power_strain(idisp),&
1011 : 3426 : & coeffs(icoeff)%terms(iterm)%strain(idisp)
1012 : : end do
1013 : 2438 : WRITE(unit_xml,'(" </term>")')
1014 : : end do
1015 : 95 : WRITE(unit_xml,'(" </coefficient>")')
1016 : : end do
1017 : 15 : WRITE(unit_xml,'("</Heff_definition>")')
1018 : : ! Close file
1019 : 15 : CLOSE(unit_xml)
1020 : : end if
1021 : :
1022 : 15 : end subroutine polynomial_coeff_writeXML
1023 : :
1024 : :
1025 : : !!***
1026 : :
1027 : : !!****f* m_polynomial_coeff/polynomial_coeff_evaluate
1028 : : !! NAME
1029 : : !! polynomial_coeff_evaluate
1030 : : !!
1031 : : !! FUNCTION
1032 : : !! Compute the energy related to the coefficients from
1033 : : !! fitted polynome
1034 : : !!
1035 : : !! INPUTS
1036 : : !! coefficients(ncoeff)<type(polynomial_coeff_type)> = list of coefficients
1037 : : !! disp(3,natom_sc) = atomics displacement between configuration and the reference
1038 : : !! natom_sc = number of atoms in the supercell
1039 : : !! natom_uc = number of atoms in the unit cell
1040 : : !! ncoeff = number of coefficients
1041 : : !! sc_size(3) = size of the supercell (2 2 2 for example)
1042 : : !! strain(6) = strain between configuration and the reference
1043 : : !! cells(ncell) = number of the cells into the supercell (1,2,3,4,5)
1044 : : !! ncell = total number of cell to treat by this cpu
1045 : : !! index_cells(3,ncell) = indexes of the cells into supercell (-1 -1 -1 ,...,1 1 1)
1046 : : !! comm=MPI communicator
1047 : : !!
1048 : : !! OUTPUT
1049 : : !! energy = contribution to the energy
1050 : : !! energy_coeff(ncoeff) = energy contribution of each anharmonic term
1051 : : !! fcart(3,natom) = contribution to the forces
1052 : : !! strten(6) = contribution to the stress tensor
1053 : : !!
1054 : : !! SOURCE
1055 : : !!
1056 : 10265 : subroutine polynomial_coeff_evaluate(coefficients,disp,energy,energy_coeff,fcart,natom_sc,natom_uc,ncoeff,sc_size,&
1057 : 10265 : & strain,strten,ncell,index_cells,comm,filename)
1058 : :
1059 : : !Arguments ------------------------------------
1060 : : ! scalar
1061 : : real(dp),intent(out):: energy
1062 : : integer, intent(in) :: ncell,ncoeff,natom_sc,natom_uc
1063 : : integer, intent(in) :: comm
1064 : : character(len=fnlen),optional,intent(in) :: filename
1065 : : ! array
1066 : : real(dp),intent(out):: strten(6)
1067 : : real(dp),intent(in) :: strain(6)
1068 : : real(dp),intent(out):: fcart(3,natom_sc)
1069 : : real(dp),intent(in) :: disp(3,natom_sc)
1070 : : real(dp),optional,intent(out):: energy_coeff(ncoeff)
1071 : : integer,intent(in) :: index_cells(4,ncell)
1072 : : integer,intent(in) :: sc_size(3)
1073 : : type(polynomial_coeff_type),intent(in) :: coefficients(ncoeff)
1074 : : !Local variables-------------------------------
1075 : : ! scalar
1076 : : integer :: i1,i2,i3,ia1,ib1,ia2,ib2,idir1,idir2,ierr,ii
1077 : : integer :: icoeff,iterm,idisp1,idisp2,idisp1_strain,idisp2_strain,icell,ndisp
1078 : : integer :: nstrain,ndisp_tot,power_disp,power_strain,unit_out
1079 : : real(dp):: coeff,disp1,disp2,tmp1,tmp2,tmp3,weight
1080 : : logical :: file_opened
1081 : : ! array
1082 : : integer :: cell_atoma1(3),cell_atoma2(3)
1083 : : integer :: cell_atomb1(3),cell_atomb2(3)
1084 : : character(len=500) :: msg
1085 : : character(len=fnlen) :: name_file
1086 : : ! *************************************************************************
1087 : :
1088 : : ! Check
1089 [ + + + - : 41060 : if (any(sc_size <= 0)) then
- + ]
1090 : 0 : write(msg,'(a,a)')' No supercell found for getEnergy'
1091 : 0 : ABI_ERROR(msg)
1092 : : end if
1093 : :
1094 [ + + ]: 10265 : if(present(filename)) name_file = filename
1095 : :
1096 : : ! Initialisation of variables
1097 : 10265 : energy = zero
1098 [ + + + + ]: 35095225 : fcart(:,:) = zero
1099 : 10265 : strten(:) = zero
1100 [ + + ]: 876307 : energy_coeff(:) = zero
1101 [ + + ]: 1758409 : do icell = 1,ncell
1102 : 1748144 : ii = index_cells(4,icell);
1103 : 1748144 : i1=index_cells(1,icell); i2=index_cells(2,icell); i3=index_cells(3,icell)
1104 : 1748144 : ia1 = 0 ; ib1 = 0
1105 : : ! Loop over coefficients
1106 [ + + ]: 173401737 : do icoeff=1,ncoeff
1107 : : ! Set the value of the coefficient
1108 : 171643328 : coeff = coefficients(icoeff)%coefficient
1109 : : ! Loop over terms of this coefficient
1110 [ + + ]: 2165241000 : do iterm=1,coefficients(icoeff)%nterm
1111 : : ! Set the weight of this term
1112 : 1991849528 : weight =coefficients(icoeff)%terms(iterm)%weight
1113 : 1991849528 : tmp1 = one
1114 : 1991849528 : ndisp = coefficients(icoeff)%terms(iterm)%ndisp
1115 : 1991849528 : nstrain = coefficients(icoeff)%terms(iterm)%nstrain
1116 : 1991849528 : ndisp_tot = ndisp + nstrain
1117 : : ! Loop over displacement and strain
1118 [ + + ]: 6148709792 : do idisp1=1,ndisp_tot
1119 : : ! Set to one the acculation of forces and strain
1120 : 4156860264 : tmp2 = one
1121 : 4156860264 : tmp3 = zero
1122 : :
1123 : : ! Strain case idisp > ndisp
1124 [ + + ]: 4156860264 : if (idisp1 > ndisp)then
1125 : 1477156968 : tmp3 = one
1126 : : ! Set the power_strain of the strain:
1127 : 1477156968 : idisp1_strain = idisp1 - ndisp
1128 : 1477156968 : power_strain = coefficients(icoeff)%terms(iterm)%power_strain(idisp1_strain)
1129 : : ! Get the direction of the displacement or strain
1130 : 1477156968 : idir1 = coefficients(icoeff)%terms(iterm)%strain(idisp1_strain)
1131 [ + + ]: 1477156968 : if(abs(strain(idir1)) > tol10)then
1132 : : ! Accumulate energy fo each displacement (\sum ((A_x-O_x)^Y(A_y-O_c)^Z))
1133 : 4049736 : tmp1 = tmp1 * (strain(idir1))**power_strain
1134 [ + + ]: 4049736 : if(power_strain > 1) then
1135 : : ! Accumulate stress for each strain (\sum (Y(eta_2)^Y-1(eta_2)^Z+...))
1136 : 523872 : tmp3 = tmp3 * power_strain*(strain(idir1))**(power_strain-1)
1137 : : end if
1138 : : else
1139 : 1473107232 : tmp1 = zero
1140 [ + + ]: 1473107232 : if(power_strain > 1) then
1141 : 1338078744 : tmp3 = zero
1142 : : end if
1143 : : end if
1144 : : else
1145 : : ! Set the power_disp of the displacement:
1146 : 2679703296 : power_disp = coefficients(icoeff)%terms(iterm)%power_disp(idisp1)
1147 : : ! Get the direction of the displacement or strain
1148 : 2679703296 : idir1 = coefficients(icoeff)%terms(iterm)%direction(idisp1)
1149 : : ! Displacement case idir = 1, 2 or 3
1150 : : ! indexes of the cell of the atom a
1151 [ + + ]: 10718813184 : cell_atoma1 = coefficients(icoeff)%terms(iterm)%cell(:,1,idisp1)
1152 [ + - + + : 2679703296 : if(cell_atoma1(1)/=0.or.cell_atoma1(2)/=0.or.cell_atoma1(3)/=0) then
+ + ]
1153 : : ! if the cell is not 0 0 0 we apply PBC:
1154 : 64800 : cell_atoma1(1) = i1 + cell_atoma1(1)
1155 : 64800 : cell_atoma1(2) = i2 + cell_atoma1(2)
1156 : 64800 : cell_atoma1(3) = i3 + cell_atoma1(3)
1157 : 64800 : call getPBCIndexes_supercell(cell_atoma1(1:3),sc_size(1:3))
1158 : : ! index of the first atom (position in the supercell if the cell is not 0 0 0)
1159 : : ia1 = (cell_atoma1(1)-1)*sc_size(2)*sc_size(3)*natom_uc+&
1160 : : & (cell_atoma1(2)-1)*sc_size(3)*natom_uc+&
1161 : : & (cell_atoma1(3)-1)*natom_uc+&
1162 : 64800 : & coefficients(icoeff)%terms(iterm)%atindx(1,idisp1)
1163 : : else
1164 : : ! index of the first atom (position in the supercell if the cell is 0 0 0)
1165 : 2679638496 : ia1 = ii + coefficients(icoeff)%terms(iterm)%atindx(1,idisp1)
1166 : : end if
1167 : :
1168 : : ! indexes of the cell of the atom b (with PBC) same as ia1
1169 [ + + ]: 10718813184 : cell_atomb1 = coefficients(icoeff)%terms(iterm)%cell(:,2,idisp1)
1170 [ + + + + : 2679703296 : if(cell_atomb1(1)/=0.or.cell_atomb1(2)/=0.or.cell_atomb1(3)/=0) then
+ + ]
1171 : 1728987576 : cell_atomb1(1) = i1 + cell_atomb1(1)
1172 : 1728987576 : cell_atomb1(2) = i2 + cell_atomb1(2)
1173 : 1728987576 : cell_atomb1(3) = i3 + cell_atomb1(3)
1174 : 1728987576 : call getPBCIndexes_supercell(cell_atomb1(1:3),sc_size(1:3))
1175 : :
1176 : : ! index of the second atom in the (position in the supercell if the cell is not 0 0 0)
1177 : : ib1 = (cell_atomb1(1)-1)*sc_size(2)*sc_size(3)*natom_uc+&
1178 : : & (cell_atomb1(2)-1)*sc_size(3)*natom_uc+&
1179 : : & (cell_atomb1(3)-1)*natom_uc+&
1180 : 1728987576 : & coefficients(icoeff)%terms(iterm)%atindx(2,idisp1)
1181 : : else
1182 : : ! index of the first atom (position in the supercell if the cell is 0 0 0)
1183 : 950715720 : ib1 = ii + coefficients(icoeff)%terms(iterm)%atindx(2,idisp1)
1184 : : end if
1185 : :
1186 : : ! Get the displacement for the both atoms
1187 : 2679703296 : disp1 = disp(idir1,ia1)
1188 : 2679703296 : disp2 = disp(idir1,ib1)
1189 : :
1190 [ + + + + ]: 2679703296 : if(abs(disp1) > tol10 .or. abs(disp2)> tol10)then
1191 : : ! Accumulate energy fo each displacement (\sum ((A_x-O_x)^Y(A_y-O_c)^Z))
1192 : 2675133608 : tmp1 = tmp1 * (disp1-disp2)**power_disp
1193 [ + + ]: 2675133608 : if(power_disp > 1) then
1194 : : ! Accumulate forces for each displacement (\sum (Y(A_x-O_x)^Y-1(A_y-O_c)^Z+...))
1195 : 2398068632 : tmp2 = tmp2 * power_disp*(disp1-disp2)**(power_disp-1)
1196 : : end if
1197 : : else
1198 : 4569688 : tmp1 = zero
1199 [ + + ]: 4569688 : if(power_disp > 1) then
1200 : 3520840 : tmp2 = zero
1201 : : end if
1202 : : end if
1203 : : end if
1204 : :
1205 [ + + ]: 13733453856 : do idisp2=1,ndisp_tot
1206 : :
1207 [ + + ]: 13733453856 : if(idisp2 /= idisp1) then
1208 [ + + ]: 5419733328 : if (idisp2 > ndisp)then
1209 : 2187237792 : idisp2_strain = idisp2 - ndisp
1210 : 2187237792 : idir2 = coefficients(icoeff)%terms(iterm)%strain(idisp2_strain)
1211 : : ! Strain case
1212 : : ! Set the power_strain of the strain:
1213 : 2187237792 : power_strain = coefficients(icoeff)%terms(iterm)%power_strain(idisp2_strain)
1214 : : ! Accumulate energy forces
1215 : 2187237792 : tmp2 = tmp2 * (strain(idir2))**power_strain
1216 : : ! Accumulate stress for each strain (\sum (Y(eta_2)^Y-1(eta_2)^Z+...))
1217 : 2187237792 : tmp3 = tmp3 * (strain(idir2))**power_strain
1218 : : else
1219 : 3232495536 : idir2 = coefficients(icoeff)%terms(iterm)%direction(idisp2)
1220 [ + + ]: 12929982144 : cell_atoma2=coefficients(icoeff)%terms(iterm)%cell(:,1,idisp2)
1221 [ + - + + : 3232495536 : if(cell_atoma2(1)/=0.or.cell_atoma2(2)/=0.or.cell_atoma2(3)/=0) then
+ + ]
1222 : 64800 : cell_atoma2(1) = i1 + cell_atoma2(1)
1223 : 64800 : cell_atoma2(2) = i2 + cell_atoma2(2)
1224 : 64800 : cell_atoma2(3) = i3 + cell_atoma2(3)
1225 : 64800 : call getPBCIndexes_supercell(cell_atoma2(1:3),sc_size(1:3))
1226 : : ! index of the first atom (position in the supercell and direction)
1227 : : ! if the cell of the atom a is not 0 0 0 (may happen)
1228 : : ia2 = (cell_atoma2(1)-1)*sc_size(2)*sc_size(3)*natom_uc+&
1229 : : & (cell_atoma2(2)-1)*sc_size(3)*natom_uc+&
1230 : : & (cell_atoma2(3)-1)*natom_uc+&
1231 : 64800 : & coefficients(icoeff)%terms(iterm)%atindx(1,idisp2)
1232 : : else
1233 : : ! index of the first atom (position in the supercell and direction)
1234 : 3232430736 : ia2 = ii + coefficients(icoeff)%terms(iterm)%atindx(1,idisp2)
1235 : : end if
1236 : :
1237 [ + + ]: 12929982144 : cell_atomb2= coefficients(icoeff)%terms(iterm)%cell(:,2,idisp2)
1238 : :
1239 [ + + + + : 3232495536 : if(cell_atomb2(1)/=0.or.cell_atomb2(2)/=0.or.cell_atomb2(3)/=0) then
+ + ]
1240 : : ! indexes of the cell2 (with PBC)
1241 : 2010815736 : cell_atomb2(1) = i1 + cell_atomb2(1)
1242 : 2010815736 : cell_atomb2(2) = i2 + cell_atomb2(2)
1243 : 2010815736 : cell_atomb2(3) = i3 + cell_atomb2(3)
1244 : 2010815736 : call getPBCIndexes_supercell(cell_atomb2(1:3),sc_size(1:3))
1245 : :
1246 : : ! index of the second atom in the (position in the supercell)
1247 : : ib2 = (cell_atomb2(1)-1)*sc_size(2)*sc_size(3)*natom_uc+&
1248 : : & (cell_atomb2(2)-1)*sc_size(3)*natom_uc+&
1249 : : & (cell_atomb2(3)-1)*natom_uc+&
1250 : 2010815736 : & coefficients(icoeff)%terms(iterm)%atindx(2,idisp2)
1251 : : else
1252 : 1221679800 : ib2 = ii + coefficients(icoeff)%terms(iterm)%atindx(2,idisp2)
1253 : : end if
1254 : :
1255 : 3232495536 : disp1 = disp(idir2,ia2)
1256 : 3232495536 : disp2 = disp(idir2,ib2)
1257 : :
1258 : : ! Set the power_disp of the displacement:
1259 : 3232495536 : power_disp = coefficients(icoeff)%terms(iterm)%power_disp(idisp2)
1260 : 3232495536 : tmp2 = tmp2 * (disp1-disp2)**power_disp
1261 : 3232495536 : tmp3 = tmp3 * (disp1-disp2)**power_disp
1262 : :
1263 : : end if
1264 : : end if
1265 : : end do
1266 : :
1267 [ + + ]: 6148709792 : if(idisp1 > ndisp)then
1268 : : ! Accumule stress tensor
1269 : 1477156968 : strten(idir1) = strten(idir1) + coeff * weight * tmp3
1270 : : else
1271 : : ! Accumule forces
1272 : 2679703296 : fcart(idir1,ia1) = fcart(idir1,ia1) + coeff * weight * tmp2
1273 : 2679703296 : fcart(idir1,ib1) = fcart(idir1,ib1) - coeff * weight * tmp2
1274 : : end if
1275 : : end do
1276 : :
1277 : 1991849528 : energy_coeff(icoeff) = energy_coeff(icoeff) + coeff * weight * tmp1
1278 : : ! accumule energy
1279 : 2163492856 : energy = energy + coeff * weight * tmp1
1280 : :
1281 : : end do
1282 : : end do
1283 : : end do
1284 : :
1285 : :
1286 : : ! MPI_SUM
1287 : 10265 : call xmpi_sum(energy, comm, ierr)
1288 : 10265 : call xmpi_sum(fcart , comm, ierr)
1289 : 10265 : call xmpi_sum(strten , comm, ierr)
1290 : :
1291 : : !Write to anharmonic_energy_terms.out ORIGINAL
1292 : 10265 : INQUIRE(FILE=name_file,OPENED=file_opened,number=unit_out)
1293 [ + + ]: 10265 : if(file_opened .eqv. .TRUE.)then
1294 [ + + ]: 6161 : do icoeff=1,ncoeff
1295 : 6060 : call xmpi_sum(energy_coeff(icoeff), comm, ierr)
1296 : : ! Marcus write energy contributions of anharmonic terms to file
1297 [ + + ]: 6161 : if(icoeff <ncoeff)then
1298 : 5959 : write(unit_out,'(A,1ES24.16)',advance='no') ' ',energy_coeff(icoeff)
1299 : : else if(icoeff==ncoeff)then
1300 : 101 : write(unit_out,'(A,1ES24.16)',advance='yes') ' ',energy_coeff(icoeff)
1301 : : end if
1302 : : enddo
1303 : : end if
1304 : :
1305 : :
1306 [ + - ]: 10265 : end subroutine polynomial_coeff_evaluate
1307 : : !!***
1308 : :
1309 : 1 : function get_crystal_cutoff(crystal) result(cutoff)
1310 : : type(crystal_t), intent(inout) :: crystal
1311 : : real(dp) :: cutoff
1312 : : integer :: ii
1313 : : ! XXX: hexu: why cutoff + sum of axis lengths?
1314 : 1 : cutoff=0.1
1315 [ + + ]: 4 : do ii=1,3
1316 : : cutoff = cutoff + sqrt(crystal%rprimd(ii,1)**2 + &
1317 : : & crystal%rprimd(ii,2)**2 + &
1318 : 4 : & crystal%rprimd(ii,3)**2)
1319 : : enddo
1320 : 1 : end function get_crystal_cutoff
1321 : :
1322 : 45 : subroutine prepare_for_getList(crystal,sc_size, dist, cell, natom, nsym, nrpt, range_ifc, symbols )
1323 : : type(crystal_t), intent(inout) :: crystal
1324 : : integer, intent(in) :: sc_size(3)
1325 : : integer, intent(inout) :: natom, nsym, nrpt
1326 : : integer, allocatable, intent(inout) :: cell(:, :)
1327 : : character(len=5),allocatable, intent(inout) :: symbols(:)
1328 : :
1329 : 45 : real(dp), allocatable:: xcart(:, :), xred(:,:), rpt(:, :)
1330 : : real(dp), allocatable:: dist(:, :, :, :)
1331 : : real(dp) :: rprimd(3,3), range_ifc(3)
1332 : : integer :: ii, ncell(3), lim1, lim2, lim3
1333 : : real(dp) :: norm
1334 : : integer :: irpt, irpt_ref, r1, r2, r3, ia, ib
1335 : :
1336 : 45 : natom = crystal%natom
1337 : 45 : nsym = crystal%nsym
1338 [ + + + + ]: 585 : rprimd = crystal%rprimd
1339 : :
1340 [ - + + - : 135 : ABI_MALLOC(xcart,(3,natom))
+ - + - -
+ - + ]
1341 [ - + + - : 90 : ABI_MALLOC(xred,(3,natom))
- + - + ]
1342 [ + + + + ]: 945 : xcart(:,:) = crystal%xcart(:,:)
1343 [ + + + + ]: 945 : xred(:,:) = crystal%xred(:,:)
1344 : :
1345 : : !Compute the max range of the ifc with respect to the trainning set
1346 : 45 : range_ifc(:) = zero
1347 [ + + ]: 180 : do ii=1,3
1348 : 135 : norm = sqrt(rprimd(ii,1)**2+ rprimd(ii,2)**2+rprimd(ii,3)**2)
1349 : 180 : range_ifc(ii) = range_ifc(ii) + norm * sc_size(ii) / 2.0
1350 : : end do
1351 : :
1352 : :
1353 : : !compute new ncell
1354 : 45 : ncell = sc_size
1355 : 45 : lim1=((ncell(1)/2)) + 1
1356 : 45 : lim2=((ncell(2)/2)) + 1
1357 : 45 : lim3=((ncell(3)/2)) + 1
1358 [ + + ]: 45 : if(mod(ncell(1),2)/=0) lim1=lim1+1
1359 [ + + ]: 45 : if(mod(ncell(2),2)/=0) lim2=lim2+1
1360 [ + + ]: 45 : if(mod(ncell(3),2)/=0) lim3=lim3+1
1361 : 45 : nrpt=(2*lim1+1)*(2*lim2+1)*(2*lim3+1)
1362 : :
1363 : 45 : ncell(1) = 2*lim1+1
1364 : 45 : ncell(2) = 2*lim2+1
1365 : 45 : ncell(3) = 2*lim3+1
1366 : :
1367 : : !Build the rpt point
1368 [ - + + - : 135 : ABI_MALLOC(rpt,(3,nrpt))
+ - + - -
+ - + ]
1369 [ - + + - : 135 : ABI_MALLOC(cell,(3,nrpt))
+ - + - -
+ - + -
+ ]
1370 : :
1371 : : !WARNING:
1372 : : !Put the reference cell into the first element
1373 : : !the code will first deal with the atoms of the first cell
1374 : 45 : irpt = 1
1375 : 45 : irpt_ref = 1
1376 [ + + ]: 180 : rpt(:,1) = zero
1377 [ + + ]: 180 : cell(:,irpt)=0
1378 : : !Fill other rpt:
1379 [ + + ]: 270 : do r1=lim1,-lim1,-1
1380 [ + + ]: 1395 : do r2=lim2,-lim2,-1
1381 [ + + ]: 6975 : do r3=lim3,-lim3,-1
1382 [ + + ]: 5625 : if(r1==0.and.r2==0.and.r3==0) then
1383 : : cycle
1384 : : end if
1385 : 5580 : irpt=irpt+1
1386 : 5580 : rpt(1,irpt)=r1*rprimd(1,1)+r2*rprimd(1,2)+r3*rprimd(1,3)
1387 : 5580 : rpt(2,irpt)=r1*rprimd(2,1)+r2*rprimd(2,2)+r3*rprimd(2,3)
1388 : 5580 : rpt(3,irpt)=r1*rprimd(3,1)+r2*rprimd(3,2)+r3*rprimd(3,3)
1389 : 6750 : cell(1,irpt)=r1;cell(2,irpt)=r2;cell(3,irpt)=r3
1390 : : end do
1391 : : end do
1392 : : end do
1393 : :
1394 [ + - ]: 45 : if (.not.(allocated(symbols))) then
1395 [ - + + - : 135 : ABI_MALLOC(symbols,(natom))
+ - + - -
+ - + -
+ ]
1396 : : endif
1397 : : call symbols_crystal(crystal%natom,crystal%ntypat,crystal%npsp,&
1398 : 45 : & symbols,crystal%typat,crystal%znucl)
1399 : :
1400 : : !Compute the distances between atoms
1401 : : !Now dist(3,ia,ib,irpt) contains the distance from atom ia to atom ib in unit cell irpt.
1402 [ - + + - : 225 : ABI_MALLOC(dist,(3,natom,natom,nrpt))
- + + - -
+ + - + -
+ - - + -
+ - + ]
1403 [ + + + + : 596295 : dist = zero
+ + + + ]
1404 [ + + ]: 270 : do ia=1,natom
1405 [ + + ]: 1395 : do ib=1,natom
1406 [ + + ]: 141975 : do irpt=1,nrpt
1407 : 140625 : dist(1,ia,ib,irpt) = xcart(1,ib)-xcart(1,ia)+rpt(1,irpt)
1408 : 140625 : dist(2,ia,ib,irpt) = xcart(2,ib)-xcart(2,ia)+rpt(2,irpt)
1409 : 141750 : dist(3,ia,ib,irpt) = xcart(3,ib)-xcart(3,ia)+rpt(3,irpt)
1410 : : end do
1411 : : end do
1412 : : end do
1413 : :
1414 : : !call polynomial_coeff_getList(cell,crystal,dist,list_symcoeff,list_symstr,&
1415 : : ! & natom,nstr_sym,ncoeff_sym,nrpt,range_ifc,cutoff,sc_size=sc_size,&
1416 : : ! & fit_iatom=fit_iatom_in)
1417 [ + - ]: 45 : ABI_SFREE(xcart)
1418 [ + - ]: 45 : ABI_SFREE(xred)
1419 [ + - ]: 45 : ABI_SFREE(rpt)
1420 : 45 : end subroutine prepare_for_getList
1421 : :
1422 : : !!****f* m_polynomial_coeff/polynomial_coeff_getList
1423 : : !!
1424 : : !! NAME
1425 : : !! polynomial_coeff_getList
1426 : : !!
1427 : : !! FUNCTION
1428 : : !! Get the list of all the possible coefficients for the polynome
1429 : : !!
1430 : : !! INPUTS
1431 : : !! cell(3,nrpt) = indexes of the cells into the supercell (-1 -1 -1, 0 0 0 ...)
1432 : : !! dist(3,natom,natom,nrpt) = distance between atoms atm1 is in the cell 0 0 0
1433 : : !! atm2 is in the nrpt cell (see cell(3,nrpt))
1434 : : !! for each component x,y and z
1435 : : !! crystal<type(crystal_t)> = datatype with all the information for the crystal
1436 : : !! natom = number of atoms in the unit cell
1437 : : !! nrpt = number of cell in the supercell
1438 : : !!
1439 : : !! OUTPUT
1440 : : !! list_symcoeff(6,ncoeff_sym,nsym) = array with the list of the coefficients,
1441 : : !! for each coefficients (ncoeff_sym), we store the symmetrics(nsym)
1442 : : !! the 6th first dimensions are :
1443 : : !! 1 = direction of the IFC
1444 : : !! 2 = index of the atom number 1 (1=>natom)
1445 : : !! 3 = index of the atom number 2 (1=>natom)
1446 : : !! 4 = indexes of the cell of the second atom
1447 : : !! (the atom number 1 is always in the cell 0 0 0)
1448 : : !! 5 = weight of the term (-1 or 1)
1449 : : !! 6 = indexes of the symmetric
1450 : : !! list_symstr(nstr_sym,nsym) = array with the list of the strain and the symmetrics
1451 : : !! nstr_sym = number of coefficient for the strain
1452 : : !! ncoeff_sym = number of coefficient for the IFC
1453 : : !! range_ifc(3) = maximum cut-off for the inter atomic forces constants in each direction
1454 : : !! sc_size(3) = optional,size of the supercell used for the fit.
1455 : : !! For example if you want to fit 2x2x2 cell the interation
1456 : : !! Sr-Ti and Sr-Ti[2 0 0] will be identical for the fit process
1457 : : !! If check_pbc is true we remove these kind of terms
1458 : : !!
1459 : : !! SOURCE
1460 : :
1461 : 29 : subroutine polynomial_coeff_getList(cell,crystal,dist,list_symcoeff,list_symstr,&
1462 : : & natom,nstr_sym,ncoeff_sym,nrpt,range_ifc,cutoff,sc_size,&
1463 : : & fit_iatom)
1464 : :
1465 : : !Arguments ------------------------------------
1466 : : !scalars
1467 : : integer,intent(in) :: natom,nrpt
1468 : : real(dp), intent(in) :: cutoff
1469 : : integer,intent(out) :: ncoeff_sym,nstr_sym
1470 : : integer,optional,intent(in):: fit_iatom
1471 : : !arrays
1472 : : integer,intent(in) :: cell(3,nrpt)
1473 : : real(dp),intent(in):: dist(3,natom,natom,nrpt)
1474 : : type(crystal_t), intent(in) :: crystal
1475 : : integer,allocatable,intent(out) :: list_symcoeff(:,:,:),list_symstr(:,:,:)
1476 : : integer,optional,intent(in) :: sc_size(3)
1477 : : real(dp),intent(in):: range_ifc(3)
1478 : : !Local variables-------------------------------
1479 : : !scalar
1480 : : integer :: ia,ib,icoeff,icoeff2,icoeff_tot,icoeff_tmp,idisy1,idisy2,ii
1481 : : integer :: ipesy1,ipesy2,isym,irpt,irpt3,irpt_ref,irpt_sym
1482 : : integer :: jj,jsym,mu,fit_iatom_in
1483 : : integer :: ncoeff,ncoeff2,ncoeff3,ncoeff_max,nu
1484 : : integer :: nsym,shift_atm1(3)
1485 : : integer :: shift_atm2(3)
1486 : : real(dp):: dist_orig,dist_sym,tolsym8
1487 : : logical :: found,check_pbc,possible
1488 : : !arrays
1489 : : integer :: isym_rec(3,3),isym_rel(3,3),sc_size_in(3)
1490 : : integer :: transl(3),min_range(3),max_range(3)
1491 : 29 : integer,allocatable :: blkval(:,:,:,:,:),list(:),list_symcoeff_tmp(:,:,:),list_symcoeff_tmp2(:,:,:)
1492 : 29 : integer,allocatable :: list_symstr_tmp(:,:,:),indsym(:,:,:) ,symrec(:,:,:),symrel(:,:,:),list_symcoeff_tmp3(:,:,:)
1493 : 29 : integer,allocatable :: index_irred(:)
1494 : 29 : real(dp),allocatable :: tnons(:,:)
1495 : 29 : real(dp),allocatable :: wkdist(:),distance(:,:,:)
1496 : : real(dp) :: difmin(3)
1497 : : real(dp) :: rprimd(3,3)
1498 : : real(dp) :: tratom(3)
1499 : : character(len=500) :: message
1500 : :
1501 : :
1502 : :
1503 : : !Initialisation of variables
1504 : 29 : irpt_sym = 0
1505 : 29 : nsym = crystal%nsym
1506 : : rprimd = crystal%rprimd
1507 : : !ABI_MALLOC(xcart,(3,natom))
1508 : : !ABI_MALLOC(xred,(3,natom))
1509 : : !xcart(:,:) = crystal%xcart(:,:)
1510 : : !xred(:,:) = crystal%xred(:,:)
1511 : 29 : ncoeff_max = nrpt*natom*natom*3*3
1512 : :
1513 : : !Found the ref cell
1514 : 29 : irpt_ref = 1
1515 [ + + ]: 3654 : do irpt=1,nrpt
1516 [ + + + + : 4553 : if(all(cell(:,irpt)==0))then
+ + ]
1517 : 29 : irpt_ref = irpt
1518 : : ! exit
1519 : : end if
1520 : : end do
1521 : :
1522 : : !Set the size of the interaction
1523 : 29 : check_pbc = .FALSE.
1524 : 29 : sc_size_in = 0
1525 : 29 : min_range = 0; max_range = 0
1526 [ + - ]: 29 : if(present(sc_size))then
1527 : 29 : sc_size_in = sc_size
1528 [ + + ]: 116 : do mu=1,3
1529 : 116 : call findBound_supercell(min_range(mu),max_range(mu),sc_size_in(mu))
1530 : : end do
1531 : : end if
1532 : :
1533 : : !Check which atom to fit, if not present do all atoms
1534 [ + + ]: 29 : if(present(fit_iatom))then
1535 : 28 : fit_iatom_in = fit_iatom
1536 : : else
1537 : : fit_iatom_in = -1
1538 : : endif
1539 [ - + ]: 29 : if(fit_iatom_in > natom)then
1540 : : write(message, '(3a)' )&
1541 : 0 : & 'fit_iatom cannot be greater than the number of atoms on the reference unit cell',ch10,&
1542 : 0 : & 'Action: Change input'
1543 : 0 : ABI_ERROR(message)
1544 : : end if
1545 : :
1546 : : !Obtain a list of rotated atom labels:
1547 [ - + + - : 116 : ABI_MALLOC(indsym,(4,nsym,natom))
- + + - +
- + - - +
- + ]
1548 [ - + + - : 87 : ABI_MALLOC(symrec,(3,3,nsym))
+ - + - -
+ - + ]
1549 [ - + + - : 58 : ABI_MALLOC(symrel,(3,3,nsym))
- + - + ]
1550 [ - + + - : 87 : ABI_MALLOC(tnons,(3,nsym))
+ - + - -
+ - + ]
1551 [ + - + - : 18154 : symrec = crystal%symrec
+ - - + -
- - - - -
- - + + +
+ + + ]
1552 [ + - + - : 18154 : symrel = crystal%symrel
+ - - + -
- - - - -
- - + + +
+ + + ]
1553 [ + - + - : 5626 : tnons = crystal%tnons
- + - - -
- - - + +
+ + ]
1554 : :
1555 : 29 : tolsym8=tol13
1556 : : call symatm(indsym,natom,nsym,symrec,tnons,&
1557 : 29 : & tolsym8,crystal%typat,crystal%xred)
1558 [ - + + - : 174 : ABI_MALLOC(blkval,(3,natom,3,natom,nrpt))
+ - - + +
- - + + -
+ - + - -
+ - + ]
1559 [ - + + - : 87 : ABI_MALLOC(list,(natom*nrpt))
+ - + - -
+ - + ]
1560 [ - + + - : 116 : ABI_MALLOC(list_symcoeff_tmp,(5,ncoeff_max,nsym))
- + + - +
- + - - +
- + ]
1561 [ - + + - : 87 : ABI_MALLOC(wkdist,(natom*nrpt))
+ - + - -
+ - + ]
1562 : :
1563 : : !1-Fill strain list
1564 [ - + + - : 116 : ABI_MALLOC(list_symstr_tmp,(6,nsym,2))
+ - + - +
- - + -
+ ]
1565 [ + + + + : 19575 : list_symstr_tmp = 1
+ + ]
1566 [ + + ]: 203 : do ia=1,6
1567 [ + - ]: 174 : if(list_symstr_tmp(ia,1,1)==0)cycle
1568 : : ! Transform the voigt notation
1569 [ + + ]: 174 : if(ia<=3)then
1570 : : mu=ia;nu=ia
1571 : : else
1572 : : select case(ia)
1573 : : case(4)
1574 : 29 : mu=2;nu=3
1575 : : case(5)
1576 : 29 : mu=1;nu=3
1577 : : case(6)
1578 [ + + + ]: 87 : mu=1;nu=2
1579 : : end select
1580 : : end if
1581 [ + + ]: 8555 : do isym=1,nsym
1582 : : ! Get the symmetry matrix
1583 [ + + + + ]: 108576 : isym_rel(:,:) = crystal%symrel(:,:,isym)
1584 [ + + ]: 33582 : do idisy1=1,3
1585 [ + + ]: 108576 : do idisy2=1,3
1586 [ + + + + ]: 100224 : if((isym_rel(mu,idisy1)/=0.and.isym_rel(nu,idisy2)/=0)) then
1587 : : ! Transform to the voig notation
1588 [ + + ]: 8352 : if(idisy1==idisy2)then
1589 : 4176 : list_symstr_tmp(ia,isym,1) = idisy1
1590 : 4176 : list_symstr_tmp(ia,isym,2) = isym_rel(mu,idisy1)
1591 : : else
1592 [ + + ]: 4176 : if(idisy1==1.or.idisy2==1)then
1593 [ + + ]: 2784 : if(idisy1==2.or.idisy2==2)then
1594 : 1392 : list_symstr_tmp(ia,isym,1) = 6
1595 : : end if
1596 [ + + ]: 2784 : if(idisy1==3.or.idisy2==3)then
1597 : 1392 : list_symstr_tmp(ia,isym,1) = 5
1598 : : end if
1599 : : else
1600 : 1392 : list_symstr_tmp(ia,isym,1) = 4
1601 : : end if
1602 : : end if
1603 : 8352 : list_symstr_tmp(ia,isym,2) = isym_rel(mu,idisy1) * isym_rel(nu,idisy2)
1604 : : end if
1605 : : end do
1606 : : end do
1607 : : ! Remove the symetric
1608 : : ! if(list_symstr_tmp(ia,isym,1) > ia) then
1609 : : ! list_symstr_tmp(list_symstr_tmp(ia,isym,1),:,1) = 0
1610 : : ! end if
1611 : : end do
1612 : : end do
1613 : :
1614 : : !Count the number of strain and transfert into the final array
1615 : 29 : nstr_sym = 0
1616 [ + + ]: 203 : do ia=1,6
1617 [ + - ]: 203 : if(list_symstr_tmp(ia,1,1)/=0) nstr_sym = nstr_sym + 1
1618 : : end do
1619 : :
1620 [ - + ]: 29 : if(allocated(list_symstr))then
1621 : 0 : ABI_SFREE(list_symstr)
1622 : : end if
1623 [ - + + - : 145 : ABI_MALLOC(list_symstr,(nstr_sym,nsym,2))
- + + - +
- + - + -
- + - + -
+ ]
1624 : :
1625 : 29 : icoeff_tmp = 1
1626 [ + + ]: 203 : do ia=1,6
1627 [ + - ]: 203 : if(list_symstr_tmp(ia,1,1)/=0) then
1628 [ + + + + ]: 17226 : list_symstr(icoeff_tmp,:,:) = list_symstr_tmp(ia,:,:)
1629 : 174 : icoeff_tmp = icoeff_tmp + 1
1630 : : end if
1631 : : end do
1632 : : !END STRAIN
1633 : :
1634 : : !Compute the distance between each atoms. Indeed the dist array contains the difference of
1635 : : !cartesian coordinate for each direction
1636 : : ! dist: vector between atom a, (0,0,0) and atom b (rpt)
1637 : :
1638 [ - + + - : 145 : ABI_MALLOC(distance,(natom,natom,nrpt))
- + + - -
+ + - + -
+ - - + -
+ ]
1639 : : ! Fortran 2008: distance = norm2(dist, dim=1)
1640 [ + + ]: 174 : do ia=1,natom
1641 [ + + ]: 899 : do ib=1,natom
1642 [ + + ]: 91495 : do irpt=1,nrpt
1643 : : distance(ia,ib,irpt) = ((dist(1,ia,ib,irpt))**2+(dist(2,ia,ib,irpt))**2+&
1644 : 91350 : & (dist(3,ia,ib,irpt))**2)**0.5
1645 : : end do
1646 : : end do
1647 : : end do
1648 : :
1649 : :
1650 : : !Set to one blkval, all the coeff have to be compute
1651 [ + + + + : 1163654 : blkval = 1
+ + + + +
+ ]
1652 : 29 : icoeff = 1
1653 : 29 : icoeff_tot = 1
1654 [ + + + + : 234901421 : list_symcoeff_tmp = 0
+ + ]
1655 : :
1656 : : !2-Fill atom list
1657 : : !Big loop over generic atom
1658 [ + + ]: 174 : do ia=1,natom
1659 [ + + ]: 290 : wkdist(:)=reshape(distance(ia,:,:),(/natom*nrpt/))
1660 [ + + ]: 90770 : do ii=1,natom*nrpt
1661 : 90770 : list(ii)=ii
1662 : : end do
1663 : : ! FIXME: hexu: I think this should be improved.
1664 : : ! It seems that it depends on the specific order of the
1665 : : ! cell, and the implementation of sort algorithm.
1666 : : ! The sort should preserve the order if two distances are the same.
1667 : : ! In the future, the order of the cell might need to be unified
1668 : : ! Also it seems to depend on the cubic cell.
1669 : 145 : call sort_dp(natom*nrpt,wkdist,list,tol8)
1670 [ + + ]: 90799 : do ii=1,natom*nrpt
1671 : : ! Get the irpt and ib
1672 : 90625 : irpt=(list(ii)-1)/natom+1
1673 : 90625 : ib=list(ii)-natom*(irpt-1)
1674 : 90625 : possible = .true.
1675 : : !Old way with the cut off
1676 [ + + ]: 90625 : if(cutoff < distance(ia,ib,irpt))then
1677 : : possible = .false.
1678 : : else
1679 : : ! in each direction jj, the d component d(jj) should be smaller than range_ifc(jj)
1680 [ + + ]: 66196 : do jj=1,3
1681 : : ! if(abs(dist(jj,ia,ib,irpt)) - range_ifc(jj) > tol10.or.&
1682 : : ! & abs(abs(dist(jj,ia,ib,irpt)) - range_ifc(jj)) < tol10)then
1683 [ + + ]: 66196 : if(abs(dist(jj,ia,ib,irpt)) - range_ifc(jj) > tol10)then
1684 : 19296 : possible = .false.
1685 : : end if
1686 : : end do
1687 : : endif
1688 : :
1689 : : ! If this distance is superior to the cutoff, we don't compute that term
1690 [ + + ]: 16549 : if(.not.possible)then
1691 [ + + + + ]: 1129908 : blkval(:,ia,:,ib,irpt)= 0
1692 : : ! remove duplication: if in reference cell, only consider ab, not ba.
1693 [ + + + + : 88116 : if(irpt==irpt_ref)blkval(:,ib,:,ia,irpt)= 0
+ + ]
1694 : : ! Stop the loop
1695 : : cycle
1696 : : end if
1697 : :
1698 : : ! If this coefficient is not possible, we cycle...
1699 [ + - + - : 3709 : if (all(blkval(:,ia,:,ib,irpt)==0)) cycle
- + + - ]
1700 : :
1701 : : ! Save the distance between the two atoms for futur checks
1702 : 3709 : dist_orig = (dist(1,ia,ib,irpt)**2+dist(2,ia,ib,irpt)**2+dist(3,ia,ib,irpt)**2)**0.5
1703 : :
1704 [ + + ]: 14981 : do mu=1,3
1705 [ + + ]: 48217 : do nu=1,3
1706 : : ! Check if : - The coefficient is not yet compute
1707 : : ! - The directions are the same
1708 : : ! - The atoms are not equivalent
1709 [ + + ]: 33381 : if (mu/=nu) then
1710 : 22254 : blkval(mu,ia,nu,ib,irpt)=0
1711 : 22254 : blkval(nu,ia,mu,ib,irpt)=0
1712 : 22254 : cycle
1713 : : end if
1714 : : ! Pass if the atoms are identical and in the ref cell
1715 [ + + ]: 11127 : if(irpt==irpt_ref.and.ia==ib) then
1716 : 435 : blkval(mu,ia,nu,ib,irpt)=0
1717 : 435 : blkval(nu,ib,mu,ia,irpt)=0
1718 : 435 : cycle
1719 : : end if
1720 [ + - ]: 10692 : if(blkval(mu,ia,nu,ib,irpt)==1)then
1721 : : ! Loop over symmetries
1722 [ + + ]: 523908 : do isym=1,nsym
1723 : : ! Get the symmetry matrix for this sym
1724 [ + + + + ]: 6671808 : isym_rec(:,:) = crystal%symrec(:,:,isym)
1725 [ + + + + ]: 6671808 : isym_rel(:,:) = crystal%symrel(:,:,isym)
1726 : : ! Get the corresponding atom and shift with the symetries
1727 : : ! For atom 1
1728 : 513216 : ipesy1 = indsym(4,isym,ia)
1729 [ + + ]: 2052864 : shift_atm1 = indsym(1:3,isym,ia)
1730 : : ! And atom 2
1731 [ + + ]: 2052864 : do jj=1,3 ! Apply transformation to original coordinates.
1732 : : tratom(jj) = dble(isym_rec(1,jj))*(crystal%xred(1,ib)+cell(1,irpt)-tnons(1,isym))&
1733 : : & +dble(isym_rec(2,jj))*(crystal%xred(2,ib)+cell(2,irpt)-tnons(2,isym))&
1734 : 2052864 : & +dble(isym_rec(3,jj))*(crystal%xred(3,ib)+cell(3,irpt)-tnons(3,isym))
1735 : :
1736 : : end do
1737 : :
1738 : : ! Find symmetrically equivalent atom
1739 : : call symchk(difmin,ipesy2,natom,tratom,transl,crystal%typat(ib),&
1740 : 513216 : & crystal%typat,crystal%xred(:,:))
1741 : :
1742 : : ! Put information into array indsym: translations and label
1743 [ + + ]: 2052864 : shift_atm2(:)= transl(:) - shift_atm1(:)
1744 : : found = .false.
1745 [ + + ]: 64665216 : do irpt3=1,nrpt
1746 : : if(cell(1,irpt3)==shift_atm2(1).and.&
1747 [ + + + + : 64152000 : & cell(2,irpt3)==shift_atm2(2).and.&
+ + ]
1748 : 513216 : & cell(3,irpt3)==shift_atm2(3))then
1749 : 64152000 : found = .true.
1750 : 64152000 : irpt_sym = irpt3
1751 : : end if
1752 : : end do
1753 : :
1754 : : ! Check the distance
1755 : : dist_sym = (dist(1,ipesy1,ipesy2,irpt_sym)**2+&
1756 : : & dist(2,ipesy1,ipesy2,irpt_sym)**2+&
1757 : 513216 : & dist(3,ipesy1,ipesy2,irpt_sym)**2)**0.5
1758 [ - + ]: 513216 : if(abs(dist_orig - dist_sym) > tol10)then
1759 : : write(message, '(a,i0,2a,I0,a,es15.8,2a,es15.8,2a)' )&
1760 : 0 : & 'The distance between the atoms for the coefficient number ',icoeff,ch10,&
1761 : 0 : & 'with the symmetry ',isym,' is ',dist_sym,ch10,'but the original distance is',&
1762 : 0 : & dist_orig,ch10,&
1763 : 0 : & 'Action: Contact abinit group'
1764 : 0 : ABI_BUG(message)
1765 : : end if
1766 : : ! Now that a symmetric perturbation has been obtained,
1767 : : ! including the expression of the symmetry matrix, see
1768 : : ! if the symmetric perturbations are available
1769 [ + + ]: 2576772 : do idisy1=1,3
1770 [ + + ]: 6671808 : do idisy2=1,3
1771 [ + + ]: 6158592 : if (idisy1/=idisy2) then
1772 : : ! Remove this term (is not computed)
1773 : : ! Also remove opposite term... (Srx-Tix) = (Tix-Srx)
1774 : 3079296 : blkval(idisy1,ipesy1,idisy2,ipesy2,irpt_sym) = 0
1775 : 3079296 : blkval(idisy2,ipesy1,idisy1,ipesy2,irpt_sym) = 0
1776 : 3079296 : cycle
1777 : : else
1778 [ + + + - ]: 1539648 : if(isym_rel(mu,idisy1)/=0.and.isym_rel(nu,idisy2)/=0)then
1779 [ + - + + : 513216 : if(.not.found.or.(irpt_sym==irpt_ref.and.ipesy1==ipesy2)) then
- + ]
1780 : : ! Remove this term (is not computed) Sr-Sr or not include in the cell
1781 : : ! Also remove oposite term... (Srx-Tix) = (Ti-Srx)
1782 : 0 : blkval(idisy1,ipesy1,idisy2,ipesy2,irpt_sym) = 0
1783 : 0 : blkval(idisy2,ipesy2,idisy1,ipesy1,irpt_sym) = 0
1784 : 0 : cycle
1785 : : else
1786 : : ! Fill the list with the coeff and symmetric (need all symmetrics)
1787 [ + + ]: 2566080 : list_symcoeff_tmp(1:4,icoeff,isym)=(/idisy1,ipesy1,ipesy2,irpt_sym/)
1788 : : ! Check the sign
1789 [ - + ]: 513216 : if(isym_rel(mu,idisy1)/=isym_rel(nu,idisy2))then
1790 : : write(message, '(a,i0,a,I0,4a)' )&
1791 : 0 : & 'The sign of coefficient number ',icoeff,' with the symmetry ',isym,ch10,&
1792 : 0 : & 'can not be found... Something is going wrong',ch10,&
1793 : 0 : & 'Action: Contact abinit group'
1794 : 0 : ABI_BUG(message)
1795 : : end if
1796 : 513216 : list_symcoeff_tmp(5,icoeff,isym)= isym_rel(nu,idisy2)
1797 : : end if
1798 : : end if
1799 : : end if
1800 : : end do
1801 : : end do
1802 : : end do ! end loop sym
1803 : 10692 : icoeff = icoeff + 1
1804 : : end if
1805 : : ! This coeff is now computed
1806 : 21819 : blkval(mu,ia,nu,ib,irpt)= 0
1807 : : end do ! end loop nu
1808 : : end do ! end loop mu
1809 : : end do ! end loop ii: iatom and irpt.
1810 : : end do ! end loop ia
1811 : :
1812 : : !Reset the output
1813 [ - + ]: 29 : ABI_SFREE(list_symcoeff)
1814 : :
1815 [ + - ]: 29 : ABI_SFREE(distance)
1816 : :
1817 : : !Transfert the final array with all the coefficients
1818 : : !With this array, we can access to all the terms presents
1819 : : !ncoeff1 + symetrics
1820 : : !first dimension is 4 (mu,ia,ib,irpt)
1821 : : !irpt is the index of the cell of the atom ib in the cell array
1822 : : !example cell(:,irpt=12) can be (-1 0 -2). The cell of ia is
1823 : : !always 0 0 0
1824 : : !Transfert the final array for the list of irreductible coeff and symetries
1825 : : !With this array, we can access to the irretuctible coefficients (ncoeff1) and
1826 : : !all the symetrics of these coefficients (nsym)
1827 : : !first dimension is 5 (mu,ia,ib,irpt,icoeff)
1828 : : !icoeff is the position of this coefficients in the list_fullcoeff array
1829 : :
1830 : : !1/ step remove the zero coeff in this array
1831 : 29 : ncoeff = 0
1832 [ + + ]: 815654 : do icoeff = 1,ncoeff_max
1833 [ + + + + : 4840319 : if(.not.(all(list_symcoeff_tmp(:,icoeff,1)==0)))then
+ + ]
1834 : 10692 : ncoeff = ncoeff + 1
1835 : : end if
1836 : : end do
1837 : :
1838 [ + + + - : 116 : ABI_MALLOC(list_symcoeff_tmp2,(6,ncoeff,nsym))
- + + - +
- + + - +
- + ]
1839 [ + + + + : 3593933 : list_symcoeff_tmp2 = 0
+ + ]
1840 : 29 : icoeff = 0
1841 [ + + ]: 815654 : do icoeff_tmp = 1,ncoeff_max
1842 [ + + + + : 4840319 : if(.not.(all(list_symcoeff_tmp(:,icoeff_tmp,1)==0)))then
+ + ]
1843 : 10692 : icoeff = icoeff + 1
1844 [ + + + + ]: 3089988 : list_symcoeff_tmp2(1:5,icoeff,:) = list_symcoeff_tmp(1:5,icoeff_tmp,:)
1845 : : end if
1846 : : end do
1847 : :
1848 : :
1849 : : !2/ set the dimension six of list_symcoeff_tmp2(6,icoeffs,1)
1850 : : ! and check is a symetric coeff is not coresspondig to an other
1851 : : ! one, in this case we set this coeff to 0
1852 : : ! ncoeff2 = zero
1853 [ + + ]: 10721 : do icoeff = 1,ncoeff
1854 : : ! found the index of each coeff in list_fullcoeff
1855 [ + + ]: 523937 : do isym = 1,nsym
1856 : : icoeff2 = getCoeffFromList(list_symcoeff_tmp2(:,:,1),&
1857 : : & list_symcoeff_tmp2(2,icoeff,isym),&
1858 : : & list_symcoeff_tmp2(3,icoeff,isym),&
1859 : : & list_symcoeff_tmp2(4,icoeff,isym),&
1860 : : & list_symcoeff_tmp2(1,icoeff,isym),&
1861 : 513216 : & ncoeff)
1862 : 523908 : list_symcoeff_tmp2(6,icoeff,isym) = icoeff2
1863 : : end do
1864 : : end do
1865 : :
1866 : :
1867 : :
1868 : :
1869 : : ! filter the list_symcoeff_tmp2 to remove the terms which the equivalent cannot be found.
1870 : : ! This can happen when the cell is too small and the equivalent is out of the cell
1871 : :
1872 : : block
1873 : 58 : integer :: list_symcoeff_tmp3(6,ncoeff,nsym)
1874 : : integer :: counter_good, counter_bad, ncoeff_orig
1875 : 58 : integer :: coeffmap(ncoeff)
1876 : 29 : ncoeff_orig=ncoeff
1877 : 29 : counter_bad=0
1878 : 29 : counter_good=0
1879 [ + + ]: 10721 : coeffmap(:)=0
1880 : : ! copy the list_symcoeff_tmp2 to list_symcoeff_tmp3
1881 [ + + + + : 3593933 : list_symcoeff_tmp3 = list_symcoeff_tmp2
+ + ]
1882 : :
1883 [ + + ]: 10721 : do icoeff = 1,ncoeff
1884 [ + + + - : 523937 : if(any(list_symcoeff_tmp2(6,icoeff,:)==0)) then
- + ]
1885 : 0 : counter_bad=counter_bad+1
1886 : : else
1887 : 10692 : counter_good=counter_good+1
1888 : 10692 : coeffmap(icoeff)=counter_good
1889 : : endif
1890 : : end do
1891 : :
1892 [ + - ]: 29 : ABI_SFREE(list_symcoeff_tmp2)
1893 : 29 : ncoeff=ncoeff-counter_bad
1894 [ + + + - : 116 : ABI_MALLOC(list_symcoeff_tmp2,(6,ncoeff,nsym))
- + + - +
- + + - +
- + - + ]
1895 : 29 : counter_good=0
1896 [ + + ]: 10721 : do icoeff = 1,ncoeff_orig
1897 [ + + + - : 523937 : if(.not. (any(list_symcoeff_tmp3(6,icoeff,:)==0))) then
+ - ]
1898 : 10692 : counter_good=counter_good+1
1899 [ + + + + ]: 3603204 : list_symcoeff_tmp2(:,counter_good,:)=list_symcoeff_tmp3(:,icoeff,:)
1900 [ + + ]: 523908 : do isym = 1,nsym
1901 [ + - ]: 523908 : if(coeffmap(list_symcoeff_tmp2(6,counter_good,isym)) /= 0) then
1902 : 513216 : list_symcoeff_tmp2(6,counter_good,isym)=coeffmap(list_symcoeff_tmp2(6,counter_good,isym))
1903 : : else
1904 : 0 : ABI_BUG("coeffmap is not correct.")
1905 : : endif
1906 : : end do
1907 : : endif
1908 : : end do
1909 : :
1910 : :
1911 : :
1912 : : end block
1913 : :
1914 : :
1915 : : !2.5/do checks
1916 [ + + ]: 10721 : do icoeff = 1,ncoeff
1917 [ + + ]: 523937 : do isym = 1,nsym
1918 [ - + ]: 523908 : if(list_symcoeff_tmp2(6,icoeff,isym)==0)then
1919 : : write(message, '(a,i0,a,I0,4a)' )&
1920 : 0 : & 'The coefficient number ',icoeff,' with the symetrie ',isym,ch10,&
1921 : 0 : & 'has no equivalent',ch10,&
1922 : 0 : & 'Action: Contact abinit group'
1923 : 0 : ABI_BUG(message)
1924 : : else
1925 [ + + ]: 513216 : if(icoeff /= list_symcoeff_tmp2(6,icoeff,isym))then
1926 [ - + ]: 487632 : if(list_symcoeff_tmp2(1,icoeff,isym)/=&
1927 : : & list_symcoeff_tmp2(1,list_symcoeff_tmp2(6,icoeff,isym),1))then
1928 : : write(message, '(a,i0,a,I0,2a,I0,4a)' )&
1929 : 0 : & 'The coefficient number ',icoeff,' with the symetrie ',isym,ch10,&
1930 : 0 : & 'does not refer to the same coefficient ',list_symcoeff_tmp2(6,icoeff,1),ch10,&
1931 : 0 : & 'because the direction is different:',ch10,&
1932 : 0 : & 'Action: Contact abinit group'
1933 : 0 : ABI_BUG(message)
1934 : : end if
1935 [ - + ]: 487632 : if(list_symcoeff_tmp2(4,icoeff,isym)/=&
1936 : : & list_symcoeff_tmp2(4,list_symcoeff_tmp2(6,icoeff,isym),1))then
1937 : : write(message, '(a,i0,a,I0,2a,I0,4a)' )&
1938 : 0 : & 'The coefficient number ',icoeff,' with the symetrie ',isym,ch10,&
1939 : 0 : & 'does not refer to the same coefficient ',list_symcoeff_tmp2(6,icoeff,1),ch10,&
1940 : 0 : & 'because the cell is different',ch10,&
1941 : 0 : & 'Action: Contact abinit group'
1942 : 0 : ABI_BUG(message)
1943 : : end if
1944 : : if((list_symcoeff_tmp2(2,icoeff,isym)/=&
1945 : : & list_symcoeff_tmp2(2,list_symcoeff_tmp2(6,icoeff,isym),1).and.&
1946 : : & list_symcoeff_tmp2(3,icoeff,isym)/=&
1947 [ - + - - : 487632 : & list_symcoeff_tmp2(3,list_symcoeff_tmp2(6,icoeff,isym),1)).and.&
- - - - ]
1948 : : & (list_symcoeff_tmp2(2,icoeff,isym)/=&
1949 : : & list_symcoeff_tmp2(3,list_symcoeff_tmp2(6,icoeff,isym),1).and.&
1950 : : & list_symcoeff_tmp2(3,icoeff,isym)/=&
1951 : : & list_symcoeff_tmp2(2,list_symcoeff_tmp2(6,icoeff,isym),1)))then
1952 : : write(message, '(a,i0,a,I0,2a,I0,4a)' )&
1953 : 0 : & 'The coefficient number ',icoeff,' with the symetrie ',isym,ch10,&
1954 : 0 : & 'does not refer to the same coefficient ',list_symcoeff_tmp2(6,icoeff,1),ch10,&
1955 : 0 : & 'because the atoms different',ch10,&
1956 : 0 : & 'Action: Contact abinit group'
1957 : 0 : ABI_BUG(message)
1958 : : end if
1959 : : end if
1960 : : end if
1961 : : end do
1962 : : end do
1963 : :
1964 : :
1965 : : !Check if the atom 2 is not (with the PBC) is in the same cell than
1966 : : !the atom 1. For example if you want to fit 2x2x2 cell the interation
1967 : : !Sr-Ti and Sr-Ti[2 0 0] will be identical for the fit process...
1968 [ + + ]: 10721 : do icoeff = 1,ncoeff
1969 [ + - - + : 10692 : if(.not.(all(list_symcoeff_tmp2(:,icoeff,1)==0)))then
+ - ]
1970 [ + + ]: 42768 : do mu=1,3
1971 : : ! out of max_range(mu)
1972 [ + - - + ]: 32076 : if( -max_range(mu) > cell(mu,list_symcoeff_tmp2(4,icoeff,1)) .or. &
1973 : 10692 : & cell(mu,list_symcoeff_tmp2(4,icoeff,1)) > max_range(mu))then
1974 : :
1975 : : ! FIXME: hexu: cell(mu, irpt). But why -1?
1976 : : ! The if above looks sufficient.
1977 : : !if(cell(mu,list_symcoeff_tmp2(4,icoeff,1)) < -1)then
1978 [ # # # # ]: 0 : list_symcoeff_tmp2(:,icoeff,:)=0
1979 : : exit
1980 : : !end if
1981 : : end if
1982 : : end do
1983 : : end if
1984 : : !MS only keep terms with ia == fit_iatom or ib == fit_iatom specified in input
1985 [ + - - + : 10721 : if(.not.(all(list_symcoeff_tmp2(:,icoeff,1)==0)) .and. fit_iatom_in > 0)then
+ + ]
1986 [ + + ]: 9390 : if (list_symcoeff_tmp2(2,icoeff,1) /= fit_iatom_in) then !& LB
1987 : : !& list_symcoeff_tmp2(3,icoeff,1) /= fit_iatom_in)then LB
1988 [ + + ]: 53214 : list_symcoeff_tmp2(:,icoeff,1) = 0
1989 : : ! else if(list_symcoeff_tmp2(2,icoeff,1) /= fit_iatom_in) then ! &
1990 : : !& .and. list_symcoeff_tmp2(3,icoeff,1) == fit_iatom_in & !LB
1991 : : !& .and. any(cell(:,list_symcoeff_tmp2(4,icoeff,1)) /= 0))then !LB
1992 : : ! list_symcoeff_tmp2(:,icoeff,1) = 0
1993 : : endif
1994 : : endif
1995 : : !For Debugging keep only terms A_x-A_x[100] etc. comment if above
1996 : : ! if(.not.(all(list_symcoeff_tmp2(:,icoeff,1)==0)))then
1997 : : ! if (list_symcoeff_tmp2(2,icoeff,1) /= 1 .or. list_symcoeff_tmp2(3,icoeff,1) /= 1)then
1998 : : ! list_symcoeff_tmp2(:,icoeff,1) = 0
1999 : : ! endif !MS only keep terms with ia == 1
2000 : : ! endif
2001 : : end do
2002 : :
2003 : :
2004 : : !3/ Remove useless terms like opposites
2005 [ + + ]: 10721 : do icoeff = 1,ncoeff
2006 [ + + + + : 60612 : if(.not.(all(list_symcoeff_tmp2(:,icoeff,1)==0)) .and. & ! valid term.
+ + + + ]
2007 : 29 : & list_symcoeff_tmp2(2,icoeff,1) /= list_symcoeff_tmp2(3,icoeff,1))then ! on-site term. ! FIXME: hexu:I don't understand if iatom=jatom but rpt/=0, it is not onsite.
2008 [ + + ]: 78008 : do isym = 1,nsym
2009 : : !icoeff2 = list_symcoeff_tmp2(6,icoeff,isym)
2010 : : !if (icoeff2> icoeff)then
2011 : : ! list_symcoeff_tmp2(:,icoeff2,1) = 0
2012 : : !end if
2013 [ + + ]: 3745976 : do jsym=1,nsym
2014 : : ! FIXME: Should the rpt be -rpt??
2015 : : ! icoeff2 = getCoeffFromList(list_symcoeff_tmp2(:,:,jsym),&
2016 : : !& list_symcoeff_tmp2(3,icoeff,isym),&
2017 : : !& list_symcoeff_tmp2(2,icoeff,isym),&
2018 : : !& list_symcoeff_tmp2(4,icoeff,isym),& !rpt
2019 : : !& list_symcoeff_tmp2(1,icoeff,isym),&
2020 : : !& ncoeff)
2021 : : icoeff2 = getCoeffFromList(list_symcoeff_tmp2(:,:,jsym),&
2022 : : & list_symcoeff_tmp2(3,icoeff,isym),&
2023 : : & list_symcoeff_tmp2(2,icoeff,isym),&
2024 : : & list_symcoeff_tmp2(4,icoeff,isym),& !rpt
2025 : : & find_opposite_irpt(cell, list_symcoeff_tmp2(1,icoeff,isym)),&
2026 : 3667968 : & ncoeff)
2027 : :
2028 [ + + ]: 3744384 : if (icoeff2> icoeff)then
2029 [ + + ]: 2091530 : list_symcoeff_tmp2(:,icoeff2,1) = 0
2030 : : end if
2031 : : end do
2032 : : end do
2033 : : end if
2034 : : end do
2035 : :
2036 : : !4/ Recount the number of coeff after step 3
2037 : 29 : ncoeff2 = 0
2038 : 29 : isym = 0
2039 [ + + ]: 10721 : do icoeff = 1,ncoeff
2040 [ + + + + : 60641 : if(.not.(all(list_symcoeff_tmp2(:,icoeff,1)==0)))then
+ + ]
2041 : 2372 : ncoeff2 = ncoeff2 + 1
2042 : : end if
2043 : : end do
2044 : :
2045 : :
2046 [ + - ]: 29 : ABI_SFREE(list_symcoeff_tmp)
2047 [ + + + - : 87 : ABI_MALLOC(list_symcoeff_tmp,(6,ncoeff,nsym))
- + + - -
+ - + -
+ ]
2048 [ - + + + : 3593962 : list_symcoeff_tmp = list_symcoeff_tmp2
+ + + + ]
2049 : :
2050 : : !4.1 Count irreducible terms
2051 [ + + ]: 10721 : do icoeff = 1,ncoeff
2052 [ + + ]: 523937 : do isym = 1,nsym
2053 : 513216 : icoeff2 = list_symcoeff_tmp2(6,icoeff,isym)
2054 [ + + ]: 523908 : if (icoeff2> icoeff)then
2055 [ + + ]: 1706712 : list_symcoeff_tmp(:,icoeff2,1) = 0
2056 : : end if
2057 : : end do
2058 : : end do
2059 : :
2060 : 29 : ncoeff3 = 0
2061 [ + + ]: 10721 : do icoeff = 1,ncoeff
2062 [ + + + + : 73937 : if(.not.(all(list_symcoeff_tmp(:,icoeff,1)==0)))then
+ + ]
2063 : 156 : ncoeff3 = ncoeff3 + 1
2064 : : end if
2065 : : end do
2066 : :
2067 : : !4.2 Put irreducible terms in front of list_symcoeff_tmp3
2068 : : ! Store index of irreducible terms
2069 [ + + + - : 116 : ABI_MALLOC(list_symcoeff_tmp3,(6,ncoeff2,nsym))
- + + - +
- + + - +
- + ]
2070 [ + + + - : 87 : ABI_MALLOC(index_irred,(ncoeff3))
+ - + + -
+ - + ]
2071 : 29 : icoeff_tmp = 0
2072 [ + + ]: 10721 : do icoeff = 1,ncoeff
2073 [ + + + + : 73937 : if(.not.(all(list_symcoeff_tmp(:,icoeff,1)==0)))then
+ + ]
2074 : 156 : icoeff_tmp = icoeff_tmp+1
2075 : 156 : index_irred(icoeff_tmp) = icoeff
2076 [ + + + + ]: 52572 : list_symcoeff_tmp3(:,icoeff_tmp,:) = list_symcoeff_tmp(:,icoeff,:)
2077 : : endif
2078 : : enddo
2079 : :
2080 : : !4.3 Put symmetric equivalents behind in list_symcoeff_tmp3
2081 : : ! Attention icoeff_tmps keeps it's value of loop before
2082 : : ! TODO for check should be equal to ncoeff2
2083 [ + + ]: 10721 : do icoeff = 1,ncoeff
2084 : : if(.not.(all(list_symcoeff_tmp2(:,icoeff,1)==0))&
2085 [ + + + + : 186342 : & .and..not. any(index_irred == icoeff))then
+ + + + +
+ ]
2086 : 2216 : icoeff_tmp = icoeff_tmp+1
2087 [ + + + + ]: 746792 : list_symcoeff_tmp3(:,icoeff_tmp,:) = list_symcoeff_tmp2(:,icoeff,:)
2088 : : endif
2089 : : enddo
2090 : :
2091 : : !4.4 A little copy round
2092 [ + - ]: 29 : ABI_SFREE(list_symcoeff_tmp)
2093 [ + + + - : 87 : ABI_MALLOC(list_symcoeff_tmp,(6,ncoeff2,nsym))
- + + - -
+ - + -
+ ]
2094 [ - + + + : 798442 : list_symcoeff_tmp = list_symcoeff_tmp3
+ + + + ]
2095 : :
2096 : : !5/ Final transfert
2097 [ + + + - : 87 : ABI_MALLOC(list_symcoeff,(6,ncoeff2,nsym))
- + + - -
+ - + -
+ ]
2098 [ + + + + : 798413 : list_symcoeff = 0
+ + ]
2099 : : icoeff = 0
2100 [ + + ]: 2401 : do icoeff = 1,ncoeff2
2101 [ + + + + ]: 799364 : list_symcoeff(1:6,icoeff,:) = list_symcoeff_tmp(1:6,icoeff,:)
2102 [ + + ]: 116257 : do isym=1,nsym
2103 : : end do
2104 : : end do
2105 : :
2106 : : !6/ reset the dimension six of list_symcoeff_tmp2(6,icoeffs,1)
2107 : : ! and check is a symetric coeff is not coresspondig to an other
2108 : : ! one, in this case we set this coeff to 0
2109 [ + + ]: 2401 : do icoeff = 1,ncoeff2
2110 : : ! found the index of each coeff in list_fullcoeff
2111 [ + + ]: 116257 : do isym = 1,nsym
2112 : : icoeff2 = getCoeffFromList(list_symcoeff(:,:,1),&
2113 : : & list_symcoeff(2,icoeff,isym),&
2114 : : & list_symcoeff(3,icoeff,isym),&
2115 : : & list_symcoeff(4,icoeff,isym),&
2116 : : & list_symcoeff(1,icoeff,isym),&
2117 : 113856 : & ncoeff2)
2118 : 116228 : list_symcoeff(6,icoeff,isym) = icoeff2
2119 : : ! list_symcoeff(6,icoeff2,isym) = icoeff
2120 : : end do
2121 : : end do
2122 : :
2123 : : !Set the max number of coeff inside list_symcoeff
2124 : 29 : ncoeff_sym = ncoeff3
2125 : :
2126 : :
2127 : : !Deallocation
2128 [ + - ]: 29 : ABI_SFREE(blkval)
2129 [ + - ]: 29 : ABI_SFREE(list)
2130 [ + - ]: 29 : ABI_SFREE(list_symcoeff_tmp)
2131 [ + - ]: 29 : ABI_SFREE(list_symcoeff_tmp2)
2132 [ + - ]: 29 : ABI_SFREE(list_symcoeff_tmp3)
2133 [ + - ]: 29 : ABI_SFREE(list_symstr_tmp)
2134 [ + - ]: 29 : ABI_SFREE(index_irred)
2135 [ + - ]: 29 : ABI_SFREE(indsym)
2136 [ + - ]: 29 : ABI_SFREE(symrec)
2137 [ + - ]: 29 : ABI_SFREE(symrel)
2138 [ + - ]: 29 : ABI_SFREE(tnons)
2139 : : !ABI_FREE(xcart)
2140 : : !ABI_FREE(xred )
2141 [ + - ]: 29 : ABI_SFREE(wkdist)
2142 : 29 : end subroutine polynomial_coeff_getList
2143 : : !!***
2144 : :
2145 : :
2146 : : !!****f* m_polynomial_coeff/polynomial_coeff_getNorder
2147 : : !!
2148 : : !! NAME
2149 : : !! polynomial_coeff_getNorder
2150 : : !!
2151 : : !! FUNCTION
2152 : : !! Compute and store into the datatype coefficients, all the possible
2153 : : !! coefficients for given orders
2154 : : !!
2155 : : !! INPUTS
2156 : : !! cutoff = cut-off for the inter atomic forces constants
2157 : : !! crystal<type(crystal_t)> = datatype with all the information for the crystal
2158 : : !! power_disps(2) = array with the minimal and maximal power_disp to be computed
2159 : : !! max_power_strain = maximum order of the strain of the strain phonon coupling
2160 : : !! option = 0 compute all terms
2161 : : !! 1 still in development
2162 : : !! sc_size(3) = size of the supercell used for the fit.
2163 : : !! For example if you want to fit 2x2x2 cell the interation
2164 : : !! Sr-Ti and Sr-Ti[2 0 0] will be identical for the fit process
2165 : : !! If check_pbc is true we remove these kind of terms
2166 : : !! comm = MPI communicator
2167 : : !! anharmstr = logical, optional : TRUE, the anharmonic strain is computed (\eta)^power_disp ...
2168 : : !! FALSE, (default) the anharmonic strain are not computed
2169 : : !! spcoupling= logical, optional : TRUE(default) the anharmonic strain-phonon coupling is computed
2170 : : !! only_odd_power = logical, optional : if TRUE return only odd power
2171 : : !! only_even_power= logical, optional : if TRUe return only even power
2172 : : !! distributed = logical, optional : True, the coefficients will be distributed on the CPU
2173 : : !! verbose = optional, flag for the verbose mode
2174 : : !!
2175 : : !! OUTPUT
2176 : : !! polynomial_coeff<(type(polynomial_coeff_type)>(ncoeff) = array of datatype with the polynomial_coeff
2177 : : !! ncoeff = number of coefficients for this CPU if distributed == true, all otherwise
2178 : : !! ncoeff_tot = total number of coefficient over the CPU
2179 : : !!
2180 : : !! SOURCE
2181 : :
2182 : 44 : subroutine polynomial_coeff_getNorder(coefficients,crystal,cutoff,ncoeff,ncoeff_tot,power_disps,&
2183 : : & max_power_strain,option,sc_size,comm,anharmstr,spcoupling,&
2184 : : & distributed,only_odd_power,only_even_power,fit_iatom,&
2185 [ - + ]: 44 : & compute_symmetric,dispterms,verbose, max_nbody)
2186 : :
2187 : : !FIXME: There is a bug that the onebody term are duplicated.
2188 : : ! e.g. (Bx-O1x)^6 and (O1x-Bx)^6 both appears
2189 : : ! We can force ia<ib if onebody.
2190 : :
2191 : : !Arguments ------------------------------------
2192 : : !scalars
2193 : : integer,intent(in) :: max_power_strain,option,comm
2194 : : integer,intent(out):: ncoeff,ncoeff_tot
2195 : : integer,optional,intent(in) :: fit_iatom
2196 : : real(dp),intent(in):: cutoff
2197 : : logical,optional,intent(in) :: anharmstr,spcoupling,distributed,verbose,dispterms
2198 : : logical,optional,intent(in) :: only_odd_power,only_even_power,compute_symmetric
2199 : : !arrays
2200 : : integer,intent(in) :: power_disps(2),sc_size(3), max_nbody(:)
2201 : : type(crystal_t), intent(inout) :: crystal
2202 : : type(polynomial_coeff_type),allocatable,intent(inout) :: coefficients(:)
2203 : : !Local variables-------------------------------
2204 : : !scalar
2205 : : integer :: icoeff,icoeff2,icoeff3,ierr,ii,iterm
2206 : : integer :: i
2207 : : integer :: master,my_rank,my_ncoeff,my_newncoeff,natom,ncombination,ncoeff_max,ncoeff_sym
2208 : : integer :: ncoeff_symsym,nirred_comb,iirred_comb,ndisp,nstrain,fit_iatom_in
2209 : : integer :: ncoeff_alone,ndisp_max,nproc,nrpt,nsym,nterm,nstr_sym,my_size
2210 : : integer :: my_icoeff,rank_to_send,rank_to_receive,rank_to_send_save
2211 : : integer :: ncombi_alone,my_ncombi_simple,my_ncombi_start,my_ncombi_end,my_ncombi,my_nirred
2212 : : logical :: iam_master,need_anharmstr,need_spcoupling,need_distributed,need_verbose
2213 : : logical :: need_only_odd_power,need_only_even_power,compute_sym,irreducible,need_compute_symmetric
2214 : : logical :: need_dispterms
2215 : : !arrays
2216 : : integer :: shape_listsymcoeff(3),shape_listsymstr(3)
2217 : 44 : integer,allocatable :: buffsize(:),buffdispl(:) !,dummylist(:) ,index_irred(:)
2218 : 44 : integer,allocatable :: offsets(:)
2219 : 44 : integer,allocatable :: cell(:,:),compatibleCoeffs(:,:)
2220 : 44 : integer,allocatable :: list_symcoeff(:,:,:),list_symstr(:,:,:),list_coeff(:),list_combination(:,:)
2221 : 44 : integer,allocatable :: list_combination_tmp(:,:)
2222 : 44 : integer,allocatable :: irank_ncombi(:),my_index_irredcomb(:)
2223 : 44 : integer,allocatable :: my_coefflist(:),my_coeffindexes(:),my_newcoeffindexes(:),my_list_combination(:,:)
2224 : 44 : type(int2d_array_type) :: my_array_combination
2225 : 44 : integer,allocatable :: my_list_combination_tmp(:,:)
2226 : : real(dp) :: rprimd(3,3),range_ifc(3)
2227 : 44 : real(dp),allocatable :: dist(:,:,:,:)
2228 : 44 : character(len=5),allocatable :: symbols(:)
2229 : : character(len=200):: name
2230 : : character(len=500) :: message
2231 : 44 : type(polynomial_coeff_type),dimension(:),allocatable :: coeffs_tmp
2232 : 44 : type(polynomial_term_type),dimension(:),allocatable :: terms
2233 : : character(len=fnlen) :: filename
2234 : : ! *************************************************************************
2235 : : !Hide filename for debugging
2236 : : ABI_UNUSED(filename)
2237 : :
2238 : 44 : call get_mpi_params()
2239 : 44 : call free_output_coefficients()
2240 : 44 : call check_and_normalize_input_params()
2241 : 44 : call getList() ! of pairs
2242 : 44 : call get_combinations_of_lists()
2243 : 44 : call get_symmetric_combinations()
2244 : 44 : call combinations_to_terms()
2245 : 44 : call distribute_coefficients_over_cpu()
2246 [ - + - + : 44 : call free_memory()
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - -
- - - - -
- - - - -
- - - - -
+ - + - -
- - - - -
- - - - -
- - - - -
- - - - +
- + - + -
+ ]
2247 : : contains
2248 : 44 : subroutine get_mpi_params()
2249 : : !MPI variables
2250 : 44 : master = 0
2251 : 44 : nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2252 : 44 : iam_master = (my_rank == master)
2253 : :
2254 : :
2255 : 44 : end subroutine get_mpi_params
2256 : :
2257 : 44 : subroutine free_output_coefficients()
2258 : : !Free the output
2259 [ + + ]: 44 : if(allocated(coefficients))then
2260 [ + + ]: 11 : do ii =1,size(coefficients)
2261 : 11 : call polynomial_coeff_free(coefficients(ii))
2262 : : end do
2263 [ + - + + : 11 : ABI_SFREE(coefficients)
- + - - -
- - - - -
- - - - -
- - - -
+ ]
2264 : : end if
2265 : 44 : end subroutine free_output_coefficients
2266 : :
2267 : 44 : subroutine check_and_normalize_input_params()
2268 : : !Check
2269 [ - + ]: 44 : if(option > power_disps(2))then
2270 : : write(message, '(3a)' )&
2271 : 0 : & 'Option can not be superior to the maximum order ',ch10,&
2272 : 0 : & 'Action: contact abinit group'
2273 : 0 : ABI_ERROR(message)
2274 : : end if
2275 : :
2276 : : !Initialisation of variables
2277 : 44 : need_anharmstr = .TRUE.
2278 [ + - ]: 44 : if(present(anharmstr)) need_anharmstr = anharmstr
2279 : 44 : need_spcoupling = .TRUE.
2280 [ + - ]: 44 : if(present(spcoupling)) need_spcoupling = spcoupling
2281 : 44 : need_distributed = .FALSE.
2282 [ + + ]: 44 : if(present(distributed)) need_distributed = distributed
2283 : 44 : need_verbose = .TRUE.
2284 [ + + ]: 44 : if(present(verbose)) need_verbose = verbose
2285 : 44 : need_only_odd_power = .FALSE.
2286 [ + - ]: 44 : if(present(only_odd_power)) need_only_odd_power = only_odd_power
2287 : 44 : need_only_even_power = .FALSE.
2288 [ + - ]: 44 : if(present(only_even_power)) need_only_even_power = only_even_power
2289 : 44 : need_compute_symmetric = .TRUE.
2290 [ + + ]: 44 : if(present(compute_symmetric)) need_compute_symmetric = compute_symmetric
2291 : 44 : need_dispterms = .TRUE.
2292 [ + + ]: 44 : if(present(dispterms)) need_dispterms = dispterms
2293 : :
2294 [ - + - - ]: 44 : if(need_only_odd_power.and.need_only_even_power)then
2295 : : write(message, '(3a)' )&
2296 : 0 : & 'need_only_odd_power and need_only_even_power are both true',ch10,&
2297 : 0 : & 'Action: contact abinit group'
2298 : 0 : ABI_ERROR(message)
2299 : : end if
2300 : :
2301 : : !Check which atom to fit, if not present do all atoms
2302 [ + + ]: 44 : if(present(fit_iatom))then
2303 : : ! FIXME: the fit_iatom as the input is the index of the atom
2304 : : ! but in the fitting subroutine, it is the irreducible atom.
2305 : 43 : fit_iatom_in = fit_iatom
2306 : : else
2307 : 1 : fit_iatom_in = -1
2308 : : endif
2309 : :
2310 : 44 : end subroutine check_and_normalize_input_params
2311 : :
2312 : 44 : subroutine getList()
2313 [ + + + + ]: 572 : rprimd = crystal%rprimd
2314 : 44 : call prepare_for_getList(crystal,sc_size, dist, cell, natom, nsym, nrpt, range_ifc , symbols)
2315 : :
2316 [ + + ]: 44 : if(iam_master)then
2317 [ + + ]: 28 : if(need_verbose)then
2318 : 22 : write(message,'(1a)')' Generation of the list of all the possible pairs of atoms within cutoff'
2319 : 22 : call wrtout(std_out,message,'COLL')
2320 : : end if
2321 : : call polynomial_coeff_getList(cell,crystal,dist,list_symcoeff,list_symstr,&
2322 : : & natom,nstr_sym,ncoeff_sym,nrpt,range_ifc,cutoff,sc_size=sc_size,&
2323 [ - + - + ]: 28 : & fit_iatom=fit_iatom_in)
2324 [ + + ]: 112 : shape_listsymcoeff = shape(list_symcoeff)
2325 [ + + ]: 112 : shape_listsymstr = shape(list_symstr)
2326 : : endif!if iam master
2327 : :
2328 : : !Broadcast Results of getList
2329 : 44 : call xmpi_bcast(shape_listsymcoeff, master, comm, ierr)
2330 : 44 : call xmpi_bcast(shape_listsymstr, master, comm, ierr)
2331 : : call xmpi_bcast(nstr_sym, master, comm, ierr)
2332 : : call xmpi_bcast(ncoeff_sym, master, comm, ierr)
2333 [ + + ]: 44 : if(.not. iam_master )then
2334 [ - + + - : 80 : ABI_MALLOC(list_symcoeff,(shape_listsymcoeff(1),shape_listsymcoeff(2),shape_listsymcoeff(3)))
- + + - -
+ + - + -
+ - - + -
+ - + ]
2335 [ - + + - : 80 : ABI_MALLOC(list_symstr,(shape_listsymstr(1),shape_listsymstr(2),shape_listsymstr(3)))
- + + - -
+ + - + -
+ - - + -
+ - + ]
2336 : : endif
2337 : : call xmpi_bcast(list_symcoeff, master, comm, ierr)
2338 : : call xmpi_bcast(list_symstr, master, comm, ierr)
2339 : 44 : call xmpi_barrier(comm)
2340 : :
2341 : 44 : ncoeff_symsym = size(list_symcoeff(1,:,1))
2342 : :
2343 : : !Compute the total number of coefficient
2344 : 44 : ncoeff_tot = ncoeff_sym+nstr_sym
2345 : :
2346 : :
2347 : : !if(iam_master)then
2348 : : !Check the distanceance bewteen coefficients and store integer:
2349 : : ! 0: the mix between these coefficient is not possible
2350 : : ! 1: the mix between these coefficient is possible
2351 [ - + + - : 176 : ABI_MALLOC(compatibleCoeffs,(ncoeff_symsym+nstr_sym,ncoeff_symsym+nstr_sym))
- + + - +
- + - - +
- + - + ]
2352 [ + + + + ]: 280808 : compatibleCoeffs(:,:) = 1
2353 : :
2354 [ + + ]: 44 : if(need_verbose)then
2355 : 38 : write(message,'(1a)')' Check the compatible coefficients with respect to the cutoff'
2356 : 38 : call wrtout(std_out,message,'COLL')
2357 : : end if
2358 : :
2359 [ + + ]: 2318 : do icoeff=1,ncoeff_symsym+nstr_sym
2360 [ + + ]: 280808 : do icoeff2=1,ncoeff_symsym+nstr_sym
2361 : : ! Select case:
2362 : : ! if both icoeff are displacement => check the distance
2363 : : ! if both icoeff are strain => check the flag
2364 : : ! Otherwise cycle (we keep the term)
2365 : : ! if(icoeff>ncoeff_sym.and.icoeff2<=ncoeff_sym)cycle
2366 : : ! if(icoeff<=ncoeff_sym.and.icoeff2>ncoeff_sym)cycle
2367 : : !Forbid anharmonic strain terms if not wanted
2368 [ + + + + : 278490 : if((icoeff>ncoeff_symsym .and. icoeff2>ncoeff_symsym).and.&
+ + ]
2369 : : & .not.need_anharmstr) then
2370 : 1440 : compatibleCoeffs(icoeff,icoeff2) = 0
2371 : 1440 : compatibleCoeffs(icoeff2,icoeff) = 0
2372 : : end if
2373 : : !Forbid strain phonon terms if not wanted
2374 [ + + + + : 278490 : if((icoeff>ncoeff_symsym.or.icoeff2>ncoeff_symsym).and.&
+ + ]
2375 : : & .not.need_spcoupling) then
2376 : 12888 : compatibleCoeffs(icoeff,icoeff2) = 0
2377 : 12888 : compatibleCoeffs(icoeff2,icoeff) = 0
2378 : : end if
2379 : :
2380 : :
2381 [ + + + + ]: 280764 : if(icoeff<=ncoeff_symsym.and.icoeff2<=ncoeff_symsym)then !Check combination of irreducible bodies and their symmetric equivalent
2382 [ - + ]: 252786 : if(list_symcoeff(2,icoeff,1)/= list_symcoeff(2,icoeff2,1)) then
2383 : 0 : ABI_BUG("The first components of the pairs in the coefficient are not equivalent.")
2384 : : end if
2385 : :
2386 : : if(abs(dist(1,list_symcoeff(2,icoeff,1),list_symcoeff(3,icoeff,1),list_symcoeff(4,icoeff,1)) & ! rx(a, b)
2387 : : & -dist(1,list_symcoeff(2,icoeff,1),list_symcoeff(3,icoeff2,1),list_symcoeff(4,icoeff2,1))) & ! rx(a, d )
2388 : : & >= (rprimd(1,1) + rprimd(1,2) + rprimd(1,3))*sc_size(1) .or. &
2389 : : abs(dist(2,list_symcoeff(2,icoeff,1),list_symcoeff(3,icoeff,1),list_symcoeff(4,icoeff,1)) & ! ry(a, b)
2390 : : & -dist(2,list_symcoeff(2,icoeff,1),list_symcoeff(3,icoeff2,1),list_symcoeff(4,icoeff2,1))) & ! ry(a, d)
2391 [ + + + + : 252786 : & >= (rprimd(2,1) + rprimd(2,2) + rprimd(2,3))*sc_size(2) .or. &
+ + ]
2392 : : abs(dist(3,list_symcoeff(2,icoeff,1),list_symcoeff(3,icoeff,1),list_symcoeff(4,icoeff,1)) & ! rz(a, b)
2393 : : & -dist(3,list_symcoeff(2,icoeff,1),list_symcoeff(3,icoeff2,1),list_symcoeff(4,icoeff2,1))) & ! rz(a, d)
2394 : : & >= (rprimd(3,1) + rprimd(3,2) + rprimd(3,3))*sc_size(3))then
2395 : 63216 : compatibleCoeffs(icoeff,icoeff2) = 0
2396 : 63216 : compatibleCoeffs(icoeff2,icoeff) = 0
2397 : : endif
2398 : : endif
2399 : : end do !end icoeff
2400 : : end do !icoeff2
2401 : :
2402 [ + - ]: 44 : ABI_SFREE(dist)
2403 : : ! Compute all the combination of coefficient up to the given order (get the number)
2404 [ + + ]: 44 : if(need_verbose)then
2405 : 38 : write(message,'(1a)')' Compute the number of possible combinations'
2406 : 38 : call wrtout(std_out,message,'COLL')
2407 : : end if
2408 : :
2409 : 88 : end subroutine getList
2410 : :
2411 : 44 : subroutine get_combinations_of_lists()
2412 [ - + - + ]: 44 : ABI_MALLOC(list_coeff,(0))
2413 [ - + - + ]: 44 : ABI_MALLOC(list_combination,(0,0))
2414 : 44 : icoeff = 1
2415 : 44 : icoeff2 = 0
2416 : 44 : nirred_comb = 0
2417 : 44 : iirred_comb = 0
2418 : : call computeCombinationFromList(cell,compatibleCoeffs,list_symcoeff,list_symstr,&
2419 : : & list_coeff,list_combination,icoeff,max_power_strain,natom,ncoeff_sym,&
2420 : : & ncoeff_symsym,iirred_comb,nirred_comb,nstr_sym,icoeff,nrpt,nsym,1,power_disps(1),power_disps(2),symbols,comm,&
2421 : : & nbody=option,compute=.false.,anharmstr=need_anharmstr,spcoupling=need_spcoupling,&
2422 : 44 : & only_odd_power=need_only_odd_power,only_even_power=need_only_even_power,disp=need_dispterms)
2423 [ + - ]: 44 : ABI_SFREE(list_coeff)
2424 [ + - ]: 44 : ABI_SFREE(list_combination)
2425 : : ! Output how much we found
2426 [ + + ]: 44 : if(need_verbose)then
2427 : 38 : write(message,"(1a,I10)") " -Number of irreducible pairs within cutoff: ", ncoeff_sym
2428 : 38 : call wrtout(std_out,message,'COLL')
2429 : 38 : write(message,"(1a,I10)") " -Number of combinations of irreducible pairs: ", nirred_comb
2430 : 38 : call wrtout(std_out,message,'COLL')
2431 : : !write(message,"(1a,I10)") " -Number of possible symmetric combinations: ", ncombination
2432 : : !call wrtout(std_out,message,'COLL')
2433 : 38 : write(message,'(a,a)') ch10,' Compute the combinations of irreducible pairs'
2434 : 38 : call wrtout(std_out,message,'COLL')
2435 : : end if
2436 [ - + - + ]: 44 : ABI_MALLOC(list_coeff,(0))
2437 [ - + + - : 176 : ABI_MALLOC(list_combination_tmp,(power_disps(2),nirred_comb))
+ + + - +
- + + - +
- + - + ]
2438 : 44 : icoeff = 1
2439 : 44 : icoeff2 = 0
2440 : 44 : iirred_comb = 0
2441 [ + + + + ]: 13689 : list_combination_tmp = 0
2442 : : ! Compute all the combination of coefficient up to the given order
2443 : : call computeCombinationFromList(cell,compatibleCoeffs,list_symcoeff,list_symstr,&
2444 : : & list_coeff,list_combination_tmp,icoeff,max_power_strain,natom,&
2445 : : & ncoeff_sym,ncoeff_symsym,iirred_comb,nirred_comb,nstr_sym,ncombination,nrpt,nsym,1,power_disps(1),&
2446 : : & power_disps(2),symbols,comm,nbody=option,compute=.true.,&
2447 : : & anharmstr=need_anharmstr,spcoupling=need_spcoupling,&
2448 : 44 : & only_odd_power=need_only_odd_power,only_even_power=need_only_even_power,disp=need_dispterms)
2449 [ + - ]: 44 : ABI_SFREE(list_coeff)
2450 : 44 : nirred_comb = size(list_combination_tmp,2)
2451 : :
2452 : 44 : end subroutine get_combinations_of_lists
2453 : :
2454 : :
2455 : 44 : subroutine get_symmetric_combinations()
2456 : 44 : type(IrreducibleCombinations_T) :: irred_combinations
2457 : : ! If we want to compute equivalent symmetric combinations go here.
2458 [ + + ]: 44 : if(need_compute_symmetric)then
2459 [ + - ]: 38 : if(need_verbose)then
2460 : 38 : write(message,'(1a)')' Distribute irreducible combinations over CPU'
2461 : 38 : call wrtout(std_out,message,'COLL')
2462 : 38 : write(message,'(1a)')' Compute symmetric combinations of combinations of irreducible pairs'
2463 : 38 : call wrtout(std_out,message,'COLL')
2464 : 38 : write(message,'(3a)')' ---> Try to match number of CPU to number combinations of irreducible pairs',ch10,&
2465 : 76 : & ' for max. speedup'
2466 : 38 : call wrtout(std_out,message,'COLL')
2467 : : endif
2468 : :
2469 : : ! COUNT IRREDUCIBLE COMBINATIONS FOR EACH PROCESSOR
2470 : 38 : ncombi_alone = mod(nirred_comb,nproc)
2471 : 38 : my_ncombi_simple = int(aint(real(nirred_comb,sp)/(nproc)))
2472 [ + + + - ]: 38 : if(ncombi_alone == 0 .and. nirred_comb >= nproc)then !ncombi > nproc and no remainder
2473 : 22 : my_ncombi_start = (my_ncombi_simple * my_rank) + 1
2474 : 22 : my_ncombi_end = my_ncombi_start + my_ncombi_simple - 1
2475 [ - + ]: 16 : else if(nirred_comb < nproc)then !ncombi smaller than nproc
2476 [ # # ]: 0 : if(my_rank + 1 <= nirred_comb)then !myrank smaller than ncombi
2477 : 0 : my_ncombi_start = my_rank + 1
2478 : 0 : my_ncombi_end = my_ncombi_start
2479 : : else
2480 : 0 : my_ncombi_start = nirred_comb + 1 !myrank bigger than ncombi
2481 : 0 : my_ncombi_end = nirred_comb + 1
2482 : : endif
2483 [ + - + - ]: 16 : else if(nirred_comb > nproc .and. ncombi_alone /= 0)then !ncombi > nproc and remainder
2484 [ + + ]: 16 : if(my_rank >= (nproc-ncombi_alone)) then
2485 : 8 : my_ncombi_start = (my_ncombi_simple * my_rank) + 1 + (my_rank - nproc + ncombi_alone)
2486 : 8 : my_ncombi_end = my_ncombi_start + my_ncombi_simple
2487 : : else
2488 : 8 : my_ncombi_start = (my_ncombi_simple * my_rank) + 1
2489 : 8 : my_ncombi_end = my_ncombi_start + my_ncombi_simple - 1
2490 : : endif
2491 : : end if
2492 : :
2493 [ + - ]: 38 : if(my_ncombi_end <= nirred_comb)then
2494 : 38 : my_nirred = my_ncombi_end+1-my_ncombi_start
2495 : : else
2496 : 0 : my_nirred = 0
2497 : : endif
2498 : : !write(std_out,*) "my_rank", my_rank, "my_ncmobi_start", my_ncombi_start, "my_ncombi_end", my_ncombi_end,'my_nirred',my_nirred
2499 : :
2500 : : !COPY IRREDUCIBLE COMBINTATIONS TO BE DONE TO EACH PROCESSOR
2501 : : !my_nirred = my_ncombi_end - my_ncombi_start + 1
2502 [ - + + - : 152 : ABI_MALLOC(my_list_combination_tmp,(power_disps(2),my_nirred))
- + + - +
- + - - +
- + - + ]
2503 [ + - + + : 9173 : if(my_nirred /= 0) my_list_combination_tmp(:,:) = list_combination_tmp(:,my_ncombi_start:my_ncombi_end)
+ + ]
2504 : :
2505 : :
2506 [ - + + - : 114 : ABI_MALLOC(my_index_irredcomb,(my_nirred))
+ - + - -
+ - + -
+ ]
2507 [ + - ]: 38 : ABI_SFREE(list_combination_tmp)
2508 : :
2509 : : !COUNT SYMMETRIC COMBINATIONS TO IRREDUCIBLE COMBINATIONS ON EACH PROCESSOR
2510 : :
2511 : : !COMPUTE SYMMETRIC COMBINATIONS
2512 : :
2513 : 38 : call irred_combinations%init()
2514 [ + + ]: 2309 : do i=1,my_nirred
2515 : 38 : associate(comb => my_list_combination_tmp(:, i))
2516 : : !ABI_MALLOC(dummylist,(0))
2517 : : ! Get number of strain and displacements for this term
2518 : 2271 : ndisp = 0
2519 : 2271 : nstrain = 0
2520 [ + + ]: 9135 : do ii = 1,power_disps(2)
2521 [ + + + + ]: 9135 : if(comb(ii) > 0 .and. comb(ii) <= ncoeff_symsym)then
2522 : 4010 : ndisp = ndisp + 1
2523 [ + + ]: 2854 : else if(comb(ii) >= ncoeff_symsym)then
2524 : 2832 : nstrain = nstrain + 1
2525 : : endif
2526 : : enddo
2527 : 2271 : compute_sym = .true.
2528 : 2271 : iterm = my_index_irredcomb(i)-1
2529 : :
2530 : : call computeSymmetricCombinations(my_array_combination,list_symcoeff,list_symstr,ndisp,nsym,&
2531 : : & comb(:ndisp+nstrain),power_disps(2),&
2532 : : & ncoeff_symsym,nstr_sym,nstrain, &
2533 : : & compatibleCoeffs,compute_sym,comm, &
2534 : : & only_even=need_only_even_power, max_nbody=max_nbody, &
2535 : 4542 : & irred_combinations=irred_combinations, cell=cell)
2536 : : end associate
2537 : : enddo
2538 : 38 : call irred_combinations%free()
2539 : :
2540 [ + - ]: 38 : ABI_SFREE(my_list_combination_tmp)
2541 : :
2542 : : ! Delete double combinations on each processor
2543 [ + - + - ]: 38 : if(need_verbose .and. my_nirred /= 0)then
2544 : 38 : write(message,'(1a,I4)')' Reduce reducible symmetric combinations on processor: ', my_rank+1
2545 : 38 : call wrtout(std_out,message,'PERS')
2546 : : endif
2547 : : !call reduce_zero_combinations(my_list_combination)
2548 : 38 : call my_array_combination%tostatic(my_list_combination, size1=power_disps(2))
2549 : 38 : call my_array_combination%finalize()
2550 : :
2551 : :
2552 : : ! Gather the Results into list_combination_tmp
2553 : 38 : my_ncombi = size(my_list_combination,2)
2554 [ - + + - : 114 : ABI_MALLOC(irank_ncombi,(nproc))
+ - + - -
+ - + -
+ ]
2555 : : call xmpi_allgather(my_ncombi,irank_ncombi,comm,ierr)
2556 [ + - ]: 38 : if(need_verbose)then
2557 : 38 : write(message,'(1a)')' Reduction on all processors finished. Gather results.'
2558 : 38 : call wrtout(std_out,message,'COLL')
2559 : : endif
2560 [ - + + - : 246 : ABI_MALLOC(list_combination_tmp,(power_disps(2),sum(irank_ncombi)))
+ + - + +
- + - + -
- + - + -
+ ]
2561 : :
2562 : :
2563 [ - + + - : 114 : ABI_MALLOC(offsets,(nproc))
+ - + - -
+ - + -
+ ]
2564 : 38 : offsets(1) = 0
2565 [ + + ]: 132 : do i=1,nproc
2566 [ + + ]: 236 : offsets(i) = sum(irank_ncombi(:i-1))*power_disps(2)
2567 : : enddo
2568 : :
2569 [ + + + + ]: 51284 : list_combination_tmp = 0
2570 [ - + + - : 114 : ABI_MALLOC(buffsize,(nproc))
+ - + - -
+ - + -
+ ]
2571 [ + + ]: 132 : do i = 1,nproc
2572 : 132 : buffsize(i) = irank_ncombi(i)*power_disps(2)
2573 : : enddo
2574 : :
2575 : :
2576 [ + + ]: 114 : call xmpi_gatherv(my_list_combination,size(my_list_combination),list_combination_tmp,buffsize,offsets,master,comm,ierr)
2577 : :
2578 : :
2579 : : !Deallocation of variables inside need_symmetric
2580 [ + - ]: 38 : ABI_SFREE(buffsize)
2581 [ + - ]: 38 : ABI_SFREE(my_list_combination)
2582 [ + - ]: 38 : ABI_SFREE(my_index_irredcomb)
2583 [ + - ]: 38 : ABI_SFREE(irank_ncombi)
2584 [ + - ]: 76 : ABI_SFREE(offsets)
2585 : : endif !compute_symmetric
2586 : :
2587 : : !Deallocation of arrays outside need_symmetric
2588 [ + - ]: 44 : ABI_SFREE(compatibleCoeffs)
2589 : :
2590 [ + + ]: 44 : if(iam_master)then
2591 : 28 : call reduce_zero_combinations(list_combination_tmp)
2592 : 28 : ncombination = size(list_combination_tmp,2)
2593 : : !ABI_MALLOC(index_irred,(1))
2594 : : !index_irred = 1
2595 : : !============ Strain phonon coupling===================
2596 [ + + ]: 28 : if(need_spcoupling)then !Check irreducibility of strain-phonon terms
2597 [ + - ]: 22 : if(need_verbose)then
2598 : 22 : write(message,'(1a)')' Reduce reducible Strain-Phonon combinations on master'
2599 : 22 : call wrtout(std_out,message,'COLL')
2600 : : endif
2601 : 22 : call irred_combinations%init()
2602 [ + + ]: 7278 : do i=1, ncombination
2603 [ + + + + : 26268 : if(any(list_combination_tmp(:,i) > ncoeff_symsym))then
+ + ]
2604 : : irreducible=irred_combinations%add_irr(list_combination_tmp(:,i), &
2605 : 2990 : & list_symcoeff, list_symstr, ncoeff_symsym, nsym, power_disps(2),cell)
2606 [ + + + + ]: 3302 : if(.not. irreducible) list_combination_tmp(:,i) = 0
2607 : : endif
2608 : : end do
2609 : 22 : call reduce_zero_combinations(list_combination_tmp)
2610 : 22 : ncombination = size(list_combination_tmp,2)
2611 : 22 : call irred_combinations%free()
2612 : : endif
2613 : :
2614 : : end if !iam_master
2615 [ + + ]: 44 : if(need_verbose)then
2616 : 38 : write(message,'(1x,I0,1a)') ncombination,' irreducible combinations generated '
2617 : 38 : call wrtout(std_out,message,'COLL')
2618 : 38 : write(message,'(1a)') ' Finished generating irreducible combinations'
2619 : 38 : call wrtout(std_out,message,'COLL')
2620 : : endif
2621 : :
2622 : : !MPI
2623 [ + + + + ]: 44 : if(need_verbose .and. nproc > 1)then
2624 : 24 : write(message,'(a,a)') ch10,' Redistribute the combinations over the CPU'
2625 : 24 : call wrtout(std_out,message,'COLL')
2626 : : end if
2627 : :
2628 : : call xmpi_bcast(ncombination, master, comm, ierr)
2629 : :
2630 : 44 : ncoeff_alone = mod(ncombination,nproc)
2631 : 44 : my_ncoeff = int(aint(real(ncombination,sp)/(nproc)))
2632 : :
2633 [ + + ]: 44 : if(my_rank >= (nproc-ncoeff_alone)) then
2634 : 6 : my_ncoeff = my_ncoeff + 1
2635 : : end if
2636 : :
2637 : : !Set the buffsize for mpi scatterv
2638 [ - + + - : 132 : ABI_MALLOC(buffsize,(nproc))
+ - + - -
+ - + -
+ ]
2639 [ - + + - : 88 : ABI_MALLOC(buffdispl,(nproc))
- + - + -
+ ]
2640 [ + + ]: 144 : do ii = 1,nproc
2641 : 100 : buffsize(ii) = int(aint(real(ncombination,sp)/(nproc))*power_disps(2))
2642 [ + + ]: 144 : if(ii > (nproc-ncoeff_alone)) then
2643 : 20 : buffsize(ii) = buffsize(ii) + power_disps(2)
2644 : : end if
2645 : : end do
2646 : :
2647 : 44 : buffdispl(1) = 0
2648 [ + + ]: 100 : do ii = 2,nproc
2649 : 100 : buffdispl(ii) = buffdispl(ii-1) + buffsize(ii-1)
2650 : : end do
2651 : :
2652 [ - + + - : 176 : ABI_MALLOC(list_combination,(power_disps(2),my_ncoeff))
+ + + - +
- + + - +
- + - + ]
2653 [ + + + + ]: 29688 : list_combination = 0
2654 : :
2655 : 44 : my_size = my_ncoeff*power_disps(2)
2656 : : call xmpi_scatterv(list_combination_tmp,buffsize,buffdispl,list_combination,my_size,master,&
2657 : : & comm,ierr)
2658 : :
2659 [ + - ]: 44 : ABI_SFREE(buffdispl)
2660 [ + - ]: 44 : ABI_SFREE(buffsize)
2661 [ + - ]: 44 : ABI_SFREE(list_combination_tmp)
2662 : :
2663 [ - + - - : 132 : end subroutine get_symmetric_combinations
- - - + ]
2664 : :
2665 : 44 : subroutine combinations_to_terms()
2666 : 44 : logical, allocatable :: reverse(:)
2667 [ + + + + ]: 44 : if(need_verbose .and. nproc > 1)then
2668 : 24 : write(message,'(1a)')' Compute the coefficients'
2669 : 24 : call wrtout(std_out,message,'COLL')
2670 : : end if
2671 [ + + + - : 7390 : ABI_MALLOC(coeffs_tmp,(my_ncoeff))
+ - + + -
+ - + - +
+ + ]
2672 : 44 : nterm = nsym
2673 : 44 : ndisp_max = power_disps(2)
2674 : 44 : ncoeff_max = my_ncoeff
2675 [ + + ]: 7302 : do ii=1,my_ncoeff
2676 [ - + + - : 370158 : ABI_MALLOC(terms,(nterm))
+ - + - -
+ - + - +
+ + ]
2677 [ - + + - : 21774 : ABI_MALLOC(reverse,(ndisp_max))
+ - + - -
+ - + -
+ ]
2678 [ + + ]: 29644 : reverse=.False.
2679 : : call generateTermsFromList(cell,list_combination(:,ii),list_symcoeff,list_symstr,ncoeff_symsym,&
2680 : 7258 : & ndisp_max,nrpt,nstr_sym,nsym,nterm,terms, reverse=reverse)
2681 : :
2682 : :
2683 : 7258 : call polynomial_coeff_init(one,nterm,coeffs_tmp(ii),terms(1:nterm), check=.true.)
2684 : : !DMSG(coeffs_tmp(ii)%debug_str)
2685 [ + - ]: 7258 : ABI_SFREE(reverse)
2686 : : ! Free the terms array
2687 [ + + ]: 355642 : do iterm=1,nterm
2688 : 355642 : call polynomial_term_free(terms(iterm))
2689 : : end do
2690 [ + - + + : 355686 : ABI_SFREE(terms)
- + - + -
+ - + - +
- + - + -
+ ]
2691 : : end do
2692 [ + - ]: 44 : ABI_SFREE(cell)
2693 [ + - ]: 44 : ABI_SFREE(list_combination)
2694 [ + - ]: 44 : ABI_SFREE(list_symcoeff)
2695 [ + - ]: 44 : ABI_SFREE(list_symstr)
2696 : :
2697 : : !Final tranfert
2698 : : !1- Count the total number of coefficient
2699 : 44 : ncoeff = 0
2700 [ + + ]: 7302 : do icoeff=1,ncoeff_max
2701 [ + + ]: 7302 : if (abs(coeffs_tmp(icoeff)%coefficient) >tol16) then
2702 : 5063 : ncoeff = ncoeff + 1
2703 : : end if
2704 : : end do
2705 : :
2706 : : !Get the total number of coefficients
2707 : : !ncoeff_max is the number of total coefficients before the symetries check
2708 : : !ncoeff_tot is the number of total coefficients after the symetries check
2709 : 44 : ncoeff_tot = ncoeff!set the output
2710 : : call xmpi_sum(ncoeff_tot,comm,ierr)
2711 : : call xmpi_sum(ncoeff_max,comm,ierr)
2712 : 44 : end subroutine combinations_to_terms
2713 : :
2714 : :
2715 : 44 : subroutine distribute_coefficients_over_cpu()
2716 : : !Need to redistribute the coefficients over the CPU
2717 : : !Get the list with the number of coeff on each CPU
2718 : : !In order to be abble to compute the my_coeffindexes array which is for example:
2719 : : ! if CPU0 has 200 Coeff and CPU1 has 203 Coeff then
2720 : : ! for CPU0:my_coeffindexes=>1-200 and for CPU1:my_coeffindexes=>201-403
2721 [ + + + + ]: 44 : if(need_verbose .and. nproc > 1)then
2722 : 24 : write(message,'(1a)')' Redistribute the coefficients over the CPU'
2723 : 24 : call wrtout(std_out,message,'COLL')
2724 : : end if
2725 : :
2726 [ - + + - : 132 : ABI_MALLOC(buffdispl,(nproc))
+ - + - -
+ - + -
+ ]
2727 [ + + ]: 144 : buffdispl = 0
2728 : 44 : buffdispl(my_rank+1) = my_ncoeff
2729 : : call xmpi_sum(buffdispl,comm,ierr)
2730 [ + + + - : 132 : ABI_MALLOC(my_coeffindexes,(my_ncoeff))
+ - + + -
+ - + -
+ ]
2731 [ + + + - : 88 : ABI_MALLOC(my_coefflist,(my_ncoeff))
- + - + -
+ ]
2732 [ + + ]: 7302 : my_coeffindexes = 0
2733 [ + + ]: 7302 : my_coefflist = 0
2734 [ + + ]: 7302 : do icoeff=1,my_ncoeff
2735 : 7258 : my_coefflist(icoeff) = icoeff
2736 [ + + ]: 7302 : if(my_rank==0) then
2737 : 5621 : my_coeffindexes(icoeff) = icoeff
2738 : : else
2739 [ + + ]: 4257 : my_coeffindexes(icoeff) = sum(buffdispl(1:my_rank)) + icoeff
2740 : : end if
2741 : : end do
2742 [ + - ]: 44 : ABI_SFREE(buffdispl)
2743 : :
2744 : : !Compute the new number of coefficient per CPU
2745 [ + + ]: 44 : if(need_distributed) then
2746 : 38 : ncoeff_alone = mod(ncoeff_tot,nproc)
2747 : 38 : my_newncoeff = int(aint(real(ncoeff_tot,sp)/(nproc)))
2748 [ - + ]: 38 : if(my_rank >= (nproc-ncoeff_alone)) then
2749 : 0 : my_newncoeff = my_newncoeff + 1
2750 : : end if
2751 : : else
2752 : 6 : my_newncoeff = ncoeff_tot
2753 : : end if
2754 : :
2755 : 44 : ncoeff = my_newncoeff ! Set the output
2756 : :
2757 : : !2:compute the number of coefficients and the list of the corresponding
2758 : : ! coefficients for each CPU.
2759 [ + + + - : 132 : ABI_MALLOC(my_newcoeffindexes,(my_newncoeff))
+ - + + -
+ - + -
+ ]
2760 [ + + ]: 44 : if(need_distributed) then
2761 [ + + ]: 4995 : do icoeff=1,my_newncoeff
2762 [ - + ]: 4995 : if(my_rank >= (nproc-ncoeff_alone))then
2763 : : my_newcoeffindexes(icoeff)=int(aint(real(ncoeff_tot,sp)/(nproc)))*(my_rank)+&
2764 : 0 : & (my_rank - (nproc-ncoeff_alone)) + icoeff
2765 : : else
2766 : 4957 : my_newcoeffindexes(icoeff)=(my_newncoeff)*(my_rank) + icoeff
2767 : : end if
2768 : : end do
2769 : : else
2770 [ + + ]: 112 : do icoeff=1,my_newncoeff
2771 : 112 : my_newcoeffindexes(icoeff) = icoeff
2772 : : end do
2773 : : end if
2774 : :
2775 : : !2- Transfer
2776 [ + + ]: 44 : if(.not.need_distributed)then
2777 [ + - ]: 6 : if(.not.allocated(coefficients))then
2778 [ + + + - : 124 : ABI_MALLOC(coefficients,(my_newncoeff))
+ - + + -
+ - + - +
+ + ]
2779 : : end if
2780 : : end if
2781 : : icoeff = 0! icoeff is the current index in the total list of coefficients
2782 : 44 : icoeff2 = 0! icoeff2 is the current index in the output coefficients array on each CPU
2783 : 44 : icoeff3 = 0! icoeff3 is the current index in total new list of coefficients
2784 : 44 : rank_to_send_save = 0
2785 : :
2786 [ + + ]: 12534 : do icoeff=1,ncoeff_max
2787 : : ! Need to send the rank with the chosen coefficient
2788 : 12490 : rank_to_send = 0
2789 : 12490 : my_icoeff = 0
2790 [ + + ]: 2331761 : do ii=1,my_ncoeff
2791 [ + + ]: 2331761 : if (my_coeffindexes(ii)==icoeff) then
2792 : 7258 : my_icoeff = ii
2793 [ + + ]: 7258 : if (abs(coeffs_tmp(my_icoeff)%coefficient) > tol16)then
2794 : 5063 : rank_to_send = my_rank
2795 : : else
2796 : 2195 : rank_to_send = -1
2797 : : ! Free the coefficient
2798 : 2195 : call polynomial_coeff_free(coeffs_tmp(ii))
2799 : : end if
2800 : : exit
2801 : : end if
2802 : : end do
2803 : : call xmpi_sum(rank_to_send, comm, ierr)
2804 : : ! This coefficient is not compute
2805 [ + + ]: 12490 : if (rank_to_send == -1) cycle
2806 : :
2807 : : ! increase icoeff3
2808 : 8759 : icoeff3 = icoeff3 + 1
2809 : :
2810 : : ! Find the receiver CPU
2811 : 8759 : rank_to_receive = 0
2812 [ + + ]: 2032208 : do ii=1,my_newncoeff
2813 [ + + ]: 2032208 : if (my_newcoeffindexes(ii)==icoeff3) then
2814 : 5063 : rank_to_receive = my_rank
2815 : : end if
2816 : : end do
2817 : : call xmpi_sum(rank_to_receive, comm, ierr)
2818 : :
2819 [ + + + + ]: 8759 : if(need_distributed.and.rank_to_send /= rank_to_send_save) then
2820 : : if(my_rank == rank_to_send_save)then
2821 : : !call polynomial_coeff_list_free(coeffs_tmp)
2822 : : !ABI_SFREE(coeffs_tmp)!Free memory if the current CPU has already distribute
2823 : : !all its own coefficients
2824 : : end if
2825 : 56 : rank_to_send_save = rank_to_send
2826 : : end if
2827 : :
2828 [ + + + + ]: 8759 : if(need_distributed.and.my_rank == rank_to_receive)then
2829 [ + + ]: 4957 : if(.not.allocated(coefficients))then
2830 [ - + + - : 5071 : ABI_MALLOC(coefficients,(my_newncoeff))
+ - + - -
+ - + - +
+ + ]
2831 : : end if
2832 : : end if
2833 : :
2834 : :
2835 [ + + ]: 30052 : if (need_distributed)then
2836 [ + + ]: 8653 : if(my_rank==rank_to_send)then
2837 [ + + + + : 774525 : if(any(my_newcoeffindexes(:)==icoeff3))then
+ + ]
2838 : 4869 : icoeff2 = icoeff2 + 1
2839 : : ! Get the name of this coefficient
2840 : 4869 : call polynomial_coeff_getName(name,coeffs_tmp(my_icoeff),symbols,recompute=.TRUE.)
2841 : : call polynomial_coeff_init(one,coeffs_tmp(my_icoeff)%nterm,coefficients(icoeff2),&
2842 : : & coeffs_tmp(my_icoeff)%terms,name=name, &
2843 : 4869 : & check=.false.)
2844 : : !DMSG(coefficients(icoeff2)%debug_str)
2845 : : else
2846 : 88 : call polynomial_coeff_MPIsend(coeffs_tmp(my_icoeff), icoeff, rank_to_receive, comm)
2847 : : end if
2848 : : ! Free the coefficient
2849 : : !call polynomial_coeff_free(coeffs_tmp(my_icoeff))
2850 : : else
2851 [ + + + + : 479557 : if(any(my_newcoeffindexes(:)==icoeff3))then
+ + ]
2852 : 88 : icoeff2 = icoeff2 + 1
2853 : 88 : call polynomial_coeff_MPIrecv(coefficients(icoeff2), icoeff, rank_to_send, comm)
2854 : 88 : call polynomial_coeff_getName(name,coefficients(icoeff2),symbols,recompute=.TRUE.)
2855 : 88 : call polynomial_coeff_SetName(name,coefficients(icoeff2))
2856 : : !DMSG(coefficients(icoeff2)%debug_str)
2857 : : end if
2858 : : end if
2859 : : else
2860 : 106 : icoeff2 = icoeff2 + 1
2861 : : ! Get the name of this coefficient
2862 [ + - ]: 106 : if(my_rank==rank_to_send)then
2863 : 106 : call polynomial_coeff_getName(name,coeffs_tmp(my_icoeff),symbols,recompute=.TRUE.)
2864 : : call polynomial_coeff_init(one,coeffs_tmp(my_icoeff)%nterm,coefficients(icoeff2),&
2865 : : & coeffs_tmp(my_icoeff)%terms,name=name, &
2866 : 106 : & check=.false.)
2867 : : !DMSG(coefficients(icoeff2)%debug_str)
2868 : : ! Free the coefficient
2869 : : !call polynomial_coeff_free(coeffs_tmp(my_icoeff))
2870 : : end if
2871 : 106 : call polynomial_coeff_broadcast(coefficients(icoeff2),rank_to_send, comm)
2872 : : end if
2873 : : end do
2874 : :
2875 : :
2876 [ + + ]: 44 : if(need_verbose)then
2877 : 38 : write(message,'(1x,I0,2a)') ncoeff_tot,' coefficients generated ',ch10
2878 : 38 : call wrtout(ab_out,message,'COLL')
2879 : 38 : call wrtout(std_out,message,'COLL')
2880 : : end if
2881 : :
2882 : :
2883 : 44 : end subroutine distribute_coefficients_over_cpu
2884 : :
2885 : 44 : subroutine free_memory()
2886 : : !Final deallocation
2887 [ + - ]: 44 : ABI_SFREE(symbols)
2888 [ + - ]: 44 : ABI_SFREE(my_coeffindexes)
2889 [ + - ]: 44 : ABI_SFREE(my_newcoeffindexes)
2890 [ + - ]: 44 : ABI_SFREE(my_coefflist)
2891 : 44 : call polynomial_coeff_list_free(coeffs_tmp)
2892 : : !ABI_SFREE(coeffs_tmp)
2893 : 44 : end subroutine free_memory
2894 : :
2895 : : end subroutine polynomial_coeff_getNorder
2896 : : !!***
2897 : :
2898 : :
2899 : : !!****f* m_polynomial_coeff/computeNorder
2900 : : !! NAME
2901 : : !! computeNorder
2902 : : !!
2903 : : !! FUNCTION
2904 : : !! Recursive routine to compute the order N of a all the possible coefficient
2905 : : !! from the list list_symcoeff and list_symstr.
2906 : : !!
2907 : : !! INPUTS
2908 : : !! cell(3,nrpt) = indexes of the cells into the supercell (-1 -1 -1, 0 0 0 ...)
2909 : : !! compatibleCoeffs(ncoeff+nstr,ncoeff+nstr) = array with the list of compatible coefficients 0 or 1
2910 : : !! list_symcoeff(6,ncoeff_sym,nsym) = array with the list of the coefficients,
2911 : : !! for each coefficients (ncoeff_sym), we store the symmetrics(nsym)
2912 : : !! the 6th first dimensions are :
2913 : : !! 1 = direction of the IFC
2914 : : !! 2 = index of the atom number 1 (1=>natom)
2915 : : !! 3 = index of the atom number 2 (1=>natom)
2916 : : !! 4 = indexes of the cell of the second atom
2917 : : !! (the atom number 1 is always in the cell 0 0 0)
2918 : : !! 5 = weight of the term (-1 or 1)
2919 : : !! 6 = indexes of the symmetric
2920 : : !! list_symstr(nstr_sym,nsym) = array with the list of the strain and the symmetrics
2921 : : !! index_coeff_in(power_disp-1) = list of previous coefficients computed (start with 0)
2922 : : !! icoeff = current indexes of the cofficients (start we 1)
2923 : : !! icoeff_tot = current number of coefficients already computed (start we 0)
2924 : : !! natom = number of atoms in the unit cell
2925 : : !! nstr = number of coefficient for related to the atomic displacment into list_symcoeff
2926 : : !! nstr = number of coefficient for related to the strain into list_symstr
2927 : : !! ncoeff_out = number of maximum coefficients
2928 : : !! nrpt = number of cell
2929 : : !! nsym = number of symmetries in the system
2930 : : !! power_disp = initial power_disp to be computed (can be < power_disp_min,
2931 : : !! this routine will skip the firts power_disp)
2932 : : !! power_disp_min = minimal power_disp to be computed
2933 : : !! power_disp_max = maximum power_disp to be computed
2934 : : !! symbols(natom) = array with the symbols of each atoms (Sr,O,Ti,...)
2935 : : !! nbody = optional, number of body for the coefficients, for example:
2936 : : !! 0 => all the terms
2937 : : !! 1 => only (Sr_x-T_y)^power_disp and (Sr_x-T_y)^power_disp\eta^power_disp ...
2938 : : !! compute = logical, optional: TRUE if we store the coefficients
2939 : : !! FALSE just to count the number of coefficient
2940 : : !! anharmstr = logical, optional : TRUE, the anharmonic strain are computed
2941 : : !! FALSE, (default) the anharmonic strain are not computed
2942 : : !! distributed = logical, optional : True, the coefficients will be distributed on the CPU
2943 : : !!
2944 : : !! OUTPUT
2945 : : !! icoeff = current indexes of the cofficients (start we 1)
2946 : : !! icoeff_tot = current number of coefficients already computed (start we 0)
2947 : : !! polynomial_coeff<(type(polynomial_coeff_type)>(ncoeff_out) = array of datatype with
2948 : : !! the polynomial_coeff
2949 : : !!
2950 : : !! SOURCE
2951 : :
2952 : : recursive subroutine computeNorder(cell,coeffs_out,compatibleCoeffs,list_coeff,list_str,&
2953 : : & index_coeff_in,icoeff,icoeff_tot,natom,ncoeff,nstr,ncoeff_out,&
2954 : : & nrpt,nsym,power_disp,power_disp_min,power_disp_max,symbols,nbody,&
2955 : : & compute,anharmstr,spcoupling,distributed)
2956 : :
2957 : : !Arguments ---------------------------------------------
2958 : : !scalar
2959 : : integer,intent(in) :: natom,ncoeff,power_disp,power_disp_min,power_disp_max,ncoeff_out,nsym,nrpt,nstr,icoeff
2960 : : integer,intent(inout) :: icoeff_tot
2961 : : logical,optional,intent(in) :: compute,anharmstr,spcoupling,distributed
2962 : : integer,optional,intent(in) :: nbody
2963 : : !arrays
2964 : : integer,intent(in) :: cell(3,nrpt),compatibleCoeffs(ncoeff+nstr,ncoeff+nstr)
2965 : : integer,intent(in) :: list_coeff(6,ncoeff,nsym),list_str(nstr,nsym,2)
2966 : : integer,intent(in) :: index_coeff_in(power_disp-1)
2967 : : type(polynomial_coeff_type),intent(inout) :: coeffs_out(ncoeff_out)
2968 : : character(len=5),intent(in) :: symbols(natom)
2969 : : !Local variables ---------------------------------------
2970 : : !scalar
2971 : : integer :: ia,ib,ii,icoeff1,icoeff_tmp
2972 : : integer :: iterm,nbody_in,ncoeff_max,pa,pb
2973 : : integer :: ndisp_max,nterm_max
2974 : : real(dp):: coefficient
2975 : : logical :: need_compute,compatible,possible,need_anharmstr,need_spcoupling,need_distributed
2976 : : !arrays
2977 : : integer,allocatable :: index_coeff(:)
2978 : : character(len=200):: name
2979 : : type(polynomial_term_type),dimension(:),allocatable :: terms
2980 : : type(polynomial_coeff_type),allocatable :: coeffs_tmp(:)
2981 : : logical, allocatable :: reverse(:)
2982 : : ! *************************************************************************
2983 : :
2984 : : !Set the inputs
2985 : : need_compute = .TRUE.
2986 : : need_anharmstr = .TRUE.
2987 : : need_spcoupling = .TRUE.
2988 : : need_distributed = .FALSE.
2989 : : nbody_in = 0 !all kind of terms
2990 : : if(present(compute)) need_compute = compute
2991 : : if(present(nbody)) nbody_in = nbody
2992 : : if(present(anharmstr)) need_anharmstr = anharmstr
2993 : : if(present(spcoupling)) need_spcoupling = spcoupling
2994 : : if(present(distributed)) need_distributed = distributed
2995 : : if(power_disp <= power_disp_max)then
2996 : :
2997 : : ! Initialisation of variables
2998 : : nterm_max = nsym
2999 : : ncoeff_max = (ncoeff+nstr)
3000 : : ndisp_max = power_disp
3001 : : icoeff_tmp = 0
3002 : : ABI_MALLOC(coeffs_tmp,(ncoeff_max))
3003 : : ABI_MALLOC(terms,(nterm_max))
3004 : : ABI_MALLOC(index_coeff,(power_disp))
3005 : :
3006 : : index_coeff(1:power_disp-1) = index_coeff_in(:)
3007 : :
3008 : : do icoeff1=icoeff,ncoeff+nstr
3009 : : ! If the distance between the 2 coefficients is superior than the cut-off,
3010 : : ! we cycle
3011 : : ! If the power_disp is one check if icoeff1 is compatible with itself
3012 : : if(power_disp==1) then
3013 : : if(icoeff1 <= ncoeff .and. compatibleCoeffs(icoeff1,icoeff1)==0)then
3014 : : cycle
3015 : : end if
3016 : : end if
3017 : : if(compatibleCoeffs(icoeff,icoeff1)==0) cycle
3018 : :
3019 : : ! Reset the flag compatible and possible
3020 : : compatible = .TRUE.
3021 : : possible = .TRUE.
3022 : :
3023 : : index_coeff(power_disp) = icoeff1
3024 : : iterm = 0
3025 : : coefficient = one
3026 : :
3027 : : if(power_disp >= power_disp_min) then
3028 : : ABI_MALLOC(reverse,(ndisp_max))
3029 : : reverse(:) = .False.
3030 : : call generateTermsFromList(cell,index_coeff,list_coeff,list_str,ncoeff,&
3031 : : & ndisp_max,nrpt,nstr,nsym,iterm,terms, reverse=reverse)
3032 : : ABI_SFREE(reverse)
3033 : :
3034 : : if(iterm > 0)then
3035 : : ! Do some checks
3036 : : ! -------------
3037 : : ! 1-Check if the coefficient is full anharmonic strain and if we need to compute it
3038 : : if(terms(1)%ndisp == 0)then
3039 : : compatible = (need_anharmstr .or. need_spcoupling)
3040 : : possible = need_anharmstr
3041 : : end if
3042 : : ! 1-Check if the coefficient is strain-coupling and if we need to compute it
3043 : : if(terms(1)%nstrain > 0.and.terms(1)%ndisp > 0)then
3044 : : possible = need_spcoupling
3045 : : compatible = need_spcoupling
3046 : : end if
3047 : : ! ------------
3048 : : ! 2-Check if this terms is compatible with nbody
3049 : : if(nbody_in > 0)then
3050 : : pa = 1 ; pb = 1
3051 : : ia = 0 ; ib = 0
3052 : : ! Count the number of terms and the power_disp
3053 : : do ii=1,terms(1)%ndisp
3054 : : if(terms(1)%nstrain > 0) then
3055 : : pb = pb*terms(1)%power_disp(ii)
3056 : : ib = ib + 1
3057 : : else
3058 : : pa = pa*terms(1)%power_disp(ii)
3059 : : ia = ia + 1
3060 : : end if
3061 : : end do
3062 : : if(ia <= nbody_in)then
3063 : : if(ia==nbody_in.and.abs(mod(pa,2)) < tol16)then
3064 : : if(ib==0)then
3065 : : compatible = .FALSE.
3066 : : possible = .TRUE.
3067 : : else if (ib==nbody_in.and.abs(mod(pb,2)) < tol16) then
3068 : : compatible = .FALSE.
3069 : : possible = .TRUE.
3070 : : else
3071 : : possible = .FALSE.
3072 : : compatible = .FALSE.
3073 : : end if
3074 : : else
3075 : : possible = .FALSE.
3076 : : compatible = .FALSE.
3077 : : end if
3078 : : else
3079 : : compatible = .FALSE.
3080 : : possible = .FALSE.
3081 : : end if
3082 : : end if
3083 : :
3084 : : if(possible)then
3085 : : ! increase coefficients and set it
3086 : : icoeff_tmp = icoeff_tmp + 1
3087 : : icoeff_tot = icoeff_tot + 1
3088 : : call polynomial_coeff_init(coefficient,iterm,coeffs_tmp(icoeff_tmp),&
3089 : : & terms(1:iterm), check=.true.)
3090 : : !DMSG(coeffs_tmp(icoeff_tmp)%debug_str)
3091 : : end if
3092 : : end if
3093 : :
3094 : : ! Deallocate the terms
3095 : : do iterm=1,nterm_max
3096 : : call polynomial_term_free(terms(iterm))
3097 : : end do
3098 : : end if!end if power_disp < power_disp_min
3099 : :
3100 : : if(compatible)then
3101 : : call computeNorder(cell,coeffs_out,compatibleCoeffs,list_coeff,list_str,index_coeff,&
3102 : : & icoeff1,icoeff_tot,natom,ncoeff,nstr,ncoeff_out,nrpt,nsym,power_disp+1,&
3103 : : & power_disp_min,power_disp_max,symbols,nbody=nbody_in,compute=need_compute,&
3104 : : & anharmstr=need_anharmstr,spcoupling=need_spcoupling)
3105 : : end if
3106 : : end do
3107 : :
3108 : : ABI_SFREE(terms)
3109 : : ABI_SFREE(index_coeff)
3110 : :
3111 : : ! Transfer in the final array
3112 : : icoeff1 = 0
3113 : : do icoeff_tmp=1,ncoeff_max
3114 : : if (abs(coeffs_tmp(icoeff_tmp)%coefficient) > tol16)then
3115 : : ! Increase icoeff and fill the coeffs_out array
3116 : : icoeff_tot = icoeff_tot + 1
3117 : : if(need_compute)then
3118 : : name = ''
3119 : : ! Get the name of this coefficient
3120 : : call polynomial_coeff_getName(name,coeffs_tmp(icoeff_tmp),symbols,recompute=.TRUE.)
3121 : : call polynomial_coeff_init(one,coeffs_tmp(icoeff_tmp)%nterm,&
3122 : : & coeffs_out(icoeff_tot),coeffs_tmp(icoeff_tmp)%terms,&
3123 : : & name=name)
3124 : : !DMSG(coeffs_out(icoeff_tot)%debug_str)
3125 : : end if
3126 : : end if
3127 : : end do
3128 : : ! Deallocation
3129 : : do icoeff1=1,ncoeff_max
3130 : : call polynomial_coeff_free(coeffs_tmp(icoeff1))
3131 : : end do
3132 : : ABI_SFREE(coeffs_tmp)
3133 : : end if
3134 : :
3135 : : end subroutine computeNorder
3136 : : !!***
3137 : :
3138 : :
3139 : : !!****f* m_polynomial_coeff/computeCombinationFromList
3140 : : !! NAME
3141 : : !! computeCombinationFromList
3142 : : !!
3143 : : !! FUNCTION
3144 : : !! Recursive routine to compute the order N of a all the possible coefficient
3145 : : !! from the list list_symcoeff and list_symstr.
3146 : : !!
3147 : : !! INPUTS
3148 : : !! cell(3,nrpt) = indexes of the cells into the supercell (-1 -1 -1, 0 0 0 ...)
3149 : : !! compatibleCoeffs(ncoeff+nstr,ncoeff+nstr) = array with the list of compatible coefficients 0 or 1
3150 : : !! list_coeff(6,ncoeff_sym,nsym) = array with the list of the coefficients,
3151 : : !! for each coefficients (ncoeff_sym), we store the symmetrics(nsym)
3152 : : !! the 6th first dimensions are :
3153 : : !! 1 = direction of the IFC
3154 : : !! 2 = index of the atom number 1 (1=>natom)
3155 : : !! 3 = index of the atom number 2 (1=>natom)
3156 : : !! 4 = indexes of the cell of the second atom
3157 : : !! (the atom number 1 is always in the cell 0 0 0)
3158 : : !! 5 = weight of the term (-1 or 1)
3159 : : !! 6 = indexes of the symmetric
3160 : : !! list_str(nstr_sym,nsym) = array with the list of the strain and the symmetrics
3161 : : !! index_coeff_in(power_disp-1) = list of previous coefficients computed (start with 0)
3162 : : !! icoeff = current indexes of the combination (start with 1)
3163 : : !! max_power_strain = maximum order of the strain of the strain phonon coupling
3164 : : !! nmodel_tot = current number of combination already computed (start with 0)
3165 : : !! natom = number of atoms in the unit cell
3166 : : !! ncoeff = number of coefficient for related to the atomic displacment into list_symcoeff
3167 : : !! nstr = number of coefficient for related to the strain into list_symstr
3168 : : !! nmodel = number of maximum models
3169 : : !! nrpt = number of cell
3170 : : !! nsym = number of symmetries in the system
3171 : : !! For example, the sum of all the term like (Sr_y-O_y)^odd, are 0 by symetrie in cubic system.
3172 : : !! Here, we build a list with: 0 this term is not allowed for odd
3173 : : !! 1 this term is allowed for odd
3174 : : !! power_disp = initial power_disp to be computed (can be < power_disp_min,
3175 : : !! this routine will skip the first power_disp)
3176 : : !! power_disp_min = minimal power_disp to be computed
3177 : : !! power_disp_max = maximum power_disp to be computed
3178 : : !! symbols(natom) = array with the symbols of each atoms (Sr,O,Ti,...)
3179 : : !! nbody = optional, number of body for the coefficients, for example:
3180 : : !! 0 => all the terms
3181 : : !! 1 => only (Sr_x-T_y)^power_disp and (Sr_x-T_y)^power_disp\eta^power_disp ...
3182 : : !! compute = logical, optional: TRUE if we store the coefficients
3183 : : !! FALSE just to count the number of coefficient
3184 : : !! anharmstr = logical, optional : TRUE, the anharmonic strain are computed
3185 : : !! FALSE, (default) the anharmonic strain are not computed
3186 : : !! distributed = logical, optional : True, the coefficients will be distributed on the CPU
3187 : : !! only_odd_power = logical, optional : if TRUE return only odd power
3188 : : !! only_even_power= logical, optional : if TRUe return only even power
3189 : : !!
3190 : : !! OUTPUT
3191 : : !! icoeff = current indexes of the cofficients (start we 1)
3192 : : !! nmodel_tot = current number of coefficients already computed (start we 0)
3193 : : !! list_combination = list of the possible combination of coefficients
3194 : : !!
3195 : : !! SOURCE
3196 : :
3197 : 18044 : recursive subroutine computeCombinationFromList(cell,compatibleCoeffs,list_coeff,list_str,&
3198 : 18044 : & index_coeff_in,list_combination,icoeff,max_power_strain,&
3199 : : & natom,ncoeff,ncoeff_sym,iirred_comb,nirred_comb,nstr,nmodel,nrpt,nsym,power_disp,power_disp_min,&
3200 : 18044 : & power_disp_max,symbols,comm,nbody,only_odd_power,only_even_power,&
3201 : : & compute,anharmstr,spcoupling,disp)
3202 : :
3203 : : !Arguments ---------------------------------------------
3204 : : !scalar
3205 : : integer,intent(in) :: natom,ncoeff,ncoeff_sym,power_disp,power_disp_min,power_disp_max
3206 : : integer,intent(in) :: max_power_strain,nmodel,nsym,nrpt,nstr,comm,icoeff
3207 : : integer,intent(inout) :: nirred_comb,iirred_comb
3208 : : logical,optional,intent(in) :: compute,anharmstr,spcoupling,disp
3209 : : integer,optional,intent(in) :: nbody
3210 : : logical,optional,intent(in) :: only_odd_power,only_even_power
3211 : : !arrays
3212 : : integer,intent(in) :: cell(3,nrpt),compatibleCoeffs(ncoeff_sym+nstr,ncoeff_sym+nstr)
3213 : : integer,intent(in) :: list_coeff(6,ncoeff_sym,nsym),list_str(nstr,nsym,2)
3214 : : integer,intent(in) :: index_coeff_in(power_disp-1)
3215 : : integer,intent(out) :: list_combination(power_disp_max,nirred_comb)
3216 : : character(len=5),intent(in) :: symbols(natom)
3217 : : !Local variables ---------------------------------------
3218 : : !scalar
3219 : : integer :: icoeff1,icoeff2,nbody_in,nbody_count, ii
3220 : : integer :: ndisp_out,nstrain
3221 : : logical :: need_compute,compatible,possible,need_anharmstr,need_spcoupling
3222 : : logical :: need_only_odd_power,need_only_even_power,compute_sym,need_disp
3223 : : !arrays
3224 : 36088 : integer :: powers(power_disp)
3225 : 18044 : integer,allocatable :: index_coeff(:)
3226 : : ! *************************************************************************
3227 : :
3228 : : !Set the inputs
3229 : 18044 : need_compute = .TRUE.
3230 : 18044 : need_anharmstr = .TRUE.
3231 : 18044 : need_spcoupling = .TRUE.
3232 : 18044 : need_disp = .TRUE.
3233 : 18044 : need_only_odd_power = .FALSE.
3234 : 18044 : need_only_even_power = .FALSE.
3235 : 18044 : compute_sym = .FALSE. !Never compute the symmetric combinations here
3236 : 18044 : nbody_in = 0 !all kind of terms
3237 [ + - ]: 18044 : if(present(compute)) need_compute = compute
3238 [ + - ]: 18044 : if(present(nbody)) nbody_in = nbody
3239 [ + - ]: 18044 : if(present(anharmstr)) need_anharmstr = anharmstr
3240 [ + - ]: 18044 : if(present(spcoupling)) need_spcoupling = spcoupling
3241 [ + - ]: 18044 : if(present(disp)) need_disp = disp
3242 [ + - ]: 18044 : if(present(only_odd_power)) need_only_odd_power = only_odd_power
3243 [ + - ]: 18044 : if(present(only_even_power)) need_only_even_power = only_even_power
3244 : :
3245 : :
3246 [ + + ]: 18044 : if(power_disp <= power_disp_max)then
3247 : :
3248 : : ! Initialisation of variables
3249 [ - + + - : 16590 : ABI_MALLOC(index_coeff,(power_disp))
+ - + - -
+ - + ]
3250 [ + + ]: 17898 : index_coeff(1:power_disp-1) = index_coeff_in(:)
3251 : : ! Loop over ncoeff+nstr
3252 [ + + ]: 27884 : do icoeff1=icoeff,ncoeff+nstr
3253 : :
3254 : : ! Reset the flag compatible and possible
3255 : 22354 : compatible = .TRUE.
3256 : 22354 : possible = .TRUE.
3257 : :
3258 : : ! If the power_disp is one, we need to set icoeff to icoeff1
3259 [ + + ]: 22354 : if(power_disp==1) then
3260 [ + + + - ]: 860 : if(icoeff1<=ncoeff .and. compatibleCoeffs(icoeff,icoeff1)==0)then
3261 : : ! is_displacement and compatible
3262 : 22354 : compatible = .FALSE.
3263 : : end if
3264 : : end if
3265 : : ! If the distance between the 2 coefficients is superior than the cut-off, we cycle.
3266 [ + + ]: 64544 : do icoeff2=1,power_disp-1
3267 : : ! write(std_out,*) "icoeff1: ", icoeff1
3268 : : ! write(std_out,*) "icoeff2: ", icoeff2, "index_icoeff2: ", index_coeff(icoeff2)
3269 [ + + - + ]: 64544 : if(icoeff1 <= ncoeff .and. index_coeff(icoeff2) <=ncoeff)then
3270 [ + - ]: 4354 : if(compatibleCoeffs(index_coeff(icoeff2),icoeff1)==0)then
3271 : 42190 : compatible = .FALSE.
3272 : : end if
3273 : : endif
3274 : : end do
3275 : :
3276 [ + - ]: 22354 : if (.not.compatible) cycle !The distance is not compatible
3277 : :
3278 : : ! Set the index of the new coeff in the list
3279 : 22354 : index_coeff(power_disp) = icoeff1
3280 : : ! Do some checks
3281 : : ! -------------
3282 : : ! 1-Check if the coefficient is full anharmonic strain and if we need to compute it
3283 [ + + + + : 49342 : if(all(index_coeff > ncoeff))then
+ + ]
3284 : 8466 : compatible = (need_anharmstr .or. need_spcoupling)
3285 : 8466 : possible = need_anharmstr
3286 : : end if
3287 : : ! 2-Check if the coefficient is strain-coupling and if we need to compute it
3288 [ + + + + : 93688 : if(any(index_coeff <= ncoeff) .and. any(index_coeff > ncoeff))then
+ + + + +
+ ]
3289 : : ! write(std_out,*) "index_coeff", index_coeff,"need_spcoupling",need_spcoupling
3290 : 10566 : possible = need_spcoupling
3291 : 10566 : compatible = need_spcoupling
3292 [ + + + + : 40446 : if(count(index_coeff > ncoeff) > max_power_strain)then
+ + ]
3293 : 4338 : possible = .false.
3294 : 4338 : compatible = .false.
3295 : : end if
3296 : : end if
3297 : : ! 3-Check if the coefficient is only disp and if we need to compute it
3298 [ + + + + : 44346 : if(all(index_coeff <= ncoeff))then
+ + ]
3299 : : ! write(std_out,*) "index_coeff", index_coeff,"need_dis",need_disp
3300 : 3322 : compatible = (need_disp .or. need_spcoupling)
3301 : 3322 : possible = need_disp
3302 : : end if
3303 : : ! 4-Count number of Strain and number of displacements for compute symmetric terms
3304 : 22354 : nstrain = 0
3305 : 22354 : ndisp_out = 0
3306 [ + + ]: 86898 : do ii=1,power_disp
3307 [ + - + + ]: 86898 : if(index_coeff(ii) > 0 .and. index_coeff(ii) <= ncoeff)then
3308 : 21992 : ndisp_out = ndisp_out + 1
3309 : : else
3310 : 42552 : nstrain = nstrain +1
3311 : 42552 : index_coeff(ii) = index_coeff(ii) - ncoeff + ncoeff_sym
3312 : : end if
3313 : : end do
3314 : :
3315 : :
3316 [ + + ]: 22354 : if(power_disp >= power_disp_min) then
3317 : : ! count the number of body
3318 : 17008 : call get_powers(index_coeff, power_disp, powers)
3319 : 17008 : nbody_count = count(powers /= 0)
3320 [ + + + + : 24752 : if(any(mod(powers(1:power_disp),2) /=0) .and. need_only_even_power) then
+ + ]
3321 : 2360 : possible = .false.
3322 : : end if
3323 [ + + + + : 47796 : if(any(mod(powers(1:power_disp),2) ==0) .and. need_only_odd_power)then
- + ]
3324 : 0 : possible = .false.
3325 : : end if
3326 : : !if(nbody_count==1) then
3327 : : ! possible = .False.
3328 : : !end if
3329 : :
3330 : : ! Check the nbody flag
3331 [ - + ]: 17008 : if(nbody_in /= 0)then
3332 [ # # # # : 0 : if(power_disp-count(powers==0) > nbody_in) then
# # ]
3333 : : possible = .false.
3334 : : compatible = .false.
3335 : : end if
3336 : : end if
3337 : :
3338 [ + + ]: 17008 : if(possible) then
3339 : : ! increase coefficients and set it
3340 : : ! nmodel_tot = nmodel_tot + 1
3341 : : ! if(need_compute)then
3342 : : ! list_combination(1:power_disp,nmodel_tot) = index_coeff
3343 : : ! end if
3344 : : !nmodel_tot_test = 0
3345 : : ! !Start from second symmetry in Symmetric Combinations
3346 : : ! isym_in_test = 2
3347 : : ! idisp_in_test = power_disp
3348 : : ! ndisp_test = power_disp
3349 : : ! index_coeff_tmp = index_coeff
3350 : : !Count anharmonic strain terms
3351 [ + + ]: 6738 : if(ndisp_out == 0 .and. nstrain > 0)then
3352 : 448 : nirred_comb = nirred_comb +1
3353 : 448 : iirred_comb = iirred_comb +1
3354 [ + + ]: 448 : if(need_compute)then
3355 [ + + ]: 1064 : list_combination(1:power_disp,iirred_comb) = index_coeff
3356 : : endif
3357 : : else !Else counst symmetric terms of atomic displacement (pure disp or disp/strain)
3358 : : !Store index for each combination of irreducible terms to later parallely compute symmetric combinations
3359 : 6290 : nirred_comb = nirred_comb +1
3360 : 6290 : iirred_comb = iirred_comb +1
3361 [ + + ]: 6290 : if(need_compute)then
3362 [ + + ]: 12581 : list_combination(:,iirred_comb) = 0
3363 [ + + ]: 12559 : list_combination(:ndisp_out+nstrain,iirred_comb) = index_coeff
3364 : : endif
3365 : : end if !ndisp_out == 0 .and.n nstrain >0
3366 : : end if!possible
3367 : : end if!end if power_disp < power_disp_min
3368 : :
3369 : : !Change back to irreducible terms ncoeff_limit
3370 [ + + ]: 86898 : do ii=1,power_disp
3371 [ + + ]: 86898 : if(index_coeff(ii) > ncoeff_sym)then
3372 : 42552 : index_coeff(ii) = index_coeff(ii) + ncoeff - ncoeff_sym
3373 : : end if
3374 : : end do
3375 : :
3376 : : ! If the model is still compatbile with the input flags, we continue.
3377 [ + + ]: 27884 : if(compatible)then
3378 : : call computeCombinationFromList(cell,compatibleCoeffs,list_coeff,list_str,&
3379 : : & index_coeff,list_combination,icoeff1,max_power_strain,&
3380 : : & natom,ncoeff,ncoeff_sym,iirred_comb,nirred_comb,nstr,nmodel,nrpt,nsym,power_disp+1,&
3381 : : & power_disp_min,power_disp_max,symbols,comm,nbody=nbody_in,&
3382 : : & compute=need_compute,anharmstr=need_anharmstr,&
3383 : : & spcoupling=need_spcoupling,only_odd_power=need_only_odd_power,&
3384 : 17956 : & only_even_power=need_only_even_power,disp=need_disp)
3385 : : end if
3386 : : end do
3387 [ + - ]: 5530 : ABI_SFREE(index_coeff)
3388 : : end if
3389 : :
3390 [ - + ]: 18044 : end subroutine computeCombinationFromList
3391 : : !!***
3392 : :
3393 : 1404 : subroutine symlist_init(self, nsym, power)
3394 : : class(symlist_t), intent(inout) :: self
3395 : : integer, intent(in) :: nsym, power
3396 [ - + ]: 1404 : if(power>7) then
3397 : 0 : ABI_ERROR("The nbody or power of terms being generate is too large.")
3398 : : end if
3399 [ + + ]: 1404 : if (power>0) then
3400 : 1256 : self%nsym = nsym
3401 : 1256 : self%power = power
3402 : 1256 : self%max=self%nsym**self%power
3403 : 1256 : self%counter=0
3404 [ - + + - : 3768 : ABI_MALLOC(self%list, (self%power))
+ - - + -
+ - + ]
3405 [ + - ]: 148 : else if (power==0) then
3406 : 148 : self%max= 1
3407 : 148 : self%nsym = nsym
3408 : 148 : self%power = power
3409 : 148 : self%counter= 0
3410 [ + - - - : 444 : ABI_MALLOC(self%list, (self%power))
+ - - + -
+ - + ]
3411 : : else
3412 : 0 : ABI_BUG("The power of the combination should be at least 1")
3413 : : end if
3414 : 1404 : end subroutine symlist_init
3415 : :
3416 : 58274260 : subroutine symlist_next(self)
3417 : : class(symlist_t), intent(inout) :: self
3418 : : integer(dp) :: j, res
3419 : : integer(dp) :: d
3420 [ + + ]: 252668692 : self%list(:) = 0
3421 : 58274260 : self%counter = self%counter +1
3422 [ - + ]: 58274260 : if(self%counter>self%max) then
3423 : 0 : ABI_BUG("The iteration exceeded the number limit in symlist_next().")
3424 : : end if
3425 [ + + ]: 58274260 : if (self%power>0) then
3426 : 58274112 : self%list(1) = 1
3427 : 58274112 : d=self%counter-1
3428 [ + + ]: 252668544 : do j=1, self%power
3429 : 194394432 : res=mod(d, self%nsym)
3430 : 194394432 : d=d/self%nsym
3431 : 252668544 : self%list(self%power-j+1) = res+1
3432 : : end do
3433 : : end if
3434 : 58274260 : end subroutine symlist_next
3435 : :
3436 : 1404 : subroutine symlist_free(self)
3437 : : class(symlist_t) :: self
3438 : 1404 : self%max=0
3439 : 1404 : self%nsym = 0
3440 : 1404 : self%power = 0
3441 : 1404 : self%counter=0
3442 [ + - ]: 1404 : ABI_SFREE(self%list)
3443 : 1404 : end subroutine symlist_free
3444 : :
3445 : :
3446 : 2271 : subroutine computeSymmetricCombinations(array_combination, &
3447 : 2271 : & list_symcoeff, list_symstr, ndisp, nsym, index_coeff_in, &
3448 : : & ndisp_max, ncoeff, nsym_str, nstrain, &
3449 [ - + ]: 2271 : & compatibleCoeffs, compute, comm, only_even, max_nbody , irred_combinations , cell)
3450 : :
3451 : : integer,intent(in) :: ndisp,nsym,ndisp_max, ncoeff,nstrain,nsym_str
3452 : : integer,intent(in) :: comm
3453 : : logical,intent(in) :: compute
3454 : : logical,optional,intent(in) :: only_even
3455 : : !scalar
3456 : : !arrays
3457 : : !integer,intent(inout) ::
3458 : : type(int2d_array_type), intent(inout) :: array_combination ! list_combination(ndisp_max, nirred*nsym**(ndisp-1))
3459 : : integer,intent(in) :: list_symcoeff(6,ncoeff,nsym),index_coeff_in(ndisp+nstrain)
3460 : : integer,intent(in) :: list_symstr(6,nsym,2),compatibleCoeffs(ncoeff+nsym_str,ncoeff+nsym_str)
3461 : : integer, intent(in) :: max_nbody(:), cell(:, :)
3462 : : type(IrreducibleCombinations_T), intent(inout) :: irred_combinations
3463 : 2271 : type(symlist_t), target :: symlist
3464 : : !Local variables-------------------------------
3465 : :
3466 : : integer :: idisp,idisp2
3467 : : logical :: irreducible, need_only_even,possible
3468 : :
3469 : : !arrays
3470 : 5946 : integer :: index_coeff_tmp(ndisp),powers(ndisp),symcoeff_found(nsym)
3471 : 4542 : integer :: comb_to_test(ndisp_max)
3472 : : !integer,allocatable :: index_isym(:)
3473 : : !integer,allocatable :: symlist(:, :)
3474 : : integer(dp) :: isymlist
3475 : : !integer, pointer :: symlist_i(:)
3476 : : integer :: nbody, totpower, nbody_strain, totpower_strain, nbody_disp, totpower_disp
3477 : 2271 : type(polyform_t) :: polyform
3478 : 2271 : integer :: max_nbody_copy(size(max_nbody))
3479 : : !Source
3480 : :
3481 : : ABI_UNUSED(compute)
3482 : : ABI_UNUSED(comm)
3483 : :
3484 : 4542 : call irred_combinations%reset_array()
3485 : 2271 : need_only_even = .FALSE.
3486 [ + - ]: 2271 : if(present(only_even))need_only_even=only_even
3487 [ + + ]: 111279 : symcoeff_found = 0
3488 : :
3489 : 2271 : irreducible = .TRUE.
3490 : :
3491 : 2271 : call get_totpower_and_nbody(index_coeff_in(ndisp+1:ndisp+nstrain), nstrain, nbody_strain, totpower_strain)
3492 : 2271 : call get_totpower_and_nbody(index_coeff_in(1:ndisp), ndisp, nbody_disp, totpower_disp)
3493 : 2271 : totpower=totpower_disp+totpower_strain
3494 : 2271 : nbody=nbody_disp+nbody_strain
3495 : :
3496 [ + + ]: 2271 : if(nbody_disp==1) then
3497 [ + + ]: 1120 : if(.not. is_right_order(index_coeff_in(1))) return
3498 : : end if
3499 : :
3500 [ + + ]: 33726 : max_nbody_copy(:)=max_nbody(:)
3501 [ - + ]: 1606 : if(max_nbody(totpower)==-1) then
3502 : 0 : max_nbody_copy(totpower)=1
3503 [ + - ]: 1606 : else if(max_nbody(totpower)>=0) then
3504 : 1606 : max_nbody_copy(totpower)=max_nbody(totpower)
3505 : : else
3506 : 0 : ABI_ERROR("max_nbody should be -1, 0, or positive.")
3507 : : end if
3508 : :
3509 : : ! skip if max_nbody=0,
3510 : : ! and skip if nbody>max_nbody
3511 [ + - + + ]: 1606 : if(max_nbody_copy(totpower)==0 .or. nbody> max_nbody_copy(totpower)) then
3512 : : return
3513 : : end if
3514 : :
3515 : :
3516 : : block
3517 : : logical:: allsym
3518 : 2808 : integer :: ibody, ind(nbody_disp)
3519 : : ! allsym: allow all combination of symmetry adapted terms.
3520 : 1404 : allsym= ( max_nbody_copy(totpower)>=totpower)
3521 [ + + ]: 3621 : ind(:)=0
3522 : :
3523 [ + + ]: 1404 : if(allsym) then
3524 : 1282 : call symlist%init(nsym, ndisp)
3525 : : else
3526 : 122 : call polyform%from_expansion(index_coeff_in(1:ndisp), ndisp)
3527 : 122 : call symlist%init(nsym, nbody_disp)
3528 : : end if
3529 [ + + ]: 58275664 : do isymlist=1, symlist%max
3530 : : ! apply symmetry to the term, and check irreducibility
3531 : 58274260 : call symlist%next()
3532 : :
3533 [ + + ]: 58274260 : if(allsym) then
3534 [ + + ]: 252428656 : do idisp=1,ndisp
3535 : 252428656 : index_coeff_tmp(idisp) = list_symcoeff(6,index_coeff_in(idisp), symlist%list(idisp))
3536 : : end do !idisp=1,ndisp
3537 : : else
3538 [ + + ]: 240036 : do ibody=1, nbody_disp
3539 : 240036 : ind(ibody)= list_symcoeff(6,polyform%ind(ibody), symlist%list(ibody))
3540 : : end do
3541 : 80868 : call expand_poly(ind, polyform%order, nbody_disp, index_coeff_tmp(:ndisp) )
3542 : : end if
3543 : :
3544 : :
3545 : 58274260 : call get_powers(index_coeff_tmp, ndisp, powers)
3546 [ + + ]: 58274260 : if(.not. allsym) then
3547 : : ! only treat the terms with nbody< nbody_max
3548 [ + + + + ]: 319908 : nbody=count(powers/=0)+nbody_strain
3549 [ + + ]: 319908 : totpower=sum(powers) + totpower_strain
3550 [ + + ]: 80868 : if(totpower==0) cycle
3551 [ + - ]: 45622 : if(nbody> max_nbody_copy(totpower)) cycle
3552 : : end if
3553 [ + + + + : 110573158 : if(any(mod(powers(1:ndisp),2) /=0) .and. need_only_even) then
- + ]
3554 [ # # ]: 0 : index_coeff_tmp(:) = 0
3555 : : end if
3556 : :
3557 : :
3558 : : !Check if symmetric combination is allowed
3559 [ + + + + : 178858900 : if(.not. any(index_coeff_tmp == 0))then ! Check if term is allowed by distance
+ + ]
3560 [ + + ]: 109592554 : do idisp=1,ndisp-1
3561 [ + + ]: 231469162 : do idisp2=idisp+1,ndisp
3562 [ - + ]: 231469162 : if(compatibleCoeffs(index_coeff_tmp(idisp),index_coeff_tmp(idisp2)) == 0) then
3563 [ # # ]: 0 : index_coeff_tmp = 0
3564 : : exit
3565 : : end if
3566 : : enddo
3567 [ + - - + : 109592554 : if(all(index_coeff_tmp == 0))exit
+ - ]
3568 : : enddo
3569 : : endif
3570 : :
3571 : :
3572 [ + + + + : 178858900 : if(any(index_coeff_tmp == 0))then ! If symmetry doesn't point to another term or isn't allowed due to distance write zeros to filter after
+ + ]
3573 : : possible = .FALSE.
3574 : : else
3575 : : possible = .TRUE.
3576 : : endif
3577 : 1404 : if(possible)then
3578 : : !loop over displacements in term
3579 [ + + ]: 142133064 : comb_to_test(:) = 0
3580 [ + + - + : 29637330 : if(.not. (all(index_coeff_tmp == 0)))then ! If symmetry doesn't point to another term or isn't allowed due to distance write zeros to filter after
+ + ]
3581 [ + + ]: 139229588 : comb_to_test(:ndisp) = index_coeff_tmp
3582 : : end if! (any(index_coeff_tmp ==0))
3583 [ + + ]: 29637330 : if(nstrain /= 0)then !If SP coupling copy strain index
3584 [ + + ]: 5080504 : comb_to_test(ndisp+1:ndisp+nstrain) = index_coeff_in(ndisp+1:ndisp+nstrain)
3585 : : end if
3586 [ + + - + : 29637330 : if(.not. (all(index_coeff_tmp == 0) .and. nstrain==0)) then
+ + + - ]
3587 : 29637330 : irreducible=irred_combinations%add_irr(comb_to_test, list_symcoeff, list_symstr, ncoeff, nsym, ndisp_max, cell)
3588 : : endif
3589 : : end if ! need compute
3590 : : end do
3591 [ + + ]: 1404 : if(.not. allsym) then
3592 : 122 : call polyform%free()
3593 : : endif
3594 : : end block
3595 : :
3596 : 1404 : call array_combination%concate(irred_combinations%array)
3597 [ - + - + : 1404 : call symlist%free()
- + - + ]
3598 : :
3599 : : contains
3600 : 1120 : function is_right_order(index) result(right)
3601 : : integer, intent(in) :: index
3602 : : logical :: right
3603 : : integer :: ia, ib, ico
3604 : 1120 : ico=list_symcoeff(6, index, 1)
3605 : 1120 : ia=list_symcoeff(1, ico, 1)
3606 : 1120 : ib=list_symcoeff(2, ico, 1)
3607 : 1120 : right=ia>=ib
3608 : 1120 : end function is_right_order
3609 : : end subroutine computeSymmetricCombinations
3610 : :
3611 : : !!****f* m_polynomial_coeff/getCoeffFromList
3612 : : !!
3613 : : !! NAME
3614 : : !! getCoeffFromList
3615 : : !!
3616 : : !! FUNCTION
3617 : : !! get the index of a coefficient into the list_coeff
3618 : : !!
3619 : : !! INPUTS
3620 : : !! list_symcoeff(6,ncoeff_sym,nsym) = array with the list of the coefficients,
3621 : : !! for each coefficients (ncoeff_sym), we store the symmetrics(nsym)
3622 : : !! the 6th first dimensions are :
3623 : : !! 1 = direction of the IFC
3624 : : !! 2 = index of the atom number 1 (1=>natom)
3625 : : !! 3 = index of the atom number 2 (1=>natom)
3626 : : !! 4 = indexes of the cell of the second atom
3627 : : !! (the atom number 1 is always in the cell 0 0 0)
3628 : : !! 5 = weight of the term (-1 or 1)
3629 : : !! 6 = indexes of the symmetric
3630 : : !! ia = index of the atom 1
3631 : : !! ib = index of the atom 1
3632 : : !! irpt = indexes of the cell of the second atom
3633 : : !! mu = direction of the IFC
3634 : : !! ncoeff = number of total coefficients in the list
3635 : : !!
3636 : : !! OUTPUT
3637 : : !! coeff = index of the coefficient
3638 : : !!
3639 : : !! SOURCE
3640 : :
3641 : 4295040 : function getCoeffFromList(list_coeff,ia,ib,irpt,mu,ncoeff) result(coeff)
3642 : :
3643 : : !Arguments ------------------------------------
3644 : : !scalar
3645 : : integer,intent(in) :: ia,ib,irpt,mu,ncoeff
3646 : : integer :: coeff
3647 : : !arrays
3648 : : integer,intent(in) :: list_coeff(6,ncoeff)
3649 : : !Local variables-------------------------------
3650 : : !scalar
3651 : : integer :: icoeff
3652 : : !arrays
3653 : :
3654 : : ! *************************************************************************
3655 : 4295040 : coeff = 0
3656 [ + + ]: 3140190700 : do icoeff = 1,ncoeff
3657 : : if(mu==list_coeff(1,icoeff).and.&
3658 : : & ia==list_coeff(2,icoeff).and.&
3659 [ + + + + : 3136923190 : & ib==list_coeff(3,icoeff).and.&
+ + + + ]
3660 : 3267510 : & irpt==list_coeff(4,icoeff))then!.and.&
3661 : : !& abs(weight-list_coeff(5,icoeff)) < tol16) then
3662 : : coeff = icoeff
3663 : : exit
3664 : : end if
3665 : : end do
3666 : :
3667 : 4295040 : end function getCoeffFromList
3668 : : !!***
3669 : :
3670 : :
3671 : : !!****f* m_polynomial_coeff/generateTermsFromList
3672 : : !!
3673 : : !! NAME
3674 : : !! generateTermsFromList
3675 : : !!
3676 : : !! FUNCTION
3677 : : !! Compute for a given list of index the correspondig set of terms
3678 : : !!
3679 : : !! INPUTS
3680 : : !! cell(3,nrpt) = indexes of the cells into the supercell (-1 -1 -1, 0 0 0 ...)
3681 : : !! index_coeff_in(ndisp) = list of coefficients to be computed
3682 : : !! list_symcoeff(6,ncoeff,nsym) = array with the list of the coefficients,
3683 : : !! for each coefficients (ncoeff_sym), we store the symmetrics(nsym)
3684 : : !! the 6th first dimensions are :
3685 : : !! 1 = direction of the IFC
3686 : : !! 2 = index of the atom number 1 (1=>natom)
3687 : : !! 3 = index of the atom number 2 (1=>natom)
3688 : : !! 4 = indexes of the cell of the second atom
3689 : : !! (the atom number 1 is always in the cell 0 0 0)
3690 : : !! 5 = weight of the term (-1 or 1)
3691 : : !! 6 = indexes of the symmetry
3692 : : !! list_symstr(nstr,nsym) = array with the list of the strain and the symmetrics
3693 : : !! ncoeff = number of maximum coefficients in the list_symcoeff
3694 : : !! ndisp = number of maximum diplacement (phonon + strain)
3695 : : !! nrpt = number of cell
3696 : : !! nsym = number of symmetries in the system
3697 : : !!
3698 : : !! OUTPUT
3699 : : !! terms<(type(polynomial_term_type)>(nterm) = list of terms
3700 : : !! nterm = number of ouput terms
3701 : : !!
3702 : : !! SOURCE
3703 : :
3704 [ + + - + : 704608 : subroutine generateTermsFromList(cell,index_coeff,list_coeff,list_str,ncoeff,ndisp_max,&
- + - + -
+ - + - +
- + + + +
- ]
3705 : 7264 : & nrpt,nstr,nsym,nterm,terms, reverse)
3706 : :
3707 : : !Arguments ------------------------------------
3708 : : !scalar
3709 : : integer,intent(in) :: ndisp_max,ncoeff,nrpt,nstr,nsym
3710 : : integer,intent(out):: nterm
3711 : : !arrays
3712 : : integer,intent(in) :: index_coeff(ndisp_max)
3713 : : integer,intent(in) :: cell(3,nrpt),list_coeff(6,ncoeff,nsym)
3714 : : integer,intent(in) :: list_str(nstr,nsym,2)
3715 : : logical, intent(in) :: reverse(ndisp_max)
3716 : : type(polynomial_term_type),intent(out) :: terms(nsym)
3717 : : !Local variables-------------------------------
3718 : : !scalar
3719 : : integer :: ia,ib,icoeff_str,idisp,irpt
3720 : : integer :: isym,ndisp,nstrain,mu
3721 : : real(dp):: weight
3722 : : !arrays
3723 : 14528 : integer :: atindx(2,ndisp_max),cells(3,2,ndisp_max),dir_int(ndisp_max),strain(ndisp_max)
3724 : 7264 : integer :: power_disps(ndisp_max),power_strain(ndisp_max)
3725 : :
3726 : : ! *************************************************************************
3727 : 7264 : nterm = 0
3728 : : !Loop over symetries
3729 [ + + ]: 355936 : do isym=1,nsym
3730 : : !Treat this coeff
3731 : 348672 : weight = 1
3732 : 348672 : ndisp = 0
3733 : 348672 : nstrain = 0
3734 [ + + ]: 1424928 : do idisp=1,ndisp_max
3735 : : ! Get index of this displacement term
3736 : : ! Check if the index is not zero
3737 [ + + ]: 1076256 : if(index_coeff(idisp)==0) cycle
3738 [ + + ]: 1421280 : if(index_coeff(idisp)<=ncoeff)then
3739 : 904800 : ndisp = ndisp + 1
3740 : 904800 : mu = list_coeff(1,index_coeff(idisp),isym)
3741 [ - + ]: 904800 : if( reverse(idisp)) then
3742 : 0 : ia = list_coeff(3,index_coeff(idisp),isym)
3743 : 0 : ib = list_coeff(2,index_coeff(idisp),isym)
3744 : 0 : irpt = list_coeff(4,index_coeff(idisp),isym)
3745 : 0 : irpt = find_opposite_irpt(cell, irpt)
3746 : 0 : weight = -weight*list_coeff(5,index_coeff(idisp),isym)
3747 : : else
3748 : 904800 : ia = list_coeff(2,index_coeff(idisp),isym)
3749 : 904800 : ib = list_coeff(3,index_coeff(idisp),isym)
3750 : 904800 : irpt = list_coeff(4,index_coeff(idisp),isym)
3751 : 904800 : weight = weight*list_coeff(5,index_coeff(idisp),isym)
3752 : : end if
3753 : : ! Fill First term arrays
3754 : 904800 : atindx(1,idisp) = ia; atindx(2,idisp) = ib;
3755 : 904800 : dir_int(idisp) = mu
3756 : 904800 : power_disps(idisp) = 1
3757 [ + + ]: 3619200 : cells(:,1,idisp) = (/0,0,0/)
3758 [ + + ]: 3619200 : cells(:,2,idisp) = cell(:,irpt)
3759 : : else
3760 : 167808 : nstrain = nstrain + 1
3761 : 167808 : icoeff_str = index_coeff(idisp)-ncoeff
3762 : 167808 : strain(nstrain) = list_str(icoeff_str,isym,1)
3763 : 167808 : power_strain(nstrain) = 1
3764 : 167808 : weight = weight*list_str(icoeff_str,isym,2)
3765 : : end if
3766 : : end do
3767 : 348672 : nterm = nterm + 1
3768 : :
3769 : : call polynomial_term_init(atindx,cells,dir_int,ndisp,nstrain,terms(nterm),power_disps,&
3770 : 355936 : & power_strain,strain,weight,check=.true., index_coeff=index_coeff)
3771 : : end do!end do sym
3772 : 7264 : end subroutine generateTermsFromList
3773 : : !!***
3774 : :
3775 : :
3776 : : !!****f* m_polynomial_coeff/polynomial_coeff_getOrder1
3777 : : !!
3778 : : !! NAME
3779 : : !! polynomial_coeff_getOrder1
3780 : : !!
3781 : : !! FUNCTION
3782 : : !! Compute the first order polynomial coefficients from the list
3783 : : !!
3784 : : !! INPUTS
3785 : : !! cell(3,nrpt) = indexes of the cells into the supercell (-1 -1 -1, 0 0 0 ...)
3786 : : !! cutoff_in = cut-off for the inter atomic forces constants
3787 : : !! list_symcoeff(6,ncoeff_sym,nsym) = array with the list of the coefficients,
3788 : : !! for each coefficients (ncoeff_sym), we store the symmetrics(nsym)
3789 : : !! the 6th first dimensions are :
3790 : : !! 1 = direction of the IFC
3791 : : !! 2 = index of the atom number 1 (1=>natom)
3792 : : !! 3 = index of the atom number 2 (1=>natom)
3793 : : !! 4 = indexes of the cell of the second atom
3794 : : !! (the atom number 1 is always in the cell 0 0 0)
3795 : : !! 5 = weight of the term (-1 or 1)
3796 : : !! 6 = indexes of the symmetric
3797 : : !! natom = number of atoms in the unit cell
3798 : : !! nrpt = number of cell
3799 : : !! nsym = number of symmetries in the system
3800 : : !! symbols(natom) = array with the symbols of each atoms (Sr,O,Ti,...)
3801 : : !! comm = MPI communicator
3802 : : !!
3803 : : !! OUTPUT
3804 : : !! polynomial_coeff<(type(polynomial_coeff_type)>(ncoeff_out) = array of datatype with
3805 : : !! the polynomial_coeff
3806 : : !! ncoeff_out = number of coefficients
3807 : : !!
3808 : : !! SOURCE
3809 : :
3810 : 0 : subroutine polynomial_coeff_getOrder1(cell,coeffs_out,list_symcoeff,&
3811 : : & natom,ncoeff_out,ncoeff,nrpt,nsym,&
3812 : 0 : & symbols)
3813 : :
3814 : : !Arguments ------------------------------------
3815 : : !scalars
3816 : : integer,intent(in) :: natom,ncoeff,nsym,nrpt
3817 : : integer,intent(out) :: ncoeff_out
3818 : : !arrays
3819 : : integer,intent(in) :: cell(3,nrpt)
3820 : : integer,intent(in) :: list_symcoeff(6,ncoeff,nsym)
3821 : : character(len=5),intent(in) :: symbols(natom)
3822 : : type(polynomial_coeff_type),allocatable,intent(inout) :: coeffs_out(:)
3823 : : !Local variables-------------------------------
3824 : : !scalar
3825 : : integer :: ia,ib,icoeff,icoeff_tmp,irpt,irpt_ref
3826 : : integer :: isym,iterm,mu,ncoeff_max,ndisp,nstrain,nterm_max
3827 : : real(dp):: coefficient,weight
3828 : : !arrays
3829 : 0 : integer,allocatable :: atindx(:,:),cells(:,:,:),dir_int(:)
3830 : 0 : integer,allocatable :: power_disps(:),power_strain(:),strain(:)
3831 : : character(len=1) :: mutodir(9) = (/"x","y","z","1","2","3","4","5","6"/)
3832 : : character(len=200):: name
3833 : : character(len=500) :: message
3834 : 0 : type(polynomial_term_type),dimension(:),allocatable :: terms
3835 : 0 : type(polynomial_coeff_type),allocatable :: coeffs_tmp(:)
3836 : : !TEST_AM
3837 : : character(len=fnlen) :: filename
3838 : : !TEST_AM
3839 : : ! *************************************************************************
3840 : :
3841 : : !Initialisation of variables
3842 : 0 : nterm_max = nsym
3843 : 0 : ncoeff_max = ncoeff
3844 : 0 : ndisp = 1
3845 : 0 : nstrain = 0
3846 [ # # # # : 0 : ABI_MALLOC(coeffs_tmp,(ncoeff_max))
# # # # #
# # # #
# ]
3847 [ # # # # : 0 : ABI_MALLOC(terms,(nterm_max))
# # # # #
# # # #
# ]
3848 : :
3849 : :
3850 : 0 : icoeff_tmp = 0
3851 [ # # ]: 0 : ABI_MALLOC(atindx,(2,ndisp))
3852 [ # # ]: 0 : ABI_MALLOC(cells,(3,2,ndisp))
3853 [ # # ]: 0 : ABI_MALLOC(dir_int,(ndisp))
3854 [ # # ]: 0 : ABI_MALLOC(power_disps,(ndisp))
3855 [ # # ]: 0 : ABI_MALLOC(power_strain,(nstrain))
3856 [ # # ]: 0 : ABI_MALLOC(strain,(nstrain))
3857 : :
3858 : : !Found the ref cell
3859 : 0 : irpt_ref = 1
3860 : 0 : do irpt=1,nrpt
3861 : : if(all(cell(:,irpt)==0))then
3862 : 0 : irpt_ref = irpt
3863 : : exit
3864 : : end if
3865 : : end do
3866 : :
3867 : 0 : write(message,'(3a)') " Irreductible coefficient and associated atom 1, atom 2 and direction:",ch10,&
3868 : 0 : & " for the 1st order"
3869 : 0 : call wrtout(std_out,message,'COLL')
3870 : :
3871 [ # # ]: 0 : do icoeff=1,ncoeff
3872 : : ! Reset counter
3873 : 0 : iterm = 0
3874 : 0 : coefficient = one
3875 [ # # ]: 0 : do isym=1,nsym
3876 : : ndisp = 1
3877 : : nstrain = 0
3878 : 0 : mu = list_symcoeff(1,icoeff,isym)
3879 : 0 : ia = list_symcoeff(2,icoeff,isym)
3880 : 0 : ib = list_symcoeff(3,icoeff,isym)
3881 : 0 : irpt = list_symcoeff(4,icoeff,isym)
3882 : 0 : weight = list_symcoeff(5,icoeff,isym)
3883 : : ! Fill First term arrays
3884 : 0 : atindx(1,1) = ia; atindx(2,1) = ib;
3885 : 0 : dir_int(1) = mu
3886 : 0 : power_disps(1) = 1
3887 [ # # ]: 0 : cells(:,1,1) = (/0,0,0/)
3888 [ # # ]: 0 : cells(:,2,1) = cell(:,irpt)
3889 : 0 : iterm = iterm + 1
3890 : : call polynomial_term_init(atindx,cells,dir_int,ndisp,nstrain,terms(iterm),&
3891 : 0 : & power_disps,power_strain,strain,weight,check=.true.)
3892 : : end do!end do sym
3893 : :
3894 [ # # ]: 0 : if(iterm > 0)then
3895 : : ! increase coefficients and set it
3896 : 0 : icoeff_tmp = icoeff_tmp + 1
3897 : : call polynomial_coeff_init(coefficient,iterm,coeffs_tmp(icoeff_tmp), &
3898 : 0 : & terms(1:iterm),check=.true.)
3899 : : !DMSG(coeffs_tmp(icoeff_tmp)%debug_str)
3900 : : end if
3901 : :
3902 : : ! Deallocate the terms
3903 [ # # ]: 0 : do iterm=1,nterm_max
3904 : 0 : call polynomial_term_free(terms(iterm))
3905 : : end do
3906 : : end do!end do coeff_sym
3907 : :
3908 [ # # # # : 0 : ABI_SFREE(terms)
# # # # #
# # # # #
# # # # ]
3909 [ # # ]: 0 : ABI_SFREE(atindx)
3910 [ # # ]: 0 : ABI_SFREE(cells)
3911 [ # # ]: 0 : ABI_SFREE(dir_int)
3912 [ # # ]: 0 : ABI_SFREE(power_disps)
3913 [ # # ]: 0 : ABI_SFREE(power_strain)
3914 [ # # ]: 0 : ABI_SFREE(strain)
3915 : :
3916 : : !Count the number of terms
3917 : 0 : ncoeff_out = 0
3918 [ # # ]: 0 : do icoeff_tmp=1,ncoeff_max
3919 [ # # ]: 0 : if (abs(coeffs_tmp(icoeff_tmp)%coefficient) > tol16)then
3920 : 0 : ncoeff_out = ncoeff_out + 1
3921 : : end if
3922 : : end do
3923 : :
3924 : : !Transfer in the final array
3925 [ # # # # : 0 : ABI_MALLOC(coeffs_out,(ncoeff_out))
# # # # #
# # # # #
# # # # ]
3926 : 0 : icoeff = 0
3927 [ # # ]: 0 : do icoeff_tmp=1,ncoeff_max
3928 [ # # ]: 0 : if (abs(coeffs_tmp(icoeff_tmp)%coefficient) > tol16)then
3929 : : ! Get the name of this coefficient
3930 : 0 : call polynomial_coeff_getName(name,coeffs_tmp(icoeff_tmp),symbols,recompute=.TRUE.)
3931 : : ! Increase icoeff and fill the coeffs_out array
3932 : 0 : icoeff = icoeff + 1
3933 : : call polynomial_coeff_init(one,coeffs_tmp(icoeff_tmp)%nterm,&
3934 : : & coeffs_out(icoeff),coeffs_tmp(icoeff_tmp)%terms,&
3935 : 0 : & name=name)
3936 : : !DMSG(coeffs_out(icoeff)%debug_str)
3937 : :
3938 : 0 : write(message,'(2a)')' ',trim(name)
3939 : 0 : call wrtout(std_out,message,'COLL')
3940 : :
3941 [ # # ]: 0 : do iterm = 1,coeffs_tmp(icoeff_tmp)%nterm
3942 : 0 : write(message,'(a,I0,a,I0,2a)') ' Atom ',coeffs_tmp(icoeff_tmp)%terms(iterm)%atindx(1,1),&
3943 : 0 : & ' and atom ',coeffs_tmp(icoeff_tmp)%terms(iterm)%atindx(2,1),&
3944 : 0 : & ' in the direction ',mutodir(coeffs_tmp(icoeff_tmp)%terms(iterm)%direction(1))
3945 [ # # # # : 0 : if(any(coeffs_tmp(icoeff_tmp)%terms(iterm)%cell(:,2,1)/=0))then
# # ]
3946 : 0 : write(message,'(2a,I0,a,I0,a,I0,a)') trim(message),' in the cell ',&
3947 : 0 : & coeffs_tmp(icoeff_tmp)%terms(iterm)%cell(1,2,1),' ',&
3948 : 0 : & coeffs_tmp(icoeff_tmp)%terms(iterm)%cell(2,2,1),' ',&
3949 : 0 : & coeffs_tmp(icoeff_tmp)%terms(iterm)%cell(3,2,1),'.'
3950 : : end if
3951 : 0 : call wrtout(std_out,message,'COLL')
3952 : : end do
3953 : : end if
3954 : : end do
3955 : :
3956 : : !TEST_AM
3957 : 0 : filename = "terms_1st_order.xml"
3958 : 0 : call polynomial_coeff_writeXML(coeffs_out,ncoeff_out,filename=filename)
3959 : : !TEST_AM
3960 : :
3961 : 0 : write(message,'(a,1x,I0,a)') ch10,&
3962 : 0 : & ncoeff_out,' fitted coefficients for the 1st order '
3963 : 0 : call wrtout(ab_out,message,'COLL')
3964 : 0 : call wrtout(std_out,message,'COLL')
3965 : :
3966 : : !Deallocation
3967 [ # # ]: 0 : do icoeff=1,ncoeff_max
3968 : 0 : call polynomial_coeff_free(coeffs_tmp(icoeff))
3969 : : end do
3970 [ # # # # : 0 : ABI_SFREE(coeffs_tmp)
# # # # #
# # # # #
# # # # #
# # # ]
3971 : :
3972 : 0 : end subroutine polynomial_coeff_getOrder1
3973 : : !!***
3974 : :
3975 : : !!****f* m_polynomial_coeff/polynomial_coeff_getEvenAnhaStrain
3976 : : !!
3977 : : !! NAME
3978 : : !! polynomial_coeff_getEvenAnhaStrain
3979 : : !!
3980 : : !! FUNCTION
3981 : : !! Get even anharmonic strain terms in defined range of order
3982 : : !!
3983 : : !! INPUTS
3984 : : !!
3985 : : !!
3986 : : !! OUTPUT
3987 : : !! polynomial_coeff<(type(polynomial_coeff_type)>(ncoeff_out) = array of datatype with
3988 : : !! the polynomial_coeff
3989 : : !! ncoeff_out = number of coefficients
3990 : : !!
3991 : : !! SOURCE
3992 [ - + ]: 1 : subroutine polynomial_coeff_getEvenAnhaStrain(strain_terms,crystal,irred_ncoeff,power_strain,comm, max_nbody)
3993 : :
3994 : : !Arguments ------------------------------------
3995 : : type(polynomial_coeff_type),allocatable,intent(inout) :: strain_terms(:)
3996 : : type(crystal_t), intent(inout) :: crystal
3997 : : integer,intent(out) :: irred_ncoeff
3998 : : integer,intent(in) :: power_strain(2), max_nbody(:)
3999 : : integer,intent(in) :: comm
4000 : : !Local variables-------------------------------
4001 : : real(dp) :: cutoff,coeff_ini
4002 : : integer :: ncoeff,ncoeff_out,power_strph,option
4003 : : integer :: icoeff1,icoeff2
4004 : : integer:: nbody, npower
4005 : 1 : type(polynomial_coeff_type),allocatable :: strain_terms_tmp(:)
4006 : : integer :: sc_size(3)
4007 : : logical :: is_duplicate
4008 : : ! *************************************************************************
4009 : :
4010 : : !Initialize empty strain_terms array
4011 : : !ABI_MALLOC(strain_terms,(0))
4012 : 1 : irred_ncoeff = 0
4013 : :
4014 : : !Initial setup for call to polynomial_coeff_getNorder
4015 : 1 : cutoff = zero
4016 : 1 : power_strph = zero
4017 : 1 : option = 0
4018 : 1 : sc_size = (/1,1,1/)
4019 : 1 : coeff_ini = 1000000
4020 : :
4021 : : ! Get the higher order terms.
4022 : : call polynomial_coeff_getNorder(strain_terms_tmp,crystal,cutoff,ncoeff,ncoeff_out,power_strain,&
4023 : : & power_strph,option,sc_size,comm,anharmstr=.true.,spcoupling=.false.,&
4024 : : & only_odd_power=.false.,only_even_power=.true.,compute_symmetric=.false.,&
4025 : 1 : verbose=.false., max_nbody=max_nbody)
4026 [ - + ]: 1 : if(allocated(strain_terms))then
4027 : 0 : ABI_ERROR("strain_terms should be unallocated")
4028 : 0 : call polynomial_coeff_list_free(strain_terms)
4029 [ # # # # : 0 : ABI_SFREE(strain_terms)
# # # # #
# # # # #
# # # # #
# # # #
# ]
4030 : : end if
4031 : :
4032 : : !Build unique list of terms by direct comparison and append
4033 [ + + ]: 57 : do icoeff1=1,ncoeff_out
4034 : : !block
4035 : : !Get nbody and power for current term
4036 : 56 : nbody=strain_terms_tmp(icoeff1)%terms(1)%get_nbody()
4037 : 56 : npower=strain_terms_tmp(icoeff1)%terms(1)%get_total_power()
4038 : :
4039 : : !Check if term is duplicate of any previously added term
4040 : 56 : is_duplicate = .false.
4041 : : ! TODO : move this to a function
4042 [ + + ]: 56 : if(allocated(strain_terms)) then
4043 [ + + ]: 115 : do icoeff2=1,size(strain_terms)
4044 [ + + ]: 115 : if(coeffs_compare(strain_terms_tmp(icoeff1),strain_terms(icoeff2))) then
4045 : : is_duplicate = .true.
4046 : : exit
4047 : : endif
4048 : : enddo
4049 : : endif
4050 : :
4051 : : !Append if unique
4052 : : !print *, nbody , npower, is_duplicate
4053 : : !print *, "nbody:", nbody, "npower:", npower, "max_nbody(npower):", max_nbody(npower)
4054 : : !print *, "is_duplicate:", is_duplicate
4055 : :
4056 : :
4057 [ + + + + ]: 57 : if((.not. is_duplicate) .and. nbody<=max_nbody(npower) ) then
4058 : 2 : irred_ncoeff = irred_ncoeff + 1
4059 : 2 : call coeffs_list_append(strain_terms, strain_terms_tmp(icoeff1), .TRUE.)
4060 : : !DMSG(strain_terms(irred_ncoeff)%debug_str)
4061 : : endif
4062 : : !end block
4063 : : enddo
4064 : :
4065 : 1 : call polynomial_coeff_list_free(strain_terms_tmp)
4066 : :
4067 : :
4068 [ - + - - : 1 : end subroutine polynomial_coeff_getEvenAnhaStrain
- - - - -
- - - - -
- - - - -
- - - ]
4069 : : !!***
4070 : :
4071 : : !!****f* m_polynomial_coeff/coeffs_compare
4072 : : !! NAME
4073 : : !! equal
4074 : : !!
4075 : : !! FUNCTION
4076 : : !!
4077 : : !! INPUTS
4078 : : !!
4079 : : !! OUTPUT
4080 : : !!
4081 : : !! SOURCE
4082 : :
4083 : 828 : function coeffs_compare(c1,c2) result (res)
4084 : : !Arguments ------------------------------------
4085 : :
4086 : : !Arguments ------------------------------------
4087 : : type(polynomial_coeff_type), intent(in) :: c1,c2
4088 : : logical :: res
4089 : : !local
4090 : : !variable
4091 : : integer :: iterm1,iterm2
4092 : : !array
4093 : : !integer,allocatable :: blkval(:,:)
4094 : 1656 : integer :: blkval(2, max(c1%nterm,c2%nterm))
4095 : : ! *************************************************************************
4096 : 828 : res = .false.
4097 [ + + + + ]: 65916 : blkval = 0
4098 [ + + ]: 18936 : do iterm1=1,c1%nterm
4099 [ + - ]: 18108 : if(blkval(1,iterm1)==1)cycle!already found
4100 [ + + ]: 384354 : do iterm2=1,c2%nterm
4101 [ + + ]: 365418 : if(blkval(2,iterm2)==1)cycle!already found
4102 [ + + ]: 383514 : if(c1%terms(iterm1)==c2%terms(iterm2)) then
4103 : 12 : blkval(1,iterm1) = 1
4104 : 12 : blkval(2,iterm2) = 1
4105 : : end if
4106 : : end do
4107 : : end do
4108 [ + + + + : 864 : if(.not.any(blkval(:,:)==0))res = .true.
+ + + + ]
4109 : :
4110 : 828 : end function coeffs_compare
4111 : : !!***
4112 : :
4113 : :
4114 : : !!****f* m_polynomial_coeff/coeffs_list_conc
4115 : : !! NAME
4116 : : !! coeff_list_conc
4117 : : !!
4118 : : !! FUNCTION
4119 : : !!
4120 : : !! Concatenate list1 and list2 of type polynomial_coeff and store it in list_out
4121 : : !!
4122 : : !! INPUTS
4123 : : !!
4124 : : !! OUTPUT
4125 : : !!
4126 : : !! SOURCE
4127 : :
4128 [ # # # # : 0 : function coeffs_list_conc(coeff_list1,coeff_list2) result (coeff_list_out)
# # ]
4129 : : !Arguments ------------------------------------
4130 : :
4131 : : !Arguments ------------------------------------
4132 : : type(polynomial_coeff_type), intent(in) :: coeff_list1(:),coeff_list2(:)
4133 : : type(polynomial_coeff_type) :: coeff_list_out(size(coeff_list1)+size(coeff_list2))
4134 : : !local
4135 : : !variable
4136 : : integer :: ncoeff1,ncoeff2,ncoeff_out,i,j
4137 : : !array
4138 : : ! *************************************************************************
4139 : :
4140 : : !Get sizes of coeff_list1/2
4141 : 0 : ncoeff1 = size(coeff_list1)
4142 : 0 : ncoeff2 = size(coeff_list2)
4143 : 0 : ncoeff_out = ncoeff1 + ncoeff2
4144 : :
4145 [ # # ]: 0 : if(ncoeff_out/= size(coeff_list_out))then
4146 : 0 : ABI_ERROR("coeff_list_out should be allocated with the size of coeff_list1+coeff_list2")
4147 : : endif
4148 : :
4149 : :
4150 [ # # ]: 0 : if (.not. allocated(coeff_list1(i)%terms)) then
4151 : 0 : ABI_ERROR("copying coeff_list1(i) but coeff_list1(i)%terms is not allocated")
4152 : : endif
4153 : :
4154 [ # # ]: 0 : do i=1,ncoeff_out
4155 [ # # ]: 0 : if(i<=ncoeff1)then
4156 : : call polynomial_coeff_init(coeff_list1(i)%coefficient,coeff_list1(i)%nterm,coeff_list_out(i),coeff_list1(i)%terms,&
4157 : 0 : & coeff_list1(i)%name, check=.TRUE.)
4158 : : !if (trim(coeff_list1(i)%debug_str)=="unintialized" &
4159 : : !& .or. trim(coeff_list1(i)%debug_str)=="freed") then
4160 : : ! ABI_ERROR("coeff_list1(i) is uninitialized or freed")
4161 : : !else
4162 : : ! coeff_list_out(i)%debug_str = "copied from" // trim(coeff_list1(i)%debug_str)
4163 : : !end if
4164 : : else
4165 : 0 : j=i-ncoeff1
4166 : : call polynomial_coeff_init(coeff_list2(j)%coefficient,coeff_list2(j)%nterm,coeff_list_out(i),coeff_list2(j)%terms,&
4167 : 0 : & coeff_list2(j)%name,check=.TRUE.)
4168 : : !if (trim(coeff_list2(j)%debug_str)=="unintialized" &
4169 : : !& .or. trim(coeff_list2(j)%debug_str)=="freed") then
4170 : : ! ABI_ERROR("coeff_list2(j) is uninitialized or freed")
4171 : : !else
4172 : : ! coeff_list_out(i)%debug_str = "copied from" // trim(coeff_list2(j)%debug_str)
4173 : : !end if
4174 : : endif
4175 : : enddo
4176 : :
4177 [ # # ]: 0 : end function coeffs_list_conc
4178 : : !!***
4179 : :
4180 : :
4181 : : !!****f* m_polynomial_coeff/coeffs_list_conc_onsite
4182 : : !! NAME
4183 : : !! coeff_list_conc_onsite
4184 : : !!
4185 : : !! FUNCTION
4186 : : !!
4187 : : !! Concatenate list1 and list2 of type polynomial_coeff and store it in list1
4188 : : !!
4189 : : !! INPUTS
4190 : : !!
4191 : : !! OUTPUT
4192 : : !!
4193 : : !! SOURCE
4194 [ # # ]: 0 : subroutine coeffs_list_conc_onsite(coeff_list1,coeff_list2)
4195 : : !Arguments ------------------------------------
4196 : :
4197 : : !Arguments ------------------------------------
4198 : : type(polynomial_coeff_type), allocatable, intent(inout) :: coeff_list1(:)
4199 : : type(polynomial_coeff_type), intent(in) ::coeff_list2(:)
4200 : : !local
4201 : : !variable
4202 : : integer :: ncoeff1,ncoeff2,ncoeff_out
4203 : : !array
4204 : : type(polynomial_coeff_type), allocatable :: coeff_list_tmp(:)
4205 : : ! *************************************************************************
4206 : :
4207 : 0 : ncoeff1=size(coeff_list1)
4208 : 0 : ncoeff2=size(coeff_list2)
4209 : 0 : ncoeff_out = ncoeff1+ncoeff2
4210 : :
4211 : : ! copy list1 to tmp
4212 [ # # # # : 0 : ABI_MALLOC(coeff_list_tmp,(ncoeff1))
# # # # #
# # # #
# ]
4213 : : !coeff_list_tmp=coeff_list1
4214 : 0 : call coeffs_list_copy(coeff_list_tmp, coeff_list1)
4215 : :
4216 : : ! allocate new list1
4217 : 0 : call polynomial_coeff_list_free(coeff_list1)
4218 [ # # # # : 0 : ABI_MALLOC(coeff_list1,(ncoeff_out))
# # # # #
# # # # #
# # # # ]
4219 : 0 : coeff_list1(1:ncoeff1) = coeff_list_tmp
4220 : 0 : coeff_list1(ncoeff1+1:ncoeff_out) = coeff_list2
4221 : 0 : call polynomial_coeff_list_free(coeff_list_tmp)
4222 : :
4223 [ # # # # : 0 : end subroutine coeffs_list_conc_onsite
# # # # #
# # # # #
# # # # #
# # # ]
4224 : : !!***
4225 : :
4226 : :
4227 : :
4228 : :
4229 : : !!****f* m_polynomial_coeff/coeffs_list_append
4230 : : !! NAME
4231 : : !! coeff_list_append
4232 : : !!
4233 : : !! FUNCTION
4234 : : !!
4235 : : !! append on coeff to the list of coeffs: coeff_list
4236 : : !!
4237 : : !! INPUTS
4238 : : !!
4239 : : !! OUTPUT
4240 : : !!
4241 : : !! SOURCE
4242 : 3 : subroutine coeffs_list_append(coeff_list,coeff, check)
4243 : : !Arguments ------------------------------------
4244 : :
4245 : : !Arguments ------------------------------------
4246 : : type(polynomial_coeff_type), allocatable, intent(inout) ::coeff_list(:)
4247 : : type(polynomial_coeff_type), intent(inout) :: coeff
4248 : : logical, intent(in) :: check
4249 : : !local variable
4250 : 3 : type(polynomial_coeff_type), allocatable :: tmp(:)
4251 : : integer :: n1, n2
4252 : : !array
4253 : : ! *************************************************************************
4254 : :
4255 : : ! check if coeff is initialized
4256 [ - + ]: 3 : if (.not.allocated(coeff%terms)) then
4257 : 0 : ABI_BUG("terms for append is not allocated.")
4258 : : endif
4259 : :
4260 [ + + ]: 3 : if (.not. allocated(coeff_list)) then
4261 : 1 : n1 = 0
4262 : 1 : n2= 1
4263 [ - + - + : 3 : ABI_MALLOC(coeff_list,(n2))
+ + - + ]
4264 : : call polynomial_coeff_init(coeff%coefficient,coeff%nterm,coeff_list(1),coeff%terms,&
4265 : 1 : & coeff%name, check=check)
4266 : : !coeff_list(n2)%debug_str = "copied from" // trim(coeff%debug_str)
4267 : : else
4268 : 2 : n1=size(coeff_list)
4269 : 2 : n2= n1+1
4270 : : ! copy list1 to tmp
4271 [ - + + - : 8 : ABI_MALLOC(tmp,(n1))
+ - + - -
+ - + +
+ ]
4272 : 2 : tmp=coeff_list
4273 : : ! allocate new list1
4274 : 2 : call polynomial_coeff_list_free(coeff_list)
4275 [ - + + - : 10 : ABI_MALLOC(coeff_list,(n2))
+ - + - -
+ - + - +
+ + + - ]
4276 : :
4277 : 2 : call coeffs_list_copy(coeff_list(:n1), tmp)
4278 : :
4279 : : call polynomial_coeff_init(coeff%coefficient,coeff%nterm,coeff_list(n2),coeff%terms,&
4280 : 2 : & coeff%name, check=check)
4281 : :
4282 : : !coeff_list(n2)%debug_str = "copied from" // trim(coeff%debug_str)
4283 : 2 : call polynomial_coeff_list_free(tmp)
4284 [ - + - - : 2 : ABI_SFREE(tmp)
- - - - -
- - - - -
- - - - -
- - - ]
4285 : : end if
4286 : :
4287 [ - + - - : 3 : end subroutine coeffs_list_append
- - - - -
- - - - -
- - - - -
- - - ]
4288 : : !!***
4289 : :
4290 : : !!****f* m_polynomial_coeff/coeffs_list_truncate
4291 : : !! NAME
4292 : : !! coeffs_list_truncate
4293 : : !!
4294 : : !! FUNCTION
4295 : : !!
4296 : : !! Truncate a list of coeffs to keep only the first n elements
4297 : : !!
4298 : : !! INPUTS
4299 : : !! coeff_list = list of coefficients to truncate
4300 : : !! n = number of coefficients to keep
4301 : : !!
4302 : : !! OUTPUT
4303 : : !! coeff_list = truncated list with only first n coefficients
4304 : : !!
4305 : : !! SOURCE
4306 : :
4307 : 2 : subroutine coeffs_list_truncate(coeff_list, n)
4308 : : !Arguments ------------------------------------
4309 : :
4310 : : !Arguments ------------------------------------
4311 : : type(polynomial_coeff_type), allocatable, intent(inout) :: coeff_list(:)
4312 : : integer, intent(in) :: n
4313 : : !local variables
4314 : 2 : type(polynomial_coeff_type), allocatable :: tmp(:)
4315 : : character(len=500) :: message
4316 : : ! *************************************************************************
4317 : :
4318 [ - + ]: 2 : if(.not. allocated(coeff_list)) then
4319 : 0 : write(message,'(a)')'The input list of polynomial_coefficients is not allocated'
4320 : 0 : ABI_ERROR(message)
4321 : : endif
4322 : :
4323 [ - + ]: 2 : if(n > size(coeff_list)) then
4324 : 0 : write(message,'(a)')'n is larger than size of input list'
4325 : 0 : ABI_ERROR(message)
4326 : : endif
4327 : :
4328 [ - + ]: 2 : if(n < 0) then
4329 : 0 : write(message,'(a)')'n must be non-negative'
4330 : 0 : ABI_ERROR(message)
4331 : : endif
4332 : :
4333 : : !Copy first n elements to temp array
4334 [ - + + - : 12 : ABI_MALLOC(tmp,(n))
+ - + - -
+ - + +
+ ]
4335 : 2 : call coeffs_list_copy(tmp, coeff_list(1:n))
4336 : :
4337 : : !Free original array
4338 : 2 : call polynomial_coeff_list_free(coeff_list)
4339 : :
4340 : : !Reallocate and copy back
4341 [ - + + - : 10 : ABI_MALLOC(coeff_list,(n))
- + - + -
+ + + +
- ]
4342 : 2 : call coeffs_list_copy(coeff_list, tmp)
4343 : :
4344 : : !Clean up
4345 : 2 : call polynomial_coeff_list_free(tmp)
4346 [ - + - - : 2 : ABI_SFREE(tmp)
- - - - -
- - - - -
- - - - -
- - - ]
4347 : :
4348 [ - + - - : 2 : end subroutine coeffs_list_truncate
- - - - -
- - - - -
- - - - -
- - - ]
4349 : : !!***
4350 : :
4351 : :
4352 : 0 : subroutine coeffs_list_reduce_duplicate(self, crystal, sc_size, fit_iatom_in, cutoff_in , ndispmax )
4353 : : type(polynomial_coeff_type), allocatable, target, intent(inout) ::self(:)
4354 : 0 : type(polynomial_coeff_type), allocatable :: tmp(:)
4355 : 0 : type(IrreducibleCombinations_T) :: irred_combinations
4356 : : type(crystal_t), target, intent(inout) :: crystal
4357 : : integer, intent(in) :: sc_size(3)
4358 : : integer, intent(in), optional :: fit_iatom_in
4359 : : real(dp), intent(in), optional :: cutoff_in
4360 : : integer , intent(in) :: ndispmax
4361 : 0 : type(sympairs_t) :: sympairs
4362 : : !integer :: list_symcoeff(:,:,:),list_symstr(:,:,:), ncoeff_symsym, nsym, power_disps(:), cell(:,:)
4363 : 0 : logical :: mask(size(self))
4364 : : integer :: i, n, counter
4365 : 0 : call SymPairs_t_init(sympairs, crystal, sc_size, fit_iatom_in, cutoff_in)
4366 [ # # ]: 0 : mask=.True.
4367 : 0 : call irred_combinations%init()
4368 : : ! Find ir terms
4369 [ # # ]: 0 : do i=1, size(self)
4370 : : mask(i)=irred_combinations%add_irr(self(i)%terms(1)%index_coeff, &
4371 : : & sympairs%list_symcoeff, sympairs%list_symstr, size(sympairs%list_symcoeff(:,1,:)), &
4372 [ # # ]: 0 : & sympairs%nsym, ndispmax, sympairs%cell)
4373 : : end do
4374 [ # # # # ]: 0 : n=count(mask)
4375 : :
4376 : : ! copy self to tmp
4377 [ # # # # : 0 : ABI_MALLOC(tmp,( size(self)))
# # # # #
# # # #
# ]
4378 [ # # ]: 0 : do i=1, size(self)
4379 : : call polynomial_coeff_init(self(i)%coefficient,self(i)%nterm,tmp(i),self(i)%terms,&
4380 : : & self(i)%name, check=.True. &
4381 : 0 : )
4382 : : end do
4383 : :
4384 : :
4385 : 0 : call polynomial_coeff_list_free(self)
4386 [ # # # # : 0 : ABI_MALLOC(self,(n))
# # # # #
# # # # #
# # ]
4387 : 0 : counter=0
4388 [ # # ]: 0 : do i=1, size(tmp)
4389 [ # # ]: 0 : if(mask(i)) then
4390 : 0 : counter =counter +1
4391 : : call polynomial_coeff_init(tmp(i)%coefficient,tmp(i)%nterm,self(counter),tmp(i)%terms,&
4392 : : & tmp(i)%name, check=.True. &
4393 : 0 : )
4394 : : end if
4395 : : end do
4396 : 0 : call polynomial_coeff_list_free(tmp)
4397 : 0 : call sympairs%free()
4398 : 0 : call irred_combinations%free()
4399 [ # # # # : 0 : end subroutine coeffs_list_reduce_duplicate
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# ]
4400 : :
4401 : : !!****f* m_polynomial_coeff/coeffs_list_copy
4402 : : !! NAME
4403 : : !! coeff_list_copy
4404 : : !!
4405 : : !! FUNCTION
4406 : : !!
4407 : : !! Copy list1 to list2 of type polynomial_coeff
4408 : : !!
4409 : : !! INPUTS
4410 : : !!
4411 : : !! OUTPUT
4412 : : !!
4413 : : !! SOURCE
4414 : :
4415 [ - + - + : 10593 : subroutine coeffs_list_copy(coeff_list_out,coeff_list_in)
+ + - + -
- - - - -
- - - - -
- - - - -
+ + + - ]
4416 : : !Arguments ------------------------------------
4417 : :
4418 : : !Arguments ------------------------------------
4419 : : type(polynomial_coeff_type), intent(in) :: coeff_list_in(:)
4420 : : type(polynomial_coeff_type),intent(out) :: coeff_list_out(:)
4421 : : !local
4422 : : !variable
4423 : : integer :: ncoeff_in,ncoeff_out,ii
4424 : : logical :: check
4425 : : character(len=500):: message
4426 : : !array
4427 : : ! *************************************************************************
4428 : :
4429 : 63 : check = .true.
4430 : :
4431 : : !Get size of coeff_lists
4432 : 63 : ncoeff_in = size(coeff_list_in)
4433 : 63 : ncoeff_out = size(coeff_list_out)
4434 : :
4435 [ - + ]: 63 : if(ncoeff_in > ncoeff_out)then
4436 : 0 : write(message,'(a,a,a)')'The input list of polynomial_coefficients is larger',ch10,&
4437 : 0 : & 'than the output list you want it assign to. Check size of lists.'
4438 : 0 : ABI_ERROR(message)
4439 : : endif
4440 : :
4441 [ + + ]: 5328 : do ii=1,ncoeff_in
4442 : :
4443 [ - + ]: 5265 : if(.not. allocated(coeff_list_in(ii)%terms)) then
4444 : 0 : ABI_BUG("terms for copy is not allocated.")
4445 : : endif
4446 : : call polynomial_coeff_init(coeff_list_in(ii)%coefficient,coeff_list_in(ii)%nterm,&
4447 : : & coeff_list_out(ii),coeff_list_in(ii)%terms, &
4448 : 5328 : & coeff_list_in(ii)%name, check)
4449 : :
4450 : : !coeff_list_out%debug_str = "copied from " // coeff_list_in(ii)%debug_str
4451 : : enddo
4452 : 63 : end subroutine coeffs_list_copy
4453 : : !!***
4454 : :
4455 : : !!****f* m_polynomial_coeff/sort_combination
4456 : : !! NAME
4457 : : !! sort_combination
4458 : : !!
4459 : : !! FUNCTION
4460 : : !!
4461 : : !! Sort a list of integer from small to large if it contains zeros will be put to highest indexe
4462 : : !!
4463 : : !! INPUTS
4464 : : !!
4465 : : !! OUTPUT
4466 : : !!
4467 : : !! SOURCE
4468 : :
4469 : 638298090 : subroutine sort_combination(combination,n_int)
4470 : : !Arguments ------------------------------------
4471 : :
4472 : : !Arguments ------------------------------------
4473 : : !scalar
4474 : : integer,intent(in) :: n_int
4475 : : !array
4476 : : integer,intent(inout) :: combination(n_int)
4477 : : !local
4478 : : !variable
4479 : : integer :: j,k,cnt,tmp_int1,tmp_int2
4480 : : !array
4481 : : ! *************************************************************************
4482 : :
4483 : 638298090 : j=2
4484 [ + + ]: 2423729328 : do while(j <= n_int)
4485 : : k = j
4486 : : cnt = 1
4487 [ + + ]: 4451760567 : do while(k >= 2 .and. cnt == 1)
4488 [ + + + + ]: 4451760567 : if(combination(k-1) > combination(k) .and. combination(k) > 0)then
4489 : 1446388836 : tmp_int1 = combination(k-1)
4490 : 1446388836 : tmp_int2 = combination(k)
4491 : 1446388836 : combination(k) = tmp_int1
4492 : 1446388836 : combination(k-1) = tmp_int2
4493 : 1446388836 : k=k-1
4494 : : else
4495 : : cnt = cnt + 1
4496 : : end if
4497 : : end do
4498 : 1785431238 : j = j+1
4499 : : end do
4500 : :
4501 : 638298090 : end subroutine sort_combination
4502 : : !!***
4503 : :
4504 : : !!****f* m_polynomial_coeff/sort_combination_list
4505 : : !! NAME
4506 : : !! sort_combination_list
4507 : : !!
4508 : : !! FUNCTION
4509 : : !!
4510 : : !! Sort a list of integer list from small to large if it contains zeros will be put to highest index
4511 : : !!
4512 : : !! INPUTS
4513 : : !!
4514 : : !! OUTPUT
4515 : : !!
4516 : : !! SOURCE
4517 : :
4518 : : subroutine sort_combination_list(combination_list,n_int,n_list)
4519 : : !Arguments ------------------------------------
4520 : :
4521 : : !Arguments ------------------------------------
4522 : : !scalar
4523 : : integer,intent(in) :: n_int,n_list
4524 : : !array
4525 : : integer,intent(inout) :: combination_list(n_int,n_list)
4526 : : !local
4527 : : !variable
4528 : : integer :: i
4529 : : !array
4530 : : ! *************************************************************************
4531 : :
4532 : : do i=1,n_list
4533 : : call sort_combination(combination_list(:,i),n_int)
4534 : : end do
4535 : :
4536 : : end subroutine sort_combination_list
4537 : : !!***
4538 : :
4539 : :
4540 : : !!****f* m_polynomial_coeff/reduce_zero_combinations
4541 : : !! NAME
4542 : : !! reduce_zero_combinations
4543 : : !!
4544 : : !! FUNCTION
4545 : : !!
4546 : : !! Sort out list of zeros in a list of integers
4547 : : !!
4548 : : !! INPUTS
4549 : : !!
4550 : : !! OUTPUT
4551 : : !!
4552 : : !! SOURCE
4553 : :
4554 : 50 : subroutine reduce_zero_combinations(combination_list)
4555 : : !Arguments ------------------------------------
4556 : :
4557 : : !Arguments ------------------------------------
4558 : : !scalar
4559 : : !array
4560 : : integer,allocatable,intent(inout) :: combination_list(:,:)
4561 : : !local
4562 : : !variable
4563 : : integer :: i,j
4564 : 50 : integer,allocatable :: combination_list_tmp(:,:)
4565 : : !array
4566 : : ! *************************************************************************
4567 : :
4568 : : !Reduce zero combinations
4569 : 50 : i = 0
4570 [ + + ]: 14668 : do j=1,size(combination_list,2)
4571 [ + + + + : 14980 : if(any(combination_list(:,j) /= 0))then
+ + ]
4572 : 14514 : i = i + 1
4573 : : endif
4574 : : enddo
4575 : :
4576 [ - + + - : 200 : ABI_MALLOC(combination_list_tmp,(size(combination_list,1),i))
+ + + - +
- + + - +
- + ]
4577 : 50 : i = 0
4578 [ + + ]: 14668 : do j=1,size(combination_list,2)
4579 [ + + + + : 14980 : if(any(combination_list(:,j) /= 0))then
+ + ]
4580 : 14514 : i = i + 1
4581 [ + + ]: 59162 : combination_list_tmp(:,i) = combination_list(:,j)
4582 : : endif
4583 : : enddo
4584 [ + - ]: 50 : ABI_SFREE(combination_list)
4585 [ - + + - : 200 : ABI_MALLOC(combination_list,(size(combination_list_tmp,1),i))
+ + + - +
- + + - +
- + - + ]
4586 [ + - + - : 59262 : combination_list = combination_list_tmp
- + - - +
+ + + ]
4587 [ + - ]: 50 : ABI_SFREE(combination_list_tmp)
4588 : :
4589 : 50 : end subroutine reduce_zero_combinations
4590 : : !!***
4591 : :
4592 [ - + - + ]: 3667970 : function find_irpt(cells, cell) result(my_irpt)
4593 : : integer, intent(in) :: cells(:, :), cell(:)
4594 : : integer :: my_irpt, i
4595 : 3667970 : my_irpt=-1
4596 [ + + ]: 462164220 : do i =1, size(cells, 2)
4597 [ + + + + : 575871290 : if (all(cell==cells(:, i))) my_irpt=i
+ + ]
4598 : : end do
4599 [ - + ]: 3667970 : if(my_irpt==-1) ABI_BUG("cell not found.")
4600 : 3667970 : end function find_irpt
4601 : :
4602 : : ! for one r in a list of rpts, find the index of -r
4603 [ - + ]: 3667968 : function find_opposite_irpt(cells, irpt) result(n)
4604 : : integer, intent(in) :: cells(:, :), irpt
4605 : : integer :: n
4606 [ + - + + ]: 14671872 : n=find_irpt(cells, -cells(:, irpt))
4607 : 3667968 : end function find_opposite_irpt
4608 : :
4609 : 1 : subroutine SymPairs_t_init(self, crystal, sc_size, fit_iatom_in, cutoff_in)
4610 : : class(SymPairs_t), intent(inout) ::self
4611 : : type(crystal_t), target, intent(inout) :: crystal
4612 : : integer, intent(in) :: sc_size(3)
4613 : : integer, intent(in), optional :: fit_iatom_in
4614 : : real(dp), intent(in), optional :: cutoff_in
4615 : 1 : self%crystal=> crystal
4616 [ + + ]: 4 : self%sc_size(:) = sc_size(:)
4617 [ - + ]: 1 : if(present(fit_iatom_in)) then
4618 : 0 : self%fit_iatom = fit_iatom_in
4619 : : else
4620 : 1 : self%fit_iatom = -1
4621 : : end if
4622 [ - + ]: 1 : if(present(cutoff_in)) then
4623 : 0 : self%cutoff=cutoff_in
4624 : : else
4625 : 1 : self%cutoff=get_crystal_cutoff(crystal)
4626 : : end if
4627 : : call prepare_for_getList(crystal, sc_size, self%dist, self%cell, &
4628 : 1 : & self%natom, self%nsym, self%nrpt, self%range_ifc, self%symbols)
4629 : : call polynomial_coeff_getList(self%cell,self%crystal,self%dist, &
4630 : : &self%list_symcoeff,self%list_symstr,&
4631 : : &self%natom,self%nstr_sym,self%ncoeff_sym,self%nrpt, &
4632 : : &self%range_ifc,self%cutoff,sc_size=self%sc_size,&
4633 [ - + - + ]: 1 : &fit_iatom=fit_iatom_in)
4634 : 1 : end subroutine SymPairs_t_init
4635 : :
4636 : 2 : subroutine SymPairs_t_free(self)
4637 : : class(SymPairs_t), intent(inout) ::self
4638 : 2 : nullify(self%crystal)
4639 [ + + ]: 2 : ABI_SFREE(self%list_symcoeff)
4640 [ + + ]: 2 : ABI_SFREE(self%list_symstr)
4641 [ + + ]: 2 : ABI_SFREE(self%cell)
4642 [ + + ]: 2 : ABI_SFREE(self%dist)
4643 [ + + ]: 2 : ABI_SFREE(self%symbols)
4644 : 2 : end subroutine SymPairs_t_free
4645 : :
4646 [ - + + + : 588 : subroutine SymPairs_t_generateTerms(self, index_coeff, power, nterm, terms, reverse)
- + - + -
+ - + - +
- + - + +
+ + - ]
4647 : : class(SymPairs_t), intent(inout) ::self
4648 : : integer,intent(in) :: index_coeff(:)
4649 : : integer, intent(in) :: power
4650 : : integer, intent(out) :: nterm
4651 : : type(polynomial_term_type),intent(out) :: terms(self%nsym)
4652 : : logical, optional, intent(in) :: reverse(power)
4653 : 6 : logical :: reverse_a(power)
4654 : : integer :: ndisp_max
4655 [ - + ]: 6 : if(present(reverse))then
4656 [ + + ]: 42 : reverse_a(:) = reverse(:)
4657 : : else
4658 [ # # ]: 0 : reverse_a(:) = .False.
4659 : : end if
4660 : 6 : ndisp_max=size(index_coeff)
4661 : : ! Note that ncoeff_sym is not the same as ncoeff_symsym
4662 : : call generateTermsFromList(self%cell,index_coeff,self%list_symcoeff, &
4663 : : &self%list_symstr,size(self%list_symcoeff, 2),power,self%nrpt,self%nstr_sym,self%nsym, &
4664 [ - + - - : 6 : &nterm,terms, reverse=reverse_a)
- - - + ]
4665 [ + - ]: 6 : end subroutine SymPairs_t_generateTerms
4666 : :
4667 : 1 : subroutine SymPairs_t_final(self)
4668 : : type(SymPairs_t), intent(inout) :: self
4669 : : ! Call the free method to ensure proper cleanup
4670 : 0 : call self%free()
4671 : 1 : end subroutine SymPairs_t_final
4672 : :
4673 : 60 : subroutine IrreducibleCombinations_init(self)
4674 : : class(IrreducibleCombinations_t), intent(inout) :: self
4675 : 60 : call self%table%init(2023)
4676 : 60 : end subroutine IrreducibleCombinations_init
4677 : :
4678 : 60 : subroutine IrreducibleCombinations_free(self)
4679 : : class(IrreducibleCombinations_t), intent(inout) :: self
4680 : 60 : call self%table%free()
4681 : 60 : call self%array%finalize()
4682 : 60 : end subroutine IrreducibleCombinations_free
4683 : :
4684 : 2271 : subroutine IrreducibleCombinations_reset_array(self)
4685 : : class(IrreducibleCombinations_t), intent(inout) :: self
4686 : 2271 : call self%array%finalize()
4687 : 0 : end subroutine IrreducibleCombinations_Reset_Array
4688 : :
4689 : :
4690 : 59280640 : function IrreducibleCombinations_add_irr(self, combination, list_symcoeff, &
4691 [ - + ]: 29640320 : & list_symstr, ncoeff_sym, nsym, ndisp, cell ) result(irreducible)
4692 : : class(IrreducibleCombinations_t), intent(inout) :: self
4693 : : integer, intent(inout) :: combination(:)
4694 : : integer,intent(in) :: ncoeff_sym,ndisp,nsym
4695 : : integer,intent(in) :: list_symcoeff(6,ncoeff_sym,nsym)
4696 : : integer,intent(in) :: list_symstr(6,nsym,2), cell(:, :)
4697 : :
4698 : 59280640 : integer :: combination_cmp_tmp(ndisp), combination_sorted(ndisp)
4699 : : logical :: irreducible
4700 : : integer :: n, isym
4701 : : ! check if the combination, or its symmetry equivalent are already in the
4702 : : ! table. If not, add it to the table.
4703 : :
4704 : : ABI_UNUSED(cell)
4705 : 29640320 : irreducible=.True.
4706 [ + + ]: 142145284 : combination_sorted(:)=combination(:)
4707 : 29640320 : call sort_combination(combination_sorted,size(combination))
4708 : 29640320 : n=size(combination)
4709 [ + + ]: 29640320 : if ( self%table%has_key_intn(combination_sorted, n)) then
4710 : 29630178 : irreducible=.False.
4711 : 29630178 : return
4712 : : else
4713 : 28127608 : isym = 2
4714 [ + + ]: 608667912 : do while(isym <= nsym)
4715 : 608657770 : call symcomb(combination_sorted, combination_cmp_tmp, isym)
4716 : 608657770 : call sort_combination(combination_cmp_tmp,size(combination_cmp_tmp))
4717 [ + + ]: 608667912 : if (self%table%has_key_intn(combination_cmp_tmp, n)) then
4718 : 28117466 : irreducible=.False.
4719 : 28117466 : return
4720 : : else
4721 : 580540304 : isym = isym + 1
4722 : : end if
4723 : : end do
4724 : : end if
4725 : 10142 : call self%table%put_intn(combination_sorted, 0.0_dp, n)
4726 : 10142 : call self%array%push(combination)
4727 : : contains
4728 : :
4729 : : function reverse(combination, isym, cell) result(rcomb)
4730 : : integer :: combination(:), isym, cell(:, :)
4731 : : integer :: rcomb(size(combination))
4732 : : integer :: ia, ib, irpt, weight, idisp
4733 : : do idisp=1,ndisp
4734 : : if(combination(idisp) /= 0 .and. combination(idisp) <= ncoeff_sym)then
4735 : : ia = list_symcoeff(3,combination(idisp),isym)
4736 : : ib = list_symcoeff(2,combination(idisp),isym)
4737 : : irpt = list_symcoeff(4,combination(idisp),isym)
4738 : : irpt = find_opposite_irpt(cell, irpt)
4739 : : rcomb(idisp)=getCoeffFromList(list_symcoeff, ia, ib, irpt, weight, ncoeff_sym)
4740 : : else
4741 : : rcomb(idisp)=combination(idisp)
4742 : : end if
4743 : : end do
4744 : : end function reverse
4745 : :
4746 : :
4747 : : ! function is_onebody_term_in_right_order(combination) result(ret)
4748 : : ! integer :: combination(:)
4749 : : ! integer :: ia, ib, irpt, weight, idisp
4750 : : ! do idisp=1,ndisp
4751 : : ! if(combination(idisp) /= 0 .and. combination(idisp) <= ncoeff_sym)then
4752 : : ! ia = list_symcoeff(2,combination(idisp),isym)
4753 : : ! ib = list_symcoeff(3,combination(idisp),isym)
4754 : : ! !irpt = list_symcoeff(4,combination(idisp),isym)
4755 : : ! !rcomb(idisp)=getCoeffFromList(list_symcoeff, ia, ib, irpt, weight, ncoeff_sym)
4756 : : ! end if
4757 : : ! end do
4758 : : ! end function is_onebody_term_in_right_order
4759 : :
4760 : :
4761 : : ! get the symmetry equivalent combination with symmetry index isym.
4762 [ - + - + ]: 608657770 : subroutine symcomb(combination, combination_cmp_tmp, isym)
4763 : : integer :: combination(:), combination_cmp_tmp(:), isym
4764 : : integer :: idisp, istrain
4765 [ + + ]: 2919882134 : do idisp=1,ndisp
4766 [ + + + + ]: 2919882134 : if(combination(idisp) /= 0 .and. combination(idisp) <= ncoeff_sym)then
4767 : 2252390476 : combination_cmp_tmp(idisp)=list_symcoeff(6,combination(idisp),isym)
4768 [ + + ]: 58833888 : else if(combination(idisp) > ncoeff_sym)then
4769 : 51841488 : istrain = combination(idisp) - ncoeff_sym
4770 : 51841488 : combination_cmp_tmp(idisp)=list_symstr(istrain,isym,1) + ncoeff_sym
4771 : : else
4772 : 6992400 : combination_cmp_tmp(idisp) = 0
4773 : : endif
4774 : : enddo
4775 : 608657770 : end subroutine symcomb
4776 : : end function IrreducibleCombinations_add_irr
4777 : :
4778 : :
4779 [ - + ]: 58295932 : subroutine get_powers(index_coeff_tmp, ndisp, powers)
4780 : : integer, intent(in) :: ndisp
4781 : : integer, intent(in) :: index_coeff_tmp(:) ! size ndisp
4782 : : integer, intent(inout) :: powers(ndisp) ! size ndisp
4783 : : integer :: ii, jj
4784 [ + + ]: 252829848 : powers(:) = 1
4785 [ + + ]: 252829848 : do ii=1,ndisp
4786 [ + + ]: 194533916 : if(index_coeff_tmp(ii)==0) then
4787 : 62684658 : powers(ii)=0
4788 : 62684658 : cycle
4789 : : end if
4790 [ + + ]: 131849258 : if (powers(ii) == 0) cycle
4791 [ + + ]: 347664292 : do jj=ii+1,ndisp
4792 [ + + ]: 169106858 : if (powers(jj) == 0) cycle
4793 [ + + ]: 296993428 : if(index_coeff_tmp(ii)==index_coeff_tmp(jj))then
4794 : 11587756 : powers(ii) = powers(ii) + 1
4795 : 11587756 : powers(jj) = 0
4796 : : end if
4797 : : end do
4798 : : end do
4799 : 58295932 : end subroutine get_powers
4800 : :
4801 : :
4802 [ - + ]: 4542 : subroutine get_totpower_and_nbody(index_coeff_tmp, ndisp, nbody, totpower)
4803 : : ! get the total power and number of body in the polynomial
4804 : : integer, intent(in) :: ndisp
4805 : : integer, intent(in) :: index_coeff_tmp(:) ! size ndisp
4806 : : integer, intent(inout) :: nbody, totpower
4807 : 9084 : integer :: powers(ndisp), i, p ! size ndisp
4808 : 4542 : nbody=0
4809 : 4542 : totpower=0
4810 : 4542 : call get_powers(index_coeff_tmp, ndisp, powers)
4811 [ + + ]: 11384 : do i=1, ndisp
4812 : 6842 : p=powers(i)
4813 [ + + ]: 11384 : if(p>0) then
4814 : 5732 : nbody=nbody+1
4815 : 5732 : totpower=totpower+p
4816 : : end if
4817 : : end do
4818 : 4542 : end subroutine get_totpower_and_nbody
4819 : :
4820 : :
4821 : 0 : subroutine polyform_init(self, n)
4822 : : class(polyform_t), intent(inout) :: self
4823 : : integer, intent(in) :: n
4824 : 0 : self%n=n
4825 [ # # # # : 0 : ABI_MALLOC(self%ind, (n))
# # # # #
# # # #
# ]
4826 [ # # # # : 0 : ABI_MALLOC(self%order, (n))
# # # # #
# ]
4827 : 0 : end subroutine polyform_init
4828 : :
4829 : 122 : subroutine polyform_free(self)
4830 : : class(polyform_t), intent(inout) :: self
4831 : 122 : self%n=0
4832 [ + - ]: 122 : ABI_SFREE(self%ind)
4833 [ + - ]: 122 : ABI_SFREE(self%order)
4834 [ - + ]: 122 : ABI_SFREE(self%expansion)
4835 : 122 : end subroutine polyform_free
4836 : :
4837 [ - + ]: 122 : subroutine polyform_from_expansion(polyform, ind, ndisp)
4838 : : class(polyform_t) , intent(inout):: polyform
4839 : : integer, intent(in) :: ind(:) ! size ndisp
4840 : : integer, intent(in) :: ndisp
4841 : 244 : integer :: ii, jj, powers(ndisp)
4842 : 122 : call get_powers(ind, ndisp, powers)
4843 [ + + + + ]: 308 : polyform%n=count(powers>0)
4844 [ + + + - : 366 : ABI_MALLOC(polyform%ind, (polyform%n))
+ - + + -
+ - + -
+ ]
4845 [ + + + - : 244 : ABI_MALLOC(polyform%order, (polyform%n))
- + - + -
+ ]
4846 : 122 : jj=1
4847 [ + + ]: 308 : do ii=1, ndisp
4848 [ + + ]: 308 : if(powers(ii)>0) then
4849 : 120 : polyform%ind(jj)=ind(ii)
4850 : 120 : polyform%order(jj)=powers(ii)
4851 : 120 : jj=jj+1
4852 : : end if
4853 : : end do
4854 : 122 : end subroutine polyform_from_expansion
4855 : :
4856 : :
4857 [ - + - + : 80868 : subroutine expand_poly(ind, order, n, list)
- + ]
4858 : : integer, intent(in) :: ind(:), order(:), n
4859 : : integer, intent(out) :: list(:)
4860 : : integer :: i, j,s, e ! s,e: start, end
4861 : 80868 : s=1
4862 [ + + ]: 240036 : do i=1, n
4863 : 159168 : e=s+order(i)-1
4864 [ + + ]: 398208 : do j=s, e
4865 : 398208 : list(j) = ind(i)
4866 : : end do
4867 : 240036 : s=e+1
4868 : : end do
4869 : 80868 : end subroutine expand_poly
4870 : :
4871 : :
4872 : 0 : subroutine polyform_expand(self, expansion)
4873 : : class(polyform_t), intent(inout) :: self
4874 : : integer, allocatable, intent(out) :: expansion(:)
4875 [ # # # # : 0 : ABI_MALLOC(expansion, (sum(self%order)))
# # # # #
# # # # #
# # ]
4876 : 0 : call expand_poly(self%ind, self%order, self%n, expansion)
4877 : 0 : end subroutine polyform_expand
4878 : :
4879 : :
4880 : : !!****f* m_polynomial_coeff/polynomial_coeff_final
4881 : : !! NAME
4882 : : !! polynomial_coeff_final
4883 : : !!
4884 : : !! FUNCTION
4885 : : !! Finalizer for polynomial_coeff_type to properly clean up allocated memory
4886 : : !!
4887 : : !! INPUTS
4888 : : !! self = The polynomial_coeff_type instance being finalized
4889 : : !!
4890 : : !! SOURCE
4891 : :
4892 : : subroutine polynomial_coeff_final(self)
4893 : : type(polynomial_coeff_type), intent(inout) :: self
4894 : : !print *, "Finalizing polynomial_coeff_type: ", trim(self%debug_str)
4895 : : call polynomial_coeff_free(self)
4896 : : end subroutine polynomial_coeff_final
4897 : : !!***
4898 : :
4899 : : subroutine starts_with(str, prefix, result)
4900 : : character(len=*), intent(in) :: str
4901 : : character(len=*), intent(in) :: prefix
4902 : : logical, intent(out) :: result
4903 : : result = str(1:len(prefix)) == prefix
4904 : : end subroutine starts_with
4905 : :
4906 : :
4907 [ - + - - : 4 : end module m_polynomial_coeff
+ - + + -
+ + - - +
- + + - -
+ - + - +
# # # # #
# ]
4908 : : !!***
|
