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
|
