commit 6b46e89eed8ddc348bec58b2b08ff1e54e04d493
Author: mitsuaki1987 <kawamitsuaki@gmail.com>
Date:   Tue Feb 13 17:10:29 2018 +0900

    Compute atom projected phonon-DOS.

diff --git a/PHonon/PH/matdyn.f90 b/PHonon/PH/matdyn.f90
index 24be505..e09e21c 100644
--- a/PHonon/PH/matdyn.f90
+++ b/PHonon/PH/matdyn.f90
@@ -144,7 +144,6 @@ PROGRAM matdyn
   USE rap_point_group,  ONLY : code_group
   USE bz_form,    ONLY : transform_label_coord
   USE parser,     ONLY : read_line
-  USE ktetra,     ONLY : tetra_dos_t
   USE rigid,       ONLY: dyndiag, nonanal, nonanal_ifc
 
   USE ifconstants, ONLY : frc, atm, zeu, tau_blk, ityp_blk, m_loc
@@ -166,7 +165,7 @@ PROGRAM matdyn
   LOGICAL :: dos, has_zstar, q_in_cryst_coord, eigen_similarity
   COMPLEX(DP), ALLOCATABLE :: dyn(:,:,:,:), dyn_blk(:,:,:,:), frc_ifc(:,:,:,:)
   COMPLEX(DP), ALLOCATABLE :: z(:,:)
-  REAL(DP), ALLOCATABLE:: tau(:,:), q(:,:), w2(:,:), freq(:,:), wq(:)
+  REAL(DP), ALLOCATABLE:: tau(:,:), q(:,:), w2(:,:), freq(:,:), wq(:), dynq(:,:,:), DOSofE(:)
   INTEGER, ALLOCATABLE:: ityp(:), itau_blk(:)
   REAL(DP) ::     omega,alat, &! cell parameters and volume
                   at_blk(3,3), bg_blk(3,3),  &! original cell
@@ -186,7 +185,7 @@ PROGRAM matdyn
   !
   LOGICAL :: readtau, la2F, xmlifc, lo_to_split, na_ifc, fd, nosym,  loto_2d 
   !
-  REAL(DP) :: qhat(3), qh, DeltaE, Emin=0._dp, Emax, E, DOSofE(2), qq
+  REAL(DP) :: qhat(3), qh, DeltaE, Emin=0._dp, Emax, E, qq
   REAL(DP) :: delta, pathL
   REAL(DP), ALLOCATABLE :: xqaux(:,:)
   INTEGER, ALLOCATABLE :: nqb(:)
@@ -210,6 +209,8 @@ PROGRAM matdyn
   CHARACTER(LEN=10) :: point_label_type
   CHARACTER(len=80) :: k_points = 'tpiba'
   !
+  REAL(DP), external       :: dos_gam
+  !
   NAMELIST /input/ flfrc, amass, asr, flfrq, flvec, fleig, at, dos,  &
        &           fldos, nk1, nk2, nk3, l1, l2, l3, ntyp, readtau, fltau, &
        &           la2F, ndos, DeltaE, q_in_band_form, q_in_cryst_coord, &
@@ -422,6 +423,7 @@ PROGRAM matdyn
         CALL mp_bcast(nq, ionode_id, world_comm)
         ALLOCATE ( q(3,nq) )
         IF (.NOT.q_in_band_form) THEN
+           ALLOCATE(wq(nq))
            DO n = 1,nq
               IF (ionode) READ (5,*) (q(i,n),i=1,3)
            END DO
@@ -530,7 +532,8 @@ PROGRAM matdyn
      END IF
 
      ALLOCATE ( dyn(3,3,nat,nat), dyn_blk(3,3,nat_blk,nat_blk) )
