Branch data Line data Source code
1 : : !!****m* ABINIT/m_gstore_converters
2 : : !! NAME
3 : : !! m_gstore_converters
4 : : !!
5 : : !! FUNCTION
6 : : !! Convert data from gstore.nc to other formats.
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2008-2026 ABINIT group (MG)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! SOURCE
15 : :
16 : : #if defined HAVE_CONFIG_H
17 : : #include "config.h"
18 : : #endif
19 : :
20 : : #include "abi_common.h"
21 : :
22 : : module m_gstore_converters
23 : :
24 : : use defs_basis
25 : : use m_abicore
26 : : use m_xmpi
27 : : use m_errors
28 : : use m_clib
29 : :
30 : : use m_io_tools, only : open_file
31 : : use m_dtset, only : dataset_type
32 : : use m_dtfil, only : datafiles_type
33 : : use m_ddb_hdr, only : ddb_hdr_type
34 : : use m_hdr, only : hdr_type
35 : : use m_fstrings, only : sjoin, itoa, strcat
36 : : use m_crystal, only : crystal_t
37 : : use m_ebands, only : ebands_t, gaps_t
38 : : use m_ifc, only : ifc_type
39 : : use m_gstore, only : gstore_t, GSTORE_GMODE_ATOM, GSTORE_GTYPE_GWPT, gstore_read_gtype
40 : :
41 : : implicit none
42 : :
43 : : private
44 : :
45 : : public :: gstore_convert
46 : : ! Convert data from gstore.nc to other formats
47 : :
48 : : !!***
49 : :
50 : : contains
51 : : !!***
52 : :
53 : : !----------------------------------------------------------------------
54 : :
55 : : !!****f* m_gstore_converters/gstore_convert
56 : : !! NAME
57 : : !! gstore_convert
58 : : !!
59 : : !! FUNCTION
60 : : !! Convert data from gstore.nc to other formats
61 : : !!
62 : : !! INPUTS
63 : : !! gstore_path=Filename of the output GSTORE.nc file
64 : : !!
65 : : !! SOURCE
66 : :
67 : 0 : subroutine gstore_convert(gstore_path, dtset, dtfil, cryst, ebands, ifc, comm)
68 : :
69 : : !Arguments ------------------------------------
70 : : !scalars
71 : : character(len=*),intent(in) :: gstore_path
72 : : type(dataset_type),target,intent(in) :: dtset
73 : : type(datafiles_type),intent(in) :: dtfil
74 : : class(crystal_t),target,intent(in) :: cryst
75 : : class(ebands_t),target,intent(in) :: ebands
76 : : class(ifc_type),target,intent(in) :: ifc
77 : : integer,intent(in) :: comm
78 : :
79 : : !Local variables-------------------------------
80 : : !scalars
81 : : integer :: nprocs, my_rank, nsppol, spin, nmodes, this_comm, unt, ib, nu, i, j, ierr
82 : : integer :: with_cplex, ik_ibz, my_is, my_ik, my_iq, iq_glob, natom, itypat, lstr_j
83 : : integer :: ibrav, idir, jdir, iat, ipert, unt_ascii, ik_glob, band_kq, band_k, mu ! jat,
84 : : logical :: with_g2dw, q_is_gamma, lborn, ascii_write
85 : : real(dp),parameter :: Ha2Ry = two
86 : : real(dp) :: weight_qq
87 : : character(len=5000) :: msg
88 : : character(len=abi_slen) :: with_gmode, gvals_name, gtype
89 : : character(len=fnlen) :: fname, fname_ascii, elphmat_dir, prefix
90 : : character(len=3) :: band_i
91 [ # # ]: 0 : type(gstore_t) :: gstore
92 : : !arrays
93 : : integer :: units(2)
94 : : real(dp) :: qpt(3), kk_bz(3), kk_ibz(3), celldm(6)
95 : 0 : real(dp),allocatable :: tau_cart(:,:)
96 : 0 : complex(dp),allocatable :: dyn_qe(:,:), g_cart(:,:,:)
97 [ # # ]: 0 : character(len=3) :: atm(cryst%ntypat)
98 : : !----------------------------------------------------------------------
99 : :
100 : : ! Only master works here as performance is not crucial.
101 : : ! In principle one could activate the q-point/spin parallelism just to distribute
102 : : ! the memory for the g's to avoid going OOM.
103 : 0 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
104 [ # # ]: 0 : if (my_rank /= 0) return
105 : :
106 [ # # ]: 0 : units = [std_out, ab_out]
107 : 0 : natom = cryst%natom; nmodes = 3 * cryst%natom; nsppol = ebands%nsppol; this_comm = xmpi_comm_self
108 : :
109 : : ! QE Bravais lattice metadata.
110 : : ! We do not classify the QE Bravais lattice from the ABINIT cell, hence ibrav = 0
111 : : ! (free lattice) and celldm(1) = alat = |a1| in Bohr (the remaining celldm are unused
112 : : ! for ibrav == 0). The full lattice is otherwise defined by the crystal structure.
113 : 0 : ibrav = 0
114 : 0 : celldm(:) = zero
115 [ # # ]: 0 : celldm(1) = sqrt(sum(cryst%rprimd(:, 1) ** 2))
116 : :
117 : : ! Work array for the phonon displacement matrix (per q-point).
118 [ # # # # : 0 : ABI_MALLOC(dyn_qe, (nmodes, nmodes))
# # # # #
# # # # #
# # ]
119 : :
120 : : ! QE tau: ionic positions in Cartesian coordinates, in units of alat (= celldm(1)).
121 [ # # # # : 0 : ABI_MALLOC(tau_cart, (3, natom))
# # # # #
# # # ]
122 [ # # ]: 0 : do iat=1,natom
123 [ # # # # : 0 : tau_cart(:, iat) = matmul(cryst%rprimd, cryst%xred(:, iat)) / celldm(1)
# # # # #
# ]
124 : : end do
125 : :
126 : : ! If ascii_write is .True. a human-readable copy of each binary elph.mat.q_* file is
127 : : ! also written, with the same content (record by record) and a ".ascii" suffix.
128 : : ! Toggle this flag (or wire it to an input variable) to disable the extra files.
129 : : ascii_write = .True.
130 : 0 : ascii_write = .False.
131 : :
132 : : ! Preliminary consistency check.
133 : 0 : call wrtout(units, sjoin(" Begin conversion GSTORE --> ", dtset%gstore_convert))
134 [ # # ]: 0 : ABI_CHECK(gstore_path /= ABI_NOFILE, sjoin("Invalid gstore_path:", gstore_path))
135 [ # # ]: 0 : ABI_CHECK(dtset%gstore_convert == "epiq", "only gstore_convert == 'epiq' is supported.")
136 : : !ABI_CHECK(.not. cryst%isalchemical(), "'epiq' format does not support alchemical pseudos.")
137 : :
138 : : ! Read g(k,q) from GSTORE and store them in gqk%my_g(nu, im_kq, my_iq, in_k, ik)
139 : : !
140 : : ! Shape of array is:
141 : : !
142 : : ! (my_npert, nb_kq, my_nq, nb_k, my_nk)
143 : : ! ( p, b1_kq, q, b2_k, k) --> <k+q, b1| D_{q,p}H |k, b2>
144 : :
145 : : ! The g's are complex and in the phonon representation.
146 : : ! Note <k+q| for the final state and |k> for the initial state. I guess epiq uses the same convention.
147 : : ! All quantities are in atomic units (Hartree and Bohr).
148 : : !
149 : : ! The gstore file produced by the GWPT code has both GWPT and KS g.
150 : : ! In this case, we convert the GWPT matrix elements but one can still select
151 : : ! the KS e-ph vertex via gstore_gname.
152 : :
153 : 0 : call gstore_read_gtype(gstore_path, gtype, this_comm)
154 : 0 : gvals_name = "gvals"
155 [ # # # # ]: 0 : if (gtype == GSTORE_GTYPE_GWPT .and. dtset%gstore_gname == "gvals_ks") gvals_name = "gvals_ks"
156 : :
157 : : ! Request the ATOM representation: gstore.nc stores g in this representation (the bare
158 : : ! deformation potential w.r.t. reduced atomic displacements), so no atom->phonon
159 : : ! conversion is performed. This avoids ephtk_gkknu_from_atm, which would otherwise zero
160 : : ! the acoustic/imaginary modes (phfrq < EPHTK_WTOL) and divide by sqrt(2*omega).
161 : 0 : with_cplex = 2; with_gmode = GSTORE_GMODE_ATOM; with_g2dw = .False.
162 : :
163 : : call gstore%from_ncpath(gstore_path, with_cplex, dtset, dtfil, cryst, ebands, ifc, &
164 : 0 : with_gmode, gvals_name, with_g2dw, this_comm)
165 : :
166 : : ! Consistency check.
167 [ # # ]: 0 : ABI_CHECK(nsppol == 1, "Don't know how to convert spin-polarized g to epiq format!")
168 : :
169 : : ! For wannierization, we need the same number of bands for m and n.
170 : : ! Also, k and q must be in the BZ without any filter.
171 : : ! Once the symmetrization of the g's has been implemented, this routine
172 : : ! will receive a gstore file in which all g(k,q) matrix elements in the BZ
173 : : ! have been reconstructed using symmetry operations.
174 [ # # ]: 0 : ABI_CHECK(gstore%same_nbands(msg), msg)
175 [ # # ]: 0 : if (gstore%check_cplex_qkzone_gmode(2, "bz", "bz", "atom", kfilter="none") /= 0) then
176 : 0 : ABI_ERROR("GSTORE.nc should have both k and q in the full BZ. See messages above.")
177 : : end if
178 : :
179 : : ! Create directory to store output files.
180 : 0 : prefix = "epiq"
181 [ # # ]: 0 : elphmat_dir = strcat(dtfil%filnam_ds(4), "_", prefix)
182 : 0 : call wrtout(units, sjoin(" Output files written to directory:", elphmat_dir))
183 : 0 : call execute_command_line(sjoin("rm -rf", elphmat_dir), exitstat=ierr)
184 : 0 : call clib_mkdir_if_needed(elphmat_dir, ierr)
185 [ # # ]: 0 : ABI_CHECK(ierr == 0, "mkdir returned ierr /= 0")
186 : :
187 : : ! Write the EPIQ input namelist (&Diff_Start_Param + KPOINTS), the analog of QE's
188 : : ! print_ph_input2epiq called inside ep_matrix_element_wannier.
189 : 0 : call write_epiq_input(ebands, dtfil, strcat(elphmat_dir, "/scf_dfpt.2epiq.in"))
190 : 0 : call wrtout(units, sjoin(" EPIQ input namelist written to:", strcat(elphmat_dir, "/scf_dfpt.2epiq.in")))
191 : :
192 : : ! Write the dynq0 file: q-mesh, number of irreducible q-points and their positions
193 : : ! in Cartesian coordinates (2pi/alat), the analog of QE's dynq0 output.
194 : 0 : call write_dynq0(gstore%ngqpt, gstore%qibz, cryst%gprimd, celldm(1), strcat(elphmat_dir, "/dynq0"))
195 : 0 : call wrtout(units, sjoin(" EPIQ dynq0 file written to:", strcat(elphmat_dir, "/dynq0")))
196 : :
197 : : ! NB: atm is character(len=3) while symbol_type returns character(len=2).
198 [ # # ]: 0 : do itypat=1, cryst%ntypat
199 : 0 : atm(itypat)(1:2) = cryst%symbol_type(itypat)
200 : 0 : atm(itypat)(3:3) = ""
201 : : end do
202 : :
203 : : ! Write one dynamical-matrix file (dynq<iq>) per irreducible q-point, in QE format
204 : : ! (same ordering as the dynq0 list).
205 [ # # ]: 0 : do iq_glob=1,gstore%nqibz
206 : : call write_dynq(cryst, ifc, gstore%qibz(:,iq_glob), celldm(1), atm, &
207 : 0 : strcat(elphmat_dir, "/dynq", itoa(iq_glob)))
208 : : end do
209 : 0 : call wrtout(units, sjoin(" EPIQ dynq<iq> dynamical-matrix files written to directory:", elphmat_dir))
210 : :
211 : : ! TODO: Need helper function to get ibrav, celldm from crystal
212 : : ! Some input from the EPIC developers would be greatly appreciated.
213 : : !call cryst%get_ibrav_celldm(ibrav, celldm)
214 : :
215 : : ! Loop over collinear spins.
216 [ # # ]: 0 : do my_is=1,gstore%my_nspins
217 : 0 : spin = gstore%my_spins(my_is)
218 : 0 : associate (gqk => gstore%gqk(my_is))
219 : :
220 : : ! Buffer for the e-ph matrix elements rotated to the Cartesian atomic-displacement basis.
221 [ # # # # : 0 : ABI_MALLOC(g_cart, (nmodes, gqk%nb_kq, gqk%nb_k))
# # # # #
# # # # #
# # # # ]
222 : :
223 : : !num_bands = gqk%nb_k
224 : :
225 : : ! Loop over q-points in the BZ.
226 [ # # ]: 0 : do my_iq=1, gqk%my_nq
227 : 0 : iq_glob = my_iq + gqk%my_qstart - 1
228 : :
229 : 0 : call gqk%myqpt(my_iq, gstore, weight_qq, qpt); q_is_gamma = sum(qpt**2) < tol14
230 : 0 : call define_band_string(iq_glob, band_i, lstr_j)
231 : :
232 [ # # ]: 0 : fname = trim(elphmat_dir)//"/"//trim(prefix)//'_elph.mat.q_'//band_i(1:lstr_j)
233 [ # # ]: 0 : if (open_file(fname, msg, newunit=unt, form="unformatted", status="unknown", action="write") /= 0) then
234 : 0 : ABI_ERROR(msg)
235 : : end if
236 : :
237 : : ! Optional human-readable companion file with the same content.
238 [ # # ]: 0 : if (ascii_write) then
239 [ # # ]: 0 : fname_ascii = trim(fname)//".ascii"
240 [ # # ]: 0 : if (open_file(fname_ascii, msg, newunit=unt_ascii, form="formatted", status="unknown", action="write") /= 0) then
241 : 0 : ABI_ERROR(msg)
242 : : end if
243 : : end if
244 : :
245 : : ! EPIQ format. See https://gitlab.com/the-epiq-team/epiq/-/blob/develop/src/io_matelem.F90
246 : : !read(unt) (xq_r(j),j=1,3)
247 : : !if (.not.fet) read(unt) noncolin, nspin, lborn !REMOVE FOR VERSION 5.1 FET
248 : : !read(unt) nel_aux
249 : : !read(unt) nbnd_min, nbnd_max, nbnd_r !nbnd_r = total number of bands in pw
250 : : !read(unt) nmodes, nk_r, nat, ntyp
251 : : !read(unt) ibrav,(celldm(j), j=1,6)
252 : : !read(unt) (atm(j),j=1,ntyp),(amass(j),j=1,ntyp), &
253 : : ! (ityp(j),j=1,nat),((tau(j,i),j=1,3),i=1,nat)
254 : : !read(unt) (w2 (nu,iqph) , nu=1,nmodes)
255 : : !read(unt) ((zz(i,j,iqph), i=1,nmodes),j=1,nmodes) !eigenvectors in the QE basis
256 : : !read(unt) ((dyn(i,j,iqph),i=1,nmodes),j=1,nmodes) ! eigenvectors divided by masses
257 : : !do k=1,num_kpts
258 : : ! read(unt) (xk_r(i,k),i=1,3)
259 : : ! read(unt) (eig(i,k),i=1,num_bands)
260 : : ! do nu=1,nmodes
261 : : ! write(unt) ((g_matrix(j, i, nu, k, iq),j=1,num_bands),i=1,num_bands)
262 : : ! end do
263 : : !end do
264 : :
265 : : ! Conventions used below (matched to QE pw/PHonon elphsum_wannier and the EPIQ reader):
266 : : ! - w2 is omega^2 (squared phonon frequency) in Ry^2, signed (< 0 for imaginary modes).
267 : : ! - lborn = .False.: BECS / dynamical quadrupoles are not exported yet.
268 : : ! - tau are atom positions in Cartesian coordinates, in units of alat (= xcart/alat).
269 : : ! - nel_aux (ebands%nelect) is the (real) number of electrons including possible doping.
270 : : ! - The e-ph matrix elements are written in the Cartesian atomic-displacement basis
271 : : ! (as in QE ep_matrix_element_wannier). gstore provides g in the atom representation
272 : : ! w.r.t. *reduced* atomic displacements; we rotate it to Cartesian directions with
273 : : ! gprimd (a per-atom 3x3 map). No mode/frequency factors are involved, so no mode is
274 : : ! zeroed. Consequently zz (the QE pattern matrix) is the identity, while dyn keeps the
275 : : ! Cartesian phonon displacements (= my_displ_cart).
276 : 0 : lborn = .False.
277 : :
278 : : ! Record 1: q-point in reduced (crystal) coordinates.
279 : 0 : write(unt) qpt
280 [ # # ]: 0 : if (ascii_write) then
281 : 0 : write(unt_ascii, '(a)') "# Record 1: q-point (reduced coordinates)"
282 : 0 : write(unt_ascii, '(3es24.15)') qpt
283 : : end if
284 : :
285 : : ! Record 2: noncolin, nspin, lborn. BECS/quadrupoles not exported yet => lborn = .False.
286 : 0 : write(unt) ebands%nspinor == 2, nsppol, lborn
287 [ # # ]: 0 : if (ascii_write) then
288 : 0 : write(unt_ascii, '(a)') "# Record 2: noncolin, nspin, lborn"
289 : 0 : write(unt_ascii, *) ebands%nspinor == 2, nsppol, lborn
290 : : end if
291 : :
292 : : ! Record 3: number of electrons (including possible doping).
293 : 0 : write(unt) ebands%nelect
294 [ # # ]: 0 : if (ascii_write) then
295 : 0 : write(unt_ascii, '(a)') "# Record 3: number of electrons"
296 : 0 : write(unt_ascii, '(es24.15)') ebands%nelect
297 : : end if
298 : :
299 : : ! Record 4: first band, last band, total number of bands in the pw calculation.
300 : 0 : write(unt) gqk%bstart_k, gqk%bstop_k, ebands%mband
301 [ # # ]: 0 : if (ascii_write) then
302 : 0 : write(unt_ascii, '(a)') "# Record 4: nbnd_min, nbnd_max, nbnd_total"
303 : 0 : write(unt_ascii, '(3i8)') gqk%bstart_k, gqk%bstop_k, ebands%mband
304 : : end if
305 : :
306 : : ! Record 5: number of modes, number of k-points, number of atoms, number of atom types.
307 : 0 : write(unt) nmodes, gqk%glob_nk, natom, cryst%ntypat
308 [ # # ]: 0 : if (ascii_write) then
309 : 0 : write(unt_ascii, '(a)') "# Record 5: nmodes, nkpt, natom, ntypat"
310 : 0 : write(unt_ascii, '(4i8)') nmodes, gqk%glob_nk, natom, cryst%ntypat
311 : : end if
312 : :
313 : : ! Record 6: Bravais lattice index and cell dimensions (see ibrav/celldm comment above).
314 : 0 : write(unt) ibrav, (celldm(j), j=1,6)
315 [ # # ]: 0 : if (ascii_write) then
316 : 0 : write(unt_ascii, '(a)') "# Record 6: ibrav, celldm(1:6)"
317 : 0 : write(unt_ascii, '(i6,6es24.15)') ibrav, (celldm(j), j=1,6)
318 : : end if
319 : :
320 : : ! Record 7: atom symbols, atomic masses (atomic mass unit), atom types and
321 : : ! Cartesian atomic positions tau in units of alat.
322 [ # # ]: 0 : write(unt) (atm(j), j=1,cryst%ntypat), (cryst%amu(j), j=1,cryst%ntypat), &
323 [ # # # # : 0 : (cryst%typat(j), j=1,natom), ((tau_cart(j,i), j=1,3), i=1,natom)
# # ]
324 [ # # ]: 0 : if (ascii_write) then
325 : 0 : write(unt_ascii, '(a)') "# Record 7: atom symbols, masses (amu), types, tau (Cartesian, alat units)"
326 : 0 : write(unt_ascii, '(*(a3,1x))') (atm(j), j=1,cryst%ntypat)
327 [ # # # # ]: 0 : write(unt_ascii, '(*(es24.15,1x))') (cryst%amu(j), j=1,cryst%ntypat)
328 [ # # # # ]: 0 : write(unt_ascii, '(*(i6,1x))') (cryst%typat(j), j=1,natom)
329 [ # # ]: 0 : do i=1,natom
330 : 0 : write(unt_ascii, '(3es24.15)') (tau_cart(j,i), j=1,3)
331 : : end do
332 : : end if
333 : :
334 : : ! Record 8: squared phonon frequencies omega^2 in Ry^2.
335 : : ! my_wnuq is the signed frequency omega in Ha (negative for imaginary modes), so
336 : : ! omega^2 = sign(omega) * (omega*Ha2Ry)^2 = omega*|omega|*Ha2Ry^2, matching QE's w2.
337 [ # # # # ]: 0 : write(unt) (gqk%my_wnuq(nu, my_iq) * abs(gqk%my_wnuq(nu, my_iq)) * Ha2Ry ** 2, nu=1,nmodes)
338 [ # # ]: 0 : if (ascii_write) then
339 : 0 : write(unt_ascii, '(a)') "# Record 8: phonon frequencies squared omega^2 (Ry^2)"
340 : : write(unt_ascii, '(*(es24.15,1x))') &
341 [ # # # # ]: 0 : (gqk%my_wnuq(nu, my_iq) * abs(gqk%my_wnuq(nu, my_iq)) * Ha2Ry ** 2, nu=1,nmodes)
342 : : end if
343 : :
344 : : ! Phonon displacements for this q-point (composite index ipert = idir + 3*(iat-1),
345 : : ! Cartesian direction fast, atom slow, matching the QE mode/perturbation ordering):
346 : : ! dyn_qe = phonon displacements (eigenvectors divided by sqrt(mass)) = my_displ_cart.
347 : : ! These are written as 'dyn' (record 10) and used below to rotate g to the Cartesian basis.
348 : : !do iat=1,natom
349 : : !do idir=1,3
350 : : !mu = (iat-1)*3+idir
351 : : !do jat=1,natom
352 : : !do jdir=1,3
353 : : !nu = (jat-1)*3+jdir
354 : : !dyn_qe(mu, nu) = cmplx(gstore%ifc_dynmat(1,idir,iat,jdir,jat),&
355 : : !gstore%ifc_dynmat(2,idir,iat,jdir,jat),kind=dp)
356 : : !end do
357 : : !end do
358 : : !end do
359 : : !end do
360 [ # # ]: 0 : do nu=1,nmodes
361 [ # # ]: 0 : do iat=1,natom
362 [ # # ]: 0 : do idir=1,3
363 : 0 : ipert = idir + 3 * (iat - 1)
364 : : dyn_qe(ipert, nu) = cmplx(gqk%my_displ_cart(1, idir, iat, nu, my_iq), &
365 : 0 : gqk%my_displ_cart(2, idir, iat, nu, my_iq), kind=dp)
366 : : end do
367 : : end do
368 : : end do
369 : :
370 : : ! Record 9: eigenvectors in the QE basis (zz).
371 : : ! Like QE's ep_matrix_element_wannier, the e-ph matrix elements are written in the
372 : : ! Cartesian (atomic-displacement) basis (see the rotation in the k-loop below), so the
373 : : ! QE pattern matrix u reduces to the identity. We therefore write the identity matrix.
374 [ # # # # : 0 : write(unt) ((cmplx(merge(one, zero, i == j), zero, kind=dp), i=1,nmodes), j=1,nmodes)
# # # # #
# ]
375 [ # # ]: 0 : if (ascii_write) then
376 : 0 : write(unt_ascii, '(a)') "# Record 9: eigenvectors in the QE basis zz = identity (Cartesian g)"
377 [ # # ]: 0 : do j=1,nmodes
378 [ # # ]: 0 : do i=1,nmodes
379 [ # # ]: 0 : write(unt_ascii, '(2i6,2es24.15)') i, j, cmplx(merge(one, zero, i == j), zero, kind=dp)
380 : : end do
381 : : end do
382 : : end if
383 : :
384 : : ! Record 10: eigenvectors divided by masses, i.e. phonon displacements (dyn).
385 : 0 : write(unt) ((dyn_qe(i, j), i=1,nmodes), j=1,nmodes)
386 [ # # ]: 0 : if (ascii_write) then
387 : 0 : write(unt_ascii, '(a)') "# Record 10: displacements dyn(component, mode)"
388 [ # # ]: 0 : do j=1,nmodes
389 [ # # ]: 0 : do i=1,nmodes
390 : 0 : write(unt_ascii, '(2i6,2es24.15)') i, j, dyn_qe(i, j)
391 : : end do
392 : : end do
393 : : end if
394 : :
395 : : ! Records 11+: per k-point, the k-point (reduced coords), band energies (Ry) and,
396 : : ! for each Cartesian atomic perturbation, the e-ph matrix elements g (Ry).
397 [ # # ]: 0 : do my_ik=1,gqk%my_nk
398 [ # # ]: 0 : kk_bz = gqk%my_kpts(:, my_ik)
399 : 0 : ik_ibz = gqk%my_k2ibz(1, my_ik)
400 : 0 : kk_ibz = ebands%kptns(:,ik_ibz)
401 : 0 : ik_glob = my_ik + gqk%my_kstart - 1
402 : :
403 : 0 : write(unt) kk_bz
404 [ # # # # ]: 0 : write(unt) (ebands%eig(ib, ik_ibz, spin) * Ha2Ry, ib=gqk%bstart_k, gqk%bstop_k)
405 [ # # ]: 0 : if (ascii_write) then
406 : 0 : write(unt_ascii, '(a,i0,a)') "# k-point ", ik_glob, " (reduced coordinates)"
407 : 0 : write(unt_ascii, '(3es24.15)') kk_bz
408 : 0 : write(unt_ascii, '(a)') "# band energies (Ry)"
409 [ # # # # ]: 0 : write(unt_ascii, '(*(es24.15,1x))') (ebands%eig(ib, ik_ibz, spin) * Ha2Ry, ib=gqk%bstart_k, gqk%bstop_k)
410 : : end if
411 : :
412 : : ! Rotate g from the atom representation (reduced atomic-displacement directions, as
413 : : ! stored in gstore) to Cartesian directions, giving the bare Cartesian deformation
414 : : ! potential expected by EPIQ:
415 : : ! d_cart(beta,kappa) = sum_alpha gprimd(beta,alpha) * g_red(alpha,kappa)
416 : : ! with the composite index mu = idir + 3*(iat-1). This is a per-atom 3x3 map: it
417 : : ! involves no frequency or mass factors, so every mode is preserved (nothing is
418 : : ! zeroed, unlike the atom->phonon->Cartesian path). The Ha->Ry factor is applied
419 : : ! on output below, hence d_cart is in Ry/Bohr there.
420 [ # # # # : 0 : g_cart = (zero, zero)
# # ]
421 [ # # ]: 0 : do iat=1,natom
422 [ # # ]: 0 : do idir=1,3 ! Cartesian direction beta
423 : 0 : mu = idir + 3 * (iat - 1)
424 [ # # ]: 0 : do jdir=1,3 ! reduced direction alpha
425 : 0 : ipert = jdir + 3 * (iat - 1)
426 : : g_cart(mu, :, :) = g_cart(mu, :, :) &
427 [ # # # # ]: 0 : + cryst%gprimd(idir, jdir) * gqk%my_g(ipert, :, my_iq, :, my_ik)
428 : : end do
429 : : end do
430 : : end do
431 : :
432 [ # # ]: 0 : do mu=1,nmodes
433 : : ! d_matrix(m, n) = <k+q, m| dV/du^cart_mu |k, n>, with the k+q (bra) band index m
434 : : ! running fastest, matching QE's el_ph_mat(jbnd, ibnd) write order.
435 : : ! The 2nd dim (nb_kq) is the k+q band, the 3rd dim (nb_k) the k band.
436 [ # # # # : 0 : write(unt) ((g_cart(mu, j, i) * Ha2Ry, j=1,gqk%nb_kq), i=1,gqk%nb_k)
# # # # ]
437 [ # # ]: 0 : if (ascii_write) then
438 : : ! Columns: k index, Cartesian perturbation, band(k+q), band(k), Re(d), Im(d) in Ry.
439 [ # # ]: 0 : do i=1,gqk%nb_k
440 : 0 : band_k = i + gqk%bstart_k - 1
441 [ # # ]: 0 : do j=1,gqk%nb_kq
442 : 0 : band_kq = j + gqk%bstart_kq - 1
443 : : write(unt_ascii, '(4i6,2es30.15)') &
444 : 0 : ik_glob, mu, band_kq, band_k, g_cart(mu, j, i) * Ha2Ry
445 : : end do
446 : : end do
447 : : end if
448 : : end do
449 : : end do ! my_ik
450 : :
451 : : ! Trailing block: symmetry operations and star of q (QE elphsum_wannier layout).
452 : 0 : call write_qe_symmetry(cryst, qpt, celldm(1), unt, ascii_write, unt_ascii)
453 : :
454 : 0 : close(unt)
455 [ # # ]: 0 : if (ascii_write) close(unt_ascii)
456 : : end do ! my_iq
457 : :
458 [ # # # # ]: 0 : ABI_FREE(g_cart)
459 : : end associate
460 : : end do ! spin
461 : :
462 : : ! TODO:
463 : : ! Output BECS, dynamical quadrupoles, dynamical matrix, group velocities
464 : : ! See m_ifc
465 : : !ifc%zeff
466 : : !fc%qdrp_cart
467 : : !ifc%eta
468 : : !ifc%rpt
469 : : !ifc%wghatm
470 : : !ifc%dynmat
471 : : !ifc%short_atmfrc
472 : :
473 [ # # ]: 0 : ABI_FREE(dyn_qe)
474 [ # # ]: 0 : ABI_FREE(tau_cart)
475 : :
476 : 0 : call gstore%free()
477 : :
478 [ # # # # : 0 : end subroutine gstore_convert
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
479 : : !!***
480 : :
481 : : !----------------------------------------------------------------------
482 : :
483 : : !!****f* m_gstore_converters/write_epiq_input
484 : : !! NAME
485 : : !! write_epiq_input
486 : : !!
487 : : !! FUNCTION
488 : : !! Write the EPIQ input file with the &Diff_Start_Param namelist and the KPOINTS
489 : : !! section, reproducing QE's print_ph_input2epiq (called from ep_matrix_element_wannier).
490 : : !! The SCF/DFPT parameters are taken from the ABINIT dataset and band structure.
491 : : !!
492 : : !! INPUTS
493 : : !! ebands<ebands_t>=band structure (only used for HOMO/LUMO of insulators).
494 : : !! dtfil<datafiles_type>=filenames; dtfil%fildvdbin (DVDB) and dtfil%filddbsin (DDB).
495 : : !! fname=name of the output file.
496 : : !!
497 : : !! NOTES
498 : : !! ALL exported parameters come from the calculation that produced the DDB/DVDB
499 : : !! (the ground-state/DFPT run), NOT from the (denser) gstore/eph run:
500 : : !! - efermi, nel_r, occopt (-> ngauss_ph), tsmear (-> sigma_ph) from the DVDB header
501 : : !! (a standard ABINIT header carrying the GS scalars including the Fermi level).
502 : : !! - the KPOINTS list from the DDB header (the GS/DFPT k-mesh).
503 : : !! Exception: HOMO/LUMO of insulators are not stored in either header (no GS
504 : : !! eigenvalues), so they are taken from ebands (gap edges, essentially mesh-independent).
505 : : !!
506 : : !! SOURCE
507 : :
508 : 0 : subroutine write_epiq_input(ebands, dtfil, fname)
509 : :
510 : : !Arguments ------------------------------------
511 : : class(ebands_t),intent(in) :: ebands
512 : : type(datafiles_type),intent(in) :: dtfil
513 : : character(len=*),intent(in) :: fname
514 : :
515 : : !Local variables-------------------------------
516 : : !scalars
517 : : integer :: unt, iunt, ik, ngauss, gap_err, fform
518 : : real(dp) :: homo, lumo, knorm
519 : : logical :: is_metal
520 : : character(len=500) :: msg
521 : : character(len=24) :: smear_label
522 : 0 : type(gaps_t) :: gaps
523 [ # # # # : 0 : type(ddb_hdr_type) :: ddb_hdr
# # # # #
# # # # #
# # # # ]
524 : 0 : type(hdr_type) :: dfpt_hdr
525 : : !----------------------------------------------------------------------
526 : :
527 [ # # ]: 0 : if (open_file(fname, msg, newunit=unt, form="formatted", status="unknown", action="write") /= 0) then
528 : 0 : ABI_ERROR(msg)
529 : : end if
530 : :
531 : : ! Read the GS/DFPT header from the DVDB: it carries the Fermi level, nelect, occopt
532 : : ! and tsmear of the run that produced the DDB/DVDB (the DDB header has no Fermi level).
533 : : ! The DVDB starts with two records (version, numv1) before the standard ABINIT header,
534 : : ! so we skip them and read the header in place (fort_read without rewind).
535 [ # # ]: 0 : if (open_file(dtfil%fildvdbin, msg, newunit=iunt, form="unformatted", status="old", action="read") /= 0) then
536 : 0 : ABI_ERROR(msg)
537 : : end if
538 : 0 : read(iunt) ! skip the DVDB version record
539 : 0 : read(iunt) ! skip the numv1 record
540 : 0 : call dfpt_hdr%fort_read(iunt, fform)
541 : 0 : close(iunt)
542 : :
543 : : ! occopt >= 3 => metallic occupation with smearing; otherwise fixed occupations (insulator).
544 : 0 : is_metal = dfpt_hdr%occopt >= 3
545 : :
546 : 0 : write(unt, '(a)') "! parameter of the SCF DFPT calculation useful for EPIq"
547 : 0 : write(unt, '(a)') "&Diff_Start_Param"
548 : 0 : write(unt, '(3x,a,f12.6,a)') "efermi=", dfpt_hdr%fermie * Ha_eV, ", ! in (eV)"
549 : 0 : write(unt, '(3x,a,f12.6,a)') "nel_r=", dfpt_hdr%nelect, ","
550 : :
551 [ # # ]: 0 : if (.not. is_metal) then
552 : : ! Insulator: report the HOMO and LUMO levels (in eV). Not in the DFPT header
553 : : ! (no GS eigenvalues), so taken from ebands (gap edges are mesh-independent).
554 : 0 : gaps = ebands%get_gaps(gap_err)
555 [ # # ]: 0 : if (gap_err == 0) then
556 : 0 : homo = gaps%vb_max(1); lumo = gaps%cb_min(1)
557 : : else
558 : : ! Could not determine a gap (semimetal?): fall back to the Fermi level.
559 : : homo = dfpt_hdr%fermie; lumo = dfpt_hdr%fermie
560 : : end if
561 : 0 : call gaps%free()
562 : 0 : write(unt, '(3x,a,f12.6,a)') "homo=", homo * Ha_eV, ", ! in (eV)"
563 : 0 : write(unt, '(3x,a,f12.6,a)') "lumo=", lumo * Ha_eV, ", ! in (eV)"
564 : : else
565 : : ! Metal: report the smearing width (Rydberg) and the QE ngauss code.
566 : : ! Map ABINIT occopt onto QE ngauss (see Modules input conventions):
567 : : ! 3 -> -99 (Fermi-Dirac), 4/5 -> -1 (cold/Marzari), 6 -> 1 (Methfessel-Paxton), 7 -> 0 (Gaussian)
568 : 0 : select case (dfpt_hdr%occopt)
569 : 0 : case (3); ngauss = -99; smear_label = "fd"
570 : 0 : case (4, 5); ngauss = -1; smear_label = "cold"
571 : 0 : case (6); ngauss = 1; smear_label = "mp"
572 : 0 : case (7); ngauss = 0; smear_label = "gauss"
573 [ # # # # : 0 : case default; ngauss = -66; smear_label = "unknown"
# ]
574 : : end select
575 : 0 : write(unt, '(3x,a,f12.6,a)') "sigma_ph=", dfpt_hdr%tsmear * two, ", ! in (Rydberg)"
576 : 0 : write(unt, '(3x,a,i3,a)') "ngauss_ph=", ngauss, ", ! "//trim(smear_label)
577 : : end if
578 : 0 : write(unt, '(a)') "/"
579 : :
580 : 0 : call dfpt_hdr%free()
581 : :
582 : : ! KPOINTS section. Use the k-mesh that produced the DDB/DVDB (the ground-state/DFPT
583 : : ! mesh), read from the DDB header, NOT the gstore/eph k-mesh in dtset/ebands.
584 : : ! The DDB header stores the explicit k-point list (no kptrlatt), so we dump it.
585 : 0 : call ddb_hdr%open_read(dtfil%filddbsin, xmpi_comm_self)
586 : 0 : call ddb_hdr%close() ! we only need the header data (k-points)
587 : :
588 [ # # ]: 0 : knorm = ddb_hdr%kptnrm; if (abs(knorm) < tol12) knorm = one
589 : :
590 : 0 : write(unt, '(/,a)') "KPOINTS"
591 : 0 : write(unt, '(a)') "crystal"
592 : 0 : write(unt, '(6x,i9)') ddb_hdr%nkpt
593 [ # # ]: 0 : do ik=1,ddb_hdr%nkpt
594 [ # # ]: 0 : write(unt, '(3x,4(es20.10,2x))') ddb_hdr%kpt(:,ik) / knorm, ddb_hdr%wtk(ik)
595 : : end do
596 : :
597 : 0 : call ddb_hdr%free()
598 : :
599 : 0 : close(unt)
600 : :
601 [ # # # # : 0 : end subroutine write_epiq_input
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
602 : : !!***
603 : :
604 : : !----------------------------------------------------------------------
605 : :
606 : : !!****f* m_gstore_converters/write_dynq0
607 : : !! NAME
608 : : !! write_dynq0
609 : : !!
610 : : !! FUNCTION
611 : : !! Write the EPIQ "dynq0" file: the q-mesh, the number of irreducible q-points
612 : : !! contained in the elph.mat files, and the list of those q-points in Cartesian
613 : : !! coordinates (units of 2pi/alat, i.e. QE tpiba). Analog of QE's dynq0 output.
614 : : !!
615 : : !! INPUTS
616 : : !! ngqpt(3)=dimensions of the q-mesh.
617 : : !! qibz(3,nqibz)=irreducible q-points in reduced (crystal) coordinates.
618 : : !! gprimd(3,3)=reciprocal lattice vectors (Bohr^-1), columns G_i/2pi (ABINIT convention).
619 : : !! alat=lattice parameter in Bohr (celldm(1)).
620 : : !! fname=name of the output file.
621 : : !!
622 : : !! SOURCE
623 : :
624 [ # # ]: 0 : subroutine write_dynq0(ngqpt, qibz, gprimd, alat, fname)
625 : :
626 : : !Arguments ------------------------------------
627 : : integer,intent(in) :: ngqpt(3)
628 : : real(dp),intent(in) :: qibz(:,:), gprimd(3,3), alat
629 : : character(len=*),intent(in) :: fname
630 : :
631 : : !Local variables-------------------------------
632 : : integer :: unt, iq, nqibz
633 : : real(dp) :: qcart(3)
634 : : character(len=500) :: msg
635 : : !----------------------------------------------------------------------
636 : :
637 : 0 : nqibz = size(qibz, 2)
638 : :
639 [ # # ]: 0 : if (open_file(fname, msg, newunit=unt, form="formatted", status="unknown", action="write") /= 0) then
640 : 0 : ABI_ERROR(msg)
641 : : end if
642 : :
643 : : ! Line 1: q-mesh. Line 2: number of irreducible q-points.
644 : 0 : write(unt, '(3i4)') ngqpt(1), ngqpt(2), ngqpt(3)
645 : 0 : write(unt, '(i4)') nqibz
646 : :
647 : : ! One line per irreducible q-point in Cartesian coordinates (2pi/alat units):
648 : : ! q_cart[tpiba] = alat * matmul(gprimd, q_red).
649 [ # # ]: 0 : do iq=1,nqibz
650 [ # # # # : 0 : qcart = alat * matmul(gprimd, qibz(:,iq))
# # ]
651 : 0 : write(unt, '(3e24.15)') qcart(1), qcart(2), qcart(3)
652 : : end do
653 : :
654 : 0 : close(unt)
655 : :
656 : 0 : end subroutine write_dynq0
657 : : !!***
658 : :
659 : : !----------------------------------------------------------------------
660 : :
661 : : !!****f* m_gstore_converters/write_dynq
662 : : !! NAME
663 : : !! write_dynq
664 : : !!
665 : : !! FUNCTION
666 : : !! Write a Quantum ESPRESSO dynamical-matrix file (dynq<iq>) for an irreducible q-point
667 : : !! AND all the q-points of its star, reproducing the layout of QE's write_dyn_on_file +
668 : : !! rotate_dvscf_star + dyndiag: a single header (cell, atoms), one "Dynamical Matrix in
669 : : !! cartesian axes" block per star member (representative q first), and a single
670 : : !! diagonalization block (frequencies and eigenvectors) for the representative q.
671 : : !!
672 : : !! The star is generated from the ABINIT crystal symmetries (q' = symrec*q, deduplicated
673 : : !! modulo a reciprocal-lattice vector), and the dynamical matrix at each star member is
674 : : !! evaluated directly with ifc%fourq (equivalent to rotating D(q) by symmetry, since
675 : : !! D(Sq) = sum_R Phi(R) exp(i Sq.R); the gauge is irrelevant as D is gauge-invariant).
676 : : !! The Cartesian dynamical matrix is rebuilt as
677 : : !! phi(ka,k'b) = sqrt(M_k M_k') * sum_nu z(ka,nu) * w2(nu) * conjg(z(k'b,nu))
678 : : !! with z the orthonormal eigenvectors, w2 = signed omega^2 in Ry^2 and M the QE
679 : : !! Rydberg atomic masses (amu * amu_emass/2). This matches QE's convention.
680 : : !!
681 : : !! INPUTS
682 : : !! cryst<crystal_t>=crystal structure.
683 : : !! ifc<ifc_type>=interatomic force constants (for Fourier interpolation at q).
684 : : !! qpt_red(3)=q-point in reduced (crystal) coordinates.
685 : : !! alat=lattice parameter in Bohr (celldm(1)).
686 : : !! atm(ntypat)=atomic symbols.
687 : : !! fname=name of the output file.
688 : : !!
689 : : !! SOURCE
690 : :
691 : 0 : subroutine write_dynq(cryst, ifc, qpt_red, alat, atm, fname)
692 : :
693 : : !Arguments ------------------------------------
694 : : class(crystal_t),intent(in) :: cryst
695 : : class(ifc_type),intent(in) :: ifc
696 : : real(dp),intent(in) :: qpt_red(3), alat
697 : : character(len=3),intent(in) :: atm(cryst%ntypat)
698 : : character(len=*),intent(in) :: fname
699 : :
700 : : !Local variables-------------------------------
701 : : !scalars
702 : : integer :: natom, nmodes, ntypat, unt, na, nb, icar, jcar, nu, it, isym, iqs, nq_star
703 : : real(dp),parameter :: Ha2Ry = two, accep = 1.0e-5_dp
704 : : real(dp) :: znorm, freq_cm, freq_thz
705 : : logical :: found
706 : : complex(dp) :: zi, zj, cs
707 : : character(len=500) :: msg
708 : : !arrays
709 : : real(dp) :: at(3,3), bg(3,3), celldm(6), aq(3), raq(3), dq(3)
710 : 0 : real(dp),allocatable :: phfrq(:), displ_cart(:,:,:,:), eigvec(:,:,:,:), w2(:), amass_qe(:), tau(:,:)
711 : 0 : real(dp),allocatable :: saq(:,:), sxq(:,:), phfrq_rep(:), eigvec_rep(:,:,:,:)
712 : 0 : complex(dp),allocatable :: phi(:,:,:,:)
713 : : !----------------------------------------------------------------------
714 : :
715 : 0 : natom = cryst%natom; nmodes = 3 * natom; ntypat = cryst%ntypat
716 : :
717 [ # # # # : 0 : ABI_MALLOC(phfrq, (nmodes))
# # # # #
# # # ]
718 [ # # # # : 0 : ABI_MALLOC(displ_cart, (2, 3, natom, nmodes))
# # # # #
# # # # #
# # ]
719 [ # # # # : 0 : ABI_MALLOC(eigvec, (2, 3, natom, nmodes))
# # # # #
# # # ]
720 [ # # # # : 0 : ABI_MALLOC(w2, (nmodes))
# # # # ]
721 [ # # # # : 0 : ABI_MALLOC(amass_qe, (ntypat))
# # # # #
# # # ]
722 [ # # # # : 0 : ABI_MALLOC(tau, (3, natom))
# # # # #
# # # ]
723 [ # # # # : 0 : ABI_MALLOC(phi, (3, 3, natom, natom))
# # # # #
# # # # #
# # ]
724 [ # # # # : 0 : ABI_MALLOC(saq, (3, cryst%nsym))
# # # # #
# # # ]
725 [ # # # # : 0 : ABI_MALLOC(sxq, (3, cryst%nsym))
# # # # ]
726 [ # # # # : 0 : ABI_MALLOC(phfrq_rep, (nmodes))
# # # # ]
727 [ # # # # : 0 : ABI_MALLOC(eigvec_rep, (2, 3, natom, nmodes))
# # # # #
# # # ]
728 : :
729 : : ! QE masses in Rydberg atomic units (amu * amu_ry, amu_ry = amu_emass/2).
730 [ # # ]: 0 : do it=1,ntypat
731 : 0 : amass_qe(it) = cryst%amu(it) * amu_emass * half
732 : : end do
733 : :
734 : : ! Cartesian cell quantities (alat / tpiba units).
735 [ # # ]: 0 : do nu=1,3
736 [ # # ]: 0 : at(:,nu) = cryst%rprimd(:,nu) / alat
737 [ # # ]: 0 : bg(:,nu) = cryst%gprimd(:,nu) * alat
738 : : end do
739 [ # # ]: 0 : do na=1,natom
740 [ # # # # : 0 : tau(:,na) = matmul(cryst%rprimd, cryst%xred(:,na)) / alat
# # # # #
# ]
741 : : end do
742 : 0 : celldm = zero; celldm(1) = alat
743 : :
744 : : ! Build the star of qpt_red from the crystal symmetries (q' = symrec*q, deduplicated
745 : : ! modulo a reciprocal-lattice vector). The representative q is stored first (member 1).
746 : 0 : aq(:) = qpt_red(:)
747 [ # # ]: 0 : nq_star = 1; saq(:,1) = aq(:)
748 [ # # ]: 0 : do isym=1,cryst%nsym
749 [ # # # # : 0 : raq = matmul(real(cryst%symrec(:,:,isym), dp), aq)
# # # # #
# ]
750 : 0 : found = .False.
751 [ # # ]: 0 : do iqs=1,nq_star
752 [ # # ]: 0 : dq = raq - saq(:,iqs)
753 [ # # # # : 0 : if (all(abs(dq - nint(dq)) < accep)) then
# # ]
754 : : found = .True.; exit
755 : : end if
756 : : end do
757 [ # # ]: 0 : if (.not. found) then
758 [ # # ]: 0 : nq_star = nq_star + 1; saq(:,nq_star) = raq(:)
759 : : end if
760 : : end do
761 [ # # ]: 0 : do iqs=1,nq_star
762 [ # # # # : 0 : sxq(:,iqs) = matmul(bg, saq(:,iqs)) ! Cartesian (2pi/alat) coordinates
# # ]
763 : : end do
764 : :
765 [ # # ]: 0 : if (open_file(fname, msg, newunit=unt, form="formatted", status="unknown", action="write") /= 0) then
766 : 0 : ABI_ERROR(msg)
767 : : end if
768 : :
769 : : ! ---- Header (written once) ----
770 : 0 : write(unt, '(a)') "Dynamical matrix file"
771 : 0 : write(unt, '(a)') "Converted from ABINIT GSTORE"
772 : : ! ntyp, nat, ibrav=0 (free lattice) followed by celldm; with ibrav=0 the basis vectors follow.
773 : 0 : write(unt, '(i3,i5,i4,6f11.7)') ntypat, natom, 0, (celldm(it), it=1,6)
774 : 0 : write(unt, '(a)') "Basis vectors"
775 [ # # ]: 0 : do nu=1,3
776 : 0 : write(unt, '(2x,3f15.9)') at(1,nu), at(2,nu), at(3,nu)
777 : : end do
778 [ # # ]: 0 : do it=1,ntypat
779 : 0 : write(unt, *) it, " '"//atm(it)//"' ", amass_qe(it)
780 : : end do
781 [ # # ]: 0 : do na=1,natom
782 : 0 : write(unt, '(2i5,3f18.10)') na, cryst%typat(na), tau(1,na), tau(2,na), tau(3,na)
783 : : end do
784 : :
785 : : ! ---- One dynamical-matrix block per q-point of the star ----
786 [ # # ]: 0 : do iqs=1,nq_star
787 : : ! Frequencies, displacements and orthonormal eigenvectors at this star member.
788 : 0 : call ifc%fourq(cryst, saq(:,iqs), phfrq, displ_cart, out_eigvec=eigvec)
789 [ # # ]: 0 : if (iqs == 1) then
790 [ # # # # : 0 : phfrq_rep = phfrq; eigvec_rep = eigvec ! keep the representative for the diag block
# # # # #
# # # #
# ]
791 : : end if
792 : :
793 : : ! Signed squared phonon frequencies in Ry^2.
794 [ # # ]: 0 : do nu=1,nmodes
795 : 0 : w2(nu) = phfrq(nu) * abs(phfrq(nu)) * Ha2Ry ** 2
796 : : end do
797 : :
798 : : ! Dynamical matrix in Cartesian axes (QE convention, see header).
799 [ # # ]: 0 : do nb=1,natom
800 [ # # ]: 0 : do na=1,natom
801 [ # # ]: 0 : do jcar=1,3
802 [ # # ]: 0 : do icar=1,3
803 : 0 : cs = czero
804 [ # # ]: 0 : do nu=1,nmodes
805 : 0 : zi = cmplx(eigvec(1,icar,na,nu), eigvec(2,icar,na,nu), kind=dp)
806 : 0 : zj = cmplx(eigvec(1,jcar,nb,nu), eigvec(2,jcar,nb,nu), kind=dp)
807 : 0 : cs = cs + zi * w2(nu) * conjg(zj)
808 : : end do
809 : 0 : phi(icar,jcar,na,nb) = sqrt(amass_qe(cryst%typat(na)) * amass_qe(cryst%typat(nb))) * cs
810 : : end do
811 : : end do
812 : : end do
813 : : end do
814 : :
815 : 0 : write(unt, '(/,5x,a)') "Dynamical Matrix in cartesian axes"
816 : 0 : write(unt, '(/,5x,a,3f14.9,a,/)') "q = ( ", sxq(1,iqs), sxq(2,iqs), sxq(3,iqs), " ) "
817 [ # # ]: 0 : do na=1,natom
818 [ # # ]: 0 : do nb=1,natom
819 : 0 : write(unt, '(2i5)') na, nb
820 [ # # ]: 0 : do icar=1,3
821 : 0 : write(unt, '(3(2f12.8,2x))') (phi(icar,jcar,na,nb), jcar=1,3)
822 : : end do
823 : : end do
824 : : end do
825 : : end do ! iqs
826 : :
827 : : ! ---- Diagonalization block (once, for the representative q = member 1) ----
828 : 0 : write(unt, '(/,5x,a)') "Diagonalizing the dynamical matrix"
829 : 0 : write(unt, '(/,5x,a,3f14.9,a,/)') "q = ( ", sxq(1,1), sxq(2,1), sxq(3,1), " ) "
830 : 0 : write(unt, '(1x,74("*"))')
831 [ # # ]: 0 : do nu=1,nmodes
832 : 0 : freq_cm = phfrq_rep(nu) * Ha_cmm1
833 : 0 : freq_thz = phfrq_rep(nu) * Ha_THz
834 : 0 : write(unt, '(5x,a,i5,a,f15.6,a,f15.6,a)') "freq (", nu, ") = ", freq_thz, " [THz] = ", freq_cm, " [cm-1]"
835 : 0 : znorm = zero
836 [ # # ]: 0 : do na=1,natom
837 [ # # ]: 0 : do icar=1,3
838 : 0 : znorm = znorm + eigvec_rep(1,icar,na,nu)**2 + eigvec_rep(2,icar,na,nu)**2
839 : : end do
840 : : end do
841 [ # # ]: 0 : znorm = sqrt(znorm); if (znorm < tol12) znorm = one
842 [ # # ]: 0 : do na=1,natom
843 : 0 : write(unt, '(1x,a,3(f10.6,1x,f10.6,3x),a)') "( ", &
844 [ # # # # ]: 0 : (eigvec_rep(1,icar,na,nu)/znorm, eigvec_rep(2,icar,na,nu)/znorm, icar=1,3), ")"
845 : : end do
846 : : end do
847 : 0 : write(unt, '(1x,74("*"))')
848 : :
849 : 0 : close(unt)
850 : :
851 [ # # ]: 0 : ABI_FREE(phfrq)
852 [ # # ]: 0 : ABI_FREE(displ_cart)
853 [ # # ]: 0 : ABI_FREE(eigvec)
854 [ # # ]: 0 : ABI_FREE(w2)
855 [ # # ]: 0 : ABI_FREE(amass_qe)
856 [ # # ]: 0 : ABI_FREE(tau)
857 [ # # ]: 0 : ABI_FREE(phi)
858 [ # # ]: 0 : ABI_FREE(saq)
859 [ # # ]: 0 : ABI_FREE(sxq)
860 [ # # ]: 0 : ABI_FREE(phfrq_rep)
861 [ # # ]: 0 : ABI_FREE(eigvec_rep)
862 : :
863 : 0 : end subroutine write_dynq
864 : : !!***
865 : :
866 : : !----------------------------------------------------------------------
867 : :
868 : : !!****f* m_gstore_converters/write_qe_symmetry
869 : : !! NAME
870 : : !! write_qe_symmetry
871 : : !!
872 : : !! FUNCTION
873 : : !! Append the symmetry-operations block and the star of q to an already open
874 : : !! elph.mat.q_* file, reproducing the layout written by Quantum ESPRESSO's
875 : : !! elphsum_wannier (PHonon/PH/elphon.f90).
876 : : !!
877 : : !! The QE quantities are reconstructed from the ABINIT crystal object using the
878 : : !! following correspondence (see PW/src/symm_base.f90, PHonon/PH/obsolete.f90 and
879 : : !! LR_Modules/star_q.f90 in the bundled QE sources):
880 : : !!
881 : : !! QE s(:,:,isym) = cryst%symrec(:,:,isym) (rotations on reduced reciprocal coords;
882 : : !! s^T = symrel^-1 acts on reduced positions)
883 : : !! QE irt(isym,na) = cryst%indsym(4,isym,na) (atom na -> atom irt under symrel^-1)
884 : : !! QE at(:,j) = cryst%rprimd(:,j) / alat (direct lattice, alat units)
885 : : !! QE bg(:,j) = cryst%gprimd(:,j) * alat (reciprocal lattice, 2pi/alat units)
886 : : !! QE tau / xau = cryst%xred (reduced atomic positions)
887 : : !!
888 : : !! invs, rtau and the star (nq, sxq, isq, imq) are then obtained by porting QE's
889 : : !! inverse_s, sgam_ph and star_q. The symmetry operations are written in the native
890 : : !! ABINIT order (1..nsym); only the order differs from QE (which sorts the small
891 : : !! group of q first), the set and the per-operation data are equivalent.
892 : : !!
893 : : !! INPUTS
894 : : !! cryst<crystal_t>=crystal structure (with symmetries).
895 : : !! qpt(3)=q-point in reduced (crystal) coordinates.
896 : : !! alat=lattice parameter in Bohr (celldm(1)).
897 : : !! unt=Fortran unit of the (open) unformatted file.
898 : : !! ascii_write=if .True. also mirror the records in the formatted file unt_ascii.
899 : : !! unt_ascii=Fortran unit of the (open) formatted file (used only if ascii_write).
900 : : !!
901 : : !! SOURCE
902 : :
903 : 0 : subroutine write_qe_symmetry(cryst, qpt, alat, unt, ascii_write, unt_ascii)
904 : :
905 : : !Arguments ------------------------------------
906 : : !scalars
907 : : class(crystal_t),intent(in) :: cryst
908 : : integer,intent(in) :: unt, unt_ascii
909 : : logical,intent(in) :: ascii_write
910 : : real(dp),intent(in) :: alat
911 : : !arrays
912 : : real(dp),intent(in) :: qpt(3)
913 : :
914 : : !Local variables-------------------------------
915 : : !scalars
916 : : integer :: nsym, natom, isym, jsym, ism1, i, j, k, na, nb, nq, imq, iq
917 : : logical :: found
918 : : real(dp),parameter :: accep = 1.0e-5_dp
919 : : !arrays
920 : : integer :: ss(3,3)
921 : : integer,parameter :: identity(3,3) = reshape([1,0,0, 0,1,0, 0,0,1], [3,3])
922 : 0 : integer,allocatable :: s(:,:,:), invs(:), irt(:,:), isq(:), nsq(:)
923 : : real(dp) :: at(3,3), bg(3,3), aq(3), raq(3), ft(3), dq(3)
924 : 0 : real(dp),allocatable :: rtau(:,:,:), sxq(:,:), saq(:,:)
925 : : !----------------------------------------------------------------------
926 : :
927 : 0 : nsym = cryst%nsym; natom = cryst%natom
928 : :
929 [ # # # # : 0 : ABI_MALLOC(s, (3, 3, nsym))
# # # # #
# # # ]
930 [ # # # # : 0 : ABI_MALLOC(invs, (nsym))
# # # # #
# # # ]
931 [ # # # # : 0 : ABI_MALLOC(irt, (nsym, natom))
# # # # #
# # # # #
# # ]
932 [ # # # # : 0 : ABI_MALLOC(isq, (nsym))
# # # # ]
933 [ # # # # : 0 : ABI_MALLOC(nsq, (nsym))
# # # # ]
934 [ # # # # : 0 : ABI_MALLOC(rtau, (3, nsym, natom))
# # # # #
# # # # #
# # ]
935 [ # # # # : 0 : ABI_MALLOC(sxq, (3, nsym))
# # # # #
# # # ]
936 [ # # # # : 0 : ABI_MALLOC(saq, (3, nsym))
# # # # ]
937 : :
938 : : ! Lattice vectors in QE units (at . bg^T = identity).
939 [ # # ]: 0 : do j=1,3
940 [ # # ]: 0 : at(:,j) = cryst%rprimd(:,j) / alat
941 [ # # ]: 0 : bg(:,j) = cryst%gprimd(:,j) * alat
942 : : end do
943 : :
944 : : ! Rotations in crystal axis and atom mapping.
945 [ # # ]: 0 : do isym=1,nsym
946 [ # # # # ]: 0 : s(:,:,isym) = cryst%symrec(:,:,isym)
947 [ # # ]: 0 : do na=1,natom
948 : 0 : irt(isym, na) = cryst%indsym(4, isym, na)
949 : : end do
950 : : end do
951 : :
952 : : ! invs(isym): index of the inverse operation (ported from QE inverse_s).
953 [ # # ]: 0 : do isym=1,nsym
954 : 0 : found = .False.
955 [ # # ]: 0 : do jsym=1,nsym
956 [ # # # # : 0 : ss = matmul(s(:,:,jsym), s(:,:,isym))
# # ]
957 [ # # # # : 0 : if (all(ss == identity)) then
# # # # ]
958 : 0 : invs(isym) = jsym; found = .True.; exit
959 : : end if
960 : : end do
961 : 0 : ABI_CHECK(found, "write_qe_symmetry: symmetry operations do not form a group.")
962 : : end do
963 : :
964 : : ! rtau(:,isym,na) = S.tau_na - tau_nb in Cartesian coords (alat units), with nb = irt.
965 : : ! Ported from QE sgam_ph using xau = reduced atomic coordinates = cryst%xred.
966 [ # # # # : 0 : rtau = zero
# # ]
967 [ # # ]: 0 : do isym=1,nsym
968 [ # # ]: 0 : do na=1,natom
969 : 0 : nb = irt(isym, na)
970 [ # # ]: 0 : do i=1,3
971 : : ft(i) = s(1,i,isym) * cryst%xred(1,na) + s(2,i,isym) * cryst%xred(2,na) &
972 : 0 : + s(3,i,isym) * cryst%xred(3,na) - cryst%xred(i,nb)
973 : : end do
974 [ # # ]: 0 : do i=1,3
975 : 0 : rtau(i, isym, na) = at(i,1) * ft(1) + at(i,2) * ft(2) + at(i,3) * ft(3)
976 : : end do
977 : : end do
978 : : end do
979 : :
980 : : ! Star of q (ported from QE star_q). aq is q in reduced (crystal) coordinates.
981 : 0 : aq(:) = qpt(:)
982 [ # # # # : 0 : nsq(:) = 0; isq(:) = 0; saq(:,:) = zero; sxq(:,:) = zero; nq = 0
# # # # #
# # # ]
983 [ # # ]: 0 : do isym=1,nsym
984 : 0 : ism1 = invs(isym)
985 [ # # ]: 0 : do i=1,3
986 : 0 : raq(i) = s(i,1,ism1) * aq(1) + s(i,2,ism1) * aq(2) + s(i,3,ism1) * aq(3)
987 : : end do
988 [ # # ]: 0 : do iq=1,nq
989 [ # # ]: 0 : dq(:) = raq(:) - saq(:,iq)
990 [ # # # # : 0 : if (all(abs(dq - nint(dq)) < accep)) then
# # ]
991 : 0 : isq(isym) = iq; nsq(iq) = nsq(iq) + 1
992 : : end if
993 : : end do
994 [ # # ]: 0 : if (isq(isym) == 0) then
995 [ # # ]: 0 : nq = nq + 1; nsq(nq) = 1; isq(isym) = nq; saq(:,nq) = raq(:)
996 [ # # ]: 0 : do i=1,3
997 : 0 : sxq(i,nq) = bg(i,1) * saq(1,nq) + bg(i,2) * saq(2,nq) + bg(i,3) * saq(3,nq)
998 : : end do
999 : : end if
1000 : : end do
1001 : :
1002 : : ! imq: index of -q in the star (0 if absent).
1003 : 0 : imq = 0
1004 [ # # ]: 0 : do iq=1,nq
1005 [ # # ]: 0 : dq(:) = -aq(:) - saq(:,iq)
1006 [ # # # # : 0 : if (all(abs(dq - nint(dq)) < accep)) imq = iq
# # ]
1007 : : end do
1008 : :
1009 : : ! Sanity check on the star degeneracy (as in QE star_q).
1010 [ # # ]: 0 : do iq=1,nq
1011 [ # # ]: 0 : if (nsq(iq) * nq /= nsym) then
1012 : 0 : ABI_WARNING(sjoin("write_qe_symmetry: unexpected star-of-q degeneracy for iq=", itoa(iq)))
1013 : : end if
1014 : : end do
1015 : :
1016 : : ! ---- Write the block in the QE elphsum_wannier order ----
1017 [ # # ]: 0 : do j=1,3
1018 : 0 : write(unt) (at(i,j), i=1,3)
1019 : : end do
1020 [ # # ]: 0 : do j=1,3
1021 : 0 : write(unt) (bg(i,j), i=1,3)
1022 : : end do
1023 : 0 : write(unt) nsym, nq, imq
1024 [ # # ]: 0 : do i=1,nsym
1025 : 0 : write(unt) i, invs(i), isq(i)
1026 [ # # ]: 0 : do j=1,3
1027 [ # # ]: 0 : do k=1,3
1028 : 0 : write(unt) k, j, s(k,j,i)
1029 : : end do
1030 : : end do
1031 [ # # ]: 0 : do j=1,natom
1032 : 0 : write(unt) j, irt(i,j)
1033 : : end do
1034 [ # # ]: 0 : do j=1,3
1035 [ # # ]: 0 : do k=1,natom
1036 : 0 : write(unt) j, i, rtau(j,i,k)
1037 : : end do
1038 : : end do
1039 [ # # ]: 0 : do j=1,3
1040 : 0 : write(unt) j, sxq(j,i)
1041 : : end do
1042 : : end do
1043 : :
1044 : : ! ---- Optional human-readable mirror ----
1045 [ # # ]: 0 : if (ascii_write) then
1046 : 0 : write(unt_ascii, '(a)') "# Symmetry: direct lattice vectors at(:,j) (alat units)"
1047 [ # # ]: 0 : do j=1,3
1048 : 0 : write(unt_ascii, '(3es24.15)') (at(i,j), i=1,3)
1049 : : end do
1050 : 0 : write(unt_ascii, '(a)') "# Symmetry: reciprocal lattice vectors bg(:,j) (2pi/alat units)"
1051 [ # # ]: 0 : do j=1,3
1052 : 0 : write(unt_ascii, '(3es24.15)') (bg(i,j), i=1,3)
1053 : : end do
1054 : 0 : write(unt_ascii, '(a)') "# Symmetry: nsym, nq, imq"
1055 : 0 : write(unt_ascii, '(3i8)') nsym, nq, imq
1056 [ # # ]: 0 : do i=1,nsym
1057 : 0 : write(unt_ascii, '(a,3i6)') "# isym, invs, isq: ", i, invs(i), isq(i)
1058 : 0 : write(unt_ascii, '(a)') "# rotation s(row k, col j)"
1059 [ # # ]: 0 : do j=1,3
1060 [ # # ]: 0 : do k=1,3
1061 : 0 : write(unt_ascii, '(3i6)') k, j, s(k,j,i)
1062 : : end do
1063 : : end do
1064 : 0 : write(unt_ascii, '(a)') "# irt(atom)"
1065 [ # # ]: 0 : do j=1,natom
1066 : 0 : write(unt_ascii, '(2i6)') j, irt(i,j)
1067 : : end do
1068 : 0 : write(unt_ascii, '(a)') "# rtau(coord, isym, atom) (alat units)"
1069 [ # # ]: 0 : do j=1,3
1070 [ # # ]: 0 : do k=1,natom
1071 : 0 : write(unt_ascii, '(2i6,es24.15)') j, i, rtau(j,i,k)
1072 : : end do
1073 : : end do
1074 : 0 : write(unt_ascii, '(a)') "# sxq(coord) (2pi/alat units)"
1075 [ # # ]: 0 : do j=1,3
1076 : 0 : write(unt_ascii, '(i6,es24.15)') j, sxq(j,i)
1077 : : end do
1078 : : end do
1079 : : end if
1080 : :
1081 [ # # ]: 0 : ABI_FREE(s)
1082 [ # # ]: 0 : ABI_FREE(invs)
1083 [ # # ]: 0 : ABI_FREE(irt)
1084 [ # # ]: 0 : ABI_FREE(isq)
1085 [ # # ]: 0 : ABI_FREE(nsq)
1086 [ # # ]: 0 : ABI_FREE(rtau)
1087 [ # # ]: 0 : ABI_FREE(sxq)
1088 [ # # ]: 0 : ABI_FREE(saq)
1089 : :
1090 : 0 : end subroutine write_qe_symmetry
1091 : : !!***
1092 : :
1093 : : !----------------------------------------------------------------------
1094 : :
1095 : : ! Helper function copied from epic/src/io_matelem.F90
1096 : 0 : subroutine define_band_string(index, string, lstr)
1097 : : integer,intent(in) :: index
1098 : : integer,intent(out) :: lstr
1099 : : character(len=3),intent(out) :: string
1100 : :
1101 : : ! here put a check on the string length
1102 : :
1103 : 0 : string=' '
1104 [ # # ]: 0 : if(index < 10) then
1105 : 0 : WRITE( string(1:1), '(I1)' ) index
1106 : 0 : lstr=1
1107 [ # # ]: 0 : elseif(index < 100) then
1108 : 0 : WRITE( string(1:2), '(I2)' ) index
1109 : 0 : lstr=2
1110 [ # # ]: 0 : elseif(index < 1000) then
1111 : 0 : WRITE( string(1:3), '(I3)' ) index
1112 : 0 : lstr=3
1113 : : endif
1114 : :
1115 : 0 : string=trim(adjustl(string))
1116 : 0 : end subroutine define_band_string
1117 : :
1118 : 0 : end module m_gstore_converters
|
