Branch data Line data Source code
1 : : !!****f* ABINIT/m_anharmonics_terms
2 : : !!
3 : : !! NAME
4 : : !! m_anharmonics_term
5 : : !!
6 : : !! FUNCTION
7 : : !! Module with datatype and tools for the anharmonics 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_anharmonics_terms
26 : :
27 : : use defs_basis
28 : : use m_errors
29 : : use m_abicore
30 : : use m_polynomial_coeff
31 : : use m_ifc, only : ifc_type
32 : : use m_supercell, only: getPBCIndexes_supercell
33 : : use m_xmpi
34 : :
35 : : implicit none
36 : :
37 : : public :: anharmonics_terms_init
38 : : public :: anharmonics_terms_free
39 : : public :: anharmonics_terms_freeCoeffs
40 : : public :: anharmonics_terms_evaluateElastic
41 : : public :: anharmonics_terms_evaluateIFCStrainCoupling
42 : : public :: anharmonics_terms_setCoeffs
43 : : public :: anharmonics_terms_setElastic3rd
44 : : public :: anharmonics_terms_setElastic4th
45 : : public :: anharmonics_terms_setElasticDispCoupling
46 : : public :: anharmonics_terms_setStrainPhononCoupling
47 : :
48 : : !!***
49 : :
50 : : !!****t* defs_abitypes/anharmonics_terms_type
51 : : !! NAME
52 : : !! anharmonics_terms_type
53 : : !!
54 : : !! FUNCTION
55 : : !! datatype for a effective potential constructed.
56 : : !!
57 : : !! SOURCE
58 : :
59 : : type, public :: anharmonics_terms_type
60 : :
61 : : integer :: ncoeff = 0
62 : : ! nterm store the number of coefficients
63 : :
64 : : logical :: has_elastic3rd
65 : : ! Flag to know if the 3rd derivatives with respect to strain is present
66 : :
67 : : logical :: has_elastic4th
68 : : ! Flag to know if the 3rd derivatives with respect to strain is present
69 : :
70 : : logical :: has_strain_coupling
71 : : ! Flag to know if the 3rd derivatives with respect to strain and 2 atom disp is present
72 : :
73 : : logical :: has_elastic_displ
74 : : ! Flag to know if the 3rd derivatives with respect to 2 strain and 3 atom disp is present
75 : :
76 : : logical :: bounded
77 : : ! True : the model is bounded
78 : :
79 : : type(polynomial_coeff_type),dimension(:),allocatable :: coefficients
80 : : ! array with all the coefficients from polynomial coefficients
81 : :
82 : : real(dp) :: elastic3rd(6,6,6)
83 : : ! elastic_constant(6,6,6)
84 : : ! Elastic tensor Hartree
85 : :
86 : : real(dp) :: elastic4th(6,6,6,6)
87 : : ! elastic_constant(6,6,6)
88 : : ! Elastic tensor Hartree
89 : :
90 : : real(dp), allocatable :: elastic_displacement(:,:,:,:)
91 : : ! elastic_displacement(6,6,3,natom)
92 : : ! internal strain tensor
93 : :
94 : : type(ifc_type),dimension(:),allocatable :: phonon_strain
95 : : ! Array of ifc with phonon_strain coupling for each strain
96 : :
97 : : end type anharmonics_terms_type
98 : : !!***
99 : :
100 : : CONTAINS !===========================================================================================
101 : :
102 : :
103 : : !!****f* m_anharmonics_terms/anharmonics_terms_init
104 : : !!
105 : : !! NAME
106 : : !! anharmonics_terms_init
107 : : !!
108 : : !! FUNCTION
109 : : !! Initialize anharmonics_terms datatype
110 : : !!
111 : : !! INPUTS
112 : : !! natom = number of atoms in primitive cell
113 : : !! ncoeff = number of coefficient for the fited polynome
114 : : !! bounded = optional, flag to now if the model in bounded
115 : : !! elastic3rd(6,6,6) = optional,3rd order of the elastic constants
116 : : !! elastic4th(6,6,6,6) = optional,4st order of the elastic constants
117 : : !! elastic_displacement(6,6,3,natom) = optional,elastic constant - force coupling
118 : : !! phonon_strain<type(ifc_type)>(6) = optional,phonon strain couling
119 : : !! coeffs<type(polynomial_coeff_type)>(ncoeff) = optional,datatype with polynomial coefficients
120 : : !!
121 : : !! OUTPUT
122 : : !! anharmonics_terms<type(anharmonics_terms_type)> = anharmonics_terms datatype to be initialized
123 : : !!
124 : : !! SOURCE
125 : :
126 [ - + - - : 76 : subroutine anharmonics_terms_init(anharmonics_terms,natom,ncoeff,&
- - - - -
- - - - -
- - - - -
- - - - +
- + - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - ]
127 : : & bounded,elastic3rd,elastic4th,elastic_displacement,&
128 : 0 : & phonon_strain,coeffs)
129 : :
130 : : !Arguments ------------------------------------
131 : : !scalars
132 : : integer, intent(in) :: natom,ncoeff
133 : : type(anharmonics_terms_type), intent(out) :: anharmonics_terms
134 : : real(dp),optional,intent(in) :: elastic_displacement(6,6,3,natom)
135 : : real(dp),optional,intent(in) :: elastic3rd(6,6,6),elastic4th(6,6,6,6)
136 : : type(polynomial_coeff_type),optional :: coeffs(ncoeff)
137 : : type(ifc_type),optional,intent(in) :: phonon_strain(6)
138 : : logical,optional,intent(in) :: bounded
139 : : !arrays
140 : : !Local variables-------------------------------
141 : : !scalar
142 : : !arrays
143 : : character(len=500) :: msg
144 : :
145 : : ! *************************************************************************
146 : :
147 : 76 : call anharmonics_terms_free(anharmonics_terms)
148 : :
149 : : ! Check the number of atoms
150 [ - + ]: 76 : if (natom < 1) then
151 : : write(msg, '(a,a,a,i10,a)' )&
152 : 0 : & 'The cell must have at least one atom.',ch10,&
153 : 0 : & 'The number of atom is ',natom,'.'
154 : 0 : ABI_BUG(msg)
155 : : end if
156 : :
157 : : !Allocation of phonon strain coupling array (3rd order)
158 [ - + ]: 76 : if(present(phonon_strain)) then
159 : 0 : call anharmonics_terms_setStrainPhononCoupling(anharmonics_terms,natom,phonon_strain)
160 : : end if
161 : :
162 : : !Set the 3rd order elastic tensor
163 [ + + + + : 19684 : anharmonics_terms%elastic3rd = zero
+ + ]
164 : 76 : anharmonics_terms%has_elastic3rd = .FALSE.
165 [ - + ]: 76 : if(present(elastic3rd))then
166 : 0 : call anharmonics_terms_setElastic3rd(anharmonics_terms,elastic3rd)
167 : : end if
168 : :
169 : : !Set the 3rd order elastic tensor
170 [ + + + + : 118180 : anharmonics_terms%elastic4th = zero
+ + + + ]
171 : 76 : anharmonics_terms%has_elastic4th = .FALSE.
172 [ - + ]: 76 : if(present(elastic4th))then
173 : 0 : call anharmonics_terms_setElastic4th(anharmonics_terms,elastic4th)
174 : : end if
175 : :
176 : : !Allocation of 3rd order with respecto to 2 strain and 1 atomic displacement
177 [ - + ]: 76 : if(present(elastic_displacement))then
178 : 0 : call anharmonics_terms_setElasticDispCoupling(anharmonics_terms,natom,elastic_displacement)
179 : : end if
180 : :
181 : 76 : anharmonics_terms%ncoeff = 0
182 : :
183 : : !Allocation of the coefficient
184 [ - + ]: 76 : if(present(coeffs))then
185 [ # # ]: 0 : if(ncoeff /= size(coeffs))then
186 : : write(msg, '(a)' )&
187 : 0 : & ' ncoeff has not the same size than coeffs array, '
188 : 0 : ABI_BUG(msg)
189 : : end if
190 : 0 : call anharmonics_terms_setCoeffs(coeffs,anharmonics_terms,ncoeff)
191 : : end if
192 : :
193 : : !Set the flag bounded
194 [ - + ]: 76 : if(present(bounded))then
195 : 0 : anharmonics_terms%bounded = bounded
196 : : else
197 : 76 : anharmonics_terms%bounded = .FALSE.
198 : : end if
199 : :
200 [ - + ]: 76 : end subroutine anharmonics_terms_init
201 : : !!***
202 : :
203 : : !****f* m_anharmonics_terms/anharmonics_terms_free
204 : : !!
205 : : !! NAME
206 : : !! anharmonics_terms_free
207 : : !!
208 : : !! FUNCTION
209 : : !! deallocate all dynamic memory for this anharmonics_terms datatype
210 : : !!
211 : : !! INPUTS
212 : : !! anharmonics_terms<type(anharmonics_terms_type)> = anharmonics_terms datatype to be free
213 : : !!
214 : : !! OUTPUT
215 : : !! anharmonics_terms<type(anharmonics_terms_type)> = anharmonics_terms datatype to be free
216 : : !!
217 : : !! SOURCE
218 : :
219 : 610 : subroutine anharmonics_terms_free(anharmonics_terms)
220 : :
221 : : !Arguments ------------------------------------
222 : : !scalars
223 : : !array
224 : : type(anharmonics_terms_type), intent(inout) :: anharmonics_terms
225 : :
226 : : !Local variables-------------------------------
227 : : !scalars
228 : : integer :: ii
229 : : !array
230 : :
231 : : ! *************************************************************************
232 : :
233 : 610 : anharmonics_terms%has_elastic3rd = .FALSE.
234 : 610 : anharmonics_terms%has_strain_coupling = .FALSE.
235 : 610 : anharmonics_terms%has_elastic_displ = .FALSE.
236 : 610 : anharmonics_terms%bounded = .FALSE.
237 : :
238 [ - + ]: 610 : if(allocated(anharmonics_terms%elastic_displacement)) then
239 [ # # # # : 0 : anharmonics_terms%elastic_displacement=zero
# # # # ]
240 : 0 : ABI_SFREE(anharmonics_terms%elastic_displacement)
241 : : end if
242 : :
243 [ + + ]: 610 : if(allocated(anharmonics_terms%phonon_strain))then
244 [ + + ]: 315 : do ii = 1,6
245 : 315 : call anharmonics_terms%phonon_strain(ii)%free()
246 : : end do
247 [ + - + + : 315 : ABI_SFREE(anharmonics_terms%phonon_strain)
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + -
+ ]
248 : : end if
249 : :
250 : 610 : call anharmonics_terms_freeCoeffs(anharmonics_terms)
251 : :
252 [ + + + + : 157990 : anharmonics_terms%elastic3rd = zero
+ + ]
253 : :
254 : :
255 : 610 : end subroutine anharmonics_terms_free
256 : : !!***
257 : :
258 : : !****f* m_anharmonics_terms/anharmonics_terms_freeCoeffs
259 : : !!
260 : : !! NAME
261 : : !! anharmonics_terms_freeCoeffs
262 : : !!
263 : : !! FUNCTION
264 : : !! deallocate all dynamic memory for the coefficients
265 : : !! of this anharmonics_terms datatype
266 : : !!
267 : : !! INPUTS
268 : : !! anharmonics_terms<type(anharmonics_terms_type)> = anharmonics_terms datatype to be free
269 : : !!
270 : : !! OUTPUT
271 : : !! anharmonics_terms<type(anharmonics_terms_type)> = anharmonics_terms datatype to be free
272 : : !!
273 : : !!
274 : : !! SOURCE
275 : :
276 : 686 : subroutine anharmonics_terms_freeCoeffs(anharmonics_terms)
277 : :
278 : : !Arguments ------------------------------------
279 : : !scalars
280 : : !array
281 : : type(anharmonics_terms_type), intent(inout) :: anharmonics_terms
282 : : !Local variables-------------------------------
283 : : !scalars
284 : : integer :: ii
285 : : !array
286 : :
287 : : ! *************************************************************************
288 : :
289 [ + + ]: 686 : if(allocated(anharmonics_terms%coefficients))then
290 [ + + ]: 2146 : do ii=1,anharmonics_terms%ncoeff
291 : 2146 : call polynomial_coeff_free(anharmonics_terms%coefficients(ii))
292 : : end do
293 [ + - + + : 2146 : ABI_SFREE(anharmonics_terms%coefficients)
- + - - -
- - - - -
- - - - -
- - - -
+ ]
294 : : end if
295 : :
296 : 686 : anharmonics_terms%ncoeff = 0
297 : :
298 : 686 : end subroutine anharmonics_terms_freeCoeffs
299 : : !!***
300 : :
301 : : !****f* m_anharmonics_terms/anharmonics_terms_setCoeffs
302 : : !!
303 : : !! NAME
304 : : !! anharmonics_terms_setCoeffs
305 : : !!
306 : : !! FUNCTION
307 : : !! Set the coefficients
308 : : !!
309 : : !! INPUTS
310 : : !! coeffs(ncoeff)<type(polynomial_coeff_type)> = array with datatype polynomial_coeff_type
311 : : !! ncoeff = number of coefficient
312 : : !!
313 : : !! OUTPUT
314 : : !! anharmonics_terms<type(anharmonics_terms_type)> = anharmonics_terms datatype
315 : : !!
316 : : !!
317 : : !! SOURCE
318 : :
319 : 192 : subroutine anharmonics_terms_setCoeffs(coeffs,anharmonics_terms,ncoeff)
320 : :
321 : : use m_polynomial_coeff
322 : :
323 : : !Arguments ------------------------------------
324 : : !scalars
325 : : integer,intent(in) :: ncoeff
326 : : !array
327 : : type(anharmonics_terms_type),intent(inout) :: anharmonics_terms
328 : : type(polynomial_coeff_type),intent(in) :: coeffs(ncoeff)
329 : : !Local variables-------------------------------
330 : : !scalar
331 : : integer :: ii
332 : : character(len=500) :: msg
333 : : !array
334 : : ! *************************************************************************
335 : :
336 [ - + ]: 96 : if(ncoeff /= size(coeffs))then
337 : : write(msg, '(a)' )&
338 : 0 : & ' ncoeff has not the same size than coeffs array, '
339 : 0 : ABI_BUG(msg)
340 : : end if
341 : :
342 : : ! 1-deallocation of the previous value
343 [ + + ]: 96 : if(allocated(anharmonics_terms%coefficients))then
344 [ + + ]: 9 : do ii=1,anharmonics_terms%ncoeff
345 : 9 : call polynomial_coeff_free(anharmonics_terms%coefficients(ii))
346 : : end do
347 [ + - + + : 9 : ABI_SFREE(anharmonics_terms%coefficients)
- + - - -
- - - - -
- - - - -
- - - -
+ ]
348 : : end if
349 : :
350 : : ! Allocation of the new array
351 : 96 : anharmonics_terms%ncoeff = ncoeff
352 [ - + + - : 2347 : ABI_MALLOC(anharmonics_terms%coefficients,(ncoeff))
+ - + - -
+ - + - +
+ + + - ]
353 [ + + ]: 2155 : do ii=1,anharmonics_terms%ncoeff
354 : : call polynomial_coeff_init(coeffs(ii)%coefficient,coeffs(ii)%nterm,&
355 : : & anharmonics_terms%coefficients(ii),&
356 : : & coeffs(ii)%terms, &
357 : 2155 : & name=coeffs(ii)%name)
358 : : end do
359 : :
360 : 96 : end subroutine anharmonics_terms_setCoeffs
361 : : !!***
362 : :
363 : : !****f* m_anharmonics_terms/anharmonics_terms_setElastic3rd
364 : : !!
365 : : !! NAME
366 : : !! anharmonics_terms_setElastic3rd
367 : : !!
368 : : !! FUNCTION
369 : : !! Set the 3rd order derivative of with respect to 3 strain
370 : : !!
371 : : !! INPUTS
372 : : !! elastics(6,6,6) = 3d order of elastics constant
373 : : !!
374 : : !! OUTPUT
375 : : !! anharmonics_terms<type(anharmonics_terms_type)> = anharmonics_terms datatype
376 : : !!
377 : : !! SOURCE
378 : :
379 : 45 : subroutine anharmonics_terms_setElastic3rd(anharmonics_terms,elastics)
380 : :
381 : : !Arguments ------------------------------------
382 : : !scalars
383 : : !array
384 : : type(anharmonics_terms_type),intent(inout) :: anharmonics_terms
385 : : real(dp),intent(in) :: elastics(6,6,6)
386 : : !Local variables-------------------------------
387 : : !scalar
388 : : !array
389 : : ! *************************************************************************
390 : :
391 : : ! 1-reinitialise the previous value
392 [ + + + + : 11655 : anharmonics_terms%elastic3rd(:,:,:) = zero
+ + ]
393 : 45 : anharmonics_terms%has_elastic3rd = .FALSE.
394 : :
395 : : ! 2-Allocation of the new array
396 [ + + + + : 11655 : anharmonics_terms%elastic3rd(:,:,:) = elastics(:,:,:)
+ + ]
397 : :
398 : : ! 3-Set the flag
399 [ + + + + : 11655 : if(any(abs(anharmonics_terms%elastic3rd)> tol15)) then
+ + + - -
+ ]
400 : 0 : anharmonics_terms%has_elastic3rd = .TRUE.
401 : : end if
402 : :
403 : 96 : end subroutine anharmonics_terms_setElastic3rd
404 : : !!***
405 : :
406 : : !****f* m_anharmonics_terms/anharmonics_terms_setElastic4th
407 : : !!
408 : : !! NAME
409 : : !! anharmonics_terms_setElastic4th
410 : : !!
411 : : !! FUNCTION
412 : : !! Set the 4th order derivative of with respect to 4 strain
413 : : !!
414 : : !! INPUTS
415 : : !! elastics = 4th order of elastics constant
416 : : !!
417 : : !! OUTPUT
418 : : !! anharmonics_terms<type(anharmonics_terms_type)> = anharmonics_terms datatype
419 : : !!
420 : : !!
421 : : !! SOURCE
422 : :
423 : 0 : subroutine anharmonics_terms_setElastic4th(anharmonics_terms,elastics)
424 : :
425 : : !Arguments ------------------------------------
426 : : !scalars
427 : : !array
428 : : type(anharmonics_terms_type),intent(inout) :: anharmonics_terms
429 : : real(dp),intent(in) :: elastics(6,6,6,6)
430 : : !Local variables-------------------------------
431 : : !scalar
432 : : !array
433 : : ! *************************************************************************
434 : :
435 : : ! 1-reinitialise the previous value
436 [ # # # # : 0 : anharmonics_terms%elastic4th(:,:,:,:) = zero
# # # # ]
437 : 0 : anharmonics_terms%has_elastic4th = .FALSE.
438 : :
439 : : ! 2-Allocation of the new array
440 [ # # # # : 0 : anharmonics_terms%elastic4th(:,:,:,:) = elastics(:,:,:,:)
# # # # ]
441 : :
442 : : ! 3-Set the flag
443 [ # # # # : 0 : if(any(abs(anharmonics_terms%elastic4th)> tol15)) then
# # # # #
# # # ]
444 : 0 : anharmonics_terms%has_elastic4th = .TRUE.
445 : : end if
446 : :
447 : 0 : end subroutine anharmonics_terms_setElastic4th
448 : : !!***
449 : :
450 : :
451 : : !****f* m_anharmonics_terms/anharmonics_terms_setStrainPhononCoupling
452 : : !!
453 : : !! NAME
454 : : !! anharmonics_terms_setStrainPhononCoupling
455 : : !!
456 : : !! FUNCTION
457 : : !! Set the strain-phonon coupling
458 : : !!
459 : : !! INPUTS
460 : : !! strain_phonon(6)<type(ifc_type) = strain-phonon coupling
461 : : !! natom = number of atoms
462 : : !!
463 : : !! OUTPUT
464 : : !! anharmonics_terms<type(anharmonics_terms_type)> = anharmonics_terms datatype
465 : : !!
466 : : !! SOURCE
467 : :
468 : 45 : subroutine anharmonics_terms_setStrainPhononCoupling(anharmonics_terms,natom,phonon_strain)
469 : :
470 : : !Arguments ------------------------------------
471 : : !scalars
472 : : integer, intent(in) :: natom
473 : : !array
474 : : type(anharmonics_terms_type),intent(inout) :: anharmonics_terms
475 : : type(ifc_type),intent(in) :: phonon_strain(6)
476 : : !Local variables-------------------------------
477 : : !scalar
478 : : integer :: ii,nrpt
479 : : character(500) :: msg
480 : : !array
481 : : ! *************************************************************************
482 : :
483 : : ! 1-Do some check
484 [ + + ]: 315 : do ii=1,6
485 : : !PROBLEM eos_gnu_13.2_openmpi
486 : : ! if(natom /= size(phonon_strain(ii)%atmfrc,2).or.&
487 [ - + ]: 270 : if( &
488 : : !ENDPROBLEM
489 : 45 : & phonon_strain(ii)%nrpt < 0)then
490 : : write(msg, '(a)' )&
491 : 0 : & ' natom or/and nrpt have not the same size than phonon_strain array. '
492 : 0 : ABI_BUG(msg)
493 : : end if
494 : : end do
495 : :
496 : : ! 1-reinitialise the previous value
497 : 45 : anharmonics_terms%has_strain_coupling = .FALSE.
498 [ - + ]: 45 : if(allocated(anharmonics_terms%phonon_strain))then
499 [ # # ]: 0 : do ii = 1,6
500 : 0 : call anharmonics_terms%phonon_strain(ii)%free()
501 : : end do
502 [ # # # # : 0 : ABI_SFREE(anharmonics_terms%phonon_strain)
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # #
# ]
503 : : end if
504 : :
505 : : ! 2-Allocation of the new array and filling
506 [ - + - + : 990 : ABI_MALLOC(anharmonics_terms%phonon_strain,(6))
+ + + + +
+ + + +
- ]
507 [ + + ]: 315 : do ii = 1,6
508 : 270 : nrpt = phonon_strain(ii)%nrpt
509 [ - + + - : 1620 : ABI_MALLOC(anharmonics_terms%phonon_strain(ii)%atmfrc,(3,natom,3,natom,nrpt))
+ - - + +
- + - - -
+ - - + -
+ - + -
+ ]
510 [ + - - - : 810 : ABI_MALLOC(anharmonics_terms%phonon_strain(ii)%cell,(3,nrpt))
+ - - + -
+ - + -
+ ]
511 : 270 : anharmonics_terms%phonon_strain(ii)%nrpt = phonon_strain(ii)%nrpt
512 [ - + - - : 270 : anharmonics_terms%phonon_strain(ii)%atmfrc(:,:,:,:,:) = phonon_strain(ii)%atmfrc(:,:,:,:,:)
- - - - -
- ]
513 [ - + - - ]: 270 : anharmonics_terms%phonon_strain(ii)%cell(:,:) = phonon_strain(ii)%cell(:,:)
514 : :
515 : : !3-Set the flag
516 [ - + - - : 315 : if(any(abs(anharmonics_terms%phonon_strain(ii)%atmfrc)> tol15)) then
- - - - -
- - - -
+ ]
517 : 0 : anharmonics_terms%has_strain_coupling = .TRUE.
518 : : ! If there is no value inside the array,
519 : : ! We don't need to store it
520 : : else
521 [ + - ]: 270 : ABI_SFREE(anharmonics_terms%phonon_strain(ii)%atmfrc)
522 [ + - ]: 270 : ABI_SFREE(anharmonics_terms%phonon_strain(ii)%cell)
523 : 270 : anharmonics_terms%phonon_strain(ii)%nrpt = 0
524 : : end if
525 : : end do
526 : :
527 : 45 : end subroutine anharmonics_terms_setStrainPhononCoupling
528 : : !!***
529 : :
530 : : !****f* m_anharmonics_terms/anharmonics_terms_setElasticDispCoupling
531 : : !!
532 : : !! NAME
533 : : !! anharmonics_terms_setElasticDispCoupling
534 : : !!
535 : : !! FUNCTION
536 : : !! Set the Elastic displacement coupling
537 : : !!
538 : : !! INPUTS
539 : : !! elastic_displacement(6,6,3,natom) = elastic displacement coupling
540 : : !! natom = number of atom
541 : : !!
542 : : !! OUTPUT
543 : : !! anharmonics_terms<type(anharmonics_terms_type)> = anharmonics_terms datatype
544 : : !!
545 : : !!
546 : : !! SOURCE
547 : :
548 : 38 : subroutine anharmonics_terms_setElasticDispCoupling(anharmonics_terms,natom,elastic_displacement)
549 : :
550 : : !Arguments ------------------------------------
551 : : !scalars
552 : : integer, intent(in) :: natom
553 : : !array
554 : : type(anharmonics_terms_type),intent(inout) :: anharmonics_terms
555 : : real(dp),intent(in) :: elastic_displacement(6,6,3,natom)
556 : : !Local variables-------------------------------
557 : : !scalar
558 : : character(500) :: msg
559 : : !array
560 : : ! *************************************************************************
561 : :
562 : : ! 1-Do some check
563 [ - + ]: 38 : if(natom /= size(elastic_displacement,4)) then
564 : : write(msg, '(a)' )&
565 : 0 : & ' natom has not the same size elastic_displacement array. '
566 : 0 : ABI_BUG(msg)
567 : : end if
568 : :
569 : : ! 1-reinitialise the previous value
570 : 38 : anharmonics_terms%has_elastic_displ = .FALSE.
571 [ - + ]: 38 : if(allocated(anharmonics_terms%elastic_displacement))then
572 : 0 : ABI_SFREE(anharmonics_terms%elastic_displacement)
573 : : end if
574 : :
575 : : ! 2-Allocation of the new array and filling
576 [ - + + - : 114 : ABI_MALLOC(anharmonics_terms%elastic_displacement,(6,6,3,natom))
+ - + - -
+ - + -
+ ]
577 [ + + + + : 26558 : anharmonics_terms%elastic_displacement(:,:,:,:) = elastic_displacement(:,:,:,:)
+ + + + ]
578 : :
579 : : ! 3-Set the flag
580 [ + + + + : 26558 : if(any(abs(anharmonics_terms%elastic_displacement)> tol15)) then
+ + + + +
- - + ]
581 : 0 : anharmonics_terms%has_elastic_displ = .TRUE.
582 : : else
583 : : ! If there is no value inside the array,
584 : : ! We don't need to store it
585 [ + - ]: 38 : ABI_SFREE(anharmonics_terms%elastic_displacement)
586 : : end if
587 : :
588 : 38 : end subroutine anharmonics_terms_setElasticDispCoupling
589 : : !!***
590 : :
591 : : !!****f* m_effective_potential/anharmonics_terms_evaluateElastic
592 : : !! NAME
593 : : !! anharmonics_terms_evaluateElastic
594 : : !!
595 : : !! FUNCTION
596 : : !! Compute the energy, stresses and forces related to the application of strain
597 : : !!
598 : : !! INPUTS
599 : : !! disp(3,natom_sc) = atomics displacement between configuration and the reference
600 : : !! natom = number of atom in the supercell
601 : : !! natom_uc = number of atom in the unit cell
602 : : !! ncell = number of cell
603 : : !! strain(6) = strain to apply
604 : : !! elastic3rd(6,6,6) = 3 order derivatives with respect to to 3 strain
605 : : !! elastic4th(6,6,66,) = 4 order derivatives with respect to to 4 strain
606 : : !! elastic_displacement(6,6,3,natom) = 3 order derivatives with respect to 2 strain and 1 Atom disp
607 : : !!
608 : : !! OUTPUT
609 : : !! energy = contribution of the ifc to the energy
610 : : !! fcart(3,natom) = contribution of the ifc to the forces
611 : : !! strten(6) = contribution to the stress tensor
612 : : !!
613 : : !! SOURCE
614 : : !!
615 : 0 : subroutine anharmonics_terms_evaluateElastic(disp,energy,fcart,natom,natom_uc,ncell,strten,strain,&
616 : : & elastic3rd,elastic4th,elastic_displacement)
617 : :
618 : : real(dp),intent(out):: energy
619 : : integer, intent(in) :: natom,natom_uc,ncell
620 : : ! array
621 : : real(dp),optional,intent(in) :: elastic3rd(6,6,6),elastic4th(6,6,6,6)
622 : : real(dp),optional,intent(in) :: elastic_displacement(6,6,3,natom)
623 : : real(dp),intent(out):: strten(6)
624 : : real(dp),intent(out):: fcart(3,natom)
625 : : real(dp),intent(in) :: disp(3,natom)
626 : : real(dp),intent(in) :: strain(6)
627 : :
628 : : !Local variables-------------------------------
629 : : ! scalar
630 : : integer :: ia,ii,mu,alpha,beta,gamma,delta,d1,d2
631 : : real(dp):: cijk
632 : : logical :: has_elastic3rd,has_elastic4th,has_elastic_displ
633 : : ! array
634 : : ! *************************************************************************
635 : :
636 : : !Reset output and flags
637 : 0 : energy = zero
638 [ # # # # ]: 0 : fcart = zero
639 : 0 : strten = zero
640 : 0 : has_elastic3rd = .FALSE.
641 : 0 : has_elastic4th = .FALSE.
642 : 0 : has_elastic_displ = .FALSE.
643 : 0 : d1=0;d2=0
644 : :
645 : : !Set the flags
646 [ # # ]: 0 : if(present(elastic3rd)) has_elastic3rd = .TRUE.
647 [ # # ]: 0 : if(present(elastic4th)) then
648 : 0 : has_elastic4th = .TRUE.
649 : 0 : d1=1;d2=6
650 : : end if
651 [ # # ]: 0 : if(present(elastic_displacement)) has_elastic_displ = .TRUE.
652 : :
653 : : !1-Treat 3rd order elastic constants
654 [ # # ]: 0 : if (has_elastic3rd.or.has_elastic4th) then
655 [ # # ]: 0 : do alpha=1,6
656 [ # # ]: 0 : do beta=1,6
657 [ # # ]: 0 : do gamma=1,6
658 : 0 : cijk = ncell*elastic3rd(alpha,beta,gamma)
659 : : ! Accumulate energy
660 : 0 : energy = energy + sixth*cijk*strain(alpha)*strain(beta)*strain(gamma)
661 : : ! Accumulate stresses contributions
662 : 0 : strten(alpha)=strten(alpha)+ half*cijk*strain(beta)*strain(gamma)
663 [ # # ]: 0 : do delta=d1,d2
664 : 0 : cijk = ncell*elastic4th(alpha,beta,gamma,delta)
665 : : ! Accumulate energy
666 : : energy = energy + (1/24.)*cijk*strain(alpha)*strain(beta)*&
667 : 0 : & strain(gamma)*strain(delta)
668 : : ! Accumulate stresses contributions
669 : : strten(alpha)=strten(alpha)+ sixth*cijk*strain(beta)*strain(gamma)*&
670 : 0 : & strain(delta)
671 : : end do
672 : : end do
673 : : end do
674 : : end do
675 : : end if
676 : :
677 : : !2-Part due to the internat strain
678 [ # # ]: 0 : if(has_elastic_displ)then
679 : : ii = 1
680 [ # # ]: 0 : do ia = 1,natom
681 [ # # ]: 0 : do mu = 1,3
682 [ # # ]: 0 : do beta=1,6
683 [ # # ]: 0 : do alpha=1,6
684 : 0 : cijk = elastic_displacement(alpha,beta,mu,ii)
685 : : ! Accumulte for this atom
686 : 0 : energy = energy + sixth*cijk*strain(alpha)*strain(beta)*disp(mu,ia)
687 : 0 : fcart(mu,ia) = fcart(mu,ia) + half*cijk*strain(alpha)*strain(beta)
688 : 0 : strten(alpha) = strten(alpha) + half*cijk*strain(beta)*disp(mu,ia)
689 : : end do
690 : : end do
691 : : end do
692 : 0 : ii = ii +1
693 : : ! Reset to 1 if the number of atoms is superior than in the initial cell
694 [ # # ]: 0 : if(ii==natom_uc+1) ii = 1
695 : : end do
696 : : end if
697 : :
698 : 0 : end subroutine anharmonics_terms_evaluateElastic
699 : : !!***
700 : :
701 : : !!****f* m_anharmonics_terms/anharmonics_terms_evaluateIFCStrainCoupling
702 : : !! NAME
703 : : !! anharmonics_terms_evaluateIFCStrainCoupling
704 : : !!
705 : : !! FUNCTION
706 : : !! This fonction compute the harmonic part of the energy
707 : : !! of the supercell in the eff_pot
708 : : !! INPUTS
709 : : !! strain_phonon(6)<type(ifc_type) = strain-phonon coupling
710 : : !! disp(3,natom_sc) = atomics displacement between configuration and the reference
711 : : !! natom = number of atoms in the supercell
712 : : !! natom_uc = number of atoms in the unit cell
713 : : !! sc_size(3) = size of the supercell
714 : : !! cells(ncell) = number of the cells into the supercell (1,2,3,4,5)
715 : : !! ncell = total number of cell to treat
716 : : !! index_cells(3,ncell) = indexes of the cells into supercell (-1 -1 -1 ,...,1 1 1)
717 : : !! comm=MPI communicator
718 : : !!
719 : : !! OUTPUT
720 : : !! energy = contribution of the ifc to the energy
721 : : !! fcart(3,natom) = contribution of the ifc to the forces
722 : : !! strten(6) = contribution to the stress tensor
723 : : !!
724 : : !! PARENT
725 : : !! effective_potential_evaluate
726 : : !!
727 : : !! SOURCE
728 : :
729 : 0 : subroutine anharmonics_terms_evaluateIFCStrainCoupling(phonon_strain,disp,energy,fcart,natom,natom_uc,&
730 : 0 : & sc_size,strain,strten,cells,ncell,&
731 : 0 : & index_cells,comm)
732 : :
733 : : !Arguments -------------------------------
734 : : ! scalars
735 : : real(dp),intent(out) :: energy
736 : : integer,intent(in) :: natom,natom_uc,ncell
737 : : integer,intent(in) :: comm
738 : : ! array
739 : : integer,intent(in) :: cells(ncell),index_cells(ncell,3)
740 : : integer,intent(in) :: sc_size(3)
741 : : type(ifc_type),intent(in) :: phonon_strain(6)
742 : : real(dp),intent(in) :: disp(3,natom)
743 : : real(dp),intent(out) :: fcart(3,natom)
744 : : real(dp),intent(out) :: strten(6)
745 : : real(dp),intent(in) :: strain(6)
746 : : !Local variables-------------------------------
747 : : ! scalar
748 : : integer :: alpha
749 : : integer :: i1,i2,i3,ia,ib,icell,ii
750 : : integer :: irpt,jj,kk,ll,mu,nu
751 : : integer :: ierr
752 : : real(dp):: ifc
753 : : ! array
754 : : integer :: cell_atom2(3)
755 : : character(500) :: msg
756 : :
757 : : ! *************************************************************************
758 : :
759 [ # # # # : 0 : if (any(sc_size <= 0)) then
# # ]
760 : 0 : write(msg,'(a,a)')' sc_size can not be inferior or equal to zero'
761 : 0 : ABI_ERROR(msg)
762 : : end if
763 : :
764 : : ! Initialisation of variables
765 : 0 : energy = zero
766 [ # # # # ]: 0 : fcart(:,:) = zero
767 : 0 : strten(:) = zero
768 : :
769 [ # # ]: 0 : do icell = 1,ncell
770 : 0 : ii = (cells(icell)-1)*natom_uc
771 : 0 : i1=index_cells(icell,1); i2=index_cells(icell,2); i3=index_cells(icell,3)
772 [ # # ]: 0 : do alpha=1,6
773 [ # # ]: 0 : do irpt = 1,phonon_strain(alpha)%nrpt
774 : : ! get the cell of atom2 (0 0 0, 0 0 1...)
775 : 0 : cell_atom2(1) = i1 + phonon_strain(alpha)%cell(1,irpt)
776 : 0 : cell_atom2(2) = i2 + phonon_strain(alpha)%cell(2,irpt)
777 : 0 : cell_atom2(3) = i3 + phonon_strain(alpha)%cell(3,irpt)
778 : 0 : call getPBCIndexes_supercell(cell_atom2(1:3),sc_size(1:3))
779 : : ! index of the second atom in the displacement array
780 : : jj = (cell_atom2(1)-1)*sc_size(2)*sc_size(3)*natom_uc+&
781 : : & (cell_atom2(2)-1)*sc_size(3)*natom_uc+&
782 : 0 : & (cell_atom2(3)-1)*natom_uc
783 [ # # ]: 0 : do ib = 1, natom_uc
784 : 0 : ll = jj + ib
785 [ # # ]: 0 : do nu=1,3
786 [ # # ]: 0 : do ia = 1, natom_uc
787 : 0 : kk = ii + ia
788 [ # # ]: 0 : do mu=1,3
789 : 0 : ifc = phonon_strain(alpha)%atmfrc(mu,ia,nu,ib,irpt)
790 : : ! accumule energy
791 : 0 : energy = energy + sixth*strain(alpha)*disp(mu,kk)*disp(nu,ll)*ifc
792 : : ! accumule forces
793 : 0 : fcart(mu,kk) = fcart(mu,kk) + half*strain(alpha)*disp(nu,ll)*ifc
794 : : ! accumule stresses
795 : 0 : strten(alpha) = strten(alpha) + half*disp(mu,kk)*disp(nu,ll)*ifc
796 : : end do
797 : : end do
798 : : end do
799 : : end do
800 : : end do
801 : : end do
802 : : end do
803 : :
804 : : ! MPI_SUM
805 : 0 : call xmpi_sum(energy, comm, ierr)
806 : 0 : call xmpi_sum(fcart , comm, ierr)
807 : 0 : call xmpi_sum(strten, comm, ierr)
808 : :
809 : 0 : end subroutine anharmonics_terms_evaluateIFCStrainCoupling
810 : : !!***
811 [ # # # # : 0 : end module m_anharmonics_terms
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # ]
812 : : !!***
|
