Line data Source code
1 : !!****f* ABINIT/m_polynomial_term
2 : !!
3 : !! NAME
4 : !! m_polynomial_term
5 : !!
6 : !! FUNCTION
7 : !! Module with the datatype polynomial terms
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2010-2026 ABINIT group (AM)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public Licence, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
15 : !!
16 : !! SOURCE
17 :
18 :
19 : #if defined HAVE_CONFIG_H
20 : #include "config.h"
21 : #endif
22 :
23 : #include "abi_common.h"
24 :
25 : module m_polynomial_term
26 :
27 : use defs_basis
28 : use m_errors
29 : use m_abicore
30 :
31 : implicit none
32 :
33 : public :: polynomial_term_init
34 : public :: polynomial_term_free
35 : public :: polynomial_term_copy
36 : public :: terms_compare
37 : public :: terms_compare_inverse
38 : !!***
39 :
40 : !!****t* m_polynomial_term/polynomial_term_type
41 : !! NAME
42 : !! polynomial_term_type
43 : !!
44 : !! FUNCTION
45 : !! Datatype for a terms (displacements or strain)
46 : !! related to a polynomial coefficient
47 : !!
48 : !! SOURCE
49 :
50 : type, public :: polynomial_term_type
51 :
52 : integer :: ndisp = 0
53 : ! Number of displacement for this terms
54 : ! 1 for (X_y-O_y)^3, 2 for (X_y-O_y)(X_x-O_y)^2...
55 :
56 : integer :: nstrain = 0
57 : ! Number of strain for this terms
58 :
59 : integer :: nindex = -1
60 : ! Number of index
61 :
62 : integer,allocatable :: atindx(:,:)
63 : ! atindx(2,ndisp)
64 : ! Indexes of the atoms a and b in the unit cell
65 :
66 : integer,allocatable :: cell(:,:,:)
67 : ! cell(3,2,ndisp)
68 : ! indexes of the cell of the atom a and b
69 :
70 : integer,allocatable :: direction(:)
71 : ! direction(ndisp)
72 : ! direction of the displacement
73 :
74 : integer,allocatable :: strain(:)
75 : ! strain(nstrain)
76 : ! strain
77 :
78 : integer,allocatable :: power_disp(:)
79 : ! power_disp(ndisp)
80 : ! power of the displacement 2 (X_z-O_z)^2 or 1 for (X_y-O_y)^1
81 :
82 : integer,allocatable :: power_strain(:)
83 : ! power_strain(nstrain)
84 : ! power of the strain 2 (\eta_1)^2 or 1 for (\eta_1)^1
85 :
86 : real(dp) :: weight = zero
87 : ! weight of the term
88 :
89 : integer,allocatable :: index_coeff(:)
90 :
91 : ! a string to identify the term
92 : character(len=500) :: debug_str = ''
93 :
94 : contains
95 : procedure :: get_nbody
96 : procedure :: get_total_power
97 : !final :: polynomial_term_finalizer
98 : end type polynomial_term_type
99 : !!***
100 :
101 :
102 : interface operator (==)
103 : module procedure terms_compare
104 : end interface
105 :
106 : CONTAINS !===========================================================================================
107 :
108 :
109 62 : function get_nbody(self) result(nbody)
110 : class(polynomial_term_type), intent(in) :: self
111 : integer :: nbody
112 62 : nbody = self%ndisp + self%nstrain
113 62 : end function get_nbody
114 :
115 :
116 56 : function get_total_power(self) result(total_power)
117 : class(polynomial_term_type), intent(in) :: self
118 : integer :: total_power
119 182 : total_power = sum(self%power_disp) + sum(self%power_strain)
120 56 : end function get_total_power
121 :
122 :
123 : !!****f* m_polynomial_term/polynomial_term_init
124 : !!
125 : !! NAME
126 : !! polynomial_term_init
127 : !!
128 : !! FUNCTION
129 : !! Initialize a polynomial_term_init for given set of displacements/strain
130 : !!
131 : !! INPUTS
132 : !! atindx(2) = Indexes of the atoms a and b in the unit cell
133 : !! cell(3,2) = Indexes of the cell of the atom a and b
134 : !! direction = direction of the perturbation => 1,2,3 for atomic displacement
135 : !! ndisp = Number of displacement for this terms
136 : !! nstrain = Number of strain for this terms
137 : !! power_disp = power_disp of the displacement 2 (X_z-O_z)^2 or 1 for (X_y-O_y)^1
138 : !! power_strain = power_strain of the strain 2 (\eta)^2 or 1 for (\eta)^1
139 : !! strain(nstrain) = index of strain, 1 2 3 4 5 or 6
140 : !! weight = Weight of the term
141 : !! check = optional,logical => if TRUE, the term will be check, same displacement/strain
142 : !! are gathered in the same displacement but with an higher
143 : !! power_disp. For example:
144 : !! ((Sr_y-O1_y)^1(Sr_y-O1_y)^1 => (Sr_y-O1_y)^2)
145 : !! if FALSE, default, do nothing
146 : !!
147 : !! OUTPUT
148 : !! polynomial_term<type(polynomial_term)> = polynomial_term datatype is now initialized
149 : !!
150 : !! SOURCE
151 :
152 1257116 : subroutine polynomial_term_init(atindx,cell,direction,ndisp,nstrain,polynomial_term,power_disp,&
153 1113818 : & power_strain,strain,weight,check, index_coeff, debug_str)
154 :
155 : implicit none
156 :
157 : !Arguments ------------------------------------
158 : !scalars
159 : integer, intent(in) :: ndisp,nstrain
160 : real(dp),intent(in) :: weight
161 : logical,optional,intent(in) :: check
162 : !arrays
163 : integer, intent(in) :: atindx(2,ndisp)
164 : integer, intent(in) :: cell(3,2,ndisp)
165 : integer, intent(in) :: direction(ndisp),power_disp(ndisp)
166 : integer, intent(in) :: strain(nstrain),power_strain(nstrain)
167 : integer, optional, intent(in) :: index_coeff(:)
168 : type(polynomial_term_type), intent(out) :: polynomial_term
169 : !Local variables-------------------------------
170 : !scalar
171 : integer :: idisp1,idisp2,ndisp_tmp,nstrain_tmp
172 : logical :: check_in
173 : !arrays
174 1257116 : integer :: power_disp_tmp(ndisp),power_strain_tmp(nstrain)
175 : character(500) :: msg
176 :
177 : character(*), optional :: debug_str
178 :
179 : ! *************************************************************************
180 :
181 1257116 : if (present(debug_str)) then
182 0 : polynomial_term%debug_str = debug_str
183 : else
184 1257116 : polynomial_term%debug_str = ' '
185 : end if
186 :
187 : !Do some checks
188 1257116 : if (size(atindx,2) /= ndisp) then
189 0 : write(msg,'(a)')' atindx and ndisp have not the same size'
190 0 : ABI_ERROR(msg)
191 : end if
192 :
193 1257116 : if (size(cell,3) /= ndisp) then
194 0 : write(msg,'(a)')' cell and ndisp have not the same size'
195 0 : ABI_ERROR(msg)
196 : end if
197 :
198 1257116 : if (size(direction) /= ndisp) then
199 0 : write(msg,'(a)')' direction and ndisp have not the same size'
200 0 : ABI_ERROR(msg)
201 : end if
202 :
203 1257116 : if (size(power_disp) /= ndisp) then
204 0 : write(msg,'(a)')' power_disp and ndisp have not the same size'
205 0 : ABI_ERROR(msg)
206 : end if
207 :
208 1257116 : if (size(power_strain) /= nstrain) then
209 0 : write(msg,'(a)')' power_strain and nstrain have not the same size'
210 0 : ABI_ERROR(msg)
211 : end if
212 :
213 1257116 : if (size(strain) /= nstrain) then
214 0 : write(msg,'(a)')' strain and nstrain have not the same size'
215 0 : ABI_ERROR(msg)
216 : end if
217 : !First free datatype before init
218 1257116 : call polynomial_term_free(polynomial_term)
219 1257116 : check_in = .false.
220 :
221 : !Copy the power array before check
222 4449860 : power_disp_tmp(:) = power_disp(:)
223 1706438 : power_strain_tmp(:) = power_strain(:)
224 :
225 1257116 : if(present(check)) check_in = check
226 :
227 1257116 : if(check_in)then
228 : !Check if displacement are identical, in this case
229 : !increase the power_disp
230 4449860 : do idisp1=1,ndisp
231 10353938 : do idisp2=idisp1,ndisp
232 : if (idisp1/=idisp2.and.&
233 : & atindx(1,idisp1) == atindx(1,idisp2).and.&
234 : & atindx(2,idisp1) == atindx(2,idisp2).and.&
235 : & direction(idisp1) == direction(idisp2).and.&
236 : & all(cell(:,1,idisp1)==cell(:,1,idisp2)).and.&
237 36267590 : & all(cell(:,2,idisp1)==cell(:,2,idisp2)).and.&
238 3192744 : & power_disp_tmp(idisp2) > 0 )then
239 60576 : power_disp_tmp(idisp1) = power_disp_tmp(idisp1) + 1
240 60576 : power_disp_tmp(idisp2) = 0
241 : end if
242 : end do
243 : end do
244 :
245 : ! Count the number of power_disp avec the previous check
246 : ! or just remove the power_disp equal to zero
247 : ndisp_tmp = 0
248 4449860 : do idisp1=1,ndisp
249 4449860 : if(power_disp_tmp(idisp1) > zero) then
250 3132168 : ndisp_tmp = ndisp_tmp + 1
251 : end if
252 : end do
253 :
254 : !Check if strain are identical, in this case
255 : !increase the power_strain
256 1706438 : do idisp1=1,nstrain
257 2220563 : do idisp2=idisp1,nstrain
258 : if (idisp1/=idisp2.and.&
259 514125 : & strain(idisp1) == strain(idisp2).and.&
260 449322 : & power_strain_tmp(idisp2) > 0 )then
261 15744 : power_strain_tmp(idisp1) = power_strain_tmp(idisp1) + 1
262 15744 : power_strain_tmp(idisp2) = 0
263 : end if
264 : end do
265 : end do
266 :
267 : ! Count the number of power_strain avec the previous check
268 : ! or just remove the power_strain equal to zero
269 : nstrain_tmp = 0
270 1706438 : do idisp1=1,nstrain
271 1706438 : if(power_strain_tmp(idisp1) > zero) then
272 433578 : nstrain_tmp = nstrain_tmp + 1
273 : end if
274 : end do
275 :
276 : else
277 : ndisp_tmp = ndisp
278 : nstrain_tmp = nstrain
279 : end if!end check
280 :
281 : !init the values
282 1257116 : polynomial_term%ndisp = ndisp_tmp
283 1257116 : polynomial_term%nstrain = nstrain_tmp
284 1257116 : polynomial_term%weight = weight
285 :
286 3771348 : ABI_MALLOC(polynomial_term%atindx,(2,polynomial_term%ndisp))
287 3771348 : ABI_MALLOC(polynomial_term%direction,(polynomial_term%ndisp))
288 3771348 : ABI_MALLOC(polynomial_term%cell,(3,2,polynomial_term%ndisp))
289 2514232 : ABI_MALLOC(polynomial_term%power_disp,(polynomial_term%ndisp))
290 3771348 : ABI_MALLOC(polynomial_term%power_strain,(polynomial_term%nstrain))
291 2514232 : ABI_MALLOC(polynomial_term%strain,(polynomial_term%nstrain))
292 :
293 : !Transfert displacement
294 1257116 : idisp2 = 0
295 4449860 : do idisp1=1,ndisp
296 4449860 : if(power_disp_tmp(idisp1) > zero)then
297 3132168 : idisp2 = idisp2 + 1
298 3132168 : polynomial_term%direction(idisp2) = direction(idisp1)
299 : polynomial_term%power_disp(idisp2) = power_disp_tmp(idisp1)
300 9396504 : polynomial_term%atindx(:,idisp2) = atindx(:,idisp1)
301 28189512 : polynomial_term%cell(:,:,idisp2) = cell(:,:,idisp1)
302 3132168 : polynomial_term%power_disp(idisp2) = power_disp_tmp(idisp1)
303 : end if
304 : end do
305 :
306 : !Transfert strain
307 : idisp2 = 0
308 1706438 : do idisp1=1,nstrain
309 1706438 : if(power_strain_tmp(idisp1) > zero)then
310 433578 : idisp2 = idisp2 + 1
311 433578 : polynomial_term%power_strain(idisp2) = power_strain_tmp(idisp1)
312 433578 : polynomial_term%strain(idisp2) = strain(idisp1)
313 : end if
314 : end do
315 :
316 1257116 : if (present(index_coeff)) then
317 1113818 : polynomial_term%nindex = size(index_coeff)
318 3341454 : ABI_MALLOC(polynomial_term%index_coeff, (polynomial_term%nindex))
319 4505326 : polynomial_term%index_coeff(:)=index_coeff(:)
320 : else
321 143298 : polynomial_term%nindex = -1
322 143298 : ABI_MALLOC(polynomial_term%index_coeff, (0))
323 : end if
324 :
325 :
326 :
327 :
328 1257116 : end subroutine polynomial_term_init
329 : !!***
330 :
331 :
332 : !!****f* m_polynomial_term/polynomial_term_free
333 : !!
334 : !! NAME
335 : !! polynomial_term_free
336 : !!
337 : !! FUNCTION
338 : !! Free polynomial_term
339 : !!
340 : !! INPUTS
341 : !! polynomial_term<type(polynomial_term)> = datatype to free
342 : !!
343 : !! OUTPUT
344 : !! polynomial_term<type(polynomial_term)> = datatype to free
345 : !!
346 : !! SOURCE
347 :
348 2564356 : subroutine polynomial_term_free(polynomial_term)
349 :
350 : implicit none
351 :
352 : !Arguments ------------------------------------
353 : !scalars
354 : !arrays
355 : type(polynomial_term_type), intent(inout) :: polynomial_term
356 : !Local variables-------------------------------
357 : !scalar
358 : !arrays
359 :
360 : ! *************************************************************************
361 :
362 2564356 : polynomial_term%ndisp = 0
363 2564356 : polynomial_term%nstrain = 0
364 2564356 : polynomial_term%weight = zero
365 :
366 2564356 : ABI_SFREE(polynomial_term%atindx)
367 2564356 : ABI_SFREE(polynomial_term%cell)
368 2564356 : ABI_SFREE(polynomial_term%direction)
369 2564356 : ABI_SFREE(polynomial_term%power_disp)
370 2564356 : ABI_SFREE(polynomial_term%power_strain)
371 2564356 : ABI_SFREE(polynomial_term%strain)
372 2564356 : ABI_SFREE(polynomial_term%index_coeff)
373 :
374 2564356 : end subroutine polynomial_term_free
375 : !!***
376 :
377 : !!****f* m_polynomial_term/polynomial_term_finalizer
378 : !!
379 : !! NAME
380 : !! polynomial_term_finalizer
381 : !!
382 : !! FUNCTION
383 : !! Finalizer procedure for polynomial_term_type to automatically free allocated memory
384 : !!
385 : !! SOURCE
386 :
387 0 : subroutine polynomial_term_finalizer(self)
388 : type(polynomial_term_type), intent(inout) :: self
389 : !print *, "Warning: polynomial_term_finalizer called. Debug str: ", self%debug_str
390 0 : call polynomial_term_free(self)
391 0 : end subroutine polynomial_term_finalizer
392 : !!***
393 :
394 : ! function polynomial_term_type_get_index_coeff(term, ndisp, nstrain) result(list)
395 : ! type(polynomial_term_type), intent(inout) :: term
396 : ! integer, intent(in) :: size
397 : ! integer :: list()
398 : ! integer :: i, counter, ip
399 : ! counter=1
400 : ! do i=1, term%ndisp
401 : ! do ip=1, term%power_disp
402 : ! list(counter) = term%
403 : ! end do
404 : ! end do
405 : ! end function polynomial_term_type_get_index_coeff
406 :
407 6 : subroutine polynomial_term_list_free(terms)
408 : type(polynomial_term_type), allocatable, intent(inout) :: terms(:)
409 : integer :: iterm
410 294 : do iterm=1, size(terms)
411 294 : call polynomial_term_free(terms(iterm))
412 : end do
413 294 : ABI_SFREE(terms)
414 6 : end subroutine polynomial_term_list_free
415 :
416 : !!****f* m_polynomial_term/terms_compare
417 : !! NAME
418 : !! equal
419 : !!
420 : !! FUNCTION
421 : !! Compare two polynomial_term_dot
422 : !!
423 : !! INPUTS
424 : !! t1<type(polynomial_term)> = datatype of the first term
425 : !! t2<type(polynomial_term)> = datatype of the second term
426 : !!
427 : !! OUTPUT
428 : !! res = logical
429 : !!
430 : !! SOURCE
431 :
432 14559351 : pure function terms_compare(t1,t2) result (res)
433 : !Arguments ------------------------------------
434 : implicit none
435 :
436 : !Arguments ------------------------------------
437 : type(polynomial_term_type), intent(in) :: t1,t2
438 : logical :: res
439 : !local
440 : !variable
441 : integer :: ia,idisp1,idisp2,mu
442 : logical :: found
443 : !array
444 29118702 : integer :: blkval(2,t1%ndisp+t1%nstrain)
445 : ! *************************************************************************
446 14559351 : res = .true.
447 137839533 : blkval(:,:) = 0
448 14559351 : if(t1%ndisp==t2%ndisp.and.t1%nstrain==t2%nstrain)then
449 : ! Check strain
450 137826825 : blkval(:,:) = 0
451 19000203 : do idisp1=1,t1%nstrain
452 4442301 : if(blkval(1,t1%ndisp+idisp1)==1)cycle!already found
453 23631924 : do idisp2=1,t2%nstrain
454 4631721 : if(blkval(2,t1%ndisp+idisp2)==1)cycle!already found
455 4588037 : found = .false.
456 4588037 : if(t1%strain(idisp1) == t2%strain(idisp2).and.&
457 : & t1%power_strain(idisp1) == t2%power_strain(idisp2))then
458 : found=.true.
459 : end if
460 4442301 : if(found)then
461 1454067 : blkval(1,t1%ndisp+idisp1) = 1
462 1454067 : blkval(2,t1%ndisp+idisp2) = 1
463 : end if
464 : end do
465 : end do
466 18854185 : if(any(blkval(:,t1%ndisp+1:t1%ndisp+t1%nstrain) == 0))then
467 14441095 : res = .false.
468 : return
469 : end if
470 : ! Check displacement
471 42396690 : do idisp1=1,t1%ndisp
472 30795972 : if(blkval(1,idisp1)==1)cycle!already found
473 128425608 : do idisp2=1,t2%ndisp
474 86028918 : if(blkval(2,idisp2)==1)cycle!already found
475 83572747 : found = .false.
476 : if(t1%atindx(1,idisp1) == t2%atindx(1,idisp2).and.&
477 : & t1%atindx(2,idisp1) == t2%atindx(2,idisp2).and.&
478 83572747 : & t1%direction(idisp1) == t2%direction(idisp2).and.&
479 30795972 : & t1%power_disp(idisp1) == t2%power_disp(idisp2))then
480 : found=.true.
481 27581955 : do ia=1,2
482 82745865 : do mu=1,3
483 73551880 : if(t1%cell(mu,ia,idisp1) /= t2%cell(mu,ia,idisp2))then
484 8930163 : found = .false.
485 8930163 : cycle
486 : end if
487 : end do
488 : end do
489 9193985 : if(found)then
490 2294476 : blkval(1,idisp1) = 1
491 2294476 : blkval(2,idisp2) = 1
492 : end if
493 : end if
494 : end do
495 : end do
496 18568671 : if(any(blkval(:,:)==0))res = .false.
497 : else
498 : res = .false.
499 : end if
500 : end function terms_compare
501 : !!***
502 :
503 :
504 0 : function terms_compare_inverse(t1, t2) result(res)
505 : type(polynomial_term_type), intent(in) :: t1,t2
506 : logical :: res
507 0 : type(polynomial_term_type) :: t3
508 0 : call polynomial_term_copy(t2, t3)
509 0 : t3%atindx(1,:) =t2%atindx(2, :)
510 0 : t3%atindx(2,:) =t2%atindx(1, :)
511 0 : t3%cell(1,:, :) = 0
512 0 : t3%cell(2, :, :) = -t2%cell(2,:, : )
513 0 : res=terms_compare(t1, t3)
514 0 : call polynomial_term_free(t3)
515 0 : end function terms_compare_inverse
516 :
517 :
518 862178 : subroutine polynomial_term_copy(in, out)
519 : type(polynomial_term_type), intent(in) :: in
520 : type(polynomial_term_type), intent(out) :: out
521 862178 : if (in%nindex>-1) then
522 : call polynomial_term_init(in%atindx,in%cell,in%direction,in%ndisp,&
523 : & in%nstrain,out,in%power_disp,&
524 : & in%power_strain,in%strain,in%weight, &
525 765146 : & check=.True., index_coeff=in%index_coeff)
526 : else
527 : call polynomial_term_init(in%atindx,in%cell,in%direction,in%ndisp,&
528 : & in%nstrain,out,in%power_disp,&
529 : & in%power_strain,in%strain,in%weight, &
530 97032 : & check=.True.)
531 : endif
532 862178 : end subroutine polynomial_term_copy
533 :
534 0 : end module m_polynomial_term
535 : !!***
|
