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 585 : 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 : !!***
|
