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