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