Branch data Line data Source code
1 : : !!****f* ABINIT/m_harmonics_terms
2 : : !!
3 : : !! NAME
4 : : !! m_harmonics_term
5 : : !!
6 : : !! FUNCTION
7 : : !! Module with datatype and tools for the harmonics 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_harmonics_terms
26 : :
27 : : use defs_basis
28 : : use m_errors
29 : : use m_abicore
30 : : use m_supercell,only: getPBCIndexes_supercell
31 : : use m_xmpi,only : xmpi_sum
32 : : use m_ifc
33 : :
34 : : implicit none
35 : :
36 : : public :: harmonics_terms_init
37 : : public :: harmonics_terms_free
38 : : public :: harmonics_terms_applySumRule
39 : : public :: harmonics_terms_evaluateIFC
40 : : public :: harmonics_terms_evaluateElastic
41 : : public :: harmonics_terms_setEffectiveCharges
42 : : public :: harmonics_terms_setDynmat
43 : : public :: harmonics_terms_setInternalStrain
44 : : !!***
45 : :
46 : : !!****t* m_harmonics_terms/harmonics_terms_type
47 : : !! NAME
48 : : !! harmonics_terms_type
49 : : !!
50 : : !! FUNCTION
51 : : !! datatype for harmonic part of effective potential.
52 : : !!
53 : : !! SOURCE
54 : :
55 : : type, public :: harmonics_terms_type
56 : :
57 : : integer :: nqpt
58 : : ! Number of qpoints
59 : :
60 : : real(dp) :: epsilon_inf(3,3)
61 : : ! epsilon_inf(3,3)
62 : : ! Dielectric tensor
63 : :
64 : : real(dp) :: elastic_constants(6,6)
65 : : ! elastic_constant(6,6)
66 : : ! Elastic tensor Hartree
67 : :
68 : : real(dp), allocatable :: strain_coupling(:,:,:)
69 : : ! strain_coupling(6,3,natom)
70 : : ! internal strain tensor
71 : :
72 : : real(dp), allocatable :: zeff(:,:,:)
73 : : ! zeff(3,3,natom) Effective charges
74 : :
75 : : type(ifc_type) :: ifcs
76 : : ! type with ifcs constants (short + ewald)
77 : : ! also contains the number of cell and the indexes
78 : :
79 : : real(dp), allocatable :: qpoints(:,:)
80 : : ! qph1l(3,nqpt)
81 : : ! List of qpoints wavevectors
82 : :
83 : : real(dp), allocatable :: dynmat(:,:,:,:,:,:)
84 : : ! dynmat(2,3,natom,3,natom,nqpt)
85 : : ! dynamical matrix for each q points
86 : :
87 : : real(dp), allocatable :: phfrq(:,:)
88 : : ! phfrq(3*natom,nqpt)
89 : : ! array with all phonons frequencies for each q points in Hartree/cm
90 : :
91 : : end type harmonics_terms_type
92 : : !!***
93 : :
94 : : CONTAINS !===========================================================================================
95 : :
96 : :
97 : : !!****f* m_harmonics_terms/harmonics_terms_init
98 : : !!
99 : : !! NAME
100 : : !! harmonics_terms_init
101 : : !!
102 : : !! FUNCTION
103 : : !! Initialize harmonics_terms datatype
104 : : !!
105 : : !! INPUTS
106 : : !! ifc<type(ifc_type)> = interatomic forces constants
107 : : !! natom = number of atoms in primitive cell
108 : : !! nrpt = number rpt (cell) in the ifc
109 : : !! dynmat(2,3,natom,3,natom,3,nqpt) = optional, dynamical matricies for each q-point
110 : : !! epsilon_inf(3,3) = optional, dielectric tensor
111 : : !! elastic_constant(6,6) = optional, elastic constant
112 : : !! strain_coupling(6,3,natom) = optional, internal strain coupling parameters
113 : : !! nqpt = optional, number of q-points
114 : : !! phfrq(3*natom,nqpt) = optional,phonons frequencies for each q points in Hartree/cm
115 : : !! qpoints(3,nqpt) = list of qpoints wavevectors
116 : : !! zeff(3,3,natom) = optional,effective charges
117 : : !!
118 : : !! OUTPUT
119 : : !! harmonics_terms<type(harmonics_terms_type)> = harmonics_terms datatype to be initialized
120 : : !!
121 : : !! SOURCE
122 : :
123 [ - + - + : 608 : subroutine harmonics_terms_init(harmonics_terms,ifcs,natom,nrpt,&
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
+ + + + ]
124 [ - + - - : 76 : & dynmat,epsilon_inf,elastic_constants,strain_coupling,&
- - ]
125 [ - + - + ]: 76 : & nqpt,phfrq,qpoints,zeff)
126 : :
127 : : implicit none
128 : :
129 : : !Arguments ------------------------------------
130 : : !scalars
131 : : integer, intent(in) :: natom,nrpt
132 : : !arrays
133 : : type(ifc_type),intent(in) :: ifcs
134 : : type(harmonics_terms_type), intent(out) :: harmonics_terms
135 : : integer, optional,intent(in) :: nqpt
136 : : real(dp),optional,intent(in) :: epsilon_inf(3,3),dynmat(:,:,:,:,:,:)
137 : : real(dp),optional,intent(in) :: elastic_constants(6,6)
138 : : real(dp),optional,intent(in) :: strain_coupling(6,3,natom),zeff(3,3,natom)
139 : : real(dp),optional,intent(in) :: phfrq(:,:),qpoints(:,:)
140 : : !Local variables-------------------------------
141 : : !scalar
142 : : !arrays
143 : : character(len=500) :: msg
144 : :
145 : : ! *************************************************************************
146 : :
147 : 76 : call harmonics_terms_free(harmonics_terms)
148 : :
149 : : ! Do some Checks
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 [ - + ]: 76 : if (nrpt < 1) then
158 : : write(msg, '(a,a,a,i10,a)' )&
159 : 0 : & 'The cell must have at least one rpt point.',ch10,&
160 : 0 : & 'The number of rpt points is ',nrpt,'.'
161 : 0 : ABI_BUG(msg)
162 : : end if
163 : :
164 [ - + ]: 76 : if (nrpt /= ifcs%nrpt) then
165 : : write(msg, '(3a,i5,a,i5,a)' )&
166 : 0 : & 'nrpt must have the same dimension as ifcs.',ch10,&
167 : 0 : & 'The number of cell is ',nrpt,' instead of ',ifcs%nrpt,'.'
168 : 0 : ABI_BUG(msg)
169 : : end if
170 : :
171 [ + - - - : 76 : if(present(nqpt).and.(.not.present(dynmat).or.&
- - ]
172 : : & .not.present(qpoints) .or.&
173 : : & .not.present(phfrq)))then
174 : : write(msg, '(a)' )&
175 : 0 : & 'nqpt is specified but dynamt,qpoints or phfrq are not.'
176 : 0 : ABI_BUG(msg)
177 : : end if
178 : :
179 [ + - - + ]: 76 : if(.not.present(nqpt).and.(present(dynmat).or.&
180 : : & present(qpoints) .or.&
181 : : & present(phfrq)))then
182 : : write(msg, '(a)' )&
183 : 0 : & ' dynamt,qpoints or phfrq are specified but nqpt is not.'
184 : 0 : ABI_BUG(msg)
185 : : end if
186 : :
187 : : !Set number of cell
188 : 76 : harmonics_terms%ifcs%nrpt = nrpt
189 : :
190 : : !Allocation of total ifc
191 [ - + + - : 456 : ABI_MALLOC(harmonics_terms%ifcs%atmfrc,(3,natom,3,natom,nrpt))
+ - - + +
- - + + -
+ - + - -
+ - + -
+ ]
192 [ + + + + : 922666 : harmonics_terms%ifcs%atmfrc(:,:,:,:,:) = ifcs%atmfrc(:,:,:,:,:)
+ + + + +
+ ]
193 : :
194 : : !Allocation of ewald part of ifc
195 [ - + + - : 304 : ABI_MALLOC(harmonics_terms%ifcs%ewald_atmfrc,(3,natom,3,natom,nrpt))
- + + - -
+ + - - +
- + - + ]
196 [ + + + + : 922666 : harmonics_terms%ifcs%ewald_atmfrc(:,:,:,:,:) = ifcs%ewald_atmfrc(:,:,:,:,:)
+ + + + +
+ ]
197 : :
198 : : !Allocation of short range part of ifc
199 [ - + + - : 304 : ABI_MALLOC(harmonics_terms%ifcs%short_atmfrc,(3,natom,3,natom,nrpt))
- + + - -
+ + - - +
- + - + ]
200 [ + + + + : 922666 : harmonics_terms%ifcs%short_atmfrc(:,:,:,:,:) = ifcs%short_atmfrc(:,:,:,:,:)
+ + + + +
+ ]
201 : :
202 : : !Allocation of cell of ifc
203 [ - + + - : 228 : ABI_MALLOC(harmonics_terms%ifcs%cell,(3,nrpt))
+ - + - -
+ - + -
+ ]
204 [ + + + + ]: 8772 : harmonics_terms%ifcs%cell(:,:) = ifcs%cell(:,:)
205 : :
206 : : !Allocation of the dynamical matrix
207 : 76 : harmonics_terms%nqpt = 0
208 [ - + - - ]: 76 : if(present(nqpt).and.present(dynmat).and.present(qpoints).and.present(phfrq))then
209 : : call harmonics_terms_setDynmat(dynmat,harmonics_terms,natom,nqpt,&
210 : 0 : & harmonics_terms%phfrq,harmonics_terms%qpoints)
211 : : end if
212 : :
213 : : !Allocation of the elastic constants
214 [ + + + + ]: 3268 : harmonics_terms%elastic_constants = zero
215 [ - + ]: 76 : if (present(elastic_constants)) then
216 [ # # # # ]: 0 : harmonics_terms%elastic_constants = elastic_constants
217 : : end if
218 : :
219 : : !Allication of the dielectric tensor
220 [ + + + + ]: 988 : harmonics_terms%epsilon_inf = zero
221 [ - + ]: 76 : if (present(epsilon_inf)) then
222 [ # # # # ]: 0 : harmonics_terms%epsilon_inf = epsilon_inf
223 : : end if
224 : :
225 : : !Allocation of Effective charges array
226 [ - + + - : 228 : ABI_MALLOC(harmonics_terms%zeff,(3,3,natom))
+ - + - -
+ - + -
+ ]
227 [ + + + + : 5198 : harmonics_terms%zeff = zero
+ + ]
228 [ - + ]: 76 : if (present(zeff)) then
229 : 0 : call harmonics_terms_setEffectiveCharges(harmonics_terms,natom,zeff)
230 [ # # # # : 0 : harmonics_terms%zeff = zeff
# # # # #
# # # # #
# # ]
231 : : end if
232 : :
233 : : !Allocation of internal strain tensor
234 [ - + + - : 228 : ABI_MALLOC(harmonics_terms%strain_coupling,(6,3,natom))
+ - + - -
+ - + -
+ ]
235 [ + + + + : 8744 : harmonics_terms%strain_coupling = zero
+ + ]
236 [ - + ]: 76 : if (present(strain_coupling)) then
237 : 0 : call harmonics_terms_setInternalStrain(harmonics_terms,natom,strain_coupling)
238 : : end if
239 : :
240 [ - + - + ]: 76 : end subroutine harmonics_terms_init
241 : : !!***
242 : :
243 : :
244 : : !****f* m_harmonics_terms/harmonics_terms_free
245 : : !!
246 : : !! NAME
247 : : !! harmonics_terms_free
248 : : !!
249 : : !! FUNCTION
250 : : !! deallocate all dynamic memory for this harmonic datatype
251 : : !!
252 : : !! INPUTS
253 : : !! harmonics_terms<type(harmonics_terms_type)> = harmonics_terms datatype to be free
254 : : !!
255 : : !! OUTPUT
256 : : !! harmonics_terms<type(harmonics_terms_type)> = harmonics_terms datatype to be free
257 : : !!
258 : : !! SOURCE
259 : :
260 : 610 : subroutine harmonics_terms_free(harmonics_terms)
261 : :
262 : : implicit none
263 : :
264 : : !Arguments ------------------------------------
265 : : !scalars
266 : : !array
267 : : type(harmonics_terms_type), intent(inout) :: harmonics_terms
268 : : !Local variables-------------------------------
269 : : !scalars
270 : : !array
271 : :
272 : : ! *************************************************************************
273 : :
274 : 610 : harmonics_terms%nqpt = 0
275 [ + + + + ]: 26230 : harmonics_terms%elastic_constants = zero
276 [ + + + + ]: 7930 : harmonics_terms%epsilon_inf = zero
277 : :
278 [ + + ]: 610 : if(allocated(harmonics_terms%zeff))then
279 [ + + + + : 7203 : harmonics_terms%zeff=zero
+ + ]
280 : 105 : ABI_FREE(harmonics_terms%zeff)
281 : : end if
282 : :
283 [ + + ]: 610 : if(allocated(harmonics_terms%strain_coupling)) then
284 [ + + + + : 12117 : harmonics_terms%strain_coupling=zero
+ + ]
285 : 105 : ABI_FREE(harmonics_terms%strain_coupling)
286 : : end if
287 : :
288 [ + + ]: 610 : if(allocated(harmonics_terms%dynmat))then
289 [ + + + + : 185742 : harmonics_terms%dynmat=zero
+ + + + +
+ + + ]
290 : 105 : ABI_FREE(harmonics_terms%dynmat)
291 : : end if
292 : :
293 [ + + ]: 610 : if(allocated(harmonics_terms%phfrq))then
294 [ + + + + ]: 3904 : harmonics_terms%phfrq=zero
295 : 105 : ABI_FREE(harmonics_terms%phfrq)
296 : : end if
297 : :
298 [ + + ]: 610 : if(allocated(harmonics_terms%qpoints))then
299 [ + + + + ]: 1165 : harmonics_terms%qpoints=zero
300 : 105 : ABI_FREE(harmonics_terms%qpoints)
301 : : end if
302 : :
303 : 610 : call harmonics_terms%ifcs%free()
304 : :
305 : 610 : end subroutine harmonics_terms_free
306 : : !!***
307 : :
308 : : !****f* m_harmonics_terms/harmonics_terms_setInternalStrain
309 : : !!
310 : : !! NAME
311 : : !! harmonics_terms_setInternalStrain
312 : : !!
313 : : !! FUNCTION
314 : : !! Set the internal strain to the harmonics_terms
315 : : !!
316 : : !! INPUTS
317 : : !! natom = number of atoms
318 : : !! strain_coupling(6,3,natom) = internal strain coupling parameters
319 : : !!
320 : : !! OUTPUT
321 : : !! harmonics_terms<type(harmonics_terms_type)> = harmonics_terms datatype
322 : : !!
323 : : !! SOURCE
324 : :
325 [ - + ]: 76 : subroutine harmonics_terms_setInternalStrain(harmonics_terms,natom,strain_coupling)
326 : :
327 : : implicit none
328 : :
329 : : !Arguments ------------------------------------
330 : : !scalars
331 : : integer,intent(in) :: natom
332 : : !array
333 : : real(dp),intent(in) :: strain_coupling(:,:,:)
334 : : type(harmonics_terms_type), intent(inout) :: harmonics_terms
335 : : !Local variables-------------------------------
336 : : !scalars
337 : : !array
338 : : character(len=500) :: msg
339 : :
340 : : ! *************************************************************************
341 : :
342 : : ! 0-Checks inputs
343 [ - + ]: 76 : if(natom /= size(strain_coupling,3)) then
344 : : write(msg, '(a)' )&
345 : 0 : & ' natom has not the same size strain_coupling array. '
346 : 0 : ABI_BUG(msg)
347 : : end if
348 : :
349 : : ! 1-deallocate old array
350 [ + - ]: 76 : if(allocated(harmonics_terms%strain_coupling))then
351 : 76 : ABI_FREE(harmonics_terms%strain_coupling)
352 : : end if
353 : :
354 : : ! 2-allocate and copy the new array
355 [ - + + - : 228 : ABI_MALLOC(harmonics_terms%strain_coupling,(6,3,natom))
+ - + - -
+ - + -
+ ]
356 [ + + + + : 8744 : harmonics_terms%strain_coupling(:,:,:) = strain_coupling(:,:,:)
+ + ]
357 : :
358 : 76 : end subroutine harmonics_terms_setInternalStrain
359 : : !!***
360 : :
361 : :
362 : : !****f* m_harmonics_terms/harmonics_terms_setEffectiveCharges
363 : : !!
364 : : !! NAME
365 : : !! harmonics_terms_setEffectiveCharges
366 : : !!
367 : : !! FUNCTION
368 : : !! Set the effectives charges to the harmonics_terms
369 : : !!
370 : : !! INPUTS
371 : : !! natom = number of atoms
372 : : !! zeff(3,natom) = effective charges
373 : : !!
374 : : !! OUTPUT
375 : : !! harmonics_terms<type(harmonics_terms_type)> = harmonics_terms datatype
376 : : !!
377 : : !! SOURCE
378 : :
379 [ - + ]: 76 : subroutine harmonics_terms_setEffectiveCharges(harmonics_terms,natom,zeff)
380 : :
381 : : implicit none
382 : :
383 : : !Arguments ------------------------------------
384 : : !scalars
385 : : integer,intent(in) :: natom
386 : : !array
387 : : real(dp),intent(in) :: zeff(:,:,:)
388 : : type(harmonics_terms_type), intent(inout) :: harmonics_terms
389 : : !Local variables-------------------------------
390 : : !scalars
391 : : !array
392 : : character(len=500) :: msg
393 : : ! *************************************************************************
394 : :
395 : : ! 0-Checks inputs
396 [ - + ]: 76 : if(natom /= size(zeff,3)) then
397 : : write(msg, '(a)' )&
398 : 0 : & ' natom has not the same size zeff array. '
399 : 0 : ABI_BUG(msg)
400 : : end if
401 : :
402 : : ! 1-deallocate old array
403 [ + - ]: 76 : if(allocated(harmonics_terms%zeff))then
404 : 76 : ABI_FREE(harmonics_terms%zeff)
405 : : end if
406 : :
407 : : ! 2-allocate and copy the new array
408 [ - + + - : 228 : ABI_MALLOC(harmonics_terms%zeff,(3,3,natom))
+ - + - -
+ - + -
+ ]
409 [ + + + + : 5198 : harmonics_terms%zeff(:,:,:) = zeff(:,:,:)
+ + ]
410 : :
411 : :
412 : 76 : end subroutine harmonics_terms_setEffectiveCharges
413 : : !!***
414 : :
415 : : !****f* m_harmonics_terms/harmonics_terms_setDynmat
416 : : !!
417 : : !! NAME
418 : : !! harmonics_terms_setDynmat
419 : : !!
420 : : !! FUNCTION
421 : : !! Set the dynamical matricies to the harmonics_terms
422 : : !!
423 : : !! INPUTS
424 : : !! natom = number of atoms
425 : : !! nqpt = number of qpoints
426 : : !! dynmat(2,3,natom,3,natom,nqpt) = dynamical matrix in cartesian coordinates
427 : : !! phfrq(3*natom,nqpt) = frequency in hartree
428 : : !! qpoints(3,nqpt) = list of qpoints
429 : : !!
430 : : !! OUTPUT
431 : : !! harmonics_terms<type(harmonics_terms_type)> = harmonics_terms datatype
432 : : !!
433 : : !! SOURCE
434 : :
435 [ - + - + : 76 : subroutine harmonics_terms_setDynmat(dynmat,harmonics_terms,natom,nqpt,phfrq,qpoints)
- + ]
436 : :
437 : : implicit none
438 : :
439 : : !Arguments ------------------------------------
440 : : !scalars
441 : : integer,intent(in) :: natom,nqpt
442 : : !array
443 : : real(dp),intent(in) :: dynmat(:,:,:,:,:,:)
444 : : real(dp),intent(in) :: qpoints(:,:)
445 : : real(dp),intent(in) :: phfrq(:,:)
446 : : type(harmonics_terms_type), intent(inout) :: harmonics_terms
447 : : !Local variables-------------------------------
448 : : !scalars
449 : : !array
450 : : character(len=500) :: msg
451 : : ! *************************************************************************
452 : :
453 : : ! 0-Checks inputs
454 [ + - - + ]: 76 : if((natom /= size(dynmat,3)).or.(natom /= size(dynmat,5))) then
455 : : write(msg, '(a)' )&
456 : 0 : & ' natom has not the same size dynmat array. '
457 : 0 : ABI_BUG(msg)
458 : : end if
459 : :
460 [ - + ]: 76 : if (nqpt /= size(dynmat,6))then
461 : : write(msg, '(a)' )&
462 : 0 : & ' nqpt has not the same size dynmat array. '
463 : 0 : ABI_BUG(msg)
464 : : end if
465 : :
466 [ - + ]: 76 : if (nqpt /= size(qpoints,2))then
467 : : write(msg, '(a)' )&
468 : 0 : & ' nqpt has not the same size qpoints array. '
469 : 0 : ABI_BUG(msg)
470 : : end if
471 : :
472 [ - + ]: 76 : if (nqpt /= size(phfrq,2))then
473 : : write(msg, '(a)' )&
474 : 0 : & ' nqpt has not the same size phfrq array. '
475 : 0 : ABI_BUG(msg)
476 : : end if
477 : :
478 : : ! 1-deallocate old array
479 [ - + ]: 76 : if(allocated(harmonics_terms%dynmat))then
480 : 0 : ABI_FREE(harmonics_terms%dynmat)
481 : : end if
482 : :
483 [ - + ]: 76 : if(allocated(harmonics_terms%phfrq))then
484 : 0 : ABI_FREE(harmonics_terms%phfrq)
485 : : end if
486 : :
487 [ - + ]: 76 : if(allocated(harmonics_terms%qpoints))then
488 : 0 : ABI_FREE(harmonics_terms%qpoints)
489 : : end if
490 : :
491 : : ! 2-allocate and copy the new array
492 : 76 : harmonics_terms%nqpt = nqpt
493 : :
494 [ - + + - : 456 : ABI_MALLOC(harmonics_terms%dynmat,(2,3,natom,3,natom,nqpt))
+ - - + +
- + + + -
+ - + + -
+ - + -
+ ]
495 [ + + + + : 126404 : harmonics_terms%dynmat(:,:,:,:,:,:) = dynmat(:,:,:,:,:,:)
+ + + + +
+ + + ]
496 : :
497 [ + - - + : 304 : ABI_MALLOC(harmonics_terms%phfrq,(3*natom,nqpt))
+ - + + +
- + - + +
- + - + -
+ ]
498 [ + + + + ]: 2662 : harmonics_terms%phfrq(:,:) = phfrq(:,:)
499 : :
500 [ + + + - : 228 : ABI_MALLOC(harmonics_terms%qpoints,(3,nqpt))
+ - + + -
+ - + -
+ ]
501 [ + + + + ]: 796 : harmonics_terms%qpoints(:,:) = qpoints(:,:)
502 : :
503 : 76 : end subroutine harmonics_terms_setDynmat
504 : : !!***
505 : : !!****f* m_harmonics_terms/harmonics_terms_evaluateIFC
506 : : !! NAME
507 : : !! harmonics_terms_evaluateIFC
508 : : !!
509 : : !! FUNCTION
510 : : !! This fonction compute the contribution of the ifc harmonic part of
511 : : !! the energy and forces.
512 : : !!
513 : : !! INPUTS
514 : : !! atmfrc(3,natom_uc,3,natom_uc,nrpt) = atomic force constants
515 : : !! disp(3,natom_sc) = atomics displacement between configuration and the reference
516 : : !! ncell = total number of cell to treat
517 : : !! nrpt = total number of rpt to treat
518 : : !! natom_sc = number of atoms in the supercell
519 : : !! natom_uc = number of atoms in the unit cell
520 : : !! nrpt = number of rpt
521 : : !! atmrpt_index(nrpt,cell) = For each cell in the supercell and each rpt,
522 : : !! give the index of the first atoms in the rpt cell
523 : : !! rpt(nrpt) = index of rpt in atmfrc (6th dimension)
524 : : !! index_cells(3,ncell) = indexes of the cells into supercell (-1 -1 -1 ,...,1 1 1)
525 : : !! comm=MPI communicator
526 : : !!
527 : : !! OUTPUT
528 : : !! energy = contribution of the ifc to the energy
529 : : !! fcart(3,natom) = contribution of the ifc to the forces
530 : : !!
531 : : !! PARENT
532 : : !! effective_potential_evaluate
533 : : !!
534 : : !! SOURCE
535 : :
536 : 27038 : subroutine harmonics_terms_evaluateIFC(atmfrc,disp,energy,fcart,natom_sc,natom_uc,&
537 : 13519 : & ncell,nrpt,atmrpt_index,index_cells,sc_size,rpt,comm)
538 : :
539 : : implicit none
540 : :
541 : : !Arguments -------------------------------
542 : : ! scalars
543 : : real(dp),intent(out) :: energy
544 : : integer,intent(in) :: natom_uc,natom_sc,ncell,nrpt
545 : : integer,intent(in) :: comm
546 : : ! array
547 : : integer,intent(in) :: sc_size(3),atmrpt_index(nrpt,ncell)
548 : : integer,intent(in) :: index_cells(4,ncell),rpt(nrpt)
549 : : real(dp),intent(in) :: atmfrc(3,natom_uc,3,natom_uc,nrpt)
550 : : real(dp),intent(in) :: disp(3,natom_sc)
551 : : real(dp),intent(out) :: fcart(3,natom_sc)
552 : :
553 : : !Local variables-------------------------------
554 : : ! scalar
555 : : integer :: i1,i2,i3,ia,ib,icell,ierr,irpt,irpt_tmp,ii,jj,kk,ll
556 : : integer :: mu,nu
557 : : real(dp):: disp1,disp2,ifc,tmp_etot1,tmp_etot2
558 : : !Variables for separation of short and dipdip ifc contribution
559 : : !real(dp):: short_ifc,ewald_ifc
560 : : !real(dp):: tmp_ewald1,tmp_ewald2,tmp_short1,tmp_short2
561 : : ! array
562 : : character(500) :: msg
563 : :
564 : : ! *************************************************************************
565 : :
566 [ + + + - : 54076 : if (any(sc_size <= 0)) then
- + ]
567 : 0 : write(msg,'(a,a)')' sc_size can not be inferior or equal to zero'
568 : 0 : ABI_ERROR(msg)
569 : : end if
570 : :
571 : : ! Initialisation of variables
572 : 13519 : energy = zero
573 [ + + + + ]: 36076719 : fcart(:,:) = zero
574 : :
575 [ + + ]: 1786383 : do icell = 1,ncell
576 : 1772864 : i1 = index_cells(1,icell)
577 : 1772864 : i2 = index_cells(2,icell)
578 : 1772864 : i3 = index_cells(3,icell)
579 : : ! index of the first atom in the current cell
580 : 1772864 : ii = index_cells(4,icell)
581 [ + + ]: 377817063 : do irpt_tmp = 1,nrpt
582 : 376030680 : irpt = rpt(irpt_tmp)
583 : : ! index of the first atom in the irpt cell
584 : 376030680 : jj = atmrpt_index(irpt_tmp,icell)
585 : : ! Loop over the atom in the cell
586 [ + + ]: 2264040144 : do ib = 1, natom_uc
587 : 1886236600 : ll = jj + ib
588 [ + + ]: 7920977080 : do nu=1,3
589 : 5658709800 : disp2 = disp(nu,ll)
590 [ + + ]: 36206405800 : do ia = 1, natom_uc
591 : 28661459400 : kk = ii + ia
592 [ + + ]: >12030*10^7 : do mu=1,3
593 : 85984378200 : disp1 = disp(mu,kk)
594 : 85984378200 : ifc = atmfrc(mu,ia,nu,ib,irpt)
595 : :
596 : : ! if(abs(ifc) > tol10)then
597 : 85984378200 : tmp_etot1 = disp2 * ifc
598 : : ! accumule energy
599 : 85984378200 : tmp_etot2 = disp1*tmp_etot1
600 : 85984378200 : energy = energy + tmp_etot2
601 : : ! accumule forces
602 : >11464*10^7 : fcart(mu,kk) = fcart(mu,kk) + tmp_etot1
603 : : ! end if
604 : : end do
605 : : end do
606 : : end do
607 : : end do
608 : : end do
609 : : end do
610 : :
611 : 13519 : energy = half * energy
612 : : ! MPI_SUM
613 : 13519 : call xmpi_sum(energy, comm, ierr)
614 : 13519 : call xmpi_sum(fcart , comm, ierr)
615 : :
616 : 13519 : end subroutine harmonics_terms_evaluateIFC
617 : : !!***
618 : :
619 : : !!****f* m_harmonics_terms/harmonics_terms_evaluateElastic
620 : : !! NAME
621 : : !! harmonics_terms_evaluateElastic
622 : : !!
623 : : !! FUNCTION
624 : : !! Compute the energy, forces and stresses related to the application of strain
625 : : !!
626 : : !! INPUTS
627 : : !! elastic_constants(6,6) = elastic constants in Hartree
628 : : !! disp(3,natom_sc) = atomics displacement between configuration and the reference
629 : : !! natom = number of atoms in the supercell
630 : : !! natom_uc = number of atoms in the unit cell
631 : : !! ncell = total number of cell
632 : : !! strain_coupling(6,3,natom) = internal strain coupling parameters
633 : : !! strain(6) = strain between configuration and the reference
634 : : !!
635 : : !! OUTPUT
636 : : !! energy = contribution to the energy
637 : : !! fcart(3,natom) = contribution to the forces
638 : : !! strten(6) = contribution to the stress tensor
639 : : !!
640 : : !! SOURCE
641 : : !!
642 : 13519 : subroutine harmonics_terms_evaluateElastic(elastic_constants,disp,energy,fcart,natom,natom_uc,ncell,&
643 : 13519 : & strain_coupling,strten,strain)
644 : :
645 : : real(dp),intent(out):: energy
646 : : integer, intent(in) :: natom,natom_uc,ncell
647 : : ! array
648 : : real(dp),intent(in) :: elastic_constants(6,6),strain_coupling(6,3,natom)
649 : : real(dp),intent(out):: strten(6)
650 : : real(dp),intent(out):: fcart(3,natom)
651 : : real(dp),intent(in) :: disp(3,natom)
652 : : real(dp),intent(in) :: strain(6)
653 : :
654 : : !Local variables-------------------------------
655 : : ! scalar
656 : : integer :: ia,ii,mu,alpha,beta
657 : : real(dp):: cij
658 : : ! array
659 : : ! *************************************************************************
660 : :
661 : 13519 : energy = zero
662 [ + + + + ]: 36076719 : fcart = zero
663 : 13519 : strten = zero
664 : :
665 : : ! write(*,*) "----- STRAIN -----"
666 : : ! write(*,*) strain
667 : :
668 : : !1- Part due to elastic constants
669 [ + + ]: 94633 : do alpha=1,6
670 [ + + ]: 581317 : do beta=1,6
671 : : ! write(*,*) "--- cij --- alpha: ", alpha, " beta: ", beta
672 : 486684 : cij = ncell*elastic_constants(alpha,beta)
673 : : ! write(*,*) cij
674 : 486684 : energy = energy + half*cij*strain(alpha)*strain(beta)
675 : 567798 : strten(alpha) = strten(alpha) + cij*strain(beta)
676 : : end do
677 : : ! write(*,*) "strten(",alpha,"): ", strten(alpha)
678 : : end do
679 : :
680 : : !2-Part due to the internal strain coupling parameters
681 : : ii = 1
682 [ + + ]: 9029319 : do ia = 1,natom
683 [ + + ]: 36063200 : do mu = 1,3
684 [ + + ]: 198347600 : do alpha=1,6
685 : 162284400 : cij = strain_coupling(alpha,mu,ii)
686 : : ! Accumulte for this atom
687 : 162284400 : energy = energy + half*cij*strain(alpha)*disp(mu,ia)
688 : 162284400 : fcart(mu,ia) = fcart(mu,ia) + half*cij*strain(alpha)
689 : 189331800 : strten(alpha) = strten(alpha) + half*cij*disp(mu,ia)
690 : : end do
691 : : end do
692 : 9015800 : ii = ii +1
693 : : ! Reset to 1 if the number of atoms is superior than in the initial cell
694 [ + + ]: 9029319 : if(ii==natom_uc+1) ii = 1
695 : : end do
696 : :
697 : : ! write(*,*) "--- STRTEN at the end --- "
698 : : ! write(*,*) strten(:)
699 : :
700 : 13519 : end subroutine harmonics_terms_evaluateElastic
701 : : !!***
702 : :
703 : : !****f* m_harmonics_terms/harmonics_terms_applySumRule
704 : : !!
705 : : !! NAME
706 : : !! harmonics_terms_applySumRule
707 : : !!
708 : : !! FUNCTION
709 : : !! Apply the acoustic sum rule on the inter-atomic force constants
710 : : !!
711 : : !! INPUTS
712 : : !! ifc<type(ifc_type)> = interatomic forces constants
713 : : !! asr = acoustic sum rule option (see anaddb help)
714 : : !! natom = number of atoms
715 : : !! option = optional if |no present asr is done on total ifc
716 : : !! |present and 1 asr is done on short part
717 : : !! |present and 2 asr is done on ewald part
718 : : !!
719 : : !! OUTPUT
720 : : !! ifc<type(ifc_type)> = interatomic forces constants
721 : : !!
722 : : !! SOURCE
723 : :
724 : 67 : subroutine harmonics_terms_applySumRule(asr,ifc,natom,option)
725 : :
726 : : implicit none
727 : :
728 : : !Arguments ------------------------------------
729 : : !scalars
730 : : integer,intent(in) :: asr
731 : : integer,intent(in) :: natom
732 : : integer,optional,intent(in) :: option
733 : : !array
734 : : type(ifc_type),target,intent(inout) :: ifc
735 : : !Local variables-------------------------------
736 : : !scalar
737 : : integer :: ia,ib,irpt,irpt_ref
738 : : integer :: mu,nu
739 : : real(dp) :: sum
740 : : character(500) :: msg
741 : : !array
742 : 67 : real(dp),pointer :: atmfrc(:,:,:,:,:)
743 : : ! *************************************************************************
744 : :
745 : 67 : irpt_ref = 0
746 : : ! Found the cell of reference
747 [ + + ]: 6494 : do irpt = 1,ifc%nrpt
748 : : if(ifc%cell(1,irpt)==0.and.&
749 [ + + + + : 6427 : & ifc%cell(2,irpt)==0.and.&
+ + ]
750 : 67 : & ifc%cell(3,irpt)==0) then
751 : 6427 : irpt_ref = irpt
752 : : cycle
753 : : end if
754 : : end do
755 : :
756 [ - + ]: 67 : if (irpt_ref<=0) then
757 : 0 : write(msg,'(a,a)')' Unable to find the cell of reference in IFC'
758 : 0 : ABI_ERROR(msg)
759 : : end if
760 : :
761 : 67 : atmfrc => ifc%atmfrc
762 [ - + ]: 67 : if (present(option)) then
763 [ # # ]: 0 : if (option == 1) then
764 : 0 : nullify(atmfrc)
765 : 0 : atmfrc => ifc%short_atmfrc
766 : 0 : write(msg,'(3a)') ch10," Impose acoustic sum rule on short range"
767 [ # # ]: 0 : else if (option == 2) then
768 : 0 : nullify(atmfrc)
769 : 0 : atmfrc => ifc%ewald_atmfrc
770 : 0 : write(msg,'(3a)') ch10," Impose acoustic sum rule on long range"
771 : : end if
772 : : else
773 : 67 : write(msg,'(3a)') ch10," Impose acoustic sum rule on total ifc"
774 : : end if
775 : 67 : call wrtout(ab_out,msg,'COLL')
776 : 67 : call wrtout(std_out,msg,'COLL')
777 : :
778 : : !impose acoustic sum rule:
779 [ + + ]: 268 : do mu=1,3
780 [ + + ]: 871 : do nu=1,3
781 [ + + ]: 4008 : do ia=1,natom
782 : : sum=zero
783 [ + + ]: 26460 : do ib=1,natom
784 : : ! Get the sum of interatomic forces acting on the atom ia,
785 : : ! either in a symmetrical manner, or an unsymmetrical one.
786 [ - + ]: 26460 : if(asr==1)then
787 [ # # ]: 0 : do irpt=1, ifc%nrpt
788 : 0 : sum=sum+atmfrc(mu,ia,nu,ib,irpt)
789 : : end do
790 [ + - ]: 23256 : else if(asr==2)then
791 [ + + ]: 1842957 : do irpt=1, ifc%nrpt
792 : : sum=sum+&
793 : : & (atmfrc(mu,ia,nu,ib,irpt)+&
794 : 1842957 : & atmfrc(nu,ia,mu,ib,irpt))/2
795 : : end do
796 : : end if
797 : : end do
798 : :
799 : : ! Correct the self-interaction in order to fulfill the ASR
800 : : atmfrc(mu,ia,nu,ia,irpt_ref)=&
801 : 3204 : & atmfrc(mu,ia,nu,ia,irpt_ref)-sum
802 [ + - ]: 3807 : if(asr==2)then
803 : : atmfrc(nu,ia,mu,ia,irpt_ref)=&
804 : 3204 : & atmfrc(mu,ia,nu,ia,irpt_ref)
805 : : end if
806 : : end do
807 : : end do
808 : : end do
809 : :
810 [ - + ]: 67 : if (present(option)) then
811 [ # # ]: 0 : if (option == 1) then
812 [ # # # # : 0 : ifc%short_atmfrc = atmfrc(:,:,:,:,:)
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
813 [ # # ]: 0 : else if (option == 2) then
814 [ # # # # : 0 : ifc%ewald_atmfrc = atmfrc(:,:,:,:,:)
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
815 : : end if
816 : : else
817 [ + - + + : 5134332 : ifc%atmfrc(:,:,:,:,:) = atmfrc(:,:,:,:,:)
+ + + + +
+ + + + +
+ + + + +
+ + + ]
818 : : end if
819 : :
820 : 67 : end subroutine harmonics_terms_applySumRule
821 : : !!***
822 : :
823 [ # # # # : 0 : end module m_harmonics_terms
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
824 : : !!***
|
