[Q-e-developers] Feature in 'dynmat.x' when using the xml-representation for dynamical matrix

[Ari P Seitsonen]

Dear Colleagues,

   We are in the middle of a Summer School here in Paris, and while 
preparing the exercises for his Day Fabio (Finocchi) pointed out a 
"feature" in the 'dynmat' code: When using the xml representation for the 
dynamical matrix from the PHonon code the unit cell gets out written out 
wrong into the .axsf file. I just reproduced the error - 'dynmat.x' works 
correctly, however, if using non-xml format for the dynamical matrix. Here 
my short analysis of the problem, using the attached example of 
eight-atom cell of silicon:

   In the xml version the unit cell is written in the format

-<clip>------------------------------------------------
...
     <CELL_DIMENSIONS type="real" size="6">
  1.026121290156904E+001
  0.000000000000000E+000
  0.000000000000000E+000
  0.000000000000000E+000
  0.000000000000000E+000
  0.000000000000000E+000
     </CELL_DIMENSIONS>
     <AT type="real" size="9" columns="3">
  1.000000000000000E+000  0.000000000000000E+000  0.000000000000000E+000
  0.000000000000000E+000  1.000000000000000E+000  0.000000000000000E+000
  0.000000000000000E+000  0.000000000000000E+000  1.000000000000000E+000
     </AT>
...
-</clip>------------------------------------------------

Thus "at" is in the units of 'celldm(1)'. In 'PHonon/PH/dynmat.f90', lines 
165-170, we have

-<clip>------------------------------------------------
          CALL read_dyn_mat_header(ntyp, nat, ibrav, nspin_mag, &
                  celldm, at, bg, omega, atm, amass_, tau, ityp, &
                  m_loc, nqs, lrigid, eps0, zstar, lraman, dchi_dtau)
          IF (nqs /= 1) CALL errore('dynmat','only q=0 matrix allowed',1)
          a0=celldm(1) ! define alat
          at = at / a0 !  bring at in units of alat
-</clip>------------------------------------------------

so now the lattice vectors in 'at' are in units of 1/a0. So it does not 
help much if in subroutine 'writexsf()' in file 
'PHonon/PH/write_eigenvectors.f90', line 217, we have

-<clip>------------------------------------------------
   write(iout,'(2(3F15.9/),3f15.9)') at(:,:)*a0*BOHR_RADIUS_ANGS
-</clip>------------------------------------------------

as the output in the .axsf file turns out to be
-<clip>------------------------------------------------
PRIMVEC
     0.529177209    0.000000000    0.000000000
     0.000000000    0.529177209    0.000000000
     0.000000000    0.000000000    0.529177209
-</clip>------------------------------------------------
no matter what the lattice constant was in reality. When trying to 
visualise the .axsf file in 'XCrysDen' the output is quite sqeezed a 
crystal.

   I am not sure about the fix, you Wiser Ones please help: The division by 
'a0' should be i) present if the file is not in XML format yet ii) absent 
if using the XML format.

   Apologies if this is already a known issue: I obtained it with the 
version 5.2.0 and revision 11709 of the SVN trunk that I just checked out. 
I performed the tests with gcc/gfortran version 4.9.2, OpenMPI 1.8.3, 
recent ATLAS library, NetLib version of ScaLAPACK 2.0.2, with command

#> mpirun -np 4 pw.x -in pw.in > pw.out && mpirun -np 4 ph.x -in ph.in > ph.out && dynmat.x -in dynmat.in > dynmat.out

     Greetings from Rainy Paris 11ème,

        apsi

-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-
   Ari Paavo Seitsonen / Ari.P.Seitsonen at iki.fi / http://www.iki.fi/~apsi/
   Ecole Normale Supérieure (ENS), Département de Chimie, Paris
   Mobile (F) : +33 789 37 24 25    (CH) : +41 79 71 90 935
-------------- next part --------------
A non-text attachment was scrubbed...
Name: QE-ph-5.2.0-dynmat_problem.tar
Type: application/x-tar
Size: 419840 bytes
Desc: 
URL: <http://lists.quantum-espresso.org/pipermail/developers/attachments/20150826/886facce/attachment.tar>