Line data Source code
1 : !!****m* ABINIT/m_ifc
2 : !! NAME
3 : !! m_ifc
4 : !!
5 : !! FUNCTION
6 : !! This module contains the declaration of data types and methods
7 : !! used to handle interatomic force constants.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2011-2026 ABINIT group (XG,MJV,EB,MG,GA)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : MODULE m_ifc
24 :
25 : use defs_basis
26 : use m_errors
27 : use m_abicore
28 : use m_xmpi
29 : use m_sort
30 : use m_cgtools
31 : use m_ddb
32 : use m_ddb_hdr
33 : use m_symkpt
34 : use netcdf
35 : use m_nctk
36 :
37 : use m_io_tools, only : open_file
38 : use m_numeric_tools, only : arth
39 : use m_fstrings, only : ktoa, int2char4, sjoin, itoa, ltoa, ftoa
40 : use m_matrix, only : matr3inv
41 : use m_special_funcs, only : abi_derfc
42 : use m_time, only : cwtime, cwtime_report, timab
43 : use m_copy, only : alloc_copy
44 : use m_pptools, only : printbxsf
45 : use m_lebedev, only : lebedev_t, lebedev_ngrids
46 : use m_ewald, only : ewald9, ewald9_2D
47 : use m_crystal, only : crystal_t
48 : use m_geometry, only : phdispl_cart2red, normv, mkrdim
49 : use m_kpts, only : kpts_ibz_from_kptrlatt, smpbz
50 : use m_bz_mesh, only : kpath_t
51 : use m_dynmat, only : canct9, dist9 , ifclo9, axial9, q0dy3_apply, q0dy3_calc, asrif9, dynmat_dq, &
52 : make_bigbox, canat9, chkrp9, ftifc_q2r, wght9, nanal9, gtdyn9, dymfz9, &
53 : massmult_and_breaksym, dfpt_phfrq, dfpt_prtph, d2cart_to_red, ftifc_r2q
54 :
55 : implicit none
56 :
57 : private
58 :
59 : public :: defwrite_nonana_terms
60 : ! Write to ncfile the phonon frequencies and displacements for q --> 0
61 : ! in the presence of non-analytical behaviour.
62 : !!***
63 :
64 : !!****t* m_ifc/ifc_type
65 : !! NAME
66 : !! ifc_type
67 : !!
68 : !! FUNCTION
69 : !! Contains the necessary data to interpolate the
70 : !! phonon bandstructure and eigenvectors in reciprocal space (ie.
71 : !! interatomic force constants and corresponding real space grid info).
72 : !!
73 : !! SOURCE
74 :
75 : type,public :: ifc_type
76 :
77 : integer :: natom = -1
78 : ! Number of atoms in the unit cell.
79 :
80 : integer :: mpert = -1
81 : ! Maximum number of ipert.
82 :
83 : integer :: asr = -1
84 : ! Option for the treatment of the Acoustic Sum Rule.
85 :
86 : integer :: brav = -1
87 : ! Option for the sampling of the BZ (anaddb input variable)
88 :
89 : real(dp) :: dielt_env = one
90 : ! Dielectric constant of environment for 2D materials electrostatics
91 :
92 : real(dp) :: dielt_thick(2) = zero
93 : ! Dielectric thickness of the 2D materials for electrostatics
94 :
95 : integer :: sys_dim = 1
96 : ! System dimmensionality for electrostatic treatment + rotational invariance
97 :
98 : integer :: dipdip = -1
99 : ! dipole dipole interaction flag.
100 :
101 : integer :: dipquad = -1
102 : ! dipole quadrupole interaction flag.
103 :
104 : integer :: quadquad = -1
105 : ! dipole quadrupole interaction flag.
106 :
107 : integer :: symdynmat = -1
108 : ! If equal to 1, the dynamical matrix is symmetrized in dfpt_phfrq before the diagonalization.
109 :
110 : integer :: nqshft = -1
111 : ! Number of shifts in the q-mesh (usually 1 since the mesh is gamma-centered!)
112 :
113 : integer :: nqibz = -1
114 : ! Number of points in the IBZ
115 :
116 : integer :: nqbz = -1
117 : ! Number of points in the full BZ
118 :
119 : integer :: nrpt = -1
120 : ! Number of real space points used to integrate IFC (for interpolation of dynamical matrices)
121 :
122 : integer :: ngqpt(3) = -1
123 : ! Number of division in the Q mesh.
124 :
125 : integer :: ewald_option = -1
126 : ! Option for the ewald sum
127 :
128 : real(dp) :: eta = -one
129 : ! eta: parameter used to split R and G-space summation. See ewald9.
130 :
131 : real(dp) :: rprim(3,3),gprim(3,3),acell(3)
132 : ! These values are used to call anaddb routines that don't use rprimd, gprimd
133 :
134 : real(dp) :: dielt(3,3)
135 : ! Dielectric tensor (Cartesian coordinates)
136 :
137 : real(dp) :: omega_minmax(2)
138 : ! Min and max frequency obtained on the initial ab-initio q-mesh (-+ 30 cmm1)
139 : ! Used to generate frequency meshes for DOSes.
140 :
141 : real(dp) :: r_inscribed_sphere
142 : ! radius of biggest sphere inscribed in the WS supercell
143 :
144 : real(dp),allocatable :: amu(:)
145 : ! amu(ntypat)
146 : ! mass of the atoms (atomic mass unit)
147 :
148 : real(dp),allocatable :: atmfrc(:,:,:,:,:)
149 : ! atmfrc(3,natom,3,natom,nrpt)
150 : ! Inter atomic forces in real space
151 :
152 : integer,allocatable :: cell(:,:)
153 : ! cell(nrpt,3)
154 : ! Give the index of the cell and irpt
155 :
156 : real(dp),allocatable :: ewald_atmfrc(:,:,:,:,:)
157 : ! Ewald_atmfrc(3,natom,3,natom,nrpt)
158 : ! Ewald Inter atomic forces in real space
159 :
160 : real(dp),allocatable :: short_atmfrc(:,:,:,:,:)
161 : ! short_atmfrc(3,natom,3,natom,nrpt)
162 : ! Short range part of Inter atomic forces in real space
163 :
164 : real(dp),allocatable :: qshft(:,:)
165 : ! qshft(3,nqshft)
166 : ! The shifts of the q-mesh
167 :
168 : real(dp), allocatable :: rpt(:,:)
169 : ! rpt(3,nrpt)
170 : ! Real space points in canonical type coordinates.
171 :
172 : real(dp),allocatable :: wghatm(:,:,:)
173 : ! wghatm(natom,natom,nrpt)
174 : ! Weights for each point and atom in the Wigner Seitz supercell in real space.
175 :
176 : real(dp),allocatable :: rcan(:,:)
177 : ! rcan(3,natom)
178 : ! Atomic position in canonical coordinates.
179 :
180 : real(dp),allocatable :: trans(:,:)
181 : ! trans(3,natom)
182 : ! Atomic translations: xred = rcan + trans
183 :
184 : real(dp),allocatable :: dyewq0(:,:,:)
185 : ! dyewq0(3,3,natom)
186 : ! Atomic electrostatic self-interaction correction to the dynamical matrix (only when dipdip = 1).
187 :
188 : real(dp),allocatable :: zeff(:,:,:)
189 : ! zeff(3,3,natom)
190 : ! Born effective charge on each atom, versus electric field and atomic displacement.
191 : ! Cartesian coordinates
192 :
193 : real(dp),allocatable :: qdrp_cart(:,:,:,:)
194 : ! qdrp_cart(3,3,3,natom)
195 : ! Quadrupole tensor on each atom
196 : ! Cartesian coordinates
197 :
198 : real(dp),allocatable :: qibz(:,:)
199 : ! qibz(3,nqibz))
200 : ! List of q-points in the IBZ
201 :
202 : real(dp),allocatable :: wtq(:)
203 : ! wtq(nqibz))
204 : ! q-point Weights.
205 :
206 : real(dp),allocatable :: qbz(:,:)
207 : ! qbz(3,nqbz))
208 : ! List of q-points in the full BZ
209 :
210 : real(dp),allocatable :: dynmat(:,:,:,:,:,:)
211 : ! dynmat(2,3,natom,3,natom,nqbz))
212 : ! dynamical matrices relative to the q points of the BZ. sampling
213 : ! Note that the long-range dip-dip part has been removed if dipdip = 1
214 : ! Moreover the array is multiplied by a phase shift in mkifc9.
215 :
216 : contains
217 :
218 : procedure :: free => ifc_free
219 : ! Release memory.
220 :
221 : procedure :: print => ifc_print
222 : ! Print info on the object.
223 :
224 : procedure :: fourq => ifc_fourq
225 : ! Use Fourier interpolation to compute interpolated frequencies w(q) and eigenvectors e(q).
226 :
227 : procedure :: get_dcdq => ifc_get_dcdq
228 : ! Compute first derivative of reciprocal-space IFCs based on IFCs.
229 :
230 : procedure :: get_dwdq => ifc_get_dwdq
231 : ! Compute phonon group velocities at an arbitrary q-point.
232 :
233 : procedure :: get_phmesh => ifc_get_phmesh
234 : ! Build linear mesh for phonons.
235 :
236 : procedure :: speedofsound => ifc_speedofsound
237 : ! Compute the speed of sound by averaging phonon group velocities.
238 :
239 : procedure :: write => ifc_write
240 : ! Print the ifc (output, netcdf and text file).
241 :
242 : procedure :: outphbtrap => ifc_outphbtrap
243 : ! Print out phonon frequencies on regular grid for BoltzTrap code.
244 :
245 : procedure :: printbxsf => ifc_printbxsf
246 : ! Output phonon isosurface in Xcrysden format.
247 :
248 : procedure :: calcnwrite_nana_terms => ifc_calcnwrite_nana_terms
249 : ! Compute phonons for q--> 0 with LO-TO.
250 :
251 : procedure :: calcnwrite_nana_terms_qpath => ifc_calcnwrite_nana_terms_qpath
252 : ! Compute phonons for q--> 0 with LO-TO from qpath_t and write results to netcdf file.
253 :
254 : procedure :: init => ifc_init
255 : ! Constructor from DDB datatype.
256 :
257 : procedure :: from_file => ifc_from_file
258 : ! Constructor from filename.
259 :
260 : procedure :: to_ddb => ifc_to_ddb
261 : ! Construct a DDB object.
262 :
263 : end type ifc_type
264 : !!***
265 :
266 : !----------------------------------------------------------------------
267 :
268 : contains !===========================================================
269 : !!***
270 :
271 : !----------------------------------------------------------------------
272 :
273 : !!****f* m_ifc/ifc_free
274 : !! NAME
275 : !! ifc_free
276 : !!
277 : !! FUNCTION
278 : !! Deallocate memory for the ifc_type structure
279 : !!
280 : !! SOURCE
281 :
282 1536 : subroutine ifc_free(ifc)
283 :
284 : !Arguments ------------------------------------
285 : class(ifc_type),intent(inout) :: ifc
286 : ! ************************************************************************
287 :
288 1536 : ABI_SFREE(ifc%amu)
289 1536 : ABI_SFREE(ifc%atmfrc)
290 1536 : ABI_SFREE(ifc%cell)
291 1536 : ABI_SFREE(ifc%ewald_atmfrc)
292 1536 : ABI_SFREE(ifc%short_atmfrc)
293 1536 : ABI_SFREE(ifc%qshft)
294 1536 : ABI_SFREE(ifc%rpt)
295 1536 : ABI_SFREE(ifc%wghatm)
296 1536 : ABI_SFREE(ifc%rcan)
297 1536 : ABI_SFREE(ifc%trans)
298 1536 : ABI_SFREE(ifc%dyewq0)
299 1536 : ABI_SFREE(ifc%qibz)
300 1536 : ABI_SFREE(ifc%wtq)
301 1536 : ABI_SFREE(ifc%qbz)
302 1536 : ABI_SFREE(ifc%zeff)
303 1536 : ABI_SFREE(ifc%qdrp_cart)
304 1536 : ABI_SFREE(ifc%dynmat)
305 :
306 1536 : end subroutine ifc_free
307 : !!***
308 :
309 : !----------------------------------------------------------------------
310 :
311 : !!****f* m_ifc/ifc_init
312 : !! NAME
313 : !! ifc_init
314 : !!
315 : !! FUNCTION
316 : !! Initialize the dynamical matrix as well as the IFCs.
317 : !! taking into account the dipole-dipole, dipole-quadrupole and quadrupole-quadrupole interaction.
318 : !!
319 : !! INPUTS
320 : !! crystal<type(crystal_t)> = Information on the crystalline structure.
321 : !! ddb<type(ddb_type)> = Database with derivatives.
322 : !! brav=bravais lattice (1 or -1=simple lattice, 2=face centered lattice, 3=centered lattice, 4=hexagonal lattice)
323 : !! asr= Option for the imposition of the ASR
324 : !! 0 => no ASR,
325 : !! 1 => modify "asymmetrically" the diagonal element
326 : !! 2 => modify "symmetrically" the diagonal element
327 : !! symdynmat=if 1, (re)symmetrize the dynamical matrix, except if Gamma wavevector with electric field added.
328 : !! dipdip=
329 : !! if 0, no dipole-dipole interaction was subtracted in atmfrc
330 : !! if 1, atmfrc has been build without dipole-dipole part
331 : !! rfmeth =
332 : !! 1 if non-stationary block
333 : !! 2 if stationary block
334 : !! 3 if third order derivatives
335 : !! dielt(3,3)=dielectric tensor.
336 : !! zeff(3,3,natom)=effective charge on each atom, versus electric field and atomic displacement
337 : !! prtsrlr: TODO: TO BE REMOVED
338 : !! enunit: TODO: TO BE REMOVED
339 : !! dielt(3,3)=dielectric tensor
340 : !! ngqpt_in = input values of ngqpt
341 : !! nqshft=Number of shifths in q-grid.
342 : !! q1shft(3,nqshft)=Shifts for q-grid
343 : !! nsphere=number of atoms to be included in the cut-off sphere for interatomic force constant.
344 : !! 0: maximum extent allowed by the grid.
345 : !! > 0: Apply cutoff
346 : !! -1: Analyze the effect of different nsphere values on the phonon spectrum, in particular the
347 : !! frequencies around gamma.
348 : !! rifcsph=radius for cutoff of IFC.
349 : !! comm=MPI communicator.
350 : !! [Ifc_coarse]=Optional.
351 : !! [dipquad] = if 1, atmfrc has been build without dipole-quadrupole part
352 : !! [quadquad] = if 1, atmfrc has been build without quadrupole-quadrupole part
353 : !! prtout=write info to the output and log files.
354 : !!
355 : !! OUTPUT
356 : !! Ifc<ifc_type>=Object containing the dynamical matrix and the IFCs.
357 : !!
358 : !! SOURCE
359 :
360 242 : subroutine ifc_init(Ifc,crystal,ddb,brav,asr,symdynmat,dipdip,&
361 242 : rfmeth,ngqpt_in,nqshft,q1shft,dielt,zeff,qdrp_cart,nsphere,rifcsph,&
362 : prtsrlr,enunit,& ! TODO: TO BE REMOVED
363 : sys_dim,comm, &
364 : Ifc_coarse,dipquad,quadquad,dielt_env,dielt_thick,prtout) ! Optional
365 :
366 : !Arguments ------------------------------------
367 : class(ifc_type),intent(inout) :: Ifc
368 : integer,intent(in) :: asr,dipdip,symdynmat,nqshft,rfmeth,nsphere,comm,sys_dim
369 : real(dp),intent(in) :: rifcsph
370 : type(crystal_t),intent(in) :: Crystal
371 : type(ddb_type),intent(in) :: ddb
372 : type(ifc_type),optional,intent(in) :: Ifc_coarse
373 : integer,optional,intent(in) :: dipquad, quadquad
374 : logical, optional, intent(in) :: prtout
375 : real(dp), optional, intent(in) :: dielt_env,dielt_thick(2)
376 :
377 : !arrays
378 : integer,intent(in) :: ngqpt_in(3)
379 : real(dp),intent(in) :: q1shft(3,nqshft)
380 : real(dp),intent(in) :: dielt(3,3),zeff(3,3,Crystal%natom)
381 : real(dp),intent(in) :: qdrp_cart(3,3,3,Crystal%natom)
382 : !anaddb variables (TO BE REMOVED)
383 : integer,intent(in) :: prtsrlr,enunit,brav
384 :
385 : !Local variables -------------------------
386 : !scalars
387 : integer,parameter :: timrev1=1,iout0=0,chksymbreak0=0
388 : integer :: mpert,iout,iqpt,mqpt,nsym,ntypat,iq_ibz,iq_bz,ii,natom
389 : integer :: nqbz,option,plus,sumg0,irpt,irpt_new, nprocs,my_rank,my_ierr,ierr
390 : logical :: prtout_
391 : real(dp),parameter :: qphnrm=one
392 : real(dp) :: cpu,wall,gflops,xval,rcut_min,r_inscribed_sphere,toldist
393 : character(len=500*4) :: msg
394 1452 : type(ifc_type) :: ifc_tmp
395 : !arrays
396 : integer :: ngqpt(9),qptrlatt(3,3)
397 484 : integer,allocatable :: qmissing(:),ibz2bz(:),bz2ibz_smap(:,:)
398 : real(dp) :: gprim(3,3),rprim(3,3),qpt(3),rprimd(3,3), gprim_tmp(3,3), rprim_tmp(3,3)
399 968 : real(dp):: rcan(3,Crystal%natom),trans(3,Crystal%natom),dyewq0(3,3,Crystal%natom)
400 484 : real(dp) :: displ_cart(2*3*Crystal%natom*3*Crystal%natom)
401 484 : real(dp) :: phfrq(3*Crystal%natom)
402 484 : real(dp) :: eigvec(2,3,Crystal%natom,3,Crystal%natom)
403 484 : real(dp),allocatable :: dyew(:,:,:,:,:),out_d2cart(:,:,:,:,:)
404 242 : real(dp),allocatable :: dynmatfull(:,:,:,:,:,:),dynmat_sr(:,:,:,:,:,:),dynmat_lr(:,:,:,:,:,:) ! for OmegaSRLR
405 242 : real(dp),allocatable :: wtq(:),wtq_folded(:),qbz(:,:)
406 : !******************************************************************
407 : ! TODO
408 : ! - Remove brav
409 : ! - Separate into ifc_init and ifc_from_ddb
410 :
411 242 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
412 242 : call cwtime(cpu, wall, gflops, "start")
413 :
414 242 : mpert = ddb%mpert
415 242 : iout = ab_out
416 242 : prtout_ = .true.
417 242 : if (present(prtout)) prtout_ = prtout
418 :
419 242 : nsym = Crystal%nsym
420 242 : natom = Crystal%natom
421 242 : ntypat = Crystal%ntypat
422 3146 : rprimd = Crystal%rprimd
423 :
424 968 : ngqpt=0; ngqpt(1:3)=ngqpt_in(1:3)
425 :
426 : ! Copy important parameters in Ifc
427 242 : Ifc%natom = natom
428 242 : Ifc%mpert = mpert
429 242 : Ifc%asr = asr
430 242 : Ifc%brav = brav
431 242 : Ifc%dielt_env=0; if (present(dielt_env)) Ifc%dielt_env = dielt_env
432 836 : Ifc%dielt_thick(:)=0; if (present(dielt_thick)) Ifc%dielt_thick = dielt_thick(1:2)
433 242 : Ifc%sys_dim = sys_dim
434 242 : Ifc%dipdip = abs(dipdip)
435 242 : Ifc%dipquad=0; if (present(dipquad)) Ifc%dipquad = dipquad
436 242 : Ifc%quadquad=0; if (present(quadquad)) Ifc%quadquad = quadquad
437 242 : Ifc%symdynmat = symdynmat
438 968 : Ifc%ngqpt = ngqpt_in(1:3)
439 242 : Ifc%nqshft = nqshft
440 242 : call alloc_copy(q1shft(:,1:Ifc%nqshft),Ifc%qshft)
441 242 : Ifc%ewald_option = 0; if (dipdip < 0) Ifc%ewald_option = 1 !HM TODO: expose this in the init?
442 :
443 :
444 : ! These two ways of defining rprim, gprim, acell should be equivalent.
445 : ! Note that the crystal object does not know about acell.
446 3146 : Ifc%rprim = ddb%rprim
447 3146 : Ifc%gprim = ddb%gprim
448 968 : Ifc%acell = ddb%acell
449 : !Ifc%rprim = crystal%rprimd
450 : !Ifc%gprim = crystal%gprimd
451 : !Ifc%acell = one
452 : !
453 : !! When brav/=1, we need specific values of rprim.
454 : !if (abs(brav)/=1) then
455 : ! if (brav==2) then
456 : ! ! Face Centered Lattice
457 : ! fac = half_sqrt2
458 : ! else if (brav==3) then
459 : ! ! Body Centered Cubic Lattice
460 : ! fac = half_sqrt3
461 : ! else if (brav==4) then
462 : ! ! Hexagonal Lattice
463 : ! fac = one
464 : ! end if
465 : ! Ifc%acell(1) = sqrt(rprimd(1,1)**2 + rprimd(2,1)**2 + rprimd(3,1)**2) / fac
466 : ! Ifc%acell(2) = sqrt(rprimd(1,2)**2 + rprimd(2,2)**2 + rprimd(3,2)**2) / fac
467 : ! Ifc%acell(3) = sqrt(rprimd(1,3)**2 + rprimd(2,3)**2 + rprimd(3,3)**2) / fac
468 : ! do ii=1,3
469 : ! Ifc%rprim(ii, 1) = rprimd(ii, 1) / Ifc%acell(1)
470 : ! Ifc%rprim(ii, 2) = rprimd(ii, 2) / Ifc%acell(2)
471 : ! Ifc%rprim(ii, 3) = rprimd(ii, 3) / Ifc%acell(3)
472 : ! Ifc%gprim(ii, 1) = crystal%gprimd(ii, 1) * Ifc%acell(1)
473 : ! Ifc%gprim(ii, 2) = crystal%gprimd(ii, 2) * Ifc%acell(2)
474 : ! Ifc%gprim(ii, 3) = crystal%gprimd(ii, 3) * Ifc%acell(3)
475 : ! end do
476 : !
477 : !end if
478 :
479 : ! Check if the rprim are coherent with the choice used in the interatomic forces generation
480 242 : call chkrp9(Ifc%brav,Ifc%rprim)
481 :
482 6050 : rprim = Ifc%rprim; gprim = Ifc%gprim
483 :
484 : ! -------------------------------------------------------------------------- !
485 : ! Compute dyewq0, the correction to be applied to the Ewald, see Eq.(71) of PRB55, 10355 (1997).
486 8549 : dyewq0 = zero
487 242 : if ((Ifc%dipdip==1.or.Ifc%dipquad==1.or.Ifc%quadquad==1).and. (Ifc%asr==1.or.Ifc%asr==2.or.Ifc%asr==6)) then
488 : ! Calculation of the non-analytical part for q=0
489 158 : sumg0=0
490 158 : qpt(:)=zero
491 790 : ABI_MALLOC(dyew,(2,3,natom,3,natom))
492 158 : if (sys_dim==1) then
493 156 : if (Ifc%dipquad==1.or.Ifc%quadquad==1) then
494 : call ewald9(Ifc%acell,dielt,dyew,Crystal%gmet,gprim,natom,qpt,Crystal%rmet,rprim,sumg0,Crystal%ucvol,&
495 3 : Crystal%xred,zeff,qdrp_cart,ifc%eta,option=ifc%ewald_option,dipquad=Ifc%dipquad,quadquad=Ifc%quadquad)
496 : else
497 : call ewald9(Ifc%acell,dielt,dyew,Crystal%gmet,gprim,natom,qpt,Crystal%rmet,rprim,sumg0,Crystal%ucvol,&
498 153 : Crystal%xred,zeff,qdrp_cart,ifc%eta, option=ifc%ewald_option)
499 : end if
500 : else
501 : ! 2D materials are embedded in a dielectric environment (typically vacuum in DFT calculations)
502 : ! which leads to a different long-range electrostatics than in 2D. The next routine allows
503 : ! to estimate it.
504 2 : call ewald9_2D(natom,ddb%acell,Crystal%xred,rprim,dielt,dyew,qpt,zeff,qdrp_cart,dielt_env,dielt_thick,sys_dim)
505 : end if
506 158 : call q0dy3_calc(natom,dyewq0,dyew,Ifc%asr)
507 158 : ABI_FREE(dyew)
508 : end if
509 :
510 : ! -------------------------------------------------------------------------- !
511 : ! Setup q-points
512 :
513 : ! Sample the Brillouin zone
514 242 : option=1
515 242 : qptrlatt = 0; qptrlatt(1,1)=ngqpt(1); qptrlatt(2,2)=ngqpt(2); qptrlatt(3,3)=ngqpt(3)
516 242 : mqpt=ngqpt(1)*ngqpt(2)*ngqpt(3)*nqshft
517 242 : if (Ifc%brav==2) mqpt=mqpt/2
518 242 : if (Ifc%brav==3) mqpt=mqpt/4
519 :
520 726 : ABI_MALLOC(qbz,(3,mqpt))
521 242 : call smpbz(Ifc%brav,ab_out,qptrlatt,mqpt,nqbz,nqshft,option,q1shft,qbz,prtout=prtout_)
522 :
523 : ! Find the irreducible zone (qibz)
524 726 : ABI_MALLOC(ibz2bz, (nqbz))
525 726 : ABI_MALLOC(wtq_folded, (nqbz))
526 484 : ABI_MALLOC(wtq, (nqbz))
527 26542 : wtq = one / nqbz ! Weights sum up to one
528 726 : ABI_MALLOC(bz2ibz_smap, (6, nqbz))
529 :
530 : ! FIXME: timrev1 should be set to 0 if TR cannot be used
531 : call symkpt(chksymbreak0,crystal%gmet,ibz2bz,iout0,qbz,nqbz,ifc%nqibz,crystal%nsym,&
532 242 : crystal%symrec,timrev1,wtq,wtq_folded, bz2ibz_smap, xmpi_comm_self)
533 :
534 242 : ABI_FREE(bz2ibz_smap)
535 :
536 726 : ABI_MALLOC(ifc%qibz, (3,ifc%nqibz))
537 726 : ABI_MALLOC(ifc%wtq, (ifc%nqibz))
538 3394 : do iq_ibz=1,ifc%nqibz
539 12608 : ifc%qibz(:,iq_ibz) = qbz(:, ibz2bz(iq_ibz))
540 3394 : ifc%wtq(iq_ibz) = wtq_folded(ibz2bz(iq_ibz))
541 : end do
542 242 : ABI_FREE(ibz2bz)
543 242 : ABI_FREE(wtq_folded)
544 242 : ABI_FREE(wtq)
545 :
546 : ! -------------------------------------------------------------------------- !
547 : ! Setup dynamical matrix and IFC matrix
548 :
549 1452 : ABI_MALLOC(Ifc%dynmat,(2,3,natom,3,natom,nqbz))
550 :
551 : ! This is needed to preserve the behavior of the old implementation with canonical coordinate.
552 242 : if (Ifc%brav == 1) then
553 2548 : gprim_tmp = Crystal%gprimd
554 196 : rprim_tmp = rprimd
555 : else
556 46 : gprim_tmp = gprim
557 46 : rprim_tmp = rprim
558 : endif
559 :
560 : ! Find symmetrical dynamical matrices
561 242 : if (.not.present(Ifc_coarse)) then
562 : ! Each q-point in the BZ mush be the symmetrical of one of the qpts in the ddb file.
563 : ! SP - gprimd and rprimd is required instead of gprim and rprim for non-diagonal supercells.
564 : call symdm9(ddb, &
565 : Ifc%dynmat,gprim_tmp,Crystal%indsym,mpert,natom,nqbz,nsym,rfmeth,rprim_tmp,qbz,&
566 241 : Crystal%symrec, Crystal%symrel, comm)
567 :
568 : else
569 : ! Symmetrize the qpts in the BZ using the q-points in the ddb.
570 : ! Then use Ifc_coarse to fill the missing entries with Fourier interpolated matrices.
571 : !
572 : ! TODO: The previous version of refineblk was hacking the DDB database to add the q-points in the **IBZ**
573 : ! Then D(q) was symmetrized in symdm9. This version avoids the symmetrization: the q-points
574 : ! in the BZ that are not in the coarse q-mesh are obtained by an explicit FT.
575 : ! This means that the final D(q) may break some symmetry in q-space if the FT does not preserve it.
576 : ! The most elegant approach would be to get D(q_ibz) via FT if q_ibz is not in the coarse mesh and then
577 : ! call symdm9 to get D(q) for each q point in the star of q_ibz.
578 1 : call wrtout(std_out,"Will fill missing qpoints in the full BZ using the coarse q-mesh")
579 :
580 : call symdm9(ddb, &
581 : Ifc%dynmat,gprim_tmp,Crystal%indsym,mpert,natom,nqbz,nsym,rfmeth,rprim_tmp,qbz,&
582 1 : Crystal%symrec,Crystal%symrel,comm, qmissing=qmissing)
583 :
584 : ! Compute dynamical matrix with Fourier interpolation on the coarse q-mesh.
585 1 : write(msg,"(a,i0,a)")"Will use Fourier interpolation to construct D(q) for ",size(qmissing)," q-points"
586 1 : call wrtout(std_out,msg)
587 :
588 4 : ABI_MALLOC(out_d2cart, (2,3,natom,3,natom))
589 97 : do ii=1,size(qmissing)
590 96 : iq_bz = qmissing(ii)
591 384 : qpt = qbz(:,iq_bz)
592 : ! TODO: check dipdip option and phase, but I think this is correct!
593 96 : call Ifc_coarse%fourq(Crystal,qpt,phfrq,displ_cart,out_d2cart=out_d2cart)
594 3361 : Ifc%dynmat(:,:,:,:,:,iq_bz) = out_d2cart
595 : end do
596 :
597 1 : ABI_FREE(qmissing)
598 1 : ABI_FREE(out_d2cart)
599 : end if
600 :
601 : ! OmegaSRLR: Store full dynamical matrix for decomposition into short- and long-range parts
602 968 : ABI_MALLOC(dynmatfull,(2,3,natom,3,natom,nqbz))
603 12537758 : dynmatfull=Ifc%dynmat
604 :
605 242 : if (Ifc%dipdip==1.or.Ifc%dipquad==1.or.Ifc%quadquad==1) then
606 : ! Take off the dipole-dipole part of the dynamical matrix
607 158 : if (prtout_) call wrtout(std_out, " Will extract the dipole-dipole part for every wavevector")
608 632 : ABI_MALLOC(dyew,(2,3,natom,3,natom))
609 :
610 11530 : do iqpt=1,nqbz
611 11372 : if (mod(iqpt, nprocs) /= my_rank) then ! mpi-parallelism
612 1988 : ifc%dynmat(:,:,:,:,:,iqpt) = zero; cycle
613 : end if
614 45472 : qpt(:)=qbz(:,iqpt)
615 11368 : sumg0=0
616 11368 : if (sys_dim==1) then
617 : ! 3D case
618 11031 : if (Ifc%dipquad==1.or.Ifc%quadquad==1) then
619 : call ewald9(Ifc%acell,dielt,dyew,Crystal%gmet,gprim,natom,qpt,Crystal%rmet,rprim,sumg0,Crystal%ucvol,&
620 : Crystal%xred,zeff,qdrp_cart, ifc%eta, &
621 48 : option=ifc%ewald_option,dipquad=Ifc%dipquad,quadquad=Ifc%quadquad)
622 : else
623 : call ewald9(Ifc%acell,dielt,dyew,Crystal%gmet,gprim,natom,qpt,Crystal%rmet,rprim,sumg0,Crystal%ucvol,&
624 : Crystal%xred,zeff,qdrp_cart, ifc%eta, &
625 10983 : option=ifc%ewald_option)
626 : end if
627 : else
628 : ! 2D case
629 337 : call ewald9_2D(natom,ddb%acell,Crystal%xred,rprim,dielt,dyew,qpt,zeff,qdrp_cart,one,dielt_thick,sys_dim)
630 : end if
631 11368 : call q0dy3_apply(natom,dyewq0,dyew)
632 11368 : plus=0
633 : ! Implement Eq.(76) of Gonze&Lee PRB 55, 10355 (1997) [[cite:Gonze1997a]], possibly generalized for quadrupoles
634 11526 : call nanal9(dyew,Ifc%dynmat,iqpt,natom,nqbz,plus)
635 : end do
636 :
637 158 : call xmpi_sum(ifc%dynmat, comm, ierr)
638 158 : ABI_FREE(dyew)
639 : !if (prtout_) call wrtout(std_out, " Done")
640 : end if
641 :
642 : ! OmegaSRLR: Store the short-range dynmat and compute long-range as difference
643 1452 : ABI_MALLOC(dynmat_sr,(2,3,natom,3,natom,nqbz))
644 968 : ABI_MALLOC(dynmat_lr,(2,3,natom,3,natom,nqbz))
645 12537758 : dynmat_sr=Ifc%dynmat
646 12537758 : dynmat_lr=dynmatfull-dynmat_sr
647 :
648 : ! Now, take care of the remaining part of the dynamical matrix
649 : ! Move to canonical normalized coordinates
650 242 : call canat9(Ifc%brav,natom,rcan,rprim,trans,Crystal%xred)
651 :
652 : ! Multiply the dynamical matrix by a phase shift
653 : option=1
654 242 : call dymfz9(Ifc%dynmat,natom,nqbz,gprim,option,qbz,trans)
655 :
656 : ! -------------------------------------------------------------------------- !
657 : ! Compute interatomic forces in real space
658 :
659 : ! Create the Big Box of R vectors in real space and compute the number of points (cells) in real space
660 242 : call make_bigbox(Ifc%brav,ifc_tmp%cell,ngqpt,nqshft,rprim,ifc_tmp%nrpt,ifc_tmp%rpt)
661 :
662 : ! Weights associated to these R points and to atomic pairs
663 1210 : ABI_MALLOC(ifc_tmp%wghatm, (natom, natom, ifc_tmp%nrpt))
664 :
665 : ! HM: this tolerance is highly dependent on the compilation/architecture
666 : ! numeric errors in the DDB text file. Try a few tolerances and check whether all the weights are found.
667 242 : toldist = tol8
668 243 : do while (toldist <= tol6)
669 : ! Note ngqpt(9) with intent(inout)!
670 : call wght9(Ifc%brav,gprim,natom,ngqpt,nqbz,nqshft,ifc_tmp%nrpt,q1shft,rcan,&
671 243 : ifc_tmp%rpt,rprimd,toldist,r_inscribed_sphere,ifc_tmp%wghatm,my_ierr)
672 243 : call xmpi_max(my_ierr, ierr, comm, ii)
673 243 : if (ierr > 0) toldist = toldist * 10
674 243 : if (ierr == 0) exit
675 : end do
676 :
677 242 : if (ierr > 0) then
678 : write(msg, '(3a,es14.4,2a,i0, 14a)' ) &
679 0 : 'The sum of the weight is not equal to nqpt.',ch10,&
680 0 : 'The sum of the weights is: ',sum(ifc_tmp%wghatm),ch10,&
681 0 : 'The number of q-points is: ',nqbz, ch10, &
682 0 : 'This might have several sources.',ch10,&
683 0 : 'If toldist is larger than 1.0e-8, the atom positions might be loose.',ch10,&
684 0 : 'and the q point weights not computed properly.',ch10,&
685 0 : 'Action: make input atomic positions more symmetric.',ch10,&
686 0 : 'Otherwise, you might increase "buffer" in m_dynmat.F90 see bigbx9 subroutine and recompile.',ch10,&
687 0 : 'Actually, this can also happen when ngqpt is 0 0 0,',ch10,&
688 0 : 'if abs(brav) /= 1, in this case you should change brav to 1. If brav is already set to 1 (default) try -1.'
689 0 : ABI_ERROR(msg)
690 : end if
691 :
692 : ! Fourier transform of the dynamical matrices (q-->R)
693 1452 : ABI_MALLOC(ifc_tmp%atmfrc, (3,natom,3,natom,ifc_tmp%nrpt))
694 242 : call ftifc_q2r(ifc_tmp%atmfrc,Ifc%dynmat,gprim,natom,nqbz,ifc_tmp%nrpt,ifc_tmp%rpt,qbz, comm)
695 :
696 : ! Eventually impose Acoustic Sum Rule on the interatomic forces
697 242 : if (Ifc%asr > 0) call asrif9(Ifc%asr,ifc_tmp%atmfrc,natom,ifc_tmp%nrpt,ifc_tmp%rpt,ifc_tmp%wghatm)
698 :
699 : ! The interatomic forces have been calculated
700 242 : write(msg, '(2a)')ch10,' The interatomic forces have been obtained '
701 634 : if (prtout_) call wrtout([std_out, ab_out], msg)
702 242 : call cwtime_report(" ifc_init1", cpu, wall, gflops)
703 :
704 : ! Apply cutoff on ifc if needed
705 242 : if (nsphere > 0 .or. abs(rifcsph) > tol10) then
706 3 : call wrtout(std_out, ' Apply cutoff on IFCs.')
707 3 : call wrtout(std_out, sjoin(" nsphere:", itoa(nsphere), ", rifcsph:", ftoa(rifcsph)))
708 3 : call wrtout(std_out, sjoin(" Radius of biggest sphere inscribed in the WS supercell: ", ftoa(r_inscribed_sphere)))
709 3 : call corsifc9(Ifc%acell,gprim,natom,ifc_tmp%nrpt,nsphere,rifcsph,rcan,rprim,ifc_tmp%rpt,rcut_min,ifc_tmp%wghatm)
710 3 : if (Ifc%asr > 0) then
711 3 : call wrtout(std_out, ' Enforcing ASR on cutoffed IFCs.')
712 3 : call asrif9(Ifc%asr,ifc_tmp%atmfrc,natom,ifc_tmp%nrpt,ifc_tmp%rpt,ifc_tmp%wghatm)
713 : end if
714 : end if
715 :
716 : ! Only conserve the necessary points in rpt: in the FT algorithm the order of the points is unimportant
717 : ! In the case of effective potential, we need to keep all the points
718 242 : Ifc%nrpt = 0
719 552916 : do irpt=1,ifc_tmp%nrpt
720 11065688 : if (sum(ifc_tmp%wghatm(:,:,irpt)) /= 0) Ifc%nrpt = Ifc%nrpt + 1
721 : end do
722 :
723 8689738 : ABI_CALLOC(Ifc%atmfrc,(3,natom,3,natom,Ifc%nrpt))
724 169430 : ABI_CALLOC(Ifc%rpt,(3,Ifc%nrpt))
725 169430 : ABI_CALLOC(Ifc%cell,(3,Ifc%nrpt))
726 851252 : ABI_CALLOC(Ifc%wghatm,(natom,natom,Ifc%nrpt))
727 8689738 : ABI_CALLOC(Ifc%short_atmfrc,(3,natom,3,natom,Ifc%nrpt))
728 8689496 : ABI_CALLOC(Ifc%ewald_atmfrc,(3,natom,3,natom,Ifc%nrpt))
729 :
730 : irpt_new = 1
731 552916 : do irpt = 1, ifc_tmp%nrpt
732 11065688 : if (sum(ifc_tmp%wghatm(:,:,irpt)) /= 0) then
733 8688528 : Ifc%atmfrc(:,:,:,:,irpt_new) = ifc_tmp%atmfrc(:,:,:,:,irpt)
734 168704 : Ifc%rpt(:,irpt_new) = ifc_tmp%rpt(:,irpt)
735 850042 : Ifc%wghatm(:,:,irpt_new) = ifc_tmp%wghatm(:,:,irpt)
736 168704 : Ifc%cell(:,irpt_new) = ifc_tmp%cell(:,irpt)
737 42176 : Ifc%r_inscribed_sphere = r_inscribed_sphere
738 42176 : irpt_new = irpt_new + 1
739 : end if
740 : end do
741 :
742 : !write(std_out,*)"nrpt before filter:", ifc_tmp%nrpt, ", after: ", ifc%nrpt
743 : !do irpt=1,ifc%nrpt
744 : ! write(std_out,*)ifc%rpt(:,irpt), (ifc%wghatm(ii,ii,irpt), ii=1,natom)
745 : !end do
746 :
747 242 : call ifc_tmp%free()
748 :
749 : ! -------------------------------------------------------------------------- !
750 :
751 : ! Copy other useful arrays.
752 3146 : Ifc%dielt = dielt
753 242 : Ifc%nqbz = nqbz
754 :
755 242 : call alloc_copy(rcan, Ifc%rcan)
756 242 : call alloc_copy(trans, Ifc%trans)
757 242 : call alloc_copy(dyewq0, Ifc%dyewq0)
758 242 : call alloc_copy(qbz(:,1:nqbz), Ifc%qbz)
759 242 : call alloc_copy(zeff, Ifc%zeff)
760 242 : call alloc_copy(qdrp_cart, Ifc%qdrp_cart)
761 242 : call alloc_copy(crystal%amu, Ifc%amu)
762 :
763 : ! Compute min/max ph frequency with ab-initio q-mesh.
764 242 : ifc%omega_minmax(1) = huge(one); ifc%omega_minmax(2) = -huge(one)
765 3394 : do iq_ibz=1,ifc%nqibz
766 3152 : if (mod(iq_ibz, nprocs) /= my_rank) cycle ! mpi-parallelism
767 3144 : call ifc%fourq(crystal, ifc%qibz(:,iq_ibz), phfrq, displ_cart)
768 38088 : ifc%omega_minmax(1) = min(ifc%omega_minmax(1), minval(phfrq))
769 38338 : ifc%omega_minmax(2) = max(ifc%omega_minmax(2), maxval(phfrq))
770 : end do
771 242 : xval = ifc%omega_minmax(1); call xmpi_min(xval, ifc%omega_minmax(1), comm, ierr)
772 242 : xval = ifc%omega_minmax(2); call xmpi_max(xval, ifc%omega_minmax(2), comm, ierr)
773 : ! Enlarge boundaries by 30 cm-1
774 242 : ifc%omega_minmax(1) = ifc%omega_minmax(1) - 30.0_dp/Ha_cmm1
775 242 : ifc%omega_minmax(2) = ifc%omega_minmax(2) + 30.0_dp/Ha_cmm1
776 :
777 : ! TODO (This is to be suppressed in a future version)
778 242 : if (prtsrlr == 1) then
779 : ! Check that the starting values are well reproduced.
780 0 : write(msg, '(2a)' )' mkifc9 : now check that the starting values ',&
781 0 : ' are reproduced after the use of interatomic forces '
782 0 : if (prtout_) call wrtout(std_out, msg)
783 0 : do iqpt=1,nqbz
784 0 : qpt(:)=Ifc%qbz(:,iqpt)
785 0 : call ifc%fourq(Crystal,qpt,phfrq,displ_cart,out_eigvec=eigvec)
786 :
787 : ! OmegaSRLR: Perform decomposition of dynamical matrix
788 : ! MG: FIXME I don't think the implementation is correct when q !=0
789 : if (prtsrlr==1) then
790 0 : call omega_decomp(Ifc%amu,natom,ntypat,Crystal%typat,dynmatfull,dynmat_sr,dynmat_lr,iqpt,nqbz,eigvec)
791 : end if
792 : ! Write the phonon frequencies (this is for checking purposes).
793 : ! Note: these phonon frequencies are not written on unit iout, only on unit std_out.
794 0 : call dfpt_prtph(displ_cart,0,enunit,-1,natom,phfrq,qphnrm,qpt)
795 : end do
796 : end if
797 :
798 : ! OmegaSRLR: deallocate memory used by dynmat decomposition
799 242 : ABI_FREE(dynmatfull)
800 242 : ABI_FREE(dynmat_sr)
801 242 : ABI_FREE(dynmat_lr)
802 242 : ABI_FREE(qbz)
803 :
804 242 : if (nsphere == -1) call ifc_autocutoff(ifc, crystal, comm)
805 :
806 242 : call cwtime_report(" ifc_init2", cpu, wall, gflops)
807 :
808 1210 : end subroutine ifc_init
809 : !!***
810 :
811 : !----------------------------------------------------------------------
812 :
813 : !!****f* m_ifc/ifc_from_file
814 : !! NAME
815 : !! ifc_from_file
816 : !!
817 : !! FUNCTION
818 : !! Initialize instance from file.
819 : !!
820 : !! INPUTS
821 : !!
822 : !! OUTPUT
823 : !!
824 : !! SOURCE
825 :
826 0 : subroutine ifc_from_file(ifc, dielt, filename,natom,ngqpt,nqshift,qshift,ucell_ddb,zeff,qdrp_cart,comm)
827 :
828 : !Arguments ------------------------------------
829 : !scalars
830 : class(ifc_type),intent(out) :: Ifc
831 : integer,intent(in) :: nqshift,comm
832 : integer,intent(inout) :: natom
833 : !arrays
834 : integer,intent(in) :: ngqpt(3)
835 : real(dp),intent(in) :: qshift(3,nqshift)
836 : character(len=*),intent(in) :: filename
837 : real(dp),intent(inout) :: dielt(3,3)
838 : real(dp),allocatable,intent(inout) :: zeff(:,:,:)
839 : real(dp),allocatable,intent(inout) :: qdrp_cart(:,:,:,:)
840 : type(crystal_t),intent(out) :: ucell_ddb
841 :
842 : !Local variables -------------------------
843 : !scalars
844 : integer :: dipdip,i,iblok,iblok_tmp
845 : logical :: file_exists
846 : character(len=500) :: msg
847 0 : type(ddb_type) :: ddb
848 0 : type(ddb_hdr_type) :: ddb_hdr
849 : !******************************************************************
850 :
851 : !check if ddb file exists
852 0 : inquire(file=filename, exist=file_exists)
853 :
854 0 : if (file_exists .eqv. .true.)then
855 : !Reading the ddb
856 0 : call ddb%from_file(filename, ddb_hdr, ucell_ddb, comm)
857 0 : call ddb_hdr%free()
858 0 : natom = ddb%natom
859 : else
860 0 : ABI_ERROR(sjoin("File:", filename, "is not present in the directory"))
861 : end if
862 :
863 : ! Get Dielectric Tensor and Effective Charges
864 : ! (initialized to one_3D and zero if the derivatives are not available in the DDB file)
865 0 : ABI_MALLOC(zeff,(3,3,natom))
866 0 : iblok = ddb%get_dielt_zeff(ucell_ddb,1,1,0,dielt,zeff)
867 :
868 : ! Try to get dielt, in case just the DDE are present
869 0 : if (iblok == 0) then
870 0 : iblok_tmp = ddb%get_dielt(1,dielt)
871 : end if
872 :
873 0 : ABI_MALLOC(qdrp_cart,(3,3,3,natom))
874 0 : iblok = ddb%get_quadrupoles(ddb_hdr%ddb_version,1,BLKTYP_d3E_xx,qdrp_cart)
875 :
876 : ! ifc to be calculated for interpolation
877 0 : write(msg, '(a,a,(80a),a,a,a,a)' ) ch10,('=',i=1,80),ch10,ch10,' Calculation of the interatomic forces ',ch10
878 0 : call wrtout([std_out, ab_out], msg)
879 0 : if ((maxval(abs(zeff)) .lt. tol10) .OR. (maxval(dielt) .gt. 100000.0)) then
880 0 : dipdip=0
881 : else
882 0 : dipdip=1
883 : end if
884 0 : call ifc%init(ucell_ddb,ddb,1,1,1,dipdip,1,ngqpt,nqshift,qshift,dielt,zeff,qdrp_cart,0,0.0_dp,0,1,1,comm)
885 :
886 : ! Free them all
887 0 : call ddb%free()
888 :
889 0 : end subroutine ifc_from_file
890 : !!***
891 :
892 : !----------------------------------------------------------------------
893 :
894 : !!****f* m_ifc/ifc_print
895 : !! NAME
896 : !! ifc_print
897 : !!
898 : !! FUNCTION
899 : !! Print info on the object
900 : !!
901 : !! INPUTS
902 : !! units=Unit numbers for output
903 : !! [header]=String to be printed as header for additional info.
904 : !! [prtvol]=Verbosity level.
905 : !!
906 : !! SOURCE
907 :
908 326 : subroutine ifc_print(ifc, units, header, prtvol)
909 :
910 : !Arguments ------------------------------------
911 : !scalars
912 : class(ifc_type),intent(in) :: ifc
913 : integer,intent(in) :: units(:)
914 : integer,optional,intent(in) :: prtvol
915 : character(len=*),optional,intent(in) :: header
916 :
917 : !Local variables-------------------------------
918 : integer :: my_prtvol,iatom,ii,idir
919 : character(len=500) :: msg
920 : ! *********************************************************************
921 :
922 163 : my_prtvol = 0; if (present(prtvol)) my_prtvol = prtvol
923 :
924 163 : msg = ' ==== Info on the interatomic force constants ==== '
925 163 : if (present(header)) msg = ' ==== '//trim(adjustl(header))//' ==== '
926 163 : call wrtout(units, msg)
927 :
928 163 : call wrtout(units,' Real(R)+Recip(G) space primitive vectors, cartesian coordinates (Bohr,Bohr^-1):')
929 652 : do ii=1,3
930 : write(msg,'(1x,a,i1,a,3f11.7,2x,a,i1,a,3f11.7)')&
931 489 : 'R(',ii,')=',ifc%rprim(:,ii),'G(',ii,')=',ifc%gprim(:,ii)
932 652 : call wrtout(units,msg)
933 : end do
934 163 : call wrtout(units, sjoin(" acell:", ltoa(ifc%acell)))
935 163 : call wrtout(units, sjoin(" Acoustic Sum Rule option (asr):", itoa(ifc%asr)))
936 163 : call wrtout(units, sjoin(" Option for the sampling of the BZ (brav):", itoa(ifc%brav)))
937 163 : call wrtout(units, sjoin(" Symmetrization flag (symdynmat):", itoa(ifc%symdynmat)))
938 163 : call wrtout(units, sjoin(" Dipole-dipole interaction flag (dipdip):", itoa(ifc%dipdip)))
939 163 : call wrtout(units, sjoin(" Dipole-quadrupole interaction flag (dipquad):", itoa(ifc%dipquad)))
940 163 : call wrtout(units, sjoin(" quadrupole-quadrupole interaction flag (quadquad):", itoa(ifc%quadquad)))
941 163 : call wrtout(units, sjoin(" Ewald option:", itoa(ifc%ewald_option)))
942 163 : call wrtout(units, sjoin(" Dielectric tensor: ", ch10, ltoa(reshape(ifc%dielt, [9]), fmt="f10.2")))
943 163 : call wrtout(units, " Effective charges:")
944 532 : do iatom=1,ifc%natom
945 532 : call wrtout(units, ltoa(reshape(ifc%zeff(:,:,iatom), [3*3]), fmt="f10.2"))
946 : end do
947 163 : call wrtout(units, " Quadrupolar terms:")
948 532 : do iatom=1,ifc%natom
949 1639 : do idir=1,3
950 1476 : call wrtout(units, ltoa(reshape(ifc%qdrp_cart(:,:,idir,iatom), [3*3]), fmt="f10.2"))
951 : end do
952 : end do
953 :
954 163 : call wrtout(units, sjoin(" Mass of the atoms (atomic mass unit): ", ltoa(ifc%amu)))
955 163 : call wrtout(units, sjoin(" Number of real-space points for IFC(R): ", itoa(ifc%nrpt)))
956 163 : call wrtout(units, sjoin(" Radius of biggest sphere inscribed in the WS supercell: ", ftoa(ifc%r_inscribed_sphere)))
957 163 : call wrtout(units, " ")
958 163 : call wrtout(units, " Q-mesh:")
959 163 : call wrtout(units, sjoin(" ngqpt:", ltoa(ifc%ngqpt),", nqshft:", itoa(ifc%nqshft)))
960 337 : do ii=1,ifc%nqshft
961 337 : call wrtout(units, sjoin(" ", ktoa(ifc%qshft(:,ii))))
962 : end do
963 :
964 163 : end subroutine ifc_print
965 : !!***
966 :
967 : !----------------------------------------------------------------------
968 :
969 : !!****f* m_ifc/ifc_fourq
970 : !! NAME
971 : !! ifc_fourq
972 : !!
973 : !! FUNCTION
974 : !! Compute the phonon frequencies and the group velocities at the specified q-point by performing
975 : !! a Fourier transform on the IFCs matrix in real space.
976 : !!
977 : !! INPUTS
978 : !! Crystal<type(crystal_t)> = Information on the crystalline structure.
979 : !! qpt(3)=q-point in reduced coordinates (unless nanaqdir is specified)
980 : !! [nanaqdir]=If present, the qpt will be treated as a vector specifying the
981 : !! direction in q-space along which the non-analytic behaviour of the dynamical
982 : !! matrix will be treated. Possible values:
983 : !! "cart" if qpt defines a direction in Cartesian coordinates
984 : !! "reduced" if qpt defines a direction in reduced coordinates
985 : !! [comm]: MPI communicator
986 : !!
987 : !! OUTPUT
988 : !! phfrq(3*natom) = Phonon frequencies in Hartree
989 : !! displ_cart(2,3,natom,3*natom) = Phonon displacement in Cartesian coordinates
990 : !! [out_d2cart(2,3,3*natom,3,3*natom)] = The (interpolated) dynamical matrix for this q-point
991 : !! [out_eigvec(2*3*natom*3*natom) = The (interpolated) eigenvectors of the dynamical matrix in Cartesian coords.
992 : !! [out_displ_red(2*3*natom*3*natom) = The (interpolated) displacement in reduced coordinates.
993 : !! [dwdq(3,3*natom)] = Group velocities i.e. d(omega(q))/dq in Cartesian coordinates.
994 : !!
995 : !! SOURCE
996 :
997 310925 : subroutine ifc_fourq(ifc, crystal, qpt, phfrq, displ_cart, &
998 : nanaqdir, comm, & ! Optional [in]
999 303185 : out_d2cart, out_eigvec, out_displ_red, dwdq) ! Optional [out]
1000 :
1001 : !Arguments ------------------------------------
1002 : !scalars
1003 : class(ifc_type),intent(in) :: Ifc
1004 : character(len=*),optional,intent(in) :: nanaqdir
1005 : type(crystal_t),intent(in) :: Crystal
1006 : integer,optional,intent(in) :: comm
1007 : !arrays
1008 : real(dp),intent(in) :: qpt(3)
1009 : real(dp),intent(out) :: displ_cart(2,3,Crystal%natom,3*Crystal%natom)
1010 : real(dp),intent(out) :: phfrq(3*Crystal%natom)
1011 : real(dp),optional,intent(out) :: out_d2cart(2,3,Crystal%natom,3,Crystal%natom)
1012 : real(dp),optional,intent(out) :: out_eigvec(2,3,Crystal%natom,3*Crystal%natom)
1013 : real(dp),optional,intent(out) :: out_displ_red(2,3,Crystal%natom,3*Crystal%natom)
1014 : real(dp),optional,intent(out) :: dwdq(3,3*crystal%natom)
1015 :
1016 : !Local variables-------------------------------
1017 : !scalars
1018 : integer :: natom, comm_
1019 : real(dp) :: qphnrm, eta
1020 : !arrays
1021 621850 : real(dp) :: my_qpt(3),eigvec(2,3,Crystal%natom,3*Crystal%natom),eigval(3*Crystal%natom)
1022 621850 : real(dp) :: d2cart(2,3,Ifc%mpert,3,Ifc%mpert),tsec(2)
1023 : ! ************************************************************************
1024 :
1025 : ! Keep track of total time spent.
1026 310925 : call timab(1748, 1, tsec)
1027 :
1028 310925 : natom = Crystal%natom
1029 : ! TODO: Rewrite and Parallelize ewald9 in gtdyn9
1030 310925 : comm_ = xmpi_comm_self; if (present(comm)) comm_ = comm
1031 :
1032 : ! Use my_qpt because dfpt_phfrq can change the q-point (very bad design)
1033 310925 : qphnrm = one; my_qpt = qpt
1034 :
1035 310925 : if (present(nanaqdir)) then
1036 : ! This will break backward compatibility because qpt is **always** in reduced coordinates.
1037 : ! while dfpt_phfrq assume cartesian coordinates !!!!!!!!!!!
1038 : ! It does not make sense to change API just to treat this particular case
1039 : ! We should **alwayse use q-points in reduced coordinates.
1040 18434 : qphnrm = zero
1041 2 : select case (nanaqdir)
1042 : case ("reduced")
1043 : ! Convert to Cartesian.
1044 26 : my_qpt = matmul(Crystal%gprimd, qpt)
1045 : case ("cart")
1046 0 : continue
1047 : case default
1048 18434 : ABI_ERROR(sjoin("Wrong value for nanaqdir:", nanaqdir))
1049 : end select
1050 : end if
1051 :
1052 : ! The dynamical matrix d2cart is calculated here:
1053 : call gtdyn9(Ifc%acell,Ifc%atmfrc,Ifc%dielt,Ifc%dipdip,Ifc%dyewq0,d2cart,Crystal%gmet,Ifc%gprim,Ifc%mpert,natom,&
1054 : Ifc%nrpt,qphnrm,my_qpt,Crystal%rmet,Ifc%rprim,Ifc%rpt,Ifc%trans,Crystal%ucvol,Ifc%wghatm,Crystal%xred,Ifc%zeff,&
1055 : Ifc%qdrp_cart,Ifc%ewald_option,eta,comm_, Ifc%sys_dim,dipquad=Ifc%dipquad,quadquad=Ifc%quadquad,&
1056 310925 : dielt_env=Ifc%dielt_env,dielt_thick=Ifc%dielt_thick)
1057 :
1058 : ! Calculate the eigenvectors and eigenvalues of the dynamical matrix
1059 : call dfpt_phfrq(Ifc%amu,displ_cart,d2cart,eigval,eigvec,Crystal%indsym,&
1060 : Ifc%mpert,Crystal%nsym,natom,Crystal%nsym,Crystal%ntypat,phfrq,qphnrm,my_qpt,&
1061 310925 : Crystal%rprimd,Ifc%symdynmat,Crystal%symrel,Crystal%symafm,Crystal%typat,Crystal%ucvol)
1062 :
1063 : ! OmegaSRLR: Perform decomposition of dynamical matrix
1064 : !if (srlr==1) call omega_decomp(amu,natom,ntypat,typat,dynmatfull,dynmatsr,dynmatlr,iqpt,nqpt,eigvec)
1065 :
1066 : ! Return the interpolated dynamical matrix and the eigenvector for this q-point
1067 10338725 : if (present(out_d2cart)) out_d2cart = d2cart(:,:,:natom,:,:natom)
1068 59527619 : if (present(out_eigvec)) out_eigvec = eigvec
1069 :
1070 : ! Return phonon displacement in reduced coordinates.
1071 310925 : if (present(out_displ_red)) call phdispl_cart2red(natom, crystal%gprimd, displ_cart, out_displ_red)
1072 :
1073 : ! Option to get vectors in reduced coordinates?
1074 : !call phdispl_cart2red(natom, crystal%gprimd, out_eigvec, out_eigvec_red)
1075 :
1076 : ! Compute group velocities.
1077 310925 : if (present(dwdq)) call ifc%get_dwdq(crystal, my_qpt, phfrq, eigvec, dwdq, comm_)
1078 :
1079 310925 : call timab(1748, 2, tsec)
1080 :
1081 341520 : end subroutine ifc_fourq
1082 : !!***
1083 :
1084 : !!****f* m_ifc/ifc_get_dcdq
1085 : !! NAME
1086 : !! ifc_get_dcdq
1087 : !!
1088 : !! FUNCTION
1089 : !! Compute the first-derivative of reciprocal space IFCs based on real-space IFCs.
1090 : !!
1091 : !! INPUTS
1092 : !! ifc<ifc_type>=Object containing the dynamical matrix and the IFCs.
1093 : !! crystal<crystal_t> = Information on the crystalline structure.
1094 : !! comm: MPI communicator
1095 : !!
1096 : !! OUTPUT
1097 : !! dcdq(2,3,natom,3,natom,3) = first derivatives of reciprocal-space IFCs in cartesian coordinates.
1098 : !!
1099 : !! NOTES
1100 : !! Using:
1101 : !!
1102 : !! dC/dq = \sum_b C(kappa alpha, kappa' beta) (0,b) (R_kappa'beta (b) - R_kappa alpha)
1103 : !!
1104 : !!
1105 : !! SOURCE
1106 :
1107 0 : subroutine ifc_get_dcdq(ifc, cryst, dcdq, dcdqdq, comm)
1108 :
1109 : !Arguments ------------------------------------
1110 : !scalars
1111 : class(ifc_type),intent(in) :: ifc
1112 : type(crystal_t),intent(in) :: cryst
1113 : integer,intent(in) :: comm!, dipdip
1114 : !arrays
1115 : real(dp), intent(inout) :: dcdq(3,cryst%natom,3,cryst%natom,3)
1116 : real(dp), intent(out) :: dcdqdq(3,cryst%natom,3,3,3)
1117 : !real(dp), intent(in) :: dyewq0(cryst%natom,3,cryst%natom,3)
1118 :
1119 : !Local variables-------------------------------
1120 : !scalars
1121 : integer :: ii,jj, mu, kk
1122 : real(dp) :: qpt(3)
1123 : !arrays
1124 0 : real(dp) :: dyntmp(2,3,cryst%natom,3,cryst%natom)
1125 0 : real(dp) :: dcdqcan(2,3,cryst%natom,3,cryst%natom,3)
1126 0 : real(dp) :: dcdqred(2,3,cryst%natom,3,cryst%natom,3)
1127 : ! ************************************************************************
1128 :
1129 : ABI_UNUSED((/comm/))
1130 0 : qpt=zero ; dyntmp = zero
1131 0 : dcdq=zero ; dcdqcan= zero ; dcdqred = zero
1132 0 : dcdqdq=zero
1133 : ! Compute the derivative based on the canonical coordinates
1134 0 : call dynmat_dq(qpt, cryst%natom, ifc%gprim, ifc%nrpt, ifc%rpt, ifc%atmfrc, ifc%wghatm, dcdqcan)
1135 : ! We have also to consider the phase factor introduced by the canonical coordinates (shift on
1136 : ! atoms coordinates in other unit cells)
1137 : ! Only a phase shift, but since we look at the derivative, this also have a contribution here
1138 0 : call ftifc_r2q(ifc%atmfrc,dyntmp, ifc%gprim, cryst%natom, 1, ifc%nrpt, ifc%rpt, qpt, ifc%wghatm, comm)
1139 : ! Move to reduced coordinates
1140 0 : do ii= 1,3
1141 0 : do jj=1,3
1142 0 : dcdqred(:,:,:,:,:,ii) = dcdqred(:,:,:,:,:,ii)+ifc%gprim(jj,ii)*dcdqcan(:,:,:,:,:,jj)
1143 : end do
1144 : end do
1145 0 : do ii=1,3
1146 0 : do mu =1,cryst%natom
1147 0 : do jj=1,3
1148 0 : do kk=1,3
1149 0 : dcdqred(2,ii,mu,jj,:,kk) = dcdqred(2,ii,mu,jj,:,kk)-dyntmp(1,ii,mu,jj,:)*(ifc%trans(kk,:)-ifc%trans(kk,mu))
1150 0 : dcdqred(2,ii,mu,jj,:,kk) = dcdqred(2,ii,mu,jj,:,kk)+dyntmp(1,ii,mu,jj,:)*(cryst%xred(kk,:)-cryst%xred(kk,mu))
1151 : end do
1152 : end do
1153 : end do
1154 : end do
1155 0 : if (ifc%dipdip==1) ABI_WARNING("2D electrostatics + IFCs derivatives from real-space IFCs not fully tested")
1156 : ! If long-range electrostatics are treated separately, we need to add here the contribution to the torque
1157 : ! Note that there is only a contribution coming from the zone-center IFCs here
1158 : ! In 3D the derivative diverges; in 2D it does not, so we can add it with finite difference
1159 : !if (dipdip==1 .and. (ifc%sys_dim==2 .or. ifc%sys_dim==3 .or. ifc%sys_dim==4)) then
1160 : !do kk=1,3
1161 : !if (kk /= ifc%sys_dim-1) then
1162 : !do jj=-1,1,2
1163 : !qtmp(:) = zero
1164 : !qtmp(kk)=jj*0.01_dp
1165 : !call ewald9_2D(cryst%natom,ifc%acell,cryst%xred,ifc%rprim,ifc%dielt,dyew_tmp,&
1166 : ! qtmp,ifc%zeff,ifc%qdrp_cart,ifc%dielt_env,ifc%dielt_thick,ifc%sys_dim)
1167 : !dcdqred(2,:,:,:,:,kk) = dcdqred(2,:,:,:,:,kk)+half*dyew_tmp(2,:,:,:,:)/qtmp(kk)
1168 : ! end do
1169 : ! end if
1170 : ! end do
1171 : !end if
1172 : ! Move to cartesian coordinates
1173 0 : do ii=1,3
1174 0 : do jj=1,3
1175 0 : dcdq(:,:,:,:,ii)=dcdq(:,:,:,:,ii)+dcdqred(2,:,:,:,:,jj)*cryst%rprimd(ii,jj)
1176 : end do
1177 : end do
1178 :
1179 0 : end subroutine ifc_get_dcdq
1180 : !!***
1181 :
1182 : !!****f* m_ifc/ifc_get_dwdq
1183 : !! NAME
1184 : !! ifc_get_dwdq
1185 : !!
1186 : !! FUNCTION
1187 : !! Compute phonon group velocities at an arbitrary q-point.
1188 : !!
1189 : !! INPUTS
1190 : !! crystal<crystal_t> = Information on the crystalline structure.
1191 : !! qpt(3)=q-point in reduced coordinates.
1192 : !! eigvec(2*3*natom*3*natom) = The eigenvectors of the dynamical matrix.
1193 : !! comm: MPI communicator
1194 : !!
1195 : !! OUTPUT
1196 : !! dwdq(3,3*natom) = Group velocities e.g. d(omega(q))/dq in Cartesian coordinates.
1197 : !!
1198 : !! NOTES
1199 : !! Using:
1200 : !!
1201 : !! D(q) u(q,nu) = w(q, nu)**2 and <u(q,nu) | u(q,nu')> = \delta_{nu, nu'}
1202 : !!
1203 : !! one can show, using the Hellman-Feynman theorem, that:
1204 : !!
1205 : !! \nabla_q w(q, nu) = 1/(2 w(q, nu)) <u(q, nu)| \nabla_q D(q) | u(q, nu)>
1206 : !!
1207 : !! SOURCE
1208 :
1209 227649 : subroutine ifc_get_dwdq(ifc, cryst, qpt, phfrq, eigvec, dwdq, comm)
1210 :
1211 : !Arguments ------------------------------------
1212 : !scalars
1213 : class(ifc_type),intent(in) :: ifc
1214 : type(crystal_t),intent(in) :: cryst
1215 : integer,intent(in) :: comm
1216 : !arrays
1217 : real(dp),intent(in) :: qpt(3)
1218 : real(dp),intent(in) :: phfrq(3*cryst%natom)
1219 : real(dp),intent(in) :: eigvec(2,3*cryst%natom,3*cryst%natom)
1220 : real(dp),intent(out) :: dwdq(3,3*cryst%natom)
1221 :
1222 : !Local variables-------------------------------
1223 : !scalars
1224 : !integer,save :: enough=0
1225 : integer,parameter :: nqpt1=1,option2=2,sumg0=0
1226 : integer :: ii,nu,natom3,jj
1227 : real(dp) :: hh, eta
1228 : !arrays
1229 455298 : real(dp) :: dddq(2,3*cryst%natom,3*cryst%natom,3),dot(2),qfd(3)
1230 455298 : real(dp) :: omat(2,3*cryst%natom,3*cryst%natom)
1231 455298 : real(dp) :: dyew(2,3*cryst%natom,3*cryst%natom)
1232 : ! ************************************************************************
1233 :
1234 : ABI_UNUSED((/comm/))
1235 227649 : natom3 = cryst%natom * 3
1236 :
1237 : ! Generate the analytical part from the interatomic forces
1238 227649 : call dynmat_dq(qpt, cryst%natom, ifc%gprim, ifc%nrpt, ifc%rpt, ifc%atmfrc, ifc%wghatm, dddq)
1239 :
1240 : ! The analytical dynamical matrix dq has been generated
1241 : ! in the normalized canonical coordinate system. Now, the
1242 : ! phase is modified, in order to recover the usual (xred) coordinate of atoms.
1243 910596 : do ii=1,3
1244 682947 : call dymfz9(dddq(:,:,:,ii), cryst%natom, nqpt1, ifc%gprim, option2, qpt, ifc%trans)
1245 124232646 : dddq(:,:,:,ii) = dddq(:,:,:,ii) * ifc%acell(ii)
1246 : end do
1247 :
1248 227649 : if (ifc%dipdip == 1.or.ifc%dipquad == 1.or.ifc%quadquad == 1) then
1249 : ! Add the gradient of the non-analytical part.
1250 : ! Note that dddq is in cartesian coordinates.
1251 : ! For the time being, the gradient is computed with finite difference and step hh.
1252 : ! TODO: should generalize ewald9 to compute dq.
1253 : !enough = enough + 1
1254 : !if (enough <= 5) ABI_WARNING("phonon velocities with dipdip==1 not yet tested.")
1255 : hh = 0.01_dp
1256 69388 : do ii=1,3
1257 331731 : do jj=-1,1,2
1258 : ! qcart --> qred
1259 104082 : qfd = zero; qfd(ii) = jj
1260 1665312 : qfd = matmul(cryst%rprimd, qfd); qfd = qfd / normv(qfd, cryst%gmet, "G")
1261 : !write(std_out,*)"normv:",normv(qfd, cryst%gmet, "G")
1262 416328 : qfd = qpt + hh * qfd
1263 :
1264 : call ewald9(ifc%acell,ifc%dielt,dyew,cryst%gmet,ifc%gprim,cryst%natom,qfd,&
1265 : cryst%rmet,ifc%rprim,sumg0,cryst%ucvol,cryst%xred,ifc%zeff,ifc%qdrp_cart, eta, &
1266 104082 : ifc%ewald_option,dipquad=ifc%dipquad,quadquad=ifc%quadquad)
1267 104082 : call q0dy3_apply(cryst%natom,ifc%dyewq0,dyew)
1268 11943351 : dddq(:,:,:,ii) = dddq(:,:,:,ii) + (jj * half / hh) * dyew
1269 : end do
1270 : end do
1271 : end if
1272 :
1273 910596 : do ii=1,3
1274 910596 : call massmult_and_breaksym(cryst%natom, cryst%ntypat, cryst%typat, ifc%amu, dddq(:,:,:,ii))
1275 : end do
1276 :
1277 : ! Compute 1/(2w(q)) <u(q)|dD(q)/dq|u(q)>
1278 910596 : do ii=1,3
1279 682947 : call zgemm('N','N',natom3,natom3,natom3,cone,dddq(:,:,:,ii),natom3,eigvec,natom3,czero,omat,natom3)
1280 5145393 : do nu=1,natom3
1281 4917744 : if (abs(phfrq(nu)) > tol12) then
1282 4234005 : dot = cg_zdotc(natom3, eigvec(1,1,nu), omat(1,1,nu))
1283 : ! abs(w) is needed to get the correct derivative if we have a purely imaginary solution.
1284 4234005 : dwdq(ii, nu) = dot(1) / (two * abs(phfrq(nu)))
1285 : else
1286 792 : dwdq(ii, nu) = zero
1287 : end if
1288 : end do
1289 : end do
1290 :
1291 227649 : end subroutine ifc_get_dwdq
1292 : !!***
1293 :
1294 : !----------------------------------------------------------------------
1295 :
1296 : !!****f* m_ifc/ifc_get_phmesh
1297 : !! NAME
1298 : !! ifc_get_phmesh
1299 : !!
1300 : !! FUNCTION
1301 : !! Build linear mesh for phonons.
1302 : !!
1303 : !! INPUTS
1304 : !!
1305 : !! OUTPUT
1306 : !!
1307 : !! SOURCE
1308 :
1309 65 : subroutine ifc_get_phmesh(ifc, ph_wstep, phmesh_size, phmesh)
1310 :
1311 : !Arguments ------------------------------------
1312 : !scalars
1313 : class(ifc_type),intent(in) :: ifc
1314 : real(dp),intent(in) :: ph_wstep
1315 : integer,intent(out) :: phmesh_size
1316 : !arrays
1317 : real(dp),allocatable,intent(out) :: phmesh(:)
1318 : !******************************************************************
1319 :
1320 65 : phmesh_size = nint((ifc%omega_minmax(2) - ifc%omega_minmax(1) ) / ph_wstep) + 1
1321 195 : ABI_MALLOC(phmesh, (phmesh_size))
1322 65 : phmesh = arth(ifc%omega_minmax(1), ph_wstep, phmesh_size)
1323 :
1324 65 : end subroutine ifc_get_phmesh
1325 : !!***
1326 :
1327 : !----------------------------------------------------------------------
1328 :
1329 : !!****f* m_ifc/ifc_speedofsound
1330 : !!
1331 : !! NAME
1332 : !! ifc_speedofsound
1333 : !!
1334 : !! FUNCTION
1335 : !! Calculate the speed of sound by averaging the phonon group velocities of the
1336 : !! three acoustic modes on a small sphere of radius qrad centered around Gamma.
1337 : !! Perform spherical integration with Lebedev-Laikov grids
1338 : !!
1339 : !! INPUTS
1340 : !! crystal<crystal_t> = Information on the crystalline structure.
1341 : !! qrad_tolkms(2):
1342 : !! qrad=Radius of the sphere in reciprocal space
1343 : !! atols_kms=Absolute tolerance in kilometer/second. The code generates spherical meshes
1344 : !! until the results are converged twice within atols_kms.
1345 : !! ncid=the id of the open NetCDF file. Use nctk_noid to disable netcdf output.
1346 : !! comm=MPI communicator.
1347 : !!
1348 : !! OUTPUT
1349 : !!
1350 : !! SOURCE
1351 :
1352 1 : subroutine ifc_speedofsound(ifc, crystal, qrad_tolkms, ncid, comm)
1353 :
1354 : !Arguments -------------------------------
1355 : !scalars
1356 : class(ifc_type),intent(in) :: ifc
1357 : integer,intent(in) :: comm,ncid
1358 : type(crystal_t),intent(in) :: crystal
1359 : !arrays
1360 : real(dp),intent(in) :: qrad_tolkms(2)
1361 :
1362 : !Local variables -------------------------
1363 : !scalars
1364 : integer,parameter :: master=0
1365 : integer :: ii,nu,igrid,my_rank,nprocs,ierr,converged,npts,num_negw,vs_ierr,ncerr
1366 : integer :: iatom,iatref,num_acoustic,isacoustic
1367 : real(dp) :: min_negw,cpu,wall,gflops, qrad,tolkms,diff
1368 : character(len=500) :: msg
1369 1 : type(lebedev_t) :: lgrid
1370 : !arrays
1371 : integer :: asnu(3)
1372 : real(dp) :: qred(3),qvers_cart(3),qvers_red(3),quad(3),prev_quad(3),vs(7,3)
1373 3 : real(dp) :: phfrqs(3*crystal%natom),dwdq(3,3*crystal%natom)
1374 2 : real(dp) :: displ_cart(2,3*crystal%natom,3*crystal%natom),eigvec(2,3*crystal%natom,3*crystal%natom)
1375 : ! *********************************************************************
1376 :
1377 1 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
1378 :
1379 1 : if (ifc%asr == 0) ABI_WARNING("Computing speed of sound with asr == 0! Use asr > 0")
1380 1 : qrad = qrad_tolkms(1); tolkms = qrad_tolkms(2)
1381 1 : ABI_CHECK(qrad > zero, "vs_qrad <= 0")
1382 1 : ABI_CHECK(tolkms > zero, "vs_tolkms <= 0")
1383 :
1384 1 : call cwtime(cpu, wall, gflops, "start")
1385 :
1386 : ! Find the index of the first acoustic modes (needed to handle systems with unstable branches at Gamma
1387 : ! In this case, indeed, we end up computing derivatives for the wrong branches, the integration becomes
1388 : ! unstable and difficult to converge.
1389 1 : qred = zero
1390 1 : call ifc%fourq(crystal, qred, phfrqs, displ_cart,out_eigvec=eigvec)
1391 : !write(std_out,*)"omega(q==Gamma): ",phfrqs
1392 :
1393 1 : num_acoustic = 0
1394 3 : do nu = 1, 3*crystal%natom
1395 : ! Check if this mode is acoustic like: scalar product of all displacement vectors are collinear
1396 12 : isacoustic = 1
1397 : ! Find reference atom with non-zero displacement
1398 12 : do iatom=1,crystal%natom
1399 93 : if(sum(displ_cart(:,(iatom-1)*3+1:(iatom-1)*3+3,nu)**2) >tol16)iatref=iatom
1400 : end do
1401 : ! Now compute scalar product, and check they are all positive
1402 12 : do iatom = 1, crystal%natom
1403 93 : if (sum(eigvec(:,(iatom-1)*3+1:(iatom-1)*3+3, nu)*eigvec(:,(iatref-1)*3+1:(iatref-1)*3+3, nu)) < tol16 ) isacoustic = 0
1404 : end do
1405 4 : if (isacoustic == 1) then
1406 3 : num_acoustic=num_acoustic+1
1407 3 : asnu(num_acoustic)=nu
1408 3 : if (num_acoustic==3) exit
1409 : end if
1410 : end do
1411 :
1412 1 : ABI_CHECK(num_acoustic == 3, sjoin("Wrong number of acoustic modes:", itoa(num_acoustic)))
1413 :
1414 1 : write(std_out,"(a,3i2,a)") "The bands with indices ",asnu(:)," will be used to calculate the sound velocities"
1415 :
1416 : ! Speed of sound along reduced directions.
1417 7 : do ii=1,6
1418 6 : qred = zero; qred(MOD(ii-1,3)+1) = one
1419 42 : if (ii >= 4 .and. ii <= 6) qred = matmul(crystal%rprimd, qred) ! Cartesian directions.
1420 : !if (ii >= 7 .and. ii <= 9) qred = matmul(crystal%rprimd, qred) ! Cartesian directions.
1421 24 : qvers_red = (qred / normv(qred, crystal%gmet, "G"))
1422 24 : qred = qrad * qvers_red
1423 : !write(std_out,*)"dir",normv(qred, crystal%gmet, "G"), qrad
1424 6 : call ifc%fourq(crystal, qred, phfrqs, displ_cart, dwdq=dwdq)
1425 :
1426 24 : do nu=1,3
1427 78 : vs(ii, nu) = sqrt(sum(dwdq(1:3,asnu(nu)) ** 2)) * Bohr_meter * 0.001_dp / Time_Sec
1428 : end do
1429 6 : write(std_out,"(a,3es12.4,a)")" ||vs(nu)||:",vs(ii,:), " [km/s]"
1430 :
1431 96 : qvers_cart = matmul(crystal%gprimd, qvers_red) * two_pi
1432 24 : do nu=1,3
1433 78 : vs(ii, nu) = dot_product(dwdq(1:3,asnu(nu)), qvers_cart) * Bohr_meter * 0.001_dp / Time_Sec
1434 : end do
1435 7 : write(std_out,"(a,3es12.4,a)")" <q|vs(nu)>:",vs(ii,:), " [km/s]"
1436 :
1437 : !do nu=1,3
1438 : ! write(std_out,"(a,3es12.4,a)")" vs(nu)_vect_red:",&
1439 : ! matmul(crystal%gprimd, dwdq(1:3,asnu(nu))) * Bohr_meter * 0.001_dp / Time_Sec, " [km/s]"
1440 : !end do
1441 : end do
1442 :
1443 : ! Spherical average with Lebedev-Laikov grids.
1444 1 : converged = 0
1445 5 : do igrid=1,lebedev_ngrids
1446 5 : call lgrid%from_index(igrid)
1447 5 : npts = lgrid%npts; quad = zero; num_negw = 0; min_negw = zero
1448 139 : do ii=1,npts
1449 134 : if (mod(ii, nprocs) /= my_rank) cycle ! mpi-parallelism
1450 :
1451 : ! Build q-point on sphere of radius qrad. qcart --> qred
1452 1742 : qred = matmul(crystal%rprimd, lgrid%versors(:, ii))
1453 536 : qred = qrad * (qred / normv(qred, crystal%gmet, "G"))
1454 : !write(std_out,*)"lebe",normv(qred, crystal%gmet, "G"), qrad
1455 134 : call ifc%fourq(crystal, qred, phfrqs, displ_cart, dwdq=dwdq)
1456 536 : if (any(phfrqs(asnu) < -tol8)) then
1457 0 : num_negw = num_negw + 1; min_negw = min(min_negw, minval(phfrqs(asnu)))
1458 : end if
1459 :
1460 541 : do nu=1,3
1461 1742 : quad(nu) = quad(nu) + lgrid%weights(ii) * sqrt(sum(dwdq(1:3,asnu(nu)) ** 2))
1462 : !quad(nu) = quad(nu) + lgrid%weights(ii) * abs(dot_product(lgrid%versors(:,ii), dwdq(:,asnu(nu))))
1463 : end do
1464 : end do
1465 :
1466 : ! Will use km/sec unit for echo purposes
1467 20 : quad = quad * Bohr_meter * 0.001_dp / Time_Sec
1468 5 : call xmpi_sum(quad, comm, ierr)
1469 5 : call xmpi_sum(num_negw, comm, ierr)
1470 5 : call lgrid%free()
1471 :
1472 : write(std_out,'(2(a,i6),a,3es12.4,a,es12.4)') &
1473 20 : " Lebedev-Laikov grid: ",igrid,", npts: ", npts, " vs_sphavg(ac_modes): ",quad, " <vs>: ",sum(quad)/3
1474 :
1475 5 : if (igrid > 1) then
1476 : diff = zero
1477 16 : do nu=1,3
1478 16 : diff = diff + abs(quad(nu) - prev_quad(nu)) / 3
1479 : end do
1480 : !if (abs(sum(quad - prev_quad)/3) < tolkms) then
1481 4 : if (diff < tolkms) then
1482 2 : converged = converged + 1
1483 : else
1484 : converged = 0
1485 : end if
1486 : end if
1487 5 : prev_quad = quad
1488 20 : vs(7, :) = quad
1489 10 : if (converged == 2) exit
1490 : end do ! igrid
1491 :
1492 1 : if (my_rank == master) then
1493 : ! vs_err: 1 if not converged, < 0 if negative freqs, == 0 if success.
1494 1 : vs_ierr = 0
1495 4 : do ii=1,3
1496 4 : write(ab_out,"(a,3es12.4,a,i1)")" Speed of sound:",vs(ii,:)," [km/s] along reduced direction: ",ii
1497 : end do
1498 : write(ab_out,'(2(a,es12.4),a,i0)') &
1499 1 : " Lebedev-Laikov integration with qradius: ", qrad, " tolkms: ",tolkms, " [km/s], npts: ", npts
1500 4 : write(ab_out,"(a,3es12.4,a,es12.4)")" Spherical average:",vs(7,:)," [km/s], ",sum(vs(7,:))/3
1501 1 : if (converged /= 2) then
1502 0 : vs_ierr = 1
1503 0 : write(msg,'(a,es12.4,a)')" WARNING: Results are not converged within: ",tolkms, " [km/s]"
1504 0 : call wrtout(ab_out, msg)
1505 0 : ABI_WARNING(msg)
1506 : end if
1507 1 : if (num_negw > 0) then
1508 0 : vs_ierr = -num_negw
1509 : write(msg,'(a,i0,a,es12.4,3a)') &
1510 0 : " WARNING: Detected ",num_negw, " negative frequencies. Minimum was: ",min_negw * Ha_meV, "[meV]",ch10,&
1511 0 : " Speed of sound could be wrong"
1512 0 : call wrtout(ab_out, msg)
1513 0 : ABI_WARNING(msg)
1514 : end if
1515 :
1516 : ! Dump results to netcdf file.
1517 1 : if (ncid /= nctk_noid) then
1518 2 : ncerr = nctk_def_arrays(ncid, [nctkarr_t("vsound", "dp", "seven, three")], defmode=.True.)
1519 1 : NCF_CHECK(ncerr)
1520 2 : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "vsound_ierr"])
1521 1 : NCF_CHECK(ncerr)
1522 3 : ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "vsound_qrad", "vsound_tolkms"])
1523 1 : NCF_CHECK(ncerr)
1524 2 : ncerr = nctk_def_arrays(ncid, [nctkarr_t("asnu", "i", "three")], defmode=.True.)
1525 1 : NCF_CHECK(ncerr)
1526 : ! Write data.
1527 1 : NCF_CHECK(nctk_set_datamode(ncid))
1528 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vsound_ierr"), vs_ierr))
1529 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vsound_qrad"), qrad))
1530 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vsound_tolkms"), tolkms))
1531 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vsound"), vs))
1532 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "asnu"), asnu))
1533 : end if
1534 : end if
1535 :
1536 1 : call cwtime_report(" ifc_speedofsound", cpu, wall, gflops)
1537 :
1538 1 : end subroutine ifc_speedofsound
1539 : !!***
1540 :
1541 : !----------------------------------------------------------------------
1542 :
1543 : !!****f* m_ifc/ifc_autocutoff
1544 : !! NAME
1545 : !! ifc_autocutoff
1546 : !!
1547 : !! FUNCTION
1548 : !! Find the value of nsphere that gives non-negative frequencies around Gamma
1549 : !! in a small sphere of radius qrad.
1550 : !! Use bisection to reduce the number of attempts although there's no guarantee
1551 : !! that the number of negative frequencies is monotonic.
1552 : !!
1553 : !! INPUTS
1554 : !! crystal<crystal_t> = Information on the crystalline structure.
1555 : !! comm=MPI communicator
1556 : !!
1557 : !! SIDE EFFECTS
1558 : !! ifc%wghatm(natom,natom,nrpt) = Weights associated to a pair of atoms and to a R vector
1559 : !! with the last cutoff found by the bisection algorithm applied.
1560 : !! ifc%atmfrc(2,3,natom,3,natom,nrpt)= ASR-imposed Interatomic Forces
1561 : !!
1562 : !! SOURCE
1563 :
1564 1 : subroutine ifc_autocutoff(ifc, crystal, comm)
1565 :
1566 : !Arguments ------------------------------------
1567 : !scalars
1568 : class(ifc_type),intent(inout) :: ifc
1569 : type(crystal_t),intent(in) :: crystal
1570 : integer,intent(in) :: comm
1571 :
1572 : !Local variables-------------------------------
1573 : !scalars
1574 : integer,parameter :: master=0
1575 : integer :: iq_ibz,ierr,my_rank,nprocs,ii,nsphere,num_negw,jl,ju,jm,natom,nrpt
1576 : real(dp),parameter :: rifcsph0 = zero
1577 : real(dp) :: adiff,qrad,min_negw,xval,rcut_min
1578 1 : type(lebedev_t) :: lgrid
1579 : !arrays
1580 2 : real(dp) :: displ_cart(2*3*ifc%natom*3*ifc%natom)
1581 2 : real(dp) :: qred(3),qred_vers(3),phfrqs(3*ifc%natom) !,dwdq(3,3*ifc%natom)
1582 1 : real(dp),allocatable :: ref_phfrq(:,:),cut_phfrq(:,:), save_wghatm(:,:,:),save_atmfrc(:,:,:,:,:)
1583 : ! *********************************************************************
1584 :
1585 1 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
1586 1 : natom = ifc%natom; nrpt = ifc%nrpt
1587 :
1588 : ! Compute frequencies on the ab-initio q-mesh without cutoff.
1589 1828 : ABI_CALLOC(ref_phfrq, (3*natom, ifc%nqibz))
1590 97 : do iq_ibz=1,ifc%nqibz
1591 96 : if (mod(iq_ibz, nprocs) /= my_rank) cycle ! mpi-parallelism
1592 97 : call ifc%fourq(crystal, ifc%qibz(:,iq_ibz), ref_phfrq(:,iq_ibz), displ_cart)
1593 : end do
1594 1 : call xmpi_sum(ref_phfrq, comm, ierr)
1595 :
1596 5 : ABI_MALLOC(save_wghatm, (natom,natom,nrpt))
1597 6 : ABI_MALLOC(save_atmfrc, (3,natom,3,natom,ifc%nrpt))
1598 512503 : save_wghatm = ifc%wghatm; save_atmfrc = ifc%atmfrc
1599 :
1600 4 : ABI_MALLOC(cut_phfrq, (3*natom, ifc%nqibz))
1601 1 : qrad = 0.01
1602 1 : call lgrid%from_index(16)
1603 :
1604 1 : if (my_rank == master) then
1605 1 : write(ab_out, "(a)")" Apply cutoff on IFCs. Using bisection algorithm to find initial guess for nsphere."
1606 1 : write(ab_out, "(a,i0)")" Maximum number of atom-centered spheres: ",natom * nrpt
1607 1 : write(ab_out, "(a,i0,a,f5.3)")" Using Lebedev-Laikov grid with npts: ",lgrid%npts, ", qrad: ",qrad
1608 1 : write(ab_out, "(/,a)")" <adiff>: Average difference between ab-initio frequencies and frequencies with cutoff."
1609 1 : write(ab_out, "(a)")" num_negw: Number of negative freqs detected in small sphere around Gamma."
1610 1 : write(ab_out, "(a)")" min_negw: Min negative frequency on the small sphere."
1611 1 : write(ab_out, "(a,/,/)")" rifcsph: Effective cutoff radius corresponding to nsphere."
1612 1 : write(ab_out, "(a)")" nsphere <adiff>[meV] num_negw min_negw[meV] rifcsph"
1613 : end if
1614 :
1615 1 : jl = 0; ju = natom * nrpt + 1 ! Initialize lower and upper limits.
1616 39 : do
1617 14 : if (ju - jl <= 1) then
1618 : exit
1619 : end if
1620 13 : jm = (ju + jl) / 2 ! Compute a midpoint
1621 13 : nsphere = jm
1622 :
1623 6662539 : ifc%wghatm = save_wghatm; ifc%atmfrc = save_atmfrc
1624 13 : call corsifc9(ifc%acell,ifc%gprim, natom, nrpt,nsphere,rifcsph0,ifc%rcan,ifc%rprim,ifc%rpt,rcut_min,ifc%wghatm)
1625 13 : if (ifc%asr > 0) call asrif9(ifc%asr,ifc%atmfrc,ifc%natom,ifc%nrpt,ifc%rpt,ifc%wghatm)
1626 :
1627 23725 : cut_phfrq = zero
1628 1261 : do iq_ibz=1,ifc%nqibz
1629 1248 : if (mod(iq_ibz, nprocs) /= my_rank) cycle ! mpi-parallelism
1630 1261 : call ifc%fourq(crystal, ifc%qibz(:,iq_ibz), cut_phfrq(:,iq_ibz), displ_cart)
1631 : !write(std_out,*)cut_phfrq(1,iq_ibz),ref_phfrq(1,iq_ibz)
1632 : end do
1633 13 : call xmpi_sum(cut_phfrq, comm, ierr)
1634 :
1635 : ! Test whether there are negative frequencies around gamma, including reciprocal lattice vectors.
1636 13 : num_negw = 0; min_negw = zero
1637 5694 : do ii=1,lgrid%npts+3
1638 5681 : if (mod(ii, nprocs) /= my_rank) cycle ! mpi-parallelism
1639 5681 : if (ii <= 3) then
1640 39 : qred = zero; qred(ii) = one
1641 : else
1642 73346 : qred = matmul(crystal%rprimd, lgrid%versors(:, ii-3))
1643 : end if
1644 22724 : qred_vers = (qred / normv(qred, crystal%gmet, "G"))
1645 22724 : qred = qrad * qred_vers
1646 5681 : call ifc%fourq(crystal, qred, phfrqs, displ_cart) !, dwdq=dwdq)
1647 96036 : if (any(phfrqs < +tol8)) then
1648 13240 : num_negw = num_negw + 1; min_negw = min(min_negw, minval(phfrqs))
1649 : end if
1650 : !do jj=1,3
1651 : ! xval = dot_product(dwdq(:,jj), matmul(crystal%gprimd, qred_vers))
1652 : ! if (xval < zero) num_negw = num_negw + 1
1653 : !end do
1654 : end do
1655 13 : call xmpi_sum(num_negw, comm, ierr)
1656 13 : xval = min_negw; call xmpi_min(xval, min_negw, comm, ierr)
1657 :
1658 23725 : adiff = sum(abs(cut_phfrq - ref_phfrq)) / (ifc%nqibz * 3 * natom)
1659 13 : if (my_rank == master) then
1660 : write(ab_out,"(a,i7,1x,es13.4,4x,i8,1x,es13.4,2x,es13.4)") &
1661 13 : "-",nsphere, adiff * Ha_meV, num_negw, min_negw * Ha_meV, rcut_min
1662 : end if
1663 :
1664 14 : if (num_negw == 0) then
1665 : jl = jm ! Replace lower limit
1666 : else
1667 9 : ju = jm ! Replace upper limit
1668 : end if
1669 : end do
1670 :
1671 1 : ABI_FREE(ref_phfrq)
1672 1 : ABI_FREE(cut_phfrq)
1673 1 : ABI_FREE(save_wghatm)
1674 1 : ABI_FREE(save_atmfrc)
1675 1 : call lgrid%free()
1676 :
1677 2 : end subroutine ifc_autocutoff
1678 : !!***
1679 :
1680 : !----------------------------------------------------------------------
1681 :
1682 : !!****f* m_ifc/corsifc9
1683 : !! NAME
1684 : !! corsifc9
1685 : !!
1686 : !! FUNCTION
1687 : !! Applies a cutoff on the ifc in real space
1688 : !!
1689 : !! INPUTS
1690 : !! acell(3)=length scales by which rprim is to be multiplied
1691 : !! gprim(3,3)=dimensionless primitive translations in reciprocal space
1692 : !! natom=number of atoms in unit cell
1693 : !! nrpt= Number of R points in the Big Box
1694 : !! rcan(3,natom)=canonical coordinates of atoms
1695 : !! rprim(3,3)=dimensionless primitive translations in real space
1696 : !! rpt(3,nrpt)=canonical coordinates of the points in the BigBox.
1697 : !! nsphere=number of atoms to be included in the cut-off sphere for interatomic
1698 : !! force constant; if = 0 : maximum extent allowed by the grid.
1699 : !! rifcsph=radius for cutoff of IFC
1700 : !! wghatm(natom,natom,nrpt) = Weights associated to a pair of atoms and to a R vector.
1701 : !!
1702 : !! OUTPUT
1703 : !! wghatm(natom,natom,nrpt) = Weights associated to a pair of atoms and to a R vector
1704 : !! with the required cutoff applied.
1705 : !! rcut_min=Effective cutoff. Defined by the minimum cutoff radius over the natom sites.
1706 : !!
1707 : !! SOURCE
1708 :
1709 16 : subroutine corsifc9(acell,gprim,natom,nrpt,nsphere,rifcsph,rcan,rprim,rpt,rcut_min,wghatm)
1710 :
1711 : !Arguments -------------------------------
1712 : !scalars
1713 : integer,intent(in) :: natom,nrpt,nsphere
1714 : real(dp),intent(in) :: rifcsph
1715 : real(dp),intent(out) :: rcut_min
1716 : !arrays
1717 : real(dp),intent(in) :: acell(3), gprim(3,3),rcan(3,natom), rprim(3,3),rpt(3,nrpt)
1718 : real(dp),intent(inout) :: wghatm(natom,natom,nrpt)
1719 :
1720 : !Local variables -------------------------
1721 : !scalars
1722 : integer :: ia,ib,ii,index,irpt
1723 : real(dp) :: rmax,rsigma,r0
1724 : !arrays
1725 16 : integer,allocatable :: list(:)
1726 16 : real(dp),allocatable :: dist(:,:,:),wkdist(:)
1727 : ! *********************************************************************
1728 :
1729 : ! Compute the distances between atoms
1730 : ! dist(ia,ib,irpt) contains the distance from atom ia to atom ib in unit cell irpt.
1731 80 : ABI_MALLOC(dist,(natom,natom,nrpt))
1732 16 : call dist9(acell,dist,gprim,natom,nrpt,rcan,rprim,rpt)
1733 :
1734 48 : ABI_MALLOC(list,(natom*nrpt))
1735 48 : ABI_MALLOC(wkdist,(natom*nrpt))
1736 :
1737 : ! loop on all generic atoms.
1738 16 : rcut_min = huge(one)
1739 99 : do ia=1,natom
1740 :
1741 166 : wkdist = reshape(dist(ia,:,:), [natom*nrpt])
1742 488916 : do ii=1,natom*nrpt
1743 488916 : list(ii)=ii
1744 : end do
1745 : ! This sorting algorithm is slow ...
1746 83 : call sort_dp(natom*nrpt,wkdist,list,tol14)
1747 83 : rmax = wkdist(natom*nrpt)
1748 :
1749 : ! zero the outside IFCs: act on wghatm
1750 :
1751 : ! fix number of spheres
1752 83 : if(nsphere/=0.and.nsphere<natom*nrpt)then
1753 80 : rcut_min = min(rcut_min, wkdist(nsphere+1))
1754 345170 : do ii=nsphere+1,natom*nrpt
1755 345090 : index=list(ii)
1756 345090 : irpt=(index-1)/natom+1
1757 345090 : ib=index-natom*(irpt-1)
1758 345170 : wghatm(ia,ib,irpt)=zero
1759 : end do
1760 : end if
1761 :
1762 : ! or fix radius of maximum ifc
1763 83 : if(rifcsph>tol10)then
1764 3648 : do ii=nsphere+1,natom*nrpt
1765 3645 : index=list(ii)
1766 : ! preserve weights for atoms inside sphere of radius rifcsph
1767 3645 : if (wkdist(ii) < rifcsph) cycle
1768 3464 : rcut_min = min(rcut_min, wkdist(ii))
1769 3464 : irpt=(index-1)/natom+1
1770 3464 : ib=index-natom*(irpt-1)
1771 3648 : wghatm(ia,ib,irpt)=zero
1772 : end do
1773 : end if
1774 :
1775 : ! filter smoothly to 0 at edge of WS cell
1776 99 : if (rifcsph < -tol10) then
1777 : ! Use different filter
1778 0 : r0 = abs(rifcsph) * rmax; rsigma = half*(rmax-r0) !one
1779 0 : rcut_min = r0 ! Set it to r0
1780 0 : do ii=nsphere+1,natom*nrpt
1781 0 : index=list(ii)
1782 0 : irpt=(index-1)/natom+1
1783 0 : ib=index-natom*(irpt-1)
1784 0 : wghatm(ia,ib,irpt) = wghatm(ia,ib,irpt) * half * abi_derfc((wkdist(ii) - r0) / rsigma)
1785 : end do
1786 : end if
1787 :
1788 : end do
1789 :
1790 16 : ABI_FREE(dist)
1791 16 : ABI_FREE(list)
1792 16 : ABI_FREE(wkdist)
1793 :
1794 16 : end subroutine corsifc9
1795 : !!***
1796 :
1797 : !----------------------------------------------------------------------
1798 :
1799 : !!****f* m_ifc/ifc_write
1800 : !! NAME
1801 : !! ifc_write
1802 : !!
1803 : !! FUNCTION
1804 : !! Adds the real-space interatomic force constants to:
1805 : !! the output file,
1806 : !! a NetCDF file which is already open on ncid
1807 : !! if prt_ifc==1, to the ifcinfo.out file
1808 : !! to a TDEP file named outfile.forceconstants_ABINIT
1809 : !!
1810 : !! INPUTS
1811 : !! ifcana= 0 => no analysis of ifc ; 1 => full analysis
1812 : !! atifcflg(natom) = atifcflg(ia) equals 1 if the analysis of ifc has to be done for atom ia; otherwise 0.
1813 : !! ifcout= Number of interatomic force constants written in the output file
1814 : !! prt_ifc = flag to print out ifc information for dynamical matrix (AI2PS)
1815 : !! ncid=the id of the open NetCDF file. Set to nctk_noid if netcdf output is not wanted.
1816 : !!
1817 : !! OUTPUT
1818 : !! written in the output file and in the NetCDF file
1819 : !!
1820 : !! NOTES
1821 : !! This routine should be executed by one processor only
1822 : !!
1823 : !! TODO:
1824 : !! 1) ifc_write should not have side effects
1825 : !!
1826 : !! 2) the code is unreadable and horrible - 3/4 different file formats for the
1827 : !! same stuff. We should make different subroutines, even if it duplicates some code
1828 : !!
1829 : !! 3) The name of the output files should respect the conventions of the code.
1830 : !! This routine should take a file prefix as argument.
1831 : !!
1832 : !! SOURCE
1833 :
1834 65 : subroutine ifc_write(Ifc,ifcana,atifcflg,ifcout,prt_ifc,ncid,prefix,&
1835 : unit_out) ! optional arguments
1836 :
1837 : !Arguments -------------------------------
1838 : !scalars
1839 : class(ifc_type),intent(inout) :: Ifc
1840 : integer,intent(in) :: ifcout,ifcana,prt_ifc,ncid
1841 : integer,optional,intent(in) :: unit_out
1842 : character(*),intent(in) :: prefix
1843 : !arrays
1844 : integer,intent(in) :: atifcflg(Ifc%natom)
1845 :
1846 : !Local variables -------------------------
1847 : !scalars
1848 : integer :: ia,ib,ii,ncerr,iatifc,ifcout1,mu,nu,iout, irpt, unit_ifc, unit_tdep
1849 : real(dp) :: detdlt, maxdist_tdep
1850 : character(len=fnlen) :: filename
1851 : character(len=500) :: msg
1852 : character(len=4) :: str1, str2
1853 : !arrays
1854 65 : integer,allocatable :: list(:),indngb(:)
1855 : real(dp) :: invdlt(3,3),ra(3),xred(3),dielt(3,3), gprimd(3,3),rprimd(3,3)
1856 65 : real(dp),allocatable :: dist(:,:,:),wkdist(:),rsiaf(:,:,:),sriaf(:,:,:),vect(:,:,:), posngb(:,:),wghia(:)
1857 : ! *********************************************************************
1858 :
1859 65 : iout = ab_out; if (present(unit_out)) iout = unit_out
1860 845 : dielt = ifc%dielt
1861 :
1862 : ! Compute the distances between atoms
1863 325 : ABI_MALLOC(dist,(Ifc%natom,Ifc%natom,Ifc%nrpt))
1864 65 : call dist9(Ifc%acell,dist,Ifc%gprim,Ifc%natom,Ifc%nrpt,Ifc%rcan,Ifc%rprim,Ifc%rpt)
1865 : ! Now dist(ia,ib,irpt) contains the distance from atom ia to atom ib in unit cell irpt.
1866 :
1867 195 : ABI_MALLOC(list,(Ifc%natom*Ifc%nrpt))
1868 195 : ABI_MALLOC(wkdist,(Ifc%natom*Ifc%nrpt))
1869 :
1870 : ! Calculating the inverse (transpose) of the dielectric tensor
1871 65 : call matr3inv(dielt,invdlt)
1872 :
1873 : ! Calculating the determinant of the dielectric tensor
1874 : detdlt=dielt(1,1)*dielt(2,2)*dielt(3,3)+dielt(1,3)*dielt(2,1)*&
1875 : dielt(3,2)+dielt(1,2)*dielt(2,3)*dielt(3,1)-dielt(1,3)*&
1876 : dielt(2,2)*dielt(3,1)-dielt(1,1)*dielt(2,3)*dielt(3,2)-&
1877 65 : dielt(1,2)*dielt(2,1)*dielt(3,3)
1878 :
1879 : ! echo to log file
1880 65 : write(std_out,'(a)' )' ifc_write: analysis of interatomic force constants '
1881 65 : call mkrdim(Ifc%acell,Ifc%rprim,rprimd)
1882 65 : call matr3inv(rprimd,gprimd)
1883 :
1884 65 : if (iout > 0) then
1885 65 : write(iout, '(/,a,/)' )' Analysis of interatomic force constants '
1886 65 : if(Ifc%dipdip==1.and.Ifc%dipquad==0.and.Ifc%quadquad==0)then
1887 17 : write(iout, '(a)' )' Are given : column(1-3), the total force constant'
1888 17 : write(iout, '(a)' )' then column(4-6), the Ewald part'
1889 17 : write(iout, '(a)' )' then column(7-9), the short-range part'
1890 17 : write(iout, '(a)' )' Column 1, 4 and 7 are related to the displacement'
1891 17 : write(iout, '(a)' )' of the generic atom along x, '
1892 17 : write(iout, '(a)' )' column 2, 5 and 8 are related to the displacement'
1893 17 : write(iout, '(a)' )' of the generic atom along y, '
1894 17 : write(iout, '(a)' )' column 3, 6 and 9 are related to the displacement'
1895 17 : write(iout, '(a)')' of the generic atom along z. '
1896 48 : else if(Ifc%dipquad==1.or.Ifc%quadquad==1)then
1897 0 : write(iout, '(a)' )' Are given : column(1-3), ONLY the short-range part!!!!'
1898 0 : write(iout, '(a)' )' column 1 is related to the displacement'
1899 0 : write(iout, '(a)' )' of the generic atom along x, '
1900 0 : write(iout, '(a)' )' column 2 is related to the displacement'
1901 0 : write(iout, '(a)' )' of the generic atom along y, '
1902 0 : write(iout, '(a)' )' column 3 is related to the displacement'
1903 0 : write(iout, '(a)' )' of the generic atom along z, '
1904 48 : else if(Ifc%dipdip==0)then
1905 48 : write(iout, '(a)' )' column 1 is related to the displacement'
1906 48 : write(iout, '(a)' )' of the generic atom along x, '
1907 48 : write(iout, '(a)' )' column 2 is related to the displacement'
1908 48 : write(iout, '(a)' )' of the generic atom along y, '
1909 48 : write(iout, '(a)' )' column 3 is related to the displacement'
1910 48 : write(iout, '(a)' )' of the generic atom along z, '
1911 : end if
1912 : end if
1913 :
1914 65 : if (ifcout>Ifc%natom*Ifc%nrpt .or. ifcout == -1) then
1915 0 : ifcout1=Ifc%natom*Ifc%nrpt
1916 : write(msg, '(3a,i0,a)' )&
1917 0 : 'The value of ifcout exceeds the number of atoms in the big box.', ch10, &
1918 0 : 'Output limited to ',Ifc%natom*Ifc%nrpt,' atoms.'
1919 0 : ABI_WARNING(msg)
1920 : else
1921 65 : ifcout1=ifcout
1922 : end if
1923 :
1924 : ! set up file for real space ifc output, if required
1925 65 : if (prt_ifc == 1) then
1926 53 : filename = trim(prefix)//'_ifcinfo.dat'
1927 53 : if (open_file(trim(filename), msg, newunit=unit_ifc, status="replace") /= 0) then
1928 0 : ABI_ERROR(msg)
1929 : end if
1930 53 : write(iout, '(a,a)' )ch10,&
1931 106 : ' NOTE: Open file _ifcinfo.dat, for the output of interatomic force constants. This is because prt_ifc==1. '
1932 :
1933 53 : filename = trim(prefix)//'_forceconstants.dat'
1934 53 : if (open_file(trim(filename), msg, newunit=unit_tdep, status="replace") /= 0) then
1935 0 : ABI_ERROR(msg)
1936 : end if
1937 53 : write(iout, '(a,a,a)' )ch10,&
1938 53 : ' NOTE: Open file _forceconstants.dat, for the output of interatomic force',&
1939 106 : ' constants in TDEP format. This is because prt_ifc==1. '
1940 : ! Print necessary stuff for TDEP
1941 53 : write(unit_tdep,"(1X,I10,15X,'How many atoms per unit cell')") Ifc%natom
1942 :
1943 : ! look at all pairs, find furthest one with weight 1
1944 : ! do ia
1945 : !Ifc%wghatm(ia,ib,irpt)
1946 53 : maxdist_tdep = Ifc%r_inscribed_sphere !maxval(dist)*0.8_dp
1947 53 : write(unit_tdep,"(1X,F20.15,5X,'Realspace cutoff (A)')") maxdist_tdep*Bohr_Ang
1948 : end if
1949 :
1950 65 : if (ncid /= nctk_noid) then
1951 : ! initialize netcdf variables
1952 : ncerr = nctk_def_dims(ncid, [nctkdim_t("natifc", SUM(atifcflg)), nctkdim_t("number_of_r_points_big_box", Ifc%nrpt), &
1953 499 : nctkdim_t("number_of_atoms_big_box", Ifc%natom*Ifc%nrpt), nctkdim_t("ifcout", ifcout1)], defmode=.True.)
1954 65 : NCF_CHECK(ncerr)
1955 :
1956 : ncerr = nctk_def_arrays(ncid, [&
1957 : nctkarr_t('ifc_atoms_indices', "i", "natifc"),&
1958 : nctkarr_t('ifc_neighbours_indices', "i", "ifcout, natifc"),&
1959 : nctkarr_t('ifc_distances', "dp", "ifcout, natifc "),&
1960 : nctkarr_t('ifc_matrix_cart_coord', "dp", "number_of_cartesian_directions,number_of_cartesian_directions, ifcout, natifc"),&
1961 : nctkarr_t('ifc_atoms_cart_coord', "dp", "number_of_cartesian_directions,ifcout, natifc"),&
1962 455 : nctkarr_t('ifc_weights', "dp", "ifcout, natifc")])
1963 65 : NCF_CHECK(ncerr)
1964 :
1965 65 : if (Ifc%dipdip==1) then
1966 : ncerr = nctk_def_arrays(ncid, [&
1967 : nctkarr_t('ifc_matrix_cart_coord_short_range', "dp", &
1968 34 : "number_of_cartesian_directions, number_of_cartesian_directions, ifcout, natifc")])
1969 17 : NCF_CHECK(ncerr)
1970 : end if
1971 :
1972 65 : if (ifcana==1) then
1973 : ncerr = nctk_def_arrays(ncid, [&
1974 128 : nctkarr_t('ifc_local_vectors', "dp", "number_of_cartesian_directions, number_of_cartesian_directions, ifcout, natifc")])
1975 64 : NCF_CHECK(ncerr)
1976 : end if
1977 :
1978 65 : NCF_CHECK(nctk_set_datamode(ncid))
1979 : end if
1980 :
1981 195 : ABI_MALLOC(rsiaf,(3,3,ifcout1))
1982 130 : ABI_MALLOC(sriaf,(3,3,ifcout1))
1983 130 : ABI_MALLOC(vect,(3,3,ifcout1))
1984 195 : ABI_MALLOC(indngb,(ifcout1))
1985 195 : ABI_MALLOC(posngb,(3,ifcout1))
1986 195 : ABI_MALLOC(wghia,(ifcout1))
1987 :
1988 65 : iatifc=0
1989 :
1990 : ! BIG loop on all generic atoms
1991 239 : do ia=1,Ifc%natom
1992 239 : if(atifcflg(ia)==1)then
1993 :
1994 141 : iatifc=iatifc+1
1995 :
1996 : ! First transform canonical coordinates to reduced coordinates
1997 564 : do ii=1,3
1998 564 : xred(ii)=Ifc%gprim(1,ii)*Ifc%rcan(1,ia)+Ifc%gprim(2,ii)*Ifc%rcan(2,ia)+Ifc%gprim(3,ii)*Ifc%rcan(3,ia)
1999 : end do
2000 : ! Then to cartesian coordinates
2001 564 : ra(:)=xred(1)*Ifc%acell(1)*Ifc%rprim(:,1)+ xred(2)*Ifc%acell(2)*Ifc%rprim(:,2)+ xred(3)*Ifc%acell(3)*Ifc%rprim(:,3)
2002 :
2003 : ! This sorting algorithm is slow ...
2004 282 : wkdist(:)=reshape(dist(ia,:,:),(/Ifc%natom*Ifc%nrpt/))
2005 383886 : do ii=1,Ifc%natom*Ifc%nrpt
2006 383886 : list(ii)=ii
2007 : end do
2008 141 : call sort_dp(Ifc%natom*Ifc%nrpt,wkdist,list,tol14)
2009 :
2010 141 : if (iout > 0) then
2011 141 : write(iout, '(a)' )
2012 141 : write(std_out,'(a,i4)' )' generic atom number',ia
2013 141 : write(iout, '(a,i4)' )' generic atom number',ia
2014 141 : write(std_out,'(a,3es16.8)' ) ' with cartesian coordinates',ra(1:3)
2015 141 : write(iout,'(a,3es16.8)' ) ' with cartesian coordinates',ra(1:3)
2016 141 : write(iout, '(a)' )
2017 : end if
2018 :
2019 : ! =========================================================================== !
2020 : ! Side effects happen here
2021 : ! The long-range part of the IFC (ewald_atmfrc) gets computed.
2022 :
2023 : ! BEGIN DEBUG
2024 : !write(*,*) 'ifc_write (1) : Writing atmfrc'
2025 : !write(*,*) 'atmfrc : ', Ifc%atmfrc
2026 : !call flush()
2027 : ! END DEBUG
2028 : call ifc_getiaf(Ifc,ifcana,ifcout1,iout,ifc%zeff,ia,ra,list,dist,invdlt,&
2029 141 : detdlt,rsiaf,sriaf,vect,indngb,posngb)
2030 :
2031 : ! =========================================================================== !
2032 :
2033 141 : if (prt_ifc == 1) then
2034 24220 : do ii=1,ifcout1
2035 24220 : if (wkdist(ii) > maxdist_tdep) exit
2036 : end do
2037 123 : ii = ii - 1
2038 123 : write(unit_tdep,"(1X,I10,15X,'How many neighbours does atom ',I3,' have')") ii, ia
2039 :
2040 24220 : do ii=1,ifcout1
2041 24097 : ib = indngb(ii)
2042 24097 : irpt = (list(ii)-1)/Ifc%natom+1
2043 24097 : wghia(ii) = Ifc%wghatm(ia,ib,irpt)
2044 : ! limit printing to maximum distance for tdep
2045 24097 : if (wkdist(ii) > maxdist_tdep) cycle
2046 :
2047 : !TDEP
2048 24097 : call int2char4(ii, str1)
2049 24097 : call int2char4(ia, str2)
2050 24097 : write(unit_tdep,"(1X,I10,15X,a,a,a,a)") ib, &
2051 24097 : & 'In the unit cell, what is the index of neighbour ', &
2052 48194 : & trim(str1), " of atom ", trim(str2)
2053 : ! The lattice vector needs to be in reduced coordinates?
2054 : ! TODO: check whether this is correct for TDEP: might need just lattice
2055 : ! vector part and not full vector, and could be in integers instead of
2056 : ! cartesian vector...
2057 313261 : write (unit_tdep,'(3es28.16)') matmul(Ifc%rpt(1:3,irpt),Ifc%gprim)
2058 :
2059 : !AI2PS
2060 24097 : write(unit_ifc,'(i6,i6)') ia,ii
2061 24097 : write(unit_ifc,'(3es28.16)') posngb(1:3,ii)
2062 96511 : do nu=1,3
2063 : !TDEp
2064 : ! And the actual short ranged forceconstant: TODO: check if
2065 : ! a transpose is needed or a swap between the nu and the mu
2066 : !write(unit_tdep,'(3f28.16)') (sriaf(nu,mu,ii)*Ha_eV/amu_emass, mu=1, 3)
2067 289164 : write(unit_tdep,'(3f28.16)') (Ifc%short_atmfrc(mu,ia,nu,ib,irpt)*Ha_eV/Bohr_Ang**2, mu=1, 3)
2068 :
2069 : !AI2PS
2070 96388 : write(unit_ifc,'(3f28.16)')(rsiaf(nu,mu,ii),mu=1,3)
2071 : end do
2072 : end do
2073 :
2074 123 : if (ncid /= nctk_noid) then
2075 246 : NCF_CHECK(nf90_put_var(ncid, vid("ifc_atoms_indices"), ia, start=[iatifc]))
2076 615 : NCF_CHECK(nf90_put_var(ncid, vid("ifc_neighbours_indices"), indngb, start=[1,iatifc], count=[ifcout1,1]))
2077 615 : NCF_CHECK(nf90_put_var(ncid, vid("ifc_distances"), wkdist(:ifcout1), start=[1,iatifc],count=[ifcout1,1]))
2078 1107 : ncerr = nf90_put_var(ncid, vid("ifc_matrix_cart_coord"), rsiaf, start=[1,1,1,iatifc], count=[3,3,ifcout1,1])
2079 123 : NCF_CHECK(ncerr)
2080 861 : NCF_CHECK(nf90_put_var(ncid, vid("ifc_atoms_cart_coord"), posngb, start=[1,1,iatifc], count=[3,ifcout1,1]))
2081 615 : NCF_CHECK(nf90_put_var(ncid, vid("ifc_weights"), wghia, start=[1,iatifc], count=[ifcout1,1]))
2082 :
2083 123 : if (Ifc%dipdip==1) then
2084 : ncerr = nf90_put_var(ncid, vid("ifc_matrix_cart_coord_short_range"), sriaf, &
2085 216 : start=[1,1,1,iatifc], count=[3,3,ifcout1,1])
2086 24 : NCF_CHECK(ncerr)
2087 : end if
2088 123 : if (ifcana==1) then
2089 1107 : ncerr = nf90_put_var(ncid, vid("ifc_local_vectors"), vect, start=[1,1,1,iatifc], count=[3,3,ifcout1,1])
2090 123 : NCF_CHECK(ncerr)
2091 : end if
2092 : end if
2093 : end if
2094 : end if ! End the condition on atifcflg
2095 : end do ! End Big loop on atoms in the unit cell, and corresponding test
2096 :
2097 :
2098 : ! NB for future use: in TDEP the following can also be provided.
2099 : ! ! Print some auxiliary information, if it is there. Such as norm of
2100 : ! ! forceconstant per shell, which shells there are and so on.
2101 : ! if ( fc%npairshells .gt. 0 .and. allocated(fc%pairshell) ) then
2102 : ! write(u,"(1X,I10,15X,'Number of irreducible coordination shells')") fc%npairshells
2103 : ! do i=1,fc%npairshells
2104 : ! write(u,"(1X,I10,1X,F16.10,1X,F16.10,15X,'number atoms in shell, radius, norm of forceconstant',I0)") fc%pairshell(i)%n,fc%pairshell(i)%rad,fc%pairshell(i)%norm,i
2105 : ! enddo
2106 : ! do i=1,fc%npairshells
2107 : ! do j=1,fc%pairshell(i)%n
2108 : ! write(u,"(1X,3(1X,F18.12),2(1X,I0))") lo_chop(matmul(p%inv_latticevectors,fc%pairshell(i)%vec(:,j)),lo_sqtol),fc%pairshell(i)%atind(j),fc%pairshell(i)%pairind(j)
2109 : ! enddo
2110 : ! enddo
2111 : ! endif
2112 :
2113 65 : if (prt_ifc == 1) then
2114 53 : close(unit_ifc)
2115 53 : close(unit_tdep)
2116 :
2117 53 : filename = trim(prefix)//'_lotosplitting.dat'
2118 53 : if (open_file(trim(filename), msg, newunit=unit_tdep, status="replace") /= 0) then
2119 0 : ABI_ERROR(msg)
2120 : end if
2121 53 : write(unit_tdep,'(3es28.16)') dielt(:,1)
2122 53 : write(unit_tdep,'(3es28.16)') dielt(:,2)
2123 53 : write(unit_tdep,'(3es28.16)') dielt(:,3)
2124 183 : do ia = 1, Ifc%natom
2125 573 : do ii = 1, 3
2126 520 : write(unit_tdep,'(3es28.16)') ifc%zeff(:,ii,ia)
2127 : end do
2128 : end do
2129 53 : close(unit_tdep)
2130 : end if
2131 :
2132 65 : ABI_FREE(rsiaf)
2133 65 : ABI_FREE(sriaf)
2134 65 : ABI_FREE(vect)
2135 65 : ABI_FREE(indngb)
2136 65 : ABI_FREE(posngb)
2137 65 : ABI_FREE(dist)
2138 65 : ABI_FREE(list)
2139 65 : ABI_FREE(wkdist)
2140 65 : ABI_FREE(wghia)
2141 :
2142 : contains
2143 885 : integer function vid(vname)
2144 : character(len=*),intent(in) :: vname
2145 885 : vid = nctk_idname(ncid, vname)
2146 : end function vid
2147 :
2148 : end subroutine ifc_write
2149 : !!***
2150 :
2151 : !----------------------------------------------------------------------
2152 :
2153 : !!****f* m_ifc/ifc_getiaf
2154 : !!
2155 : !! NAME
2156 : !! ifc_getiaf
2157 : !!
2158 : !! FUNCTION
2159 : !! Extracts the IFCs needed for the output for one atom in the
2160 : !! unit cell. Accumulates the results for writing in the NetCDF file.
2161 : !! Prints to the output file
2162 : !!
2163 : !! INPUTS
2164 : !! ifcana= 0 => no analysis of ifc ; 1 => full analysis
2165 : !! ifcout= Number of interatomic force constants written in the output file
2166 : !! iout=unit number for nice output
2167 : !! zeff(3,3,natom)=effective charge on each atom, versus electric field and atomic displacement
2168 : !! ia=index of the atom in the unit cell for which the IFCs are being analyzed
2169 : !! ra(3)=position of atom ia in cartesian coordinates
2170 : !! list(ifcout)=index of permutation for distances from atom ia in ascending order
2171 : !! dist(natom,natom,nrpt)=distance from atom ia to atom ib in unit cell irpt.
2172 : !! invdlt(3,3)=inverse (transpose) of the dielectric tensor
2173 : !! detdlt=determinant of the dielectric tensor
2174 : !!
2175 : !! OUTPUT
2176 : !! rsiaf(3,3,ifcout)=list of real space IFCs
2177 : !! sriaf(3,3,ifcout)=list of the short range part of the real space IFCs
2178 : !! vect(3,3,ifcout)=base vectors for local coordinates (longitudinal/transverse), if ifc_getiaf is able to find
2179 : !! a third atom not aligned with the two atoms characterizing the IFC. If no, the second and third vectors are set to zero.
2180 : !! indngb(ifcout)=indices in the unit cell of the neighbouring atoms
2181 : !! posngb(3,ifcout)=position of the neighbouring atoms in cartesian coordinates
2182 : !! output file
2183 : !!
2184 : !! SIDE EFFECTS
2185 : !! The long-range part of the IFC (ewald_atmfrc) gets computed.
2186 : !!
2187 : !! NOTES
2188 : !! This routine should be executed by one processor only
2189 : !!
2190 : !! SOURCE
2191 :
2192 141 : subroutine ifc_getiaf(Ifc,ifcana,ifcout,iout,zeff,ia,ra,list,&
2193 141 : dist,invdlt,detdlt,rsiaf,sriaf,vect,indngb,posngb)
2194 :
2195 : !Arguments -------------------------------
2196 : !scalars
2197 : class(ifc_type),intent(inout) :: Ifc
2198 : integer,intent(in) :: ia,ifcana,ifcout,iout
2199 : real(dp), intent(in) :: detdlt
2200 : !arrays
2201 : real(dp),intent(in) :: invdlt(3,3),ra(3)
2202 : real(dp),intent(in) :: dist(Ifc%natom,Ifc%natom,Ifc%nrpt)
2203 : real(dp),intent(in) :: zeff(3,3,Ifc%natom)
2204 : integer,intent(in) :: list(Ifc%natom*Ifc%nrpt)
2205 : integer,intent(out) :: indngb(ifcout)
2206 : real(dp),intent(out) :: rsiaf(3,3,ifcout),sriaf(3,3,ifcout),vect(3,3,ifcout),posngb(3,ifcout)
2207 :
2208 : !Local variables -------------------------
2209 : !scalars
2210 : integer :: flag,ib,ii,index,jj,kk,mu,nu,irpt
2211 : real(dp) :: ew1,rsq,scprod,trace1,trace2,trace3, yy,dist1
2212 : character(len=500) :: msg
2213 : !arrays
2214 : real(dp) :: ewiaf0(3,3),ewiaf1(3,3),ewloc(3,3),ifcloc(3,3)
2215 : real(dp) :: rcart(3),rdiff(3),rsloc(3,3), srloc(3,3),vect1(3),vect2(3),vect3(3),work(3),xx(3)
2216 : ! *********************************************************************
2217 :
2218 141 : if(ifcana==1)then
2219 : ! Generate the local coordinate system for the atom ia
2220 138 : index=list(2)
2221 138 : write(std_out,*)index
2222 138 : call canct9(Ifc%acell,Ifc%gprim,ib,index,irpt,Ifc%natom,Ifc%nrpt,Ifc%rcan,rcart,Ifc%rprim,Ifc%rpt)
2223 138 : vect2(1)=rcart(1)-ra(1)
2224 138 : vect2(2)=rcart(2)-ra(2)
2225 138 : vect2(3)=rcart(3)-ra(3)
2226 138 : flag=0
2227 160 : do ii=3,Ifc%natom*Ifc%nrpt
2228 160 : index=list(ii)
2229 160 : call canct9(Ifc%acell,Ifc%gprim,ib,index,irpt,Ifc%natom,Ifc%nrpt,Ifc%rcan,rcart,Ifc%rprim,Ifc%rpt)
2230 160 : vect1(1)=(rcart(1)-ra(1))-vect2(1)
2231 160 : vect1(2)=(rcart(2)-ra(2))-vect2(2)
2232 160 : vect1(3)=(rcart(3)-ra(3))-vect2(3)
2233 160 : scprod=0.0_dp
2234 640 : do jj=1,3
2235 640 : scprod=scprod+vect1(jj)**2
2236 : end do
2237 640 : do jj=1,3
2238 640 : vect1(jj)=vect1(jj)/scprod**0.5
2239 : end do
2240 : scprod=0.0_dp
2241 640 : do jj=1,3
2242 640 : scprod=scprod+vect2(jj)*vect1(jj)
2243 : end do
2244 640 : do jj=1,3
2245 640 : work(jj)=vect2(jj)-vect1(jj)*scprod
2246 : end do
2247 : scprod=0.0_dp
2248 640 : do jj=1,3
2249 640 : scprod=scprod+work(jj)**2
2250 : end do
2251 160 : if(scprod>1.0d-10)then
2252 : flag=1
2253 : end if
2254 0 : if(flag==1)exit
2255 : end do
2256 138 : if(flag==0)then
2257 : write(msg, '(3a)' )&
2258 0 : 'Unable to find a third atom not aligned with the two selected ones.',ch10,&
2259 0 : 'The local analysis (longitudinal/transverse) will not be done. The two transverse vectors are set to zero.'
2260 0 : ABI_WARNING(msg)
2261 0 : vect2(:)=zero ; vect3(:)=zero
2262 : else
2263 138 : vect2(1)=work(1)/scprod**0.5
2264 138 : vect2(2)=work(2)/scprod**0.5
2265 138 : vect2(3)=work(3)/scprod**0.5
2266 138 : vect3(1)=vect1(2)*vect2(3)-vect1(3)*vect2(2)
2267 138 : vect3(2)=vect1(3)*vect2(1)-vect1(1)*vect2(3)
2268 138 : vect3(3)=vect1(1)*vect2(2)-vect1(2)*vect2(1)
2269 138 : if (iout > 0) then
2270 138 : write(iout, '(a)' )' Third atom defining local coordinates : '
2271 138 : write(iout, '(a,i4,a,i4)' )' ib = ',ib,' irpt = ',irpt
2272 : end if
2273 : endif
2274 : end if
2275 :
2276 : ! Analysis and output of force constants, ordered with respect to the distance from atom ia
2277 24596 : do ii=1,ifcout
2278 24455 : index=list(ii)
2279 24455 : call canct9(Ifc%acell,Ifc%gprim,ib,index,irpt,Ifc%natom,Ifc%nrpt,Ifc%rcan,posngb(:,ii),Ifc%rprim,Ifc%rpt)
2280 24455 : indngb(ii)=ib
2281 24455 : dist1=dist(ia,ib,irpt)
2282 24455 : if (iout > 0) then
2283 24455 : write(iout, '(a)' )
2284 24455 : write(iout, '(i4,a,i6,a,i8)' )ii,' interaction with atom',ib,' cell',irpt
2285 97820 : write(iout, '(a,3es16.6)' )' with coordinates ',posngb(1:3,ii)*(one+tol8)
2286 24455 : write(iout, '(a,es16.6)' )' and distance ',dist1
2287 : end if
2288 :
2289 24455 : if(ifcana==1.and.ii/=1)then
2290 24257 : vect1(1)=(posngb(1,ii)-ra(1))/dist1
2291 24257 : vect1(2)=(posngb(2,ii)-ra(2))/dist1
2292 24257 : vect1(3)=(posngb(3,ii)-ra(3))/dist1
2293 : end if
2294 :
2295 24596 : if(Ifc%dipdip==0.or.Ifc%dipquad==1.or.Ifc%quadquad==1)then
2296 : ! Get the "total" force constants (=real space FC)
2297 : ! without taking into account the dipole-dipole interaction
2298 79140 : do mu=1,3
2299 257205 : do nu=1,3
2300 237420 : rsiaf(mu,nu,ii)=Ifc%atmfrc(mu,ia,nu,ib,irpt) * Ifc%wghatm(ia,ib,irpt)
2301 : end do
2302 : end do
2303 : ! Output of the ifcs in cartesian coordinates
2304 19785 : if (iout > 0) then
2305 79140 : do nu=1,3
2306 237420 : write(iout, '(1x,3f9.5)' )(rsiaf(mu,nu,ii)+tol10,mu=1,3)
2307 : ! transfer short range and long range
2308 257205 : do mu=1,3
2309 237420 : Ifc%short_atmfrc(mu,ia,nu,ib,irpt) = rsiaf(mu,nu,ii) + tol10
2310 : end do
2311 :
2312 : end do
2313 : end if
2314 :
2315 19785 : if(ifcana==1)then
2316 : ! Further analysis
2317 19785 : trace1=rsiaf(1,1,ii)+rsiaf(2,2,ii)+rsiaf(3,3,ii)
2318 19785 : if (iout > 0) then
2319 19785 : write(iout, '(a,f9.5)' ) ' Trace ',trace1+tol10
2320 : end if
2321 19785 : if(flag==1)then
2322 19785 : if(ii/=1)then
2323 19680 : call axial9(rsiaf(:,:,ii),vect1,vect2,vect3)
2324 : end if
2325 19785 : if (iout > 0) then
2326 19785 : write(iout, '(a)' )' Transformation to local coordinates '
2327 19785 : write(iout, '(a,3f16.6)' ) ' First local vector :',vect1
2328 19785 : write(iout, '(a,3f16.6)' ) ' Second local vector :',vect2
2329 19785 : write(iout, '(a,3f16.6)' ) ' Third local vector :',vect3
2330 : end if
2331 19785 : call ifclo9(rsiaf(:,:,ii),ifcloc,vect1,vect2,vect3)
2332 19785 : if (iout > 0) then
2333 79140 : do nu=1,3
2334 257205 : write(iout, '(1x,3f9.5)' )(ifcloc(mu,nu)+tol10,mu=1,3)
2335 : end do
2336 : end if
2337 : endif ! flag==1
2338 :
2339 79140 : vect(:,1,ii) = vect1
2340 79140 : vect(:,2,ii) = vect2
2341 79140 : vect(:,3,ii) = vect3
2342 :
2343 : end if ! Further analysis finished
2344 :
2345 4670 : else if(Ifc%dipdip==1)then
2346 :
2347 : !write(iout,'(a)')
2348 : !write(iout,'(a)')' Enter dipdip section, for debugging'
2349 : !write(iout,'(a)')
2350 :
2351 : ! Get the Coulomb part
2352 18680 : do jj=1,3
2353 18680 : rdiff(jj)=ra(jj)-posngb(jj,ii)
2354 : end do
2355 4670 : rsq=0.0_dp
2356 4670 : xx(1:3)=0.0_dp
2357 18680 : do jj=1,3
2358 60710 : do kk=1,3
2359 42030 : ewiaf0(jj,kk)=0.0_dp
2360 42030 : rsq=rsq+rdiff(jj)*invdlt(kk,jj)*rdiff(kk)
2361 56040 : xx(kk)=xx(kk)+invdlt(kk,jj)*rdiff(jj)
2362 : end do
2363 : end do
2364 4670 : yy=sqrt(rsq)
2365 : ! Avoid zero denominators in term:
2366 4670 : if (sqrt(rsq)>=tol12) then
2367 18536 : do mu=1,3
2368 60242 : do nu=1,3
2369 55608 : ewiaf0(mu,nu)=(-3*xx(nu)*xx(mu)+invdlt(nu,mu)*yy**2)/yy**5/sqrt(detdlt)
2370 : end do
2371 : end do
2372 : else
2373 36 : if (ia/=ib)then
2374 : write(msg, '(a,a,a,a,a,i5,a,i5,a)' )&
2375 0 : 'The distance between two atoms vanishes.',ch10,&
2376 0 : 'This is not allowed.',ch10,&
2377 0 : 'Action: check the input for the atoms number',ia,' and',ib,'.'
2378 0 : ABI_ERROR(msg)
2379 : end if
2380 : end if
2381 :
2382 : ! Take into account the effective charge tensor
2383 18680 : do mu=1,3
2384 60710 : do nu=1,3
2385 42030 : ew1=zero
2386 42030 : if(ii==1)then
2387 324 : ew1=-Ifc%dyewq0(mu,nu,ia)
2388 : end if
2389 168120 : do jj=1,3
2390 546390 : do kk=1,3
2391 504360 : ew1=ew1+zeff(jj,mu,ia)*(zeff(kk,nu,ib)* ewiaf0(jj,kk))
2392 : end do
2393 : end do
2394 56040 : ewiaf1(mu,nu)=ew1
2395 : end do
2396 : end do
2397 : ! Get the short-range force constants and the
2398 : ! "total" force constants (=real space FC)
2399 18680 : do mu=1,3
2400 60710 : do nu=1,3
2401 42030 : sriaf(mu,nu,ii)=Ifc%atmfrc(mu,ia,nu,ib,irpt)* Ifc%wghatm(ia,ib,irpt)
2402 56040 : rsiaf(mu,nu,ii)=ewiaf1(mu,nu)+sriaf(mu,nu,ii)
2403 : end do
2404 : end do
2405 :
2406 : ! Output of the results
2407 4670 : if (iout > 0) then
2408 18680 : do nu=1,3
2409 : write(iout, '(1x,3(3f9.5,1x))' )&
2410 56040 : (rsiaf(mu,nu,ii) +tol10,mu=1,3),&
2411 56040 : (ewiaf1(mu,nu)+tol10,mu=1,3),&
2412 168120 : (sriaf(mu,nu,ii) +tol10,mu=1,3)
2413 :
2414 : ! transfer short range and long range
2415 60710 : do mu=1,3
2416 42030 : Ifc%short_atmfrc(mu,ia,nu,ib,irpt) = sriaf(mu,nu,ii) + tol10
2417 56040 : Ifc%ewald_atmfrc(mu,ia,nu,ib,irpt) = ewiaf1(mu,nu) + tol10
2418 : end do
2419 : end do
2420 : end if
2421 :
2422 4670 : if(ifcana==1)then
2423 : ! Further analysis
2424 4610 : if (iout > 0) then
2425 4610 : write(iout, '(a)' )' Traces (and ratios) :'
2426 : end if
2427 4610 : trace1=rsiaf(1,1,ii)+rsiaf(2,2,ii)+rsiaf(3,3,ii)
2428 4610 : trace2=ewiaf1(1,1)+ewiaf1(2,2)+ewiaf1(3,3)
2429 4610 : trace3=sriaf(1,1,ii)+sriaf(2,2,ii)+sriaf(3,3,ii)
2430 4610 : if (iout > 0) then
2431 4610 : write(iout,'(3(f9.5,17x))')trace1+tol10,trace2+tol10,trace3+tol10
2432 4610 : write(iout,'(3(f9.5,17x))')1.0,(trace2+tol10)/(trace1+tol10),(trace3+tol10)/(trace1+tol10) !
2433 : end if
2434 :
2435 4610 : if(flag==1)then
2436 4610 : if(ii/=1)then
2437 4577 : call axial9(rsiaf(:,:,ii),vect1,vect2,vect3)
2438 : end if
2439 4610 : if (iout > 0) then
2440 4610 : write(iout, '(a)' )' Transformation to local coordinates '
2441 4610 : write(iout, '(a,3f16.6)' )' First local vector :',vect1
2442 4610 : write(iout, '(a,3f16.6)' )' Second local vector :',vect2
2443 4610 : write(iout, '(a,3f16.6)' )' Third local vector :',vect3
2444 : end if
2445 4610 : call ifclo9(rsiaf(:,:,ii),rsloc,vect1,vect2,vect3)
2446 4610 : call ifclo9(ewiaf1,ewloc,vect1,vect2,vect3)
2447 4610 : call ifclo9(sriaf(:,:,ii),srloc,vect1,vect2,vect3)
2448 4610 : if (iout > 0) then
2449 18440 : do nu=1,3
2450 : write(iout, '(1x,3(3f9.5,1x))' )&
2451 55320 : & (rsloc(mu,nu)+tol10,mu=1,3),&
2452 55320 : & (ewloc(mu,nu)+tol10,mu=1,3),&
2453 170570 : & (srloc(mu,nu)+tol10,mu=1,3)
2454 : end do
2455 4610 : if(ii/=1)then
2456 4577 : write(iout, '(a)' )' Ratio with respect to the longitudinal ifc'
2457 : else
2458 33 : write(iout, '(a)' )' Ratio with respect to the (1,1) element'
2459 : end if
2460 18440 : do nu=1,3
2461 : write(iout, '(1x,3(3f9.5,1x))' )&
2462 55320 : & (rsloc(mu,nu)/rsloc(1,1)+tol10,mu=1,3),&
2463 55320 : & (ewloc(mu,nu)/rsloc(1,1)+tol10,mu=1,3),&
2464 170570 : & (srloc(mu,nu)/rsloc(1,1)+tol10,mu=1,3)
2465 : end do
2466 : end if
2467 : endif ! flag==1
2468 :
2469 18440 : vect(:,1,ii) = vect1
2470 18440 : vect(:,2,ii) = vect2
2471 18440 : vect(:,3,ii) = vect3
2472 :
2473 : end if ! Further analysis finished
2474 : end if ! End the condition on dipdip
2475 : end do ! End loop over all atoms in BigBox:
2476 :
2477 141 : end subroutine ifc_getiaf
2478 : !!***
2479 :
2480 : !----------------------------------------------------------------------
2481 :
2482 : !!****f* m_ifc/omega_decomp
2483 : !!
2484 : !! NAME
2485 : !! omega_decomp
2486 : !!
2487 : !! FUNCTION
2488 : !! Compute and return the eigenvalues (frequencies) of the short-range and
2489 : !! long-range part of the dynamical matrix See Europhys. Lett. 33 p.713 (1996) for details.
2490 : !! (included by U. Aschauer and EB)
2491 : !!
2492 : !! INPUTS
2493 : !!
2494 : !! OUTPUT
2495 : !!
2496 : !! SOURCE
2497 :
2498 0 : subroutine omega_decomp(amu,natom,ntypat,typat,dynmatfl,dynmatsr,dynmatlr,iqpt,nqpt,eigenvec)
2499 :
2500 : !Arguments -------------------------------
2501 : !scalars
2502 : integer,intent(in) :: natom,ntypat
2503 : integer,intent(in) :: iqpt,nqpt
2504 : !arrays
2505 : integer,intent(in) :: typat(natom)
2506 : real(dp),intent(in) :: amu(ntypat)
2507 : real(dp),intent(inout) :: dynmatfl(2,3,natom,3,natom,nqpt)
2508 : real(dp),intent(inout) :: dynmatsr(2,3,natom,3,natom,nqpt)
2509 : real(dp),intent(inout) :: dynmatlr(2,3,natom,3,natom,nqpt)
2510 : real(dp),intent(in) :: eigenvec(2*3*natom*3*natom)
2511 :
2512 : !Local variables -------------------------
2513 : !scalars
2514 : integer :: i1,i2,idir1,idir2,imode,ipert1,ipert2,index1,index2
2515 : real(dp),parameter :: break_symm=1.0d-12
2516 : real(dp) :: fac
2517 : !arrays
2518 : real(dp) :: omegafl, omegasr, omegalr, sumfl,sumlr,sumsr,asr
2519 : real(dp) :: nearidentity(3,3)
2520 : ! *********************************************************************
2521 :
2522 : !write(ab_out,*)''
2523 : !write(std_out,*) 'SR/LR decomposition: enter for wavevector number :',iqpt
2524 :
2525 : !apply asr (note the lr part is already asred by construction in mkifc9)
2526 0 : do ipert1=1,natom
2527 0 : do idir1=1,3
2528 0 : do idir2=1,3
2529 : asr=0.0d0
2530 0 : do ipert2=1,natom
2531 0 : asr=asr+dynmatfl(1,idir1,ipert1,idir2,ipert2,iqpt)
2532 : end do
2533 0 : dynmatfl(1,idir1,ipert1,idir2,ipert1,iqpt)=dynmatfl(1,idir1,ipert1,idir2,ipert1,iqpt)-asr
2534 0 : dynmatsr(1,idir1,ipert1,idir2,ipert1,iqpt)=dynmatsr(1,idir1,ipert1,idir2,ipert1,iqpt)-asr
2535 : end do
2536 : end do
2537 : end do
2538 :
2539 : !This slight breaking of the symmetry allows the results to be more portable between machines
2540 0 : nearidentity(:,:)=1.0
2541 0 : nearidentity(1,1)=1.0+break_symm
2542 0 : nearidentity(3,3)=1.0-break_symm
2543 :
2544 : !Include Mass
2545 0 : do ipert1=1,natom
2546 0 : do ipert2=1,natom
2547 :
2548 0 : fac=1.0d0/sqrt(amu(typat(ipert1))*amu(typat(ipert2)))/amu_emass
2549 :
2550 0 : do idir1=1,3
2551 0 : do idir2=1,3
2552 :
2553 : dynmatfl(1,idir1,ipert1,idir2,ipert2,iqpt)=&
2554 : & dynmatfl(1,idir1,ipert1,idir2,ipert2,iqpt)*&
2555 0 : & fac*nearidentity(idir1,idir2)
2556 :
2557 : dynmatsr(1,idir1,ipert1,idir2,ipert2,iqpt)=&
2558 : & dynmatsr(1,idir1,ipert1,idir2,ipert2,iqpt)*&
2559 0 : & fac*nearidentity(idir1,idir2)
2560 :
2561 : dynmatlr(1,idir1,ipert1,idir2,ipert2,iqpt)=&
2562 : & dynmatlr(1,idir1,ipert1,idir2,ipert2,iqpt)*&
2563 0 : & fac*nearidentity(idir1,idir2)
2564 :
2565 : ! This is to break slightly the translation invariance, and make
2566 : ! the automatic tests more portable
2567 0 : if(ipert1==ipert2 .and. idir1==idir2)then
2568 : dynmatfl(1,idir1,ipert1,idir2,ipert2,iqpt)=&
2569 : & dynmatfl(1,idir1,ipert1,idir2,ipert2,iqpt)+&
2570 0 : & break_symm*natom/amu_emass/idir1*0.01d0
2571 :
2572 : dynmatsr(1,idir1,ipert1,idir2,ipert2,iqpt)=&
2573 : & dynmatsr(1,idir1,ipert1,idir2,ipert2,iqpt)+&
2574 0 : & break_symm*natom/amu_emass/idir1*0.01d0
2575 :
2576 : dynmatlr(1,idir1,ipert1,idir2,ipert2,iqpt)=&
2577 : & dynmatlr(1,idir1,ipert1,idir2,ipert2,iqpt)+&
2578 0 : & break_symm*natom/amu_emass/idir1*0.01d0
2579 : end if
2580 :
2581 : end do
2582 : end do
2583 : end do
2584 : end do
2585 :
2586 : !Calculation of <eigvec|Dyn_tot,Dyn_SR,Dyn_LR|eigenvec>=omega**2
2587 :
2588 : !write(ab_out,*)''
2589 : !write(ab_out,*)'==============================================================================='
2590 0 : write(ab_out,*)''
2591 0 : write(ab_out,*) 'Long-Range/Short-Range decomposed phonon freq. (cm-1)**2'
2592 0 : write(ab_out,*) 'at wavevector number:',iqpt
2593 0 : write(ab_out,*)''
2594 0 : write(ab_out,'(a13,1x,a16,2x,a16,2x,a16)') ' Mode number.','tot**2','SR**2','LR**2'
2595 0 : write(std_out,'(a13,1x,a16,2x,a16,2x,a16)') ' Mode number.','tot**2','SR**2','LR**2'
2596 : !write(ab_out,'(a12,2x,a10,2x,a10,2x,a10,2x,a16,2x,a16,2x,a16)') 'Mode number.','tot','SR','LR','tot**2','SR**2','LR**2'
2597 : !write(std_out,'(a12,2x,a10,2x,a10,2x,a10,2x,a16,2x,a16,2x,a16)') 'Mode number.','tot','SR','LR','tot**2','SR**2','LR**2'
2598 :
2599 0 : do imode=1,3*natom
2600 0 : sumfl=zero; sumlr=zero; sumsr=zero
2601 :
2602 0 : do ipert1=1,natom
2603 0 : do ipert2=1,natom
2604 0 : do i1=1,3
2605 0 : do i2=1,3
2606 :
2607 0 : index1=i1+(ipert1-1)*3+3*natom*(imode-1)
2608 0 : index2=i2+(ipert2-1)*3+3*natom*(imode-1)
2609 : ! MG FIXME: I don't think these expressions are correct when q != 0
2610 : ! We should also include the imaginary part
2611 :
2612 0 : sumfl = sumfl + eigenvec(2*index1-1) * dynmatfl(1,i1,ipert1,i2,ipert2,iqpt) * eigenvec(2*index2-1)
2613 0 : sumlr = sumlr + eigenvec(2*index1-1) * dynmatlr(1,i1,ipert1,i2,ipert2,iqpt) * eigenvec(2*index2-1)
2614 0 : sumsr = sumsr + eigenvec(2*index1-1) * dynmatsr(1,i1,ipert1,i2,ipert2,iqpt) * eigenvec(2*index2-1)
2615 : end do
2616 : end do
2617 : end do
2618 : end do
2619 :
2620 0 : sumfl = sumfl * Ha_cmm1 * Ha_cmm1
2621 0 : sumsr = sumsr * Ha_cmm1 * Ha_cmm1
2622 0 : sumlr = sumlr * Ha_cmm1 * Ha_cmm1
2623 :
2624 : ! Compute omega=sqrt(omega**2)
2625 : if(sumfl>=1.0d-16)then
2626 : omegafl=sqrt(sumfl)
2627 : else if(sumfl>=-1.0d-16)then
2628 : omegafl=zero
2629 : else
2630 : omegafl=-sqrt(-sumfl)
2631 : end if
2632 :
2633 : if(sumsr>=1.0d-16)then
2634 : omegasr=sqrt(sumsr)
2635 : else if(sumsr>=-1.0d-16)then
2636 : omegasr=zero
2637 : else
2638 : omegasr=-sqrt(-sumsr)
2639 : end if
2640 :
2641 : if(sumlr>=1.0d-16)then
2642 : omegalr=sqrt(sumlr)
2643 : else if(sumlr>=-1.0d-16)then
2644 : omegalr=zero
2645 : else
2646 : omegalr=-sqrt(-sumlr)
2647 : end if
2648 :
2649 : ! Output
2650 0 : write(ab_out,'(i4,10x,s,f16.4,2x,f16.4,2x,f16.4)') imode,sumfl,sumsr,sumlr !vz_d
2651 0 : write(std_out,'(i4,10x,s,f16.4,2x,f16.4,2x,f16.4)') imode,sumfl,sumsr,sumlr !vz_d
2652 : ! write(ab_out,'(i4,8x,f10.4,2x,f10.4,2x,f10.4,2x,s,f16.6,2x,f16.6,2x,f16.6)') imode,omegafl,omegasr,omegalr,sumfl,sumsr,sumlr
2653 : ! write(std_out,'(i4,8x,f10.4,2x,f10.4,2x,f10.4,2x,s,f16.6,2x,f16.6,2x,f16.6)') imode,omegafl,omegasr,omegalr,sumfl,sumsr,sumlr
2654 : end do
2655 :
2656 0 : end subroutine omega_decomp
2657 : !!***
2658 :
2659 : !----------------------------------------------------------------------
2660 :
2661 : !!****f* m_ifc/ifc_outphbtrap
2662 : !! NAME
2663 : !! ifc_outphbtrap
2664 : !!
2665 : !! FUNCTION
2666 : !! Print out phonon frequencies on regular grid for BoltzTrap
2667 : !! Flag in input file is outboltztrap=1
2668 : !!
2669 : !! INPUTS
2670 : !! Crystal<crystal_t>=Info on the crystal structure
2671 : !! basename = file name for output to disk
2672 : !! ngqpt(3)=Divisions of the q-mesh
2673 : !! nqshft=Number of shifts
2674 : !! qshft(3,nqshft)=Shifts of the q-mesh.
2675 : !!
2676 : !! OUTPUT
2677 : !! only write to file. This routine should be called by a single processor.
2678 : !!
2679 : !! SOURCE
2680 :
2681 2 : subroutine ifc_outphbtrap(ifc, cryst, ngqpt, nqshft, qshft, basename)
2682 :
2683 : !Arguments -------------------------------
2684 : !scalars
2685 : class(ifc_type),intent(in) :: ifc
2686 : integer,intent(in) :: nqshft
2687 : character(len=*),intent(in) :: basename
2688 : type(crystal_t),intent(in) :: cryst
2689 : !arrays
2690 : integer,intent(in) :: ngqpt(3)
2691 : real(dp),intent(in) :: qshft(3,nqshft)
2692 :
2693 : !Local variables -------------------------
2694 : !scalars
2695 : integer,parameter :: qptopt1=1
2696 : integer :: natom,imode,iq_ibz,nqbz,nqibz, nreals,unit_btrap,iatom,idir
2697 : character(len=500) :: msg,format_nreals,format_line_btrap
2698 : character(len=fnlen) :: outfile
2699 : !arrays
2700 : integer :: qptrlatt(3,3)
2701 4 : real(dp) :: d2cart(2,3,cryst%natom,3,cryst%natom),displ(2*3*cryst%natom*3*cryst%natom)
2702 4 : real(dp) :: phfrq(3*cryst%natom),qphon(3)
2703 2 : real(dp),allocatable :: qbz(:,:),qibz(:,:),wtq(:)
2704 : ! *********************************************************************
2705 :
2706 : DBG_ENTER("COLL")
2707 2 : natom = cryst%natom
2708 :
2709 : ! Setup IBZ, weights and BZ. Always use q --> -q symmetry for phonons even in systems wo inversion
2710 2 : qptrlatt = 0; qptrlatt(1,1) = ngqpt(1); qptrlatt(2,2) = ngqpt(2); qptrlatt(3,3) = ngqpt(3)
2711 2 : call kpts_ibz_from_kptrlatt(cryst, qptrlatt, qptopt1, nqshft, qshft, nqibz, qibz, wtq, nqbz, qbz)
2712 :
2713 2 : outfile = trim(basename) // '_BTRAP'
2714 2 : write(msg, '(3a)')ch10,' Will write phonon FREQS in BoltzTrap format to file ',trim(outfile)
2715 6 : call wrtout([std_out, ab_out], msg)
2716 :
2717 2 : if (open_file(outfile,msg,newunit=unit_btrap,status="replace") /= 0) then
2718 0 : ABI_ERROR(msg)
2719 : end if
2720 :
2721 2 : write(unit_btrap,'(a)') '#'
2722 2 : write(unit_btrap,'(a)') '# ABINIT package : Boltztrap phonon file. With old BT versions remove this header before feeding to BT'
2723 2 : write(unit_btrap,'(a)') '# for compatibility with PHON output the freq are in Ry (before the square)'
2724 2 : write(unit_btrap,'(a)') '#'
2725 2 : write(unit_btrap,'(a)') '# nq, nband '
2726 2 : write(unit_btrap,'(a)') '# qx, qy, qz '
2727 2 : write(unit_btrap,'(a)') '# qpt weight '
2728 2 : write(unit_btrap,'(a)') '# freq_1^2, dynmat column for mode 1 '
2729 2 : write(unit_btrap,'(a)') '# etc for mode 2,3,4... qpt 2,3,4... '
2730 2 : write(unit_btrap,'(2I6)') nqibz, 3*natom
2731 :
2732 : ! Loop over irreducible q-points
2733 197 : do iq_ibz=1,nqibz
2734 780 : qphon(:)=qibz(:,iq_ibz)
2735 :
2736 195 : call ifc%fourq(cryst, qphon, phfrq, displ, out_d2cart=d2cart)
2737 :
2738 195 : write (unit_btrap,'(3E20.10)') qphon
2739 195 : write (unit_btrap,'(E20.10)') wtq(iq_ibz)
2740 195 : nreals=1+2*3*natom
2741 195 : call appdig(nreals,'(',format_nreals)
2742 195 : format_line_btrap=trim(format_nreals)//'E20.10)'
2743 442 : do iatom = 1, natom
2744 1175 : do idir = 1, 3
2745 735 : imode = idir + 3*(iatom-1)
2746 : ! factor two for Ry output - this may change in definitive BT and abinit formats
2747 980 : write (unit_btrap,trim(format_line_btrap))phfrq(imode)*two,d2cart(1:2,1:3,1:natom,idir,iatom)
2748 : end do
2749 : end do
2750 :
2751 : end do !irred q-points
2752 2 : close (unit_btrap)
2753 :
2754 2 : ABI_FREE(qibz)
2755 2 : ABI_FREE(qbz)
2756 2 : ABI_FREE(wtq)
2757 :
2758 : DBG_EXIT("COLL")
2759 :
2760 2 : end subroutine ifc_outphbtrap
2761 : !!***
2762 :
2763 : !----------------------------------------------------------------------
2764 :
2765 : !!****f* m_ifc/ifc_printbxsf
2766 : !! NAME
2767 : !! ifc_printbxsf
2768 : !!
2769 : !! FUNCTION
2770 : !! Output phonon isosurface in Xcrysden format.
2771 : !!
2772 : !! INPUTS
2773 : !! crystal<crystal_t>=Info on the crystal structure
2774 : !! ngqpt(3)=Divisions of the q-mesh
2775 : !! nqshft=Number of shifts
2776 : !! qshft(3,nqshft)=Shifts of the q-mesh.
2777 : !! path=File name for output to disk
2778 : !! comm=MPI communicator.
2779 : !!
2780 : !! OUTPUT
2781 : !! Only write to file
2782 : !!
2783 : !! SOURCE
2784 :
2785 1 : subroutine ifc_printbxsf(ifc, cryst, ngqpt, nqshft, qshft, path, comm)
2786 :
2787 : !Arguments -------------------------------
2788 : !scalars
2789 : class(ifc_type),intent(in) :: ifc
2790 : integer,intent(in) :: nqshft,comm
2791 : character(len=*),intent(in) :: path
2792 : type(crystal_t),intent(in) :: cryst
2793 : !arrays
2794 : integer,intent(in) :: ngqpt(3)
2795 : real(dp),intent(in) :: qshft(3,nqshft)
2796 :
2797 : !Local variables -------------------------
2798 : !scalars
2799 : integer,parameter :: nsppol1=1,master=0,qptopt1=1
2800 : integer :: my_rank,nprocs,iq_ibz,nqibz,nqbz,ierr
2801 : character(len=500) :: msg
2802 : !arrays
2803 2 : integer :: qptrlatt(3,3),dummy_symafm(cryst%nsym)
2804 2 : real(dp) :: displ_cart(2,3*cryst%natom,3*cryst%natom)
2805 1 : real(dp),allocatable :: qibz(:,:),wtq(:),qbz(:,:),freqs_qibz(:,:)
2806 : ! *********************************************************************
2807 :
2808 1 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
2809 :
2810 : ! Setup IBZ, weights and BZ. Always use q --> -q symmetry for phonons even in systems wo inversion
2811 1 : qptrlatt = 0; qptrlatt(1,1) = ngqpt(1); qptrlatt(2,2) = ngqpt(2); qptrlatt(3,3) = ngqpt(3)
2812 1 : call kpts_ibz_from_kptrlatt(cryst, qptrlatt, qptopt1, nqshft, qshft, nqibz, qibz, wtq, nqbz, qbz)
2813 1 : ABI_FREE(qbz)
2814 1 : ABI_FREE(wtq)
2815 :
2816 : ! Compute phonon frequencies in the irreducible wedge.
2817 585 : ABI_CALLOC(freqs_qibz, (3*cryst%natom, nqibz))
2818 :
2819 146 : do iq_ibz=1,nqibz
2820 145 : if (mod(iq_ibz, nprocs) /= my_rank) cycle ! mpi parallelism.
2821 146 : call ifc%fourq(cryst, qibz(:,iq_ibz), freqs_qibz(:,iq_ibz), displ_cart)
2822 : end do
2823 1 : call xmpi_sum(freqs_qibz, comm, ierr)
2824 :
2825 : ! Output phonon isosurface.
2826 1 : if (my_rank == master) then
2827 49 : dummy_symafm = 1
2828 : call printbxsf(freqs_qibz, zero, zero, cryst%gprimd, qptrlatt, 3*cryst%natom,&
2829 1 : nqibz, qibz, cryst%nsym, .False., cryst%symrec, dummy_symafm, .True., nsppol1, qshft, nqshft, path, ierr)
2830 1 : if (ierr /=0) then
2831 0 : msg = "Cannot produce BXSF file with phonon isosurface, see log file for more info"
2832 0 : ABI_WARNING(msg)
2833 0 : call wrtout(ab_out, msg)
2834 : end if
2835 : end if
2836 :
2837 1 : ABI_FREE(freqs_qibz)
2838 1 : ABI_FREE(qibz)
2839 :
2840 1 : end subroutine ifc_printbxsf
2841 : !!***
2842 :
2843 : !----------------------------------------------------------------------
2844 :
2845 : !!****f* m_ifc/ifc_calcnwrite_nana_terms
2846 : !! NAME
2847 : !! ifc_calcnwrite_nana_terms
2848 : !!
2849 : !! FUNCTION
2850 : !! Compute frequencies and phonon displacement for q-->0 in the presence of non-analytical behaviour.
2851 : !!
2852 : !! INPUTS
2853 : !! nph2l=Number of qpoints.
2854 : !! qph2l(3,nph2l)=List of phonon wavevector directions along which the non-analytical correction
2855 : !! to the Gamma-point phonon frequencies will be calculated
2856 : !! The direction is in CARTESIAN COORDINATES
2857 : !! qnrml2(nph2l)=Normalization factor.
2858 : !!
2859 : !! OUTPUT
2860 : !! (Optional)
2861 : !! phfrq2l(3*crystal%natom,nph2l)=List of phonon frequencies
2862 : !! polarity2l(3,3*crystal%natom,nph2l)=List of mode-polarities
2863 : !! (see Eq.(41) of Veithen et al, PRB71, 125107 (2005) [[cite:Veithen2005]])
2864 : !!
2865 : !! NOTES:
2866 : !! This routine should be called by master node and when ifcflag == 1.
2867 : !!
2868 : !! SOURCE
2869 :
2870 32 : subroutine ifc_calcnwrite_nana_terms(ifc, crystal, nph2l, qph2l, &
2871 6 : qnrml2, ncid, phfrq2l, polarity2l) ! optional arguments
2872 :
2873 : !Arguments ------------------------------------
2874 : class(ifc_type),intent(in) :: ifc
2875 : integer,intent(in) :: nph2l
2876 : integer,optional,intent(in) :: ncid
2877 : type(crystal_t),intent(in) :: crystal
2878 : !arrays
2879 : real(dp),intent(in) :: qph2l(3, nph2l)
2880 : real(dp),optional,intent(in) :: qnrml2(nph2l)
2881 : real(dp),optional,intent(out) :: phfrq2l(3*crystal%natom,nph2l), polarity2l(3,3*crystal%natom,nph2l)
2882 :
2883 : !Local variables-------------------------------
2884 : !scalars
2885 : integer :: iatom,idir,imode,iphl2, ncerr
2886 : real(dp) :: eta
2887 : !arrays
2888 : real(dp) :: qphnrm(3),qphon(3,3)
2889 32 : real(dp),allocatable :: displ_cart(:,:,:),phfrq(:),d2cart(:,:,:),eigvec(:,:,:),eigval(:)
2890 : ! ************************************************************************
2891 :
2892 32 : if (nph2l == 0) return
2893 :
2894 : ! Now treat the second list of vectors (only at the Gamma point, but can include non-analyticities)
2895 96 : ABI_MALLOC(phfrq, (3*crystal%natom))
2896 128 : ABI_MALLOC(displ_cart, (2, 3*crystal%natom, 3*crystal%natom))
2897 128 : ABI_MALLOC(d2cart, (2, 3*ifc%mpert, 3*ifc%mpert))
2898 96 : ABI_MALLOC(eigvec, (2, 3*crystal%natom, 3*crystal%natom))
2899 64 : ABI_MALLOC(eigval, (3*crystal%natom))
2900 :
2901 : ! Before examining every direction or the dielectric tensor, generates the dynamical matrix at gamma
2902 128 : qphon(:,1)=zero; qphnrm(1)=zero
2903 :
2904 : ! Generation of the dynamical matrix in cartesian coordinates
2905 : ! Get d2cart using the interatomic forces and the long-range coulomb interaction through Ewald summation
2906 : call gtdyn9(ifc%acell,ifc%atmfrc,ifc%dielt,ifc%dipdip, &
2907 : ifc%dyewq0,d2cart,crystal%gmet,ifc%gprim,ifc%mpert,crystal%natom, &
2908 : ifc%nrpt,qphnrm(1),qphon,crystal%rmet,ifc%rprim,ifc%rpt, &
2909 32 : ifc%trans,crystal%ucvol,ifc%wghatm,crystal%xred,ifc%zeff,ifc%qdrp_cart,ifc%ewald_option,eta,xmpi_comm_self,ifc%sys_dim)
2910 :
2911 32 : if (present(ncid)) then
2912 29 : iphl2 = 0
2913 29 : call defwrite_nonana_terms(ncid, iphl2, nph2l, qph2l, crystal%natom, phfrq, displ_cart, mode="define")
2914 : ! Add epsinf, Born effective charges and some useful metadata.
2915 : ncerr = nctk_def_arrays(ncid, [ &
2916 : nctkarr_t('emacro_cart', "dp", 'number_of_cartesian_directions, number_of_cartesian_directions'), &
2917 : nctkarr_t('becs_cart', "dp", "number_of_cartesian_directions, number_of_cartesian_directions, number_of_atoms")], &
2918 87 : defmode=.True.)
2919 29 : NCF_CHECK(ncerr)
2920 : ! TODO chneut is missing
2921 : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: &
2922 116 : "asr", "dipdip", "symdynmat"])
2923 29 : NCF_CHECK(ncerr)
2924 29 : NCF_CHECK(nctk_set_datamode(ncid))
2925 29 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'emacro_cart'), ifc%dielt))
2926 29 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'becs_cart'), ifc%zeff))
2927 : ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: &
2928 : "asr", "dipdip", "symdynmat"], &
2929 203 : [ifc%asr, ifc%dipdip, ifc%symdynmat])
2930 29 : NCF_CHECK(ncerr)
2931 : end if
2932 :
2933 : ! Examine every wavevector of this list
2934 60069 : do iphl2=1,nph2l
2935 : ! Initialisation of the phonon wavevector
2936 240148 : qphon(:,1) = qph2l(:,iphl2)
2937 60037 : qphnrm(1)=zero
2938 60037 : if(present(qnrml2)) qphnrm(1) = qnrml2(iphl2)
2939 :
2940 : ! Calculation of the eigenvectors and eigenvalues of the dynamical matrix
2941 : call dfpt_phfrq(ifc%amu,displ_cart,d2cart,eigval,eigvec,crystal%indsym, &
2942 : ifc%mpert,crystal%nsym,crystal%natom,crystal%nsym,crystal%ntypat,phfrq,qphnrm(1),qphon, &
2943 60037 : crystal%rprimd,ifc%symdynmat,crystal%symrel,crystal%symafm,crystal%typat,crystal%ucvol)
2944 :
2945 : ! Write the phonon frequencies
2946 : !call dfpt_prtph(displ_cart,inp%eivec,inp%enunit,ab_out,natom,phfrq,qphnrm(1),qphon)
2947 :
2948 480037 : if(present(phfrq2l)) phfrq2l(:,iphl2)=phfrq(:)
2949 :
2950 60037 : if(present(polarity2l))then
2951 1740000 : polarity2l(:,:,iphl2)=zero
2952 480000 : do imode=1,3*crystal%natom
2953 1500000 : do iatom=1,crystal%natom
2954 4500000 : do idir=1,3
2955 13260000 : polarity2l(:,imode,iphl2)=polarity2l(:,imode,iphl2)+ifc%zeff(:,idir,iatom)*displ_cart(1,idir+(iatom-1)*3,imode)
2956 : enddo
2957 : enddo
2958 : enddo
2959 : endif
2960 :
2961 60069 : if(present(ncid))then
2962 : ! Loop is not MPI-parallelized --> no need for MPI-IO API.
2963 37 : call defwrite_nonana_terms(ncid, iphl2, nph2l, qph2l, crystal%natom, phfrq, displ_cart, mode="write")
2964 : endif
2965 : end do ! iphl2
2966 :
2967 32 : ABI_FREE(phfrq)
2968 32 : ABI_FREE(displ_cart)
2969 32 : ABI_FREE(d2cart)
2970 32 : ABI_FREE(eigvec)
2971 32 : ABI_FREE(eigval)
2972 :
2973 35 : end subroutine ifc_calcnwrite_nana_terms
2974 : !!***
2975 :
2976 : !----------------------------------------------------------------------
2977 :
2978 1 : subroutine ifc_calcnwrite_nana_terms_qpath(ifc, qpath, cryst, ncid, units)
2979 :
2980 : !Arguments ------------------------------------
2981 : class(ifc_type),intent(in) :: ifc
2982 : integer,intent(in) :: ncid, units(:)
2983 : type(kpath_t),intent(in) :: qpath
2984 : type(crystal_t),intent(in) :: cryst
2985 :
2986 : !Local variables-------------------------------
2987 : !scalars
2988 : integer :: iq, nq_vers, natom3, nu
2989 : character(len=500) :: msg
2990 : !arrays
2991 2 : real(dp) :: phfreqs(3*cryst%natom), displ_cart(2,3,cryst%natom,3*cryst%natom)
2992 1 : real(dp),allocatable :: qvers_red(:,:), qvers_cart(:,:)
2993 : ! ************************************************************************
2994 :
2995 1 : natom3 = 3 * cryst%natom
2996 1 : call qpath%get_versors(nq_vers, qvers_red, qvers_cart)
2997 :
2998 3 : do iq=1,nq_vers
2999 2 : call ifc%fourq(cryst, qvers_red(:,iq), phfreqs, displ_cart, nanaqdir="reduced")
3000 :
3001 2 : call wrtout(units, sjoin("Phonon frequencies in meV along reduced direction:", ktoa(qvers_red(:,iq)), "(reciprocal space)"))
3002 2 : if (iq == 1) then
3003 1 : call defwrite_nonana_terms(ncid, iq, nq_vers, qvers_cart, cryst%natom, phfreqs, displ_cart, mode="define")
3004 : end if
3005 2 : call defwrite_nonana_terms(ncid, iq, nq_vers, qvers_cart, cryst%natom, phfreqs, displ_cart, mode="write")
3006 :
3007 15 : do nu=1,natom3
3008 12 : write(msg, "(1x,i0, es16.6)") nu, phfreqs(nu) * Ha_meV
3009 14 : call wrtout(units, msg)
3010 : end do
3011 : end do
3012 :
3013 1 : ABI_FREE(qvers_red)
3014 1 : ABI_FREE(qvers_cart)
3015 :
3016 1 : end subroutine ifc_calcnwrite_nana_terms_qpath
3017 : !!***
3018 :
3019 : !----------------------------------------------------------------------
3020 :
3021 : !!****f* m_ifc/ifc_to_ddb
3022 : !! NAME
3023 : !! ifc_to_ddb
3024 : !!
3025 : !! FUNCTION
3026 : !! Construct a DDB object from the IFC object.
3027 : !!
3028 : !! INPUTS
3029 : !! crystal = Crystal object
3030 : !!
3031 : !! OUTPUT
3032 : !! ddb = DDB object
3033 : !!
3034 : !! SOURCE
3035 :
3036 44 : subroutine ifc_to_ddb(ifc, ddb, crystal)
3037 :
3038 : !Arguments -----------------------------------
3039 : !scalars
3040 : class(ifc_type),intent(inout) :: ifc
3041 : type(ddb_type),intent(out) :: ddb
3042 : type(crystal_t),intent(in) :: crystal
3043 :
3044 : !Local variables-------------------------------
3045 : !scalars
3046 : integer,parameter :: qptopt=1, nqshft=1
3047 : integer :: jj,iqpt,idir1,idir2,ipert1,ipert2, mpert,msize, nqibz,nqbz
3048 : real(dp) :: qj,qptnrm
3049 : !arrays
3050 : integer :: qptrlatt(3,3)
3051 44 : integer,allocatable :: flg(:,:,:,:), flg_gamma(:,:,:,:)
3052 : real(dp) :: qpt(3), qshft(1,3)
3053 44 : real(dp),allocatable :: qibz(:,:),qbz(:,:), wtq(:), d2cart(:,:,:,:,:),d2red(:,:,:,:,:)
3054 : ! *********************************************************************
3055 :
3056 : ! Initialize DDB with minimal info
3057 : ! ================================
3058 :
3059 : ! Crystal info
3060 44 : ddb%natom = crystal%natom
3061 44 : ddb%ntypat = crystal%ntypat
3062 572 : ddb%rprim = ifc%rprim
3063 572 : ddb%gprim = ifc%gprim
3064 176 : ddb%acell = ifc%acell
3065 44 : call alloc_copy(crystal%amu, ddb%amu)
3066 :
3067 : ! Block size
3068 : ! GA: TODO need to compute mpert in a more transparent way.
3069 : ! Also, mpert should be consistent with ifc_init.
3070 44 : mpert = ddb%natom + 6
3071 44 : msize = 3*mpert*3*mpert
3072 44 : ddb%mpert = mpert
3073 44 : ddb%msize = msize
3074 :
3075 : ! Generate q-points
3076 44 : qshft = zero
3077 44 : qptrlatt = zero
3078 44 : qptrlatt(1,1)=ifc%ngqpt(1)
3079 44 : qptrlatt(2,2)=ifc%ngqpt(2)
3080 44 : qptrlatt(3,3)=ifc%ngqpt(3)
3081 :
3082 44 : call kpts_ibz_from_kptrlatt(crystal,qptrlatt,qptopt,nqshft,qshft, nqibz,qibz,wtq,nqbz,qbz)
3083 :
3084 44 : ddb%nblok = nqibz
3085 2166087 : ABI_MALLOC(ddb%flg,(msize,nqibz)) ; ddb%flg = one
3086 11911 : ABI_MALLOC(ddb%nrm,(3,nqibz)) ; ddb%nrm = zero ; ddb%nrm(1,:) = one
3087 23602 : ABI_MALLOC(ddb%qpt,(9,nqibz)) ; ddb%qpt = zero
3088 9476 : ABI_MALLOC(ddb%omega,(3,nqibz)) ; ddb%omega = zero
3089 6493215 : ABI_MALLOC(ddb%val,(2,msize,nqibz)); ddb%val = zero
3090 132 : ABI_MALLOC(ddb%typ,(nqibz))
3091 :
3092 2391 : do iqpt=1,nqibz
3093 2347 : ddb%typ(iqpt) = BLKTYP_d2E_ns
3094 9432 : do jj = 1,3
3095 7041 : qj = qibz(jj,iqpt)
3096 7041 : if (abs(qj).lt.tol8) qj = zero
3097 9388 : ddb%qpt(jj,iqpt) = qj
3098 : end do
3099 : end do
3100 :
3101 44 : ABI_FREE(wtq)
3102 44 : ABI_FREE(qibz)
3103 44 : ABI_FREE(qbz)
3104 :
3105 : ! Compute DDB by Fourier transform of the IFC
3106 : ! ===========================================
3107 :
3108 44 : qptnrm = one
3109 220 : ABI_MALLOC(d2cart,(2,3,mpert,3,mpert))
3110 132 : ABI_MALLOC(d2red,(2,3,mpert,3,mpert))
3111 :
3112 : ! Set up the flags
3113 42184 : ABI_CALLOC(flg,(3,mpert,3,mpert))
3114 42096 : ABI_CALLOC(flg_gamma,(3,mpert,3,mpert))
3115 150 : do ipert1=1,ddb%natom
3116 664 : do ipert2=1,ddb%natom
3117 2162 : do idir1=1,3
3118 6682 : do idir2=1,3
3119 4626 : flg(idir1,ipert1,idir2,ipert2) = one
3120 6168 : flg_gamma(idir1,ipert1,idir2,ipert2) = one
3121 : end do
3122 : end do
3123 : end do
3124 : end do
3125 :
3126 : ! Activate dielt and zeff flag if we have them
3127 44 : if (ifc%dipdip > 0) then
3128 :
3129 : ! Dielectric tensor
3130 5 : ipert1 = ddb%natom + 2; ipert2 = ddb%natom + 2;
3131 20 : do idir1=1,3
3132 65 : do idir2=1,3
3133 60 : flg_gamma(idir1,ipert1,idir2,ipert2) = one
3134 : end do
3135 : end do
3136 :
3137 : ! Born effective charges
3138 16 : ipert1 = ddb%natom + 2
3139 16 : do ipert2=1,ddb%natom
3140 49 : do idir1=1,3
3141 143 : do idir2=1,3
3142 99 : flg_gamma(idir1,ipert1,idir2,ipert2) = one
3143 132 : flg_gamma(idir1,ipert2,idir2,ipert1) = one
3144 : end do
3145 : end do
3146 : end do
3147 : end if
3148 :
3149 2391 : do iqpt=1,DDB%nblok
3150 :
3151 9388 : qpt(:) = DDB%qpt(1:3,iqpt)
3152 :
3153 : ! Get d2cart using the interatomic forces and the
3154 : ! long-range coulomb interaction through Ewald summation
3155 : call gtdyn9(ddb%acell,ifc%atmfrc,ifc%dielt,ifc%dipdip,ifc%dyewq0,d2cart, &
3156 : crystal%gmet,ddb%gprim,ddb%mpert,crystal%natom,ifc%nrpt,qptnrm,qpt,&
3157 : crystal%rmet,ddb%rprim,ifc%rpt,ifc%trans,crystal%ucvol, &
3158 : ifc%wghatm,crystal%xred,ifc%zeff,ifc%qdrp_cart,ifc%ewald_option, ifc%eta,&
3159 2347 : xmpi_comm_self,ifc%sys_dim)
3160 :
3161 : ! Impose the acoustic sum rule
3162 : !asrq0 = ddb%get_asrq0(1,1,crystal%xcart)
3163 : !call asrq0%apply(ddb%natom, ddb%mpert, ddb%msize, crystal%xcart, d2cart)
3164 :
3165 : ! Transform d2cart into reduced coordinates.
3166 : call d2cart_to_red(d2cart,d2red,crystal%gprimd,crystal%rprimd,ddb%mpert, &
3167 2347 : crystal%natom,crystal%ntypat,crystal%typat,crystal%ucvol,crystal%zion)
3168 :
3169 : ! Check if this q-point is Gamma
3170 2347 : call gamma9(jj, qpt, qptnrm, DDB_QTOL)
3171 :
3172 2391 : if (jj == 1) then
3173 44 : call ddb%set_d2matr(iqpt, d2red, flg_gamma)
3174 : else
3175 2303 : call ddb%set_d2matr(iqpt, d2red, flg)
3176 : end if
3177 :
3178 : end do
3179 :
3180 44 : ABI_FREE(d2cart)
3181 44 : ABI_FREE(d2red)
3182 44 : ABI_FREE(flg)
3183 44 : ABI_FREE(flg_gamma)
3184 :
3185 44 : end subroutine ifc_to_ddb
3186 : !!***
3187 :
3188 : !!****f* m_ifc/defwrite_nonana_terms
3189 : !! NAME
3190 : !! defwrite_nonana_terms
3191 : !!
3192 : !! FUNCTION
3193 : !! Write to ncfile the phonon frequencies and displacements for q --> 0 in the presence of non-analytical behaviour.
3194 : !!
3195 : !! INPUTS
3196 : !! ncid=netcdf file id.
3197 : !! iq_dir=Index of the q-point to be written to file
3198 : !! ndirs=Number of qpoints.
3199 : !! qdirs_cart(3,ndirs)=List of phonon wavevector directions along which the non-analytical correction
3200 : !! to the Gamma-point phonon frequencies will be calculated. The direction is in CARTESIAN COORDINATES
3201 : !! natom=Number of atoms
3202 : !! phfrq(3*natom)=Phonon frequencies in Ha
3203 : !! cart_displ(2,3*natom,3*natom)=displacements in CARTESIAN coordinates.
3204 : !!
3205 : !! OUTPUT
3206 : !! Only writing.
3207 : !!
3208 : !! SOURCE
3209 :
3210 130 : subroutine defwrite_nonana_terms(ncid, iq_dir, ndirs, qdirs_cart, natom, phfrq, cart_displ, mode)
3211 :
3212 : !Arguments ------------------------------------
3213 : !scalars
3214 : integer,intent(in) :: ncid,iq_dir,ndirs,natom
3215 : character(len=*),intent(in) :: mode
3216 : !arrays
3217 : real(dp),intent(in) :: qdirs_cart(3, ndirs), phfrq(3*natom), cart_displ(2,3*natom,3*natom)
3218 :
3219 : !Local variables-------------------------------
3220 : integer :: ncerr, na_phmodes_varid, na_phdispl_varid
3221 : ! *************************************************************************
3222 :
3223 54 : select case (mode)
3224 : case ("define")
3225 : !NCF_CHECK(nctk_def_basedims(ncid, defmode=.True.))
3226 108 : ncerr = nctk_def_dims(ncid, [nctkdim_t("number_of_non_analytical_directions", ndirs)], defmode=.True.)
3227 54 : NCF_CHECK(ncerr)
3228 :
3229 : ncerr = nctk_def_arrays(ncid, [&
3230 : nctkarr_t('non_analytical_directions', "dp", "number_of_cartesian_directions, number_of_non_analytical_directions"),&
3231 : nctkarr_t('non_analytical_phonon_modes', "dp", "number_of_phonon_modes, number_of_non_analytical_directions"),&
3232 : nctkarr_t('non_analytical_phdispl_cart', "dp", &
3233 216 : "two, number_of_phonon_modes, number_of_phonon_modes, number_of_non_analytical_directions")])
3234 54 : NCF_CHECK(ncerr)
3235 :
3236 54 : NCF_CHECK(nctk_set_datamode(ncid))
3237 54 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "non_analytical_directions"), qdirs_cart))
3238 :
3239 : case ("write")
3240 76 : NCF_CHECK(nf90_inq_varid(ncid, "non_analytical_phonon_modes", na_phmodes_varid))
3241 1154 : NCF_CHECK(nf90_put_var(ncid,na_phmodes_varid,phfrq*Ha_eV,start=[1, iq_dir], count=[3*natom, 1]))
3242 76 : NCF_CHECK(nf90_inq_varid(ncid, "non_analytical_phdispl_cart", na_phdispl_varid))
3243 : ncerr = nf90_put_var(ncid,na_phdispl_varid,cart_displ*Bohr_Ang,&
3244 33426 : start=[1,1,1,iq_dir], count=[2,3*natom,3*natom, 1])
3245 76 : NCF_CHECK(ncerr)
3246 :
3247 : case default
3248 130 : ABI_ERROR(sjoin("Wrong value for mode", mode))
3249 : end select
3250 :
3251 130 : end subroutine defwrite_nonana_terms
3252 : !!***
3253 :
3254 : !----------------------------------------------------------------------
3255 :
3256 0 : end module m_ifc
3257 : !!***
|