-     ALLOCATE ( z(3*nat,3*nat), w2(3*nat,nq), f_of_q(3,3,nat,nat) )
+     ALLOCATE ( z(3*nat,3*nat), w2(3*nat,nq), f_of_q(3,3,nat,nat), &
+     &          dynq(3*nat,nq,nat), DosofE(nat) )
      ALLOCATE ( tmp_w2(3*nat), abs_similarity(3*nat,3*nat), mask(3*nat) )
 
      if(la2F.and.ionode) open(unit=300,file='dyna2F',status='unknown')
@@ -615,6 +618,14 @@ PROGRAM matdyn
         
 
         CALL dyndiag(nat,ntyp,amass,ityp,dyn,w2(1,n),z)
+        ! Atom projection of dynamical matrix
+        DO i = 1, 3*nat
+           DO na = 1, nat
+              dynq(i, n, na) = DOT_PRODUCT(z(3*(na-1)+1:3*na, i), &
+              &                            z(3*(na-1)+1:3*na, i)  ) &
+              &              * amu_ry * amass(ityp(na))
+           END DO
+        END DO
         IF (ionode.and.iout_eig.ne.0) &
          & CALL write_eigenvectors(nat,ntyp,amass,ityp,q(1,n),w2(1,n),z,iout_eig)
         !
@@ -728,16 +739,26 @@ PROGRAM matdyn
            ndos = NINT ( (Emax - Emin) / DeltaE + 1.51d0 )
         end if
         IF (ionode) OPEN (unit=2,file=fldos,status='unknown',form='formatted')
+        IF (ionode) WRITE (2, *) "# Frequency[cm^-1] DOS PDOS"
+        DO na = 1, nat
+           dynq(1:3*nat,1:nq,na) = dynq(1:3*nat,1:nq,na) * freq(1:3*nat,1:nq)
+        END DO
         DO n= 1, ndos
            E = Emin + (n - 1) * DeltaE
-           CALL tetra_dos_t(freq, 1, 3*nat, nq, E, DOSofE)
+           DO na = 1, nat
+              DOSofE(na) = 0d0
+              DO i = 1, 3*nat
+                 DOSofE(na) = DOSofE(na) &
+                 & + dos_gam(3*nat, nq, i, dynq(1:3*nat,1:nq,na), freq, E)
+              END DO
+           END DO
            !
            ! The factor 0.5 corrects for the factor 2 in dos_t,
            ! that accounts for the spin in the electron DOS.
            !
            !WRITE (2, '(F15.10,F15.2,F15.6,F20.5)') &
            !     E, E*RY_TO_CMM1, E*RY_TO_THZ, 0.5d0*DOSofE(1)
-           IF (ionode) WRITE (2, '(ES12.4,ES12.4)') E, 0.5d0*DOSofE(1)
+           IF (ionode) WRITE (2, '(ES12.4,1000ES12.4)') E, SUM(DOSofE(1:nat)), DOSofE(1:nat)
         END DO
         IF (ionode) CLOSE(unit=2)
      END IF  !dos
@@ -2044,7 +2065,6 @@ SUBROUTINE a2Fdos &
   USE mp_images,   ONLY : intra_image_comm
   USE ifconstants, ONLY : zeu, tau_blk
   USE constants,  ONLY : pi, RY_TO_THZ
-  USE ktetra,     ONLY : tetra_init
   USE constants, ONLY : K_BOLTZMANN_RY
   !
   IMPLICIT NONE
@@ -2198,8 +2218,7 @@ SUBROUTINE a2Fdos &
               !
            enddo
            lambda = lambda + 2.d0 * dos_tot/E * DeltaE
-           write (ifn, '(3X,2F12.6)') E, dos_tot
-           write (ifn, '(6F16.8)') (dos_a2F(j),j=1,nmodes)
+           write (ifn, '(3X,1000E16.6)') E, dos_tot, dos_a2F(1:nmodes)
         enddo  !ndos
         write(ifn,*) " lambda =",lambda,'   Delta = ',DeltaE
         close (ifn)