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 32 : 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 53304 : 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 26652 : check_in = .false.
248 26652 : if(present(check)) check_in = check
249 :
250 26652 : if (present(debug_str)) then
251 1770 : polynomial_coeff%debug_str = debug_str
252 : else
253 24882 : polynomial_coeff%debug_str = "initialized"
254 : end if
255 :
256 26652 : 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 7130 : nterm_tmp = nterm
287 7130 : coefficient_tmp = coefficient
288 249512 : weights(:) = terms(:)%weight
289 : end if!end Check
290 :
291 26652 : if(present(name))then
292 17610 : name_tmp = name
293 : else
294 9042 : name_tmp = ""
295 : end if
296 :
297 : !Initilisation
298 26652 : polynomial_coeff%name = name_tmp
299 26652 : polynomial_coeff%nterm = nterm_tmp
300 26652 : polynomial_coeff%coefficient = coefficient_tmp
301 942134 : ABI_MALLOC(polynomial_coeff%terms,(polynomial_coeff%nterm))
302 26652 : iterm1 = 0
303 1050335 : do ii = 1,nterm
304 1050335 : if(abs(weights(ii)) > tol16)then
305 862178 : iterm1 = iterm1 + 1
306 862178 : call polynomial_term_copy(terms(ii), polynomial_coeff%terms(iterm1))
307 862178 : if(sum(polynomial_coeff%terms(iterm1)%power_strain)>1) then
308 : endif
309 : end if
310 : end do
311 26652 : 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 : !!***
|
