The ABINIT Users Mailing List ( CLOSED )

["Matteo Giantomassi" <Matteo.Giantomassi@uclouvain.be>]

>>  I had a look at my implementation and it seems that I took into
>>  account the possibility of calculating atom-projected PHDOS
>>   The array pjdos is indeed defined as
>>
>> pjdos(nomega,3,natom)
>>
>> so it contains the partial projected PHDOS for each atom and for
>> each cartesian direction. One can easily extract from this array
>> the contribution to the PHDOS given by a particular atom.
>>
> ----------------------------------------
> Zhenhua: It seem not like that. I found the  PJDOSs are given type by
> type, no atom by atom.

 Dear Zhenhua,

 The values reported in the PHDOS file are relative to atom types
 since the routine writes on file the quantities stored in
pjdos_typ(nomega,ntypat):

do io=1,nomega
 write(unt,'(3f11.5)',advance='NO')omega(io),phdos(io),phdos_int(io)
 do itype=1,ntypat
  write(unt,frmt,advance='NO')pjdos_typ(io,itype),pjdos_typ_int(io,itype)
 end do
 write(unt,*)
end do

pjdos_typ is obtained by just summing over the atoms of the same type
and over the three Cartesian directions (few lines above the previous
piece of code)

 do iat=1,natom
  itype=typat(iat)
  do io=1,nomega
   pjdos_xyz_typ(io,:,itype)=pjdos_xyz_typ(io,:,itype)+pjdos(io,:,iat)
   pjdos_typ(io,itype)=pjdos_typ(io,itype)+sum(pjdos(io,:,iat))
  end do
  if (anaddb_dtset%prtdos==2) then
   do io=1,nomega
    pjdos_typ_int(io,itype)=pjdos_typ_int(io,itype)+SUM(pjdos_int(io,:,iat))
   end do
  end if
 end do

I decided to use this verbose approach since it's easy to extract
selected contributions from the two arrays pjdos_xyz_type and pjdos

If you need the atom-type projected contribution to the total DOS,
you can easily obtain this information from pjdos(nomega,3,natom)
by summing over the three Cartesian components.

pjdos_xyz_typ(nomega,3,ntypat) gives, for each Cartesian direction, the
contribution
to the PHDOS given by a particular atom type.


> If I want to get PJDOS atom by atom, I have to trick the routine by
> setting 'ntypat 2  znucl  13  13'
> I use V5.6.4 now. I don't know whether the problem is improved in the
> version going on.

 Well I have no idea of the system you are studying but, to be on-the-safe
side,
 I suggest avoiding dirty tricks unless you know what you are doing.
 Everything you need is already calculated inside mkphdos.
 You just have to add a print inside the routine at the point
 where the quantity you are interested in is available.

>
> I found there is a file 'out_PHFRQ' generated when done PHDOS calcualtion,
> whose second line is
> 'phonon frequencies (in Ha) on qph1l list of qpoints'
> However, the just one line date in this file. So I want to know whether
> the output of this file is ok.
> BTW: It seem the purpose of this file is to include the projected PHDOS of
> particular Qpoint, which is also very useful, interesting, and original
> drive of this post.


 I haven't written this part but, looking at the code (anaddb.F90),
 it seems that out_PHFRQ is used to store the phonon frequencies
 for the list of q-points specified by qph1l.

 Likely your input file doesn't specify qph1l thus it's not strange that
 the file is empty. Please have a look at the meaning of nph1l and qph1l

 Be careful that out_PHFRQ doesn't report the atomic eigendisplacements
 thus it cannot be used to analyze the projected PHDOS.
 If I remember well it's possible to plot phonon bands showing the
contribution
 given by atom types using eivec==4. The only way to have projected PHDOS
is, as far as I know, using mkphdos.F90.

>
>
>>
>>>
>>>  No example and reference output files are provided yet (I know, I
>>> know,
>>> shame on me!).
>>>  Digging into my files I found an anaddb input that can be used to
>>> calculate the PHDOS,
>>>  The flags that have to be added in the input file are:
>>>
>>> ifcflag  1  #IFC's are needed
>>>
>>> #DOS related variables
>>> prtdos  1              # 1 for gaussian, 2 for tetrahedrons
>>> ng2qpt  35 35 7        # Dense q-mesh for Fourier interpolation
>>> q2shft  0.0 0.0 0.0    # Shift for the dense mesh
>>> dosdeltae  9.112670E-7 # 0.1 cm-1  frequency step for PHDOS
>>> dosmear    1.366900E-5 # 3 cm-1    broadening (only for Gaussian
>>> method)
>>>
>>> At the end of the calculation the results are stored in the PHDOS file
>>> There are several columns:
>>>
>>>  # omega     PHDOS    IPHDOS   PJDOS[1]  IPJDOS[1] ..
>>>
>>>  PHDOS refers to the total phonon DOS, IPHDOS is the integrated DOS,
>>> then
>>> we have the
>>>  contribution due to all the atoms of the first type as well as the
>>> corresponding integrated dos.
>>>
>>> 3 Comments:
>>>
>>> 1) To call mkphdos you have to modify a bit 17ddb/invars9
>>>   commenting out the call to leave_new in the following piece of code
>>>
>>>  if(anaddb_dtset%prtdos/=0 .and. anaddb_dtset%ifcflag/=1) then
>>>   write(message, '(5a)' )&
>>> &  ' invars9 : ERROR -',ch10,&
>>> &  '  ifcflag must be 1 when the calculation of the phonon DOS is
>>> required
>>> ',ch10,&
>>> &  '  Action : correct ifcflag in your input file.'
>>>   call wrtout(6,message,'COLL')
>>>   call leave_new('COLL')
>>>  end if
>>>
>>>
>>> 2) If I remember well in old versions of abinit there's a small bug
>>> (just
>>> a missing factor)
>>>    in the calculation of the integrated type-projected PHDOS. So don't
>>> rely on the columns
>>>    reporting the type-projected IPHDOS. If you need, you can easily
>>> integrate the curve using
>>>    a plotting software such as xmgrace.

this problem has been solved in 5.7. The values obtained with the gaussian
technique
were correct, there was a small bug only in case of tetrahedrons. Anyway
the tetrahedron technique is not available in versions < 5.7


 Hope it helps.
 Best Regards,
 Matteo Giantomassi