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