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

[Andrea Dal Corso]

Probably you are right. The division by a0 is a bug. In any case that
instruction is within a

IF (xmldyn) THEN
...
ELSE
ENDIF

so it will not be used in the non-xml case.

Andrea

On Wed, 2015-08-26 at 22:32 +0200, Ari P Seitsonen wrote:
> 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
> _______________________________________________ Q-e-developers mailing list Q-e-developers at qe-forge.org http://qe-forge.org/mailman/listinfo/q-e-developers

-- 
Andrea Dal Corso                    Tel. 0039-040-3787428
SISSA, Via Bonomea 265              Fax. 0039-040-3787249
I-34136 Trieste (Italy)             e-mail: dalcorso at sissa.it