1!
2! Copyright (C) 2010 Quantum ESPRESSO group
3! This file is distributed under the terms of the
4! GNU General Public License. See the file `License'
5! in the root directory of the present distribution,
6! or http://www.gnu.org/copyleft/gpl.txt .
7!
8!
9!----------------------------------------------------------------------------
10MODULE io_dyn_mat
11  !----------------------------------------------------------------------------
12  !
13  ! ... this module contains methods to read and write the dynamical
14  !     matrix and the interatomic force constants files in xml format.
15  !
16  USE kinds,     ONLY : DP
17  USE io_global, ONLY : ionode, ionode_id
18  USE mp_images, ONLY : intra_image_comm
19  USE mp,        ONLY : mp_bcast
20  USE xmltools
21  !
22  IMPLICIT NONE
23  !
24  SAVE
25  !
26  PRIVATE
27  !
28  PUBLIC :: write_dyn_mat_header, write_dyn_mat, write_dyn_mat_tail, &
29            write_ifc, read_dyn_mat_param, read_dyn_mat_header, read_dyn_mat, &
30            read_dyn_mat_tail, read_ifc, read_ifc_param
31  !
32  INTEGER, PRIVATE :: iunout
33  !
34  CONTAINS
35    !
36    SUBROUTINE write_dyn_mat_header( fildyn, ntyp, nat, ibrav, nspin_mag,  &
37               celldm, at, bg, omega, atm, amass, tau, ityp, m_loc, &
38               nqs, epsil, zstareu, lraman, ramtns)
39    !
40  USE constants,     ONLY : FPI, BOHR_RADIUS_ANGS
41
42    INTEGER, INTENT(IN) :: ntyp, nat, ibrav, nspin_mag, nqs
43    CHARACTER(LEN=256), INTENT(IN) :: fildyn
44    CHARACTER(LEN=3), INTENT(IN) :: atm(ntyp)
45    REAL(DP), INTENT(IN) :: celldm(6)
46    REAL(DP), INTENT(IN) :: at(3,3)
47    REAL(DP), INTENT(IN) :: bg(3,3)
48    REAL(DP), INTENT(IN) :: omega
49    REAL(DP), INTENT(IN) :: amass(ntyp)
50    REAL(DP), INTENT(IN) :: tau(3,nat)
51    REAL(DP), INTENT(IN) :: m_loc(3,nat)
52    REAL(DP), INTENT(IN), OPTIONAL :: epsil(3,3)
53    REAL(DP), INTENT(IN), OPTIONAL :: zstareu(3,3,nat)
54    LOGICAL, INTENT(IN), OPTIONAL  :: lraman
55    REAL(DP), INTENT(IN), OPTIONAL :: ramtns(3,3,3,nat)
56
57    INTEGER, INTENT(IN) :: ityp(nat)
58
59    LOGICAL :: epsil_,raman_, zstareu_
60
61    INTEGER :: na, nt, kc
62    REAL(DP) :: aux(3,3)
63    REAL (DP), PARAMETER ::   convfact = BOHR_RADIUS_ANGS**2
64
65    IF ( ionode ) THEN
66       !
67       ! ... open XML descriptor
68       !
69       iunout = xml_openfile (TRIM( fildyn ) // '.xml' )
70       !
71    ENDIF
72    CALL mp_bcast( iunout, ionode_id, intra_image_comm )
73    !
74    IF ( iunout == -1 ) CALL errore( 'write_dyn_mat_header', &
75         'error opening the dyn mat file ', 1 )
76    !
77    IF (ionode) THEN
78       !
79       call add_attr( 'version','1.0')
80       call add_attr( 'encoding','UTF-8')
81       CALL xmlw_writetag ( 'xml', '?' )
82       CALL xmlw_opentag ( 'Root' )
83       !
84       CALL xmlw_opentag("GEOMETRY_INFO" )
85       !
86       CALL xmlw_writetag ( "NUMBER_OF_TYPES", ntyp )
87       CALL xmlw_writetag( "NUMBER_OF_ATOMS", nat )
88       CALL xmlw_writetag( "BRAVAIS_LATTICE_INDEX", ibrav )
89       CALL xmlw_writetag( "SPIN_COMPONENTS", nspin_mag )
90       CALL xmlw_writetag( "CELL_DIMENSIONS", celldm )
91       CALL xmlw_writetag( "AT", at )
92       CALL xmlw_writetag( "BG", bg )
93       CALL xmlw_writetag( "UNIT_CELL_VOLUME_AU", omega )
94       DO nt=1, ntyp
95          CALL xmlw_writetag( "TYPE_NAME."//i2c(nt),atm(nt))
96          CALL xmlw_writetag( "MASS." // i2c(nt),amass(nt))
97       ENDDO
98       DO na=1,nat
99          CALL add_attr( "SPECIES", atm(ityp(na)) )
100          CALL add_attr( "INDEX", ityp(na) )
101          CALL add_attr( "TAU", &
102               r2c(tau(1,na)) //' '// r2c(tau(2,na)) //' '// r2c(tau(3,na)) )
103          CALL xmlw_writetag( "ATOM." // i2c(na), '' )
104          IF (nspin_mag==4) &
105             CALL xmlw_writetag( "STARTING_MAG_."//i2c(na), m_loc(:,na) )
106       END DO
107       CALL xmlw_writetag( "NUMBER_OF_Q",nqs )
108       CALL xmlw_closetag( )
109       !
110       epsil_=.false.
111       IF (present(epsil)) epsil_=.true.
112       zstareu_=.false.
113       IF (present(zstareu)) zstareu_=.true.
114       raman_=.false.
115       IF (PRESENT(lraman)) raman_=.true.
116       !
117       CALL add_attr ( "epsil", epsil_)
118       CALL add_attr ( "zstar", zstareu_)
119       CALL add_attr ( "raman", raman_)
120       CALL xmlw_opentag( "DIELECTRIC_PROPERTIES" )
121       IF ( epsil_ ) THEN
122          CALL xmlw_writetag( "EPSILON",epsil)
123          IF ( zstareu_ ) THEN
124             CALL xmlw_opentag( "ZSTAR" )
125             DO na=1, nat
126                CALL xmlw_writetag( "Z_AT_."//i2c(na), zstareu(:,:,na) )
127             ENDDO
128             CALL xmlw_closetag( )
129          ENDIF
130          IF ( raman_) THEN
131             CALL xmlw_opentag( "RAMAN_TENSOR_A2")
132             DO na = 1, nat
133                DO kc = 1, 3
134                   aux(:,:) = ramtns(:, :, kc, na)*omega/fpi*convfact
135                   CALL xmlw_writetag( "RAMAN_S_ALPHA."//i2c(na)//'.'//i2c(kc),&
136                           aux )
137                ENDDO
138             ENDDO
139             CALL xmlw_closetag( )
140          ENDIF
141       ENDIF
142       CALL xmlw_closetag( )
143    ENDIF
144    !
145    RETURN
146    END SUBROUTINE write_dyn_mat_header
147
148    SUBROUTINE write_dyn_mat(nat,iq,xq,phi)
149
150    INTEGER, INTENT(IN) :: nat, iq
151    REAL(DP), INTENT(IN) :: xq(3)
152    COMPLEX(DP), INTENT(IN) :: phi(3,3,nat,nat)
153
154    INTEGER :: na, nb
155
156    IF (.NOT.ionode) RETURN
157
158    CALL xmlw_opentag( "DYNAMICAL_MAT_."//i2c(iq) )
159
160    CALL xmlw_writetag( "Q_POINT", xq )
161
162    DO na=1, nat
163       DO nb=1,nat
164          CALL xmlw_writetag( "PHI."//i2c(na)//'.'//i2c(nb),phi(:,:,na,nb) )
165       ENDDO
166    ENDDO
167
168    CALL xmlw_closetag( )
169
170    RETURN
171    END SUBROUTINE write_dyn_mat
172
173    SUBROUTINE write_dyn_mat_tail(nat,omega2,u)
174
175    USE constants, ONLY : RY_TO_THZ, RY_TO_CMM1
176
177    INTEGER, INTENT(IN) :: nat
178    REAL(DP), INTENT(IN) :: omega2(3*nat)
179    COMPLEX(DP), INTENT(IN) :: u(3*nat,3*nat)
180
181    REAL(DP) :: omega(2), om
182    INTEGER :: mu
183
184    IF (.NOT. ionode) RETURN
185
186    CALL xmlw_opentag( "FREQUENCIES_THZ_CMM1" )
187
188    DO mu=1,3*nat
189       om = SIGN( SQRT( ABS(omega2(mu)) ),  omega2(mu) )
190       omega(1) = om * RY_TO_THZ
191       omega(2) = om * RY_TO_CMM1
192       CALL xmlw_writetag( "OMEGA."//i2c(mu), omega )
193       CALL xmlw_writetag( "DISPLACEMENT."//i2c(mu), u(:,mu) )
194    END DO
195    CALL xmlw_closetag( )
196    !
197    CALL xmlw_closetag( ) ! Root
198    CALL xml_closefile( )
199    !
200    RETURN
201    END SUBROUTINE write_dyn_mat_tail
202
203    SUBROUTINE write_ifc( nr1, nr2, nr3, nat, phid)
204
205    INTEGER, INTENT(IN) :: nr1, nr2, nr3, nat
206    COMPLEX(DP), INTENT(IN) :: phid(nr1*nr2*nr3,3,3,nat,nat)
207
208    INTEGER :: meshfft(3)
209    INTEGER :: na, nb, nn, m1, m2, m3
210    REAL(DP) :: aux(3,3)
211
212    IF (.NOT.ionode) RETURN
213
214    meshfft(1)=nr1
215    meshfft(2)=nr2
216    meshfft(3)=nr3
217    CALL xmlw_opentag( "INTERATOMIC_FORCE_CONSTANTS" )
218    CALL xmlw_writetag( "MESH_NQ1_NQ2_NQ3", meshfft )
219
220    DO na=1,nat
221       DO nb=1,nat
222          nn=0
223          DO m3=1,nr3
224             DO m2=1,nr2
225                DO m1=1,nr1
226                   nn=nn+1
227                   CALL xmlw_opentag( "s_s1_m1_m2_m3." // i2c(na) //'.'// &
228                      i2c(nb) //'.'// i2c(m1) //'.'// i2c(m2) //'.'// i2c(m3))
229                   aux(:,:)=DBLE(phid(nn,:,:,na,nb))
230                   CALL xmlw_writetag( 'IFC', aux )
231                   CALL xmlw_closetag( )
232                ENDDO
233             ENDDO
234          ENDDO
235       ENDDO
236    ENDDO
237    CALL xmlw_closetag()
238    !
239    CALL xmlw_closetag() ! Root
240    CALL xml_closefile()
241    !
242    !----------------------------------------------------------------------------
243    END SUBROUTINE write_ifc
244    !----------------------------------------------------------------------------
245    !
246    !----------------------------------------------------------------------------
247    SUBROUTINE read_dyn_mat_param(fildyn, ntyp, nat)
248    !----------------------------------------------------------------------------
249    !!
250    !! Read paramters from the dynamical matrix
251    !!
252    USE io_global,   ONLY : ionode
253    !
254    IMPLICIT NONE
255    !
256    CHARACTER(LEN = 256), INTENT(in) :: fildyn
257    !! Name of the file to read
258    INTEGER, INTENT(out) :: ntyp
259    !! Number of type of atoms
260    INTEGER, INTENT(out) :: nat
261    !! Number of atoms
262    !
263    ! Open XML descriptor
264    !
265    IF (ionode) iunout = xml_openfile( TRIM(fildyn) // '.xml')
266    !
267    CALL mp_bcast(iunout, ionode_id, intra_image_comm)
268    IF ( iunout == -1 ) &
269         CALL errore('read_dyn_mat_param', 'error opening the dyn mat file ',1)
270    !
271    IF (ionode) THEN
272      CALL xmlr_opentag( "GEOMETRY_INFO")
273      CALL xmlr_readtag( "NUMBER_OF_TYPES", ntyp)
274      CALL xmlr_readtag( "NUMBER_OF_ATOMS", nat)
275      CALL xmlr_closetag() ! GEOMETRY_INFO
276      REWIND(iunout)
277    ENDIF
278    !
279    CALL mp_bcast(ntyp, ionode_id, intra_image_comm)
280    CALL mp_bcast(nat, ionode_id, intra_image_comm)
281    !
282    RETURN
283    !----------------------------------------------------------------------------
284    END SUBROUTINE read_dyn_mat_param
285    !----------------------------------------------------------------------------
286    !
287    !----------------------------------------------------------------------------
288    SUBROUTINE read_dyn_mat_header(ntyp, nat, ibrav, nspin_mag,     &
289               celldm, at, bg, omega, atm, amass, tau, ityp, m_loc, &
290               nqs, lrigid, epsil, zstareu, lraman, ramtns)
291    !----------------------------------------------------------------------------
292    !!
293    !! Read the dynamical matrix
294    !!
295    USE kinds,       ONLY : DP
296    USE io_global,   ONLY : ionode
297    USE xmltools
298    !
299    IMPLICIT NONE
300    !
301    CHARACTER(LEN = 3), INTENT(out) :: atm(ntyp)
302    !! Atom
303    LOGICAL, INTENT(out), OPTIONAL :: lrigid
304    !!
305    LOGICAL, INTENT(out), OPTIONAL :: lraman
306    !! Raman
307    INTEGER, INTENT(in) :: ntyp
308    !! Number of type of atoms
309    INTEGER, INTENT(in) :: nat
310    !! Number of atoms
311    INTEGER, INTENT(out) :: ibrav
312    !! Bravais lattice
313    INTEGER, INTENT(out) :: nspin_mag
314    !!
315    INTEGER, INTENT(out) :: nqs
316    !!
317    INTEGER,  INTENT(out) :: ityp(nat)
318    !! Atom type
319    REAL(KIND = DP), INTENT(out) :: celldm(6)
320    !! Celldm
321    REAL(KIND = DP), INTENT(out) :: at(3, 3)
322    !! Real-space lattice
323    REAL(KIND = DP), INTENT(out) :: bg(3, 3)
324    !! Reciprocal-space latrice
325    REAL(KIND = DP), INTENT(out) :: omega
326    !! Volume of primitive cell
327    REAL(KIND = DP), INTENT(out) :: amass(ntyp)
328    !! Atom mass
329    REAL(KIND = DP), INTENT(out) :: tau(3, nat)
330    !! Atom position
331    REAL(KIND = DP), INTENT(out) :: m_loc(3, nat)
332    !!
333    REAL(KIND = DP), INTENT(out), OPTIONAL :: epsil(3, 3)
334    !! Dielectric cst
335    REAL(KIND = DP), INTENT(out), OPTIONAL :: zstareu(3, 3, nat)
336    !!
337    REAL(KIND = DP), INTENT(out), OPTIONAL :: ramtns(3, 3, 3, nat)
338    !!
339    CHARACTER(LEN = 80) :: dummy
340    !!
341    LOGICAL :: found_z
342    !!
343    LOGICAL :: lrigid_
344    !!
345    LOGICAL :: raman_
346    !! Is Raman present
347    INTEGER :: nt
348    !! Type of atoms
349    INTEGER :: na
350    !! Number of atoms
351    INTEGER :: kc
352    !! Cartesian direction
353    INTEGER :: ierr
354    REAL(KIND = DP) :: aux(3, 3)
355    !! Auxiliary
356    !
357    IF (ionode) THEN
358      CALL xmlr_opentag("GEOMETRY_INFO")
359      CALL xmlr_readtag("BRAVAIS_LATTICE_INDEX", ibrav)
360      CALL xmlr_readtag("SPIN_COMPONENTS", nspin_mag)
361      CALL xmlr_readtag("CELL_DIMENSIONS", celldm)
362      CALL xmlr_readtag("AT", at)
363      CALL xmlr_readtag("BG", bg)
364      CALL xmlr_readtag("UNIT_CELL_VOLUME_AU", omega)
365      DO nt = 1, ntyp
366        CALL xmlr_readtag("TYPE_NAME."//i2c(nt), atm(nt))
367        CALL xmlr_readtag("MASS." // i2c(nt), amass(nt))
368      ENDDO
369      DO na = 1, nat
370        CALL xmlr_readtag("ATOM." // i2c(na), dummy)
371        CALL get_attr("INDEX",  ityp(na))
372        CALL get_attr("TAU", dummy )
373        READ(dummy,*) tau(1, na), tau(2, na), tau(3, na)
374        IF (nspin_mag == 4) THEN
375          CALL xmlr_readtag("STARTING_MAG_."//i2c(na), m_loc(:, na))
376        ENDIF
377      ENDDO
378      CALL xmlr_readtag("NUMBER_OF_Q", nqs)
379      CALL xmlr_closetag() ! GEOMETRY_INFO
380      !
381      IF (PRESENT(epsil)) THEN
382        CALL xmlr_opentag("DIELECTRIC_PROPERTIES", ierr)
383        IF (ierr == -1) THEN
384          IF (PRESENT(lrigid))  lrigid = .false.
385          IF (PRESENT(lraman))  lraman = .false.
386          epsil = 0.0_dp
387          IF (PRESENT(zstareu)) zstareu = 0.0_DP
388          IF (PRESENT(ramtns))  ramtns = 0.0_DP
389          GOTO 10
390        ENDIF
391        CALL get_attr("epsil", lrigid_)
392        IF (PRESENT(lrigid)) lrigid = lrigid_
393        CALL get_attr("zstar", found_z)
394        CALL get_attr("raman", raman_)
395        IF (PRESENT(lraman)) lraman = raman_
396        IF (lrigid_) THEN
397          CALL xmlr_readtag( "EPSILON", epsil)
398          IF (found_z) THEN
399            CALL xmlr_opentag( "ZSTAR" )
400            DO na = 1, nat
401              CALL xmlr_readtag( "Z_AT_."//i2c(na), aux(:, :))
402              IF (PRESENT(zstareu)) zstareu(:, :, na) = aux
403            ENDDO
404            CALL xmlr_closetag() ! ZSTAR
405          ELSE
406            IF (PRESENT(zstareu)) zstareu = 0.0_DP
407          ENDIF
408          IF (raman_) THEN
409            CALL xmlr_opentag("RAMAN_TENSOR_A2" )
410            IF (PRESENT(ramtns)) THEN
411              DO na = 1, nat
412                DO kc = 1, 3
413                  CALL xmlr_readtag("RAMAN_S_ALPHA."//i2c(na)//'.'//i2c(kc), aux)
414                  IF (PRESENT(ramtns)) ramtns(:, :, kc, na) = aux(:, :)
415                ENDDO
416              ENDDO
417            ELSE
418              IF (PRESENT(ramtns)) ramtns = 0.0_DP
419            ENDIF
420            CALL xmlr_closetag() ! RAMAN_TENSOR_A2
421          ENDIF
422        ELSE
423           IF (PRESENT(epsil)) epsil = 0.0_DP
424           IF (PRESENT(zstareu)) zstareu = 0.0_DP
425           IF (PRESENT(ramtns)) ramtns = 0.0_DP
426        ENDIF ! lrigid
427        CALL xmlr_closetag() ! DIELECTRIC_PROPERTIES
428      ENDIF ! epsil
429      10 CONTINUE
430    ENDIF ! ionode
431    CALL mp_bcast(ibrav, ionode_id, intra_image_comm)
432    CALL mp_bcast(nspin_mag, ionode_id, intra_image_comm)
433    CALL mp_bcast(celldm, ionode_id, intra_image_comm)
434    CALL mp_bcast(at, ionode_id, intra_image_comm)
435    CALL mp_bcast(bg, ionode_id, intra_image_comm)
436    CALL mp_bcast(omega, ionode_id, intra_image_comm)
437    CALL mp_bcast(atm, ionode_id, intra_image_comm)
438    CALL mp_bcast(amass, ionode_id, intra_image_comm)
439    CALL mp_bcast(ityp, ionode_id, intra_image_comm)
440    CALL mp_bcast(tau, ionode_id, intra_image_comm)
441    CALL mp_bcast(m_loc, ionode_id, intra_image_comm)
442    CALL mp_bcast(nqs, ionode_id, intra_image_comm)
443    IF (PRESENT(lrigid)) CALL mp_bcast(lrigid, ionode_id, intra_image_comm)
444    IF (PRESENT(epsil)) CALL mp_bcast(epsil, ionode_id, intra_image_comm)
445    IF (PRESENT(zstareu)) CALL mp_bcast(zstareu, ionode_id, intra_image_comm)
446    IF (PRESENT(lraman)) CALL mp_bcast(lraman, ionode_id, intra_image_comm)
447    IF (PRESENT(ramtns)) CALL mp_bcast(ramtns, ionode_id, intra_image_comm)
448    RETURN
449    !----------------------------------------------------------------------------
450    END SUBROUTINE read_dyn_mat_header
451    !----------------------------------------------------------------------------
452    !
453    !----------------------------------------------------------------------------
454    SUBROUTINE read_dyn_mat(nat, iq, xq, dyn)
455    !----------------------------------------------------------------------------
456    !!
457    !! This routine reads the dynamical matrix file. The file is assumed to
458    !! be already opened. iq is the number of the dynamical matrix to read.
459    !!
460    USE kinds,       ONLY : DP
461    USE io_global,   ONLY : ionode
462    !
463    IMPLICIT NONE
464    !
465    INTEGER, INTENT(in) :: nat
466    !! Number of atoms
467    INTEGER, INTENT(in) :: iq
468    !! Q-point index
469    REAL(KIND = DP), INTENT(out) :: xq(3)
470    !! Q-point value
471    COMPLEX(KIND = DP), INTENT(out) :: dyn(3, 3, nat, nat)
472    !! Dynamical matrix
473    !
474    ! Local variables
475    INTEGER :: na, nb
476    !! Number of atoms
477    !
478    IF (ionode) THEN
479      CALL xmlr_opentag("DYNAMICAL_MAT_."//i2c(iq))
480      CALL xmlr_readtag("Q_POINT", xq)
481      DO na = 1, nat
482        DO nb = 1,nat
483          CALL xmlr_readtag( "PHI."//i2c(na)//'.'//i2c(nb), dyn(:, :, na, nb))
484        ENDDO
485      ENDDO
486      CALL xmlr_closetag() ! DYNAMICAL_MAT_.
487    ENDIF
488    CALL mp_bcast(xq, ionode_id, intra_image_comm)
489    CALL mp_bcast(dyn, ionode_id, intra_image_comm)
490    RETURN
491    !----------------------------------------------------------------------------
492    END SUBROUTINE read_dyn_mat
493    !----------------------------------------------------------------------------
494    !
495    !----------------------------------------------------------------------------
496    SUBROUTINE read_dyn_mat_tail(nat, omega, u)
497    !----------------------------------------------------------------------------
498    !!
499    !! The output of the routine in a.u.
500    !!
501    USE kinds,     ONLY : DP
502    USE constants, ONLY : RY_TO_THZ
503    !
504    INTEGER, INTENT(in) :: nat
505    !! Number of atoms
506    REAL(KIND = DP), INTENT(out), OPTIONAL :: omega(3 * nat)
507    !! Phonon freq.
508    COMPLEX(KIND = DP), INTENT(out), OPTIONAL :: u(3 * nat, 3 * nat)
509    !! Eigen displacement vectors
510    !
511    ! Local variables
512    REAL(KIND = DP) :: omega_(2)
513    !! Phonon freq
514    INTEGER :: mu
515    !!
516    !
517    IF (PRESENT(u) .AND. .NOT. PRESENT(omega)) &
518       CALL errore('read_dyn_mat_tail','omega must be present to read u',1)
519
520    IF (ionode) THEN
521      IF (PRESENT(omega)) THEN
522        CALL xmlr_opentag("FREQUENCIES_THZ_CMM1")
523        DO mu = 1, 3 * nat
524          CALL xmlr_readtag("OMEGA."//i2c(mu), omega_)
525          omega(mu) = omega_(1) / RY_TO_THZ
526          IF (PRESENT(u)) CALL xmlr_readtag("DISPLACEMENT."//i2c(mu),u(:,mu))
527        END DO
528        CALL xmlr_closetag() ! FREQUENCIES_THZ_CMM1
529      ENDIF
530      CALL xml_closefile()
531    ENDIF
532    IF (PRESENT(omega)) CALL mp_bcast(omega, ionode_id, intra_image_comm)
533    IF (PRESENT(u)) CALL mp_bcast(u, ionode_id, intra_image_comm)
534    !
535    RETURN
536    !----------------------------------------------------------------------------
537    END SUBROUTINE read_dyn_mat_tail
538    !----------------------------------------------------------------------------
539    !
540    !----------------------------------------------------------------------------
541    SUBROUTINE read_ifc_param(nr1, nr2, nr3)
542    !----------------------------------------------------------------------------
543    !!
544    !! Read IFC parameters
545    !!
546    !! The following sequence should be used:
547    !! read_dyn_mat_param
548    !! read_dyn_mat_header
549    !! read_ifc_param
550    !! read_ifc
551    !
552    IMPLICIT NONE
553    !
554    INTEGER, INTENT(out) :: nr1, nr2, nr3
555    !! Grid size
556    ! Local varialbes
557    INTEGER :: meshfft(3)
558    !! Mesh
559    IF (ionode) THEN
560      CALL xmlr_opentag( "INTERATOMIC_FORCE_CONSTANTS")
561      CALL xmlr_readtag( "MESH_NQ1_NQ2_NQ3", meshfft)
562      nr1 = meshfft(1)
563      nr2 = meshfft(2)
564      nr3 = meshfft(3)
565      CALL xmlr_closetag( )
566    ENDIF
567    CALL mp_bcast(nr1, ionode_id, intra_image_comm)
568    CALL mp_bcast(nr2, ionode_id, intra_image_comm)
569    CALL mp_bcast(nr3, ionode_id, intra_image_comm)
570    RETURN
571    !----------------------------------------------------------------------------
572    END SUBROUTINE read_ifc_param
573    !----------------------------------------------------------------------------
574    !
575    !----------------------------------------------------------------------------
576    SUBROUTINE read_ifc(nr1, nr2, nr3, nat, phid)
577    !----------------------------------------------------------------------------
578    !!
579    !! Read IFC in XML format
580    !!
581    USE kinds,       ONLY : DP
582    USE io_global,   ONLY : ionode
583    !
584    IMPLICIT NONE
585    !
586    INTEGER, INTENT(in) :: nr1, nr2, nr3
587    !! Grid size
588    INTEGER, INTENT(in) :: nat
589    !! Number of atoms
590    REAL(KIND = DP), INTENT(out) :: phid(nr1 * nr2 * nr3, 3, 3, nat, nat)
591    !!
592    ! Local variables
593    INTEGER :: na, nb
594    !! Atoms
595    INTEGER :: nn, ierr
596    !!
597    INTEGER :: m1, m2, m3
598    !! nr dimension
599    REAL(KIND = DP) :: aux(3, 3)
600    !! Auxiliary
601    !
602    IF (ionode) THEN
603      CALL xmlr_opentag( "INTERATOMIC_FORCE_CONSTANTS", ierr)
604      DO na = 1, nat
605        DO nb = 1, nat
606          nn = 0
607          DO m3 = 1, nr3
608            DO m2 = 1, nr2
609              DO m1 = 1, nr1
610                nn = nn + 1
611                CALL xmlr_opentag( "s_s1_m1_m2_m3." // i2c(na) //'.'// &
612                      i2c(nb) //'.'// i2c(m1) //'.'// i2c(m2) //'.'// i2c(m3))
613                CALL xmlr_readtag( 'IFC', aux)
614                phid(nn, :, :, na, nb) = aux(:, :)
615                CALL xmlr_closetag( )
616              ENDDO ! m1
617            ENDDO ! m2
618          ENDDO ! m3
619        ENDDO ! nb
620      ENDDO ! na
621      CALL xmlr_closetag( )
622      CALL xmlr_closetag( ) ! Root
623      CALL xml_closefile( )
624    ENDIF
625    CALL mp_bcast(phid, ionode_id, intra_image_comm)
626    RETURN
627    !----------------------------------------------------------------------------
628    END SUBROUTINE read_ifc
629    !----------------------------------------------------------------------------
630  !----------------------------------------------------------------------------
631  END MODULE io_dyn_mat
632  !----------------------------------------------------------------------------
633