TABLE OF CONTENTS


ABINIT/m_phonons [ Modules ]

[ Top ] [ Modules ]

NAME

 m_phonons

FUNCTION

 Module for the phonon density of states.
 Container type is defined, and destruction, print subroutines
 as well as the central mkphdos

COPYRIGHT

 Copyright (C) 1999-2018 ABINIT group (XG, MG, MJV, GMR)
 This file is distributed under the terms of the
 GNU General Public Licence, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

SOURCE

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_phonons
27 
28  use defs_basis
29  use m_errors
30  use m_xmpi
31  use m_abicore
32  use m_tetrahedron
33  use m_nctk
34  use iso_c_binding
35  use m_crystal_io
36  use m_atprj
37  use m_sortph
38  use m_ddb
39 #ifdef HAVE_NETCDF
40  use netcdf
41 #endif
42  use m_supercell
43 
44  use m_fstrings,        only : itoa, ftoa, sjoin, ktoa, strcat, basename, replace
45  use m_numeric_tools,   only : simpson_int, wrap2_pmhalf
46  use m_symtk,           only : matr3inv
47  use m_time,            only : cwtime
48  use m_io_tools,        only : open_file
49  use defs_abitypes,     only : dataset_type
50  use m_geometry,        only : mkrdim, symredcart
51  use m_dynmat,          only : gtdyn9, dfpt_phfrq, dfpt_prtph
52  use m_crystal,         only : crystal_t
53  use m_bz_mesh,         only : isamek, make_path, kpath_t, kpath_new, kpath_free
54  use m_ifc,             only : ifc_type, ifc_fourq, ifc_calcnwrite_nana_terms
55  use m_anaddb_dataset,  only : anaddb_dataset_type
56  use m_kpts,            only : kpts_ibz_from_kptrlatt, get_full_kgrid
57  use m_special_funcs,   only : bose_einstein
58  use m_sort,            only : sort_dp
59 
60  implicit none
61 
62  private
63 
64 ! TODO Write object to store the bands
65  public :: mkphbs                        ! Compute phonon band structure
66  public :: phonons_write_xmgrace         ! Write phonons bands in Xmgrace format.
67  public :: phonons_write_gnuplot         ! Write phonons bands in gnuplot format.
68  public :: ifc_mkphbs                    ! Compute the phonon band structure from the IFC and write data to file(s)
69  public :: dfpt_symph                    ! Determine the symmetry character of the different phonon modes at Gamma
70 
71  public :: zacharias_supercell_make
72  public :: zacharias_supercell_print
73  public :: thermal_supercell_make
74  public :: thermal_supercell_free
75  public :: thermal_supercell_print

m_phonons/dfpt_symph [ Functions ]

[ Top ] [ m_phonons ] [ Functions ]

NAME

 dfpt_symph

FUNCTION

 Determine the symmetry character of the different phonon modes.

INPUTS

 acell(3)=length scales of primitive translations (bohr)
 eigvec(2*3*natom*3*natom)=eigenvectors of the dynamical matrix
 indsym(4,nsym,natom)=indirect indexing array : for each
   isym,iatom, fourth element is label of atom into which iatom is sent by
   INVERSE of symmetry operation isym; first three elements are the primitive
   translations which must be subtracted after the transformation to get back
   to the original unit cell.
 iout=unit number to which output is written
 natom=number of atoms in unit cell
 nsym=number of space group symmetries
 phfrq(3*natom)=phonon frequencies (Hartree units)
 rprim(3,3)=dimensionless primitive translations in real space
 symrel(3,3,nsym)=matrices of the group symmetries (real space)

OUTPUT

PARENTS

      anaddb,m_phonons

CHILDREN

      matr3inv,mkrdim,wrtout

SOURCE

3573 subroutine dfpt_symph(iout,acell,eigvec,indsym,natom,nsym,phfrq,rprim,symrel)
3574 
3575 
3576 !This section has been created automatically by the script Abilint (TD).
3577 !Do not modify the following lines by hand.
3578 #undef ABI_FUNC
3579 #define ABI_FUNC 'dfpt_symph'
3580 !End of the abilint section
3581 
3582  implicit none
3583 
3584 !Arguments ------------------------------------
3585 !scalars
3586  integer,intent(in) :: iout,natom,nsym
3587 !arrays
3588  integer,intent(in) :: indsym(4,nsym,natom),symrel(3,3,nsym)
3589  real(dp),intent(in) :: acell(3),eigvec(2*3*natom*3*natom),phfrq(3*natom)
3590  real(dp),intent(in) :: rprim(3,3)
3591 
3592 !Local variables -------------------------
3593 !scalars
3594  integer :: iad1,iad2,iad3,iatom,idir,ii1,ii2,ii3,imode,isym,itol,jad,jatom,jj
3595  integer :: jmode,kk,ntol
3596  character(len=500) :: message
3597 !arrays
3598  integer,allocatable :: degeneracy(:),integer_characters(:),symind(:,:)
3599  real(dp) :: gprimd(3,3),rprimd(3,3)
3600  real(dp),allocatable :: eigvtr(:),redvec(:),redvtr(:),symph(:,:)
3601 
3602 !******************************************************************
3603 
3604 !Compute dimensional primitive translations rprimd and its inverse gprimd
3605  call mkrdim(acell,rprim,rprimd)
3606  call matr3inv(rprimd,gprimd)
3607 
3608 !Build the symmetry index (inverse of indsym(4,:,:))
3609  ABI_ALLOCATE(symind,(nsym,natom))
3610  do isym=1,nsym
3611    do iatom=1,natom
3612      symind(isym,indsym(4,isym,iatom))=iatom
3613    end do
3614  end do
3615 
3616  ABI_ALLOCATE(symph,(nsym,3*natom))
3617  ABI_ALLOCATE(redvec,(2*3*natom))
3618  ABI_ALLOCATE(redvtr,(2*3*natom))
3619  ABI_ALLOCATE(eigvtr,(2*3*natom))
3620 
3621 !Loop over the vibration modes
3622  do imode=1,3*natom
3623 
3624 !  Compute eigvec for this mode in reduced coordinates redvec
3625    do iatom=1,natom
3626      iad1=3*(iatom-1)+1
3627      ii1=2*3*natom*(imode-1)+2*(iad1-1)+1
3628      iad2=3*(iatom-1)+2
3629      ii2=2*3*natom*(imode-1)+2*(iad2-1)+1
3630      iad3=3*(iatom-1)+3
3631      ii3=2*3*natom*(imode-1)+2*(iad3-1)+1
3632      do idir=1,3
3633        jad=3*(iatom-1)+idir
3634        jj=2*(jad-1)+1
3635        redvec(jj)=gprimd(1,idir)*eigvec(ii1)+&
3636 &       gprimd(2,idir)*eigvec(ii2)+&
3637 &       gprimd(3,idir)*eigvec(ii3)
3638        redvec(jj+1)=gprimd(1,idir)*eigvec(ii1+1)+&
3639 &       gprimd(2,idir)*eigvec(ii2+1)+&
3640 &       gprimd(3,idir)*eigvec(ii3+1)
3641      end do !idir
3642    end do !iatom
3643 
3644 !  Apply each transformation to redvec and store at the correct location in redvtr (iatom -> jatom)
3645    do isym=1,nsym
3646      do iatom=1,natom
3647        jatom=symind(isym,iatom)
3648        iad1=3*(iatom-1)+1
3649        ii1=2*(iad1-1)+1
3650        iad2=3*(iatom-1)+2
3651        ii2=2*(iad2-1)+1
3652        iad3=3*(iatom-1)+3
3653        ii3=2*(iad3-1)+1
3654        do idir=1,3
3655          jad=3*(jatom-1)+idir
3656          jj=2*(jad-1)+1
3657          redvtr(jj)=dble(symrel(idir,1,isym))*redvec(ii1)+&
3658 &         dble(symrel(idir,2,isym))*redvec(ii2)+&
3659 &         dble(symrel(idir,3,isym))*redvec(ii3)
3660          redvtr(jj+1)=dble(symrel(idir,1,isym))*redvec(ii1+1)+&
3661 &         dble(symrel(idir,2,isym))*redvec(ii2+1)+&
3662 &         dble(symrel(idir,3,isym))*redvec(ii3+1)
3663 
3664        end do !idir
3665      end do !iatom
3666 
3667 !    Compute redvtr in cartesian coordinates eigvtr
3668      do iatom=1,natom
3669        iad1=3*(iatom-1)+1
3670        ii1=2*(iad1-1)+1
3671        iad2=3*(iatom-1)+2
3672        ii2=2*(iad2-1)+1
3673        iad3=3*(iatom-1)+3
3674        ii3=2*(iad3-1)+1
3675        do idir=1,3
3676          jad=3*(iatom-1)+idir
3677          jj=2*(jad-1)+1
3678          eigvtr(jj)=rprimd(idir,1)*redvtr(ii1)+&
3679 &         rprimd(idir,2)*redvtr(ii2)+&
3680 &         rprimd(idir,3)*redvtr(ii3)
3681          eigvtr(jj+1)=rprimd(idir,1)*redvtr(ii1+1)+&
3682 &         rprimd(idir,2)*redvtr(ii2+1)+&
3683 &         rprimd(idir,3)*redvtr(ii3+1)
3684        end do !idir
3685      end do !iatom
3686 
3687 !    Compute scalar product...
3688      symph(isym,imode)=0.0_dp
3689      do jad=1,3*natom
3690        jj=2*(jad-1)+1
3691        kk=2*3*natom*(imode-1)+2*(jad-1)+1
3692        symph(isym,imode)=symph(isym,imode)+eigvtr(jj)*eigvec(kk)+eigvtr(jj+1)*eigvec(kk+1)
3693      end do
3694 
3695    end do !isym
3696  end do !imode
3697 
3698 !Treat degeneracies (different tolerances will be tried)
3699 !Compute the order of the degeneracy, and
3700 !attribute it to the lowest of the degenerate modes
3701 !Also attribute the characters to the lowest mode
3702 !When all the characters are integers, consider that the
3703 !mode is non-degenerate. The maximum difference in frequency
3704 !that is tolerated is on the order of 4cm-1 (which is large...)
3705  ABI_ALLOCATE(degeneracy,(3*natom))
3706  ABI_ALLOCATE(integer_characters,(3*natom))
3707  degeneracy(:)=1
3708  integer_characters(:)=0
3709  do itol=1,20
3710    ntol=itol
3711    do imode=3*natom,2,-1
3712      if(integer_characters(imode)==0)then
3713        do jmode=imode-1,1,-1
3714          if(integer_characters(jmode)==0)then
3715            if(abs(phfrq(imode)-phfrq(jmode))<itol*tol6)then
3716              degeneracy(jmode)=degeneracy(jmode)+degeneracy(imode)
3717              degeneracy(imode)=0
3718              symph(:,jmode)=symph(:,jmode)+symph(:,imode)
3719              symph(:,imode)=0.0_dp
3720            end if
3721          end if !integer_characters(jmode)==0
3722        end do !jmode
3723      end if !integer_characters(imode)==0
3724    end do !imode
3725    do imode=1,3*natom
3726      if(maxval(abs( symph(:,imode)-nint(symph(:,imode)) ))<0.05_dp)then
3727        integer_characters(imode)=1
3728      end if
3729    end do
3730    if(sum(integer_characters(:))==3*natom)exit
3731  end do !itol
3732 
3733 !DEBUG
3734 !write(std_out,*)' dfpt_symph : degeneracy=',degeneracy(:)
3735 !ENDDEBUG
3736 
3737  write(message,'(a,a,es8.2,a)')ch10,&
3738 & ' Analysis of degeneracies and characters (maximum tolerance=',ntol*tol6,' a.u.)'
3739  call wrtout(iout,message,'COLL')
3740  call wrtout(std_out,message,'COLL')
3741 
3742  do imode=1,3*natom
3743    if(degeneracy(imode)/=0)then
3744      write(message,'(a,i4)') ' Symmetry characters of vibration mode #',imode
3745      call wrtout(iout,message,'COLL')
3746      call wrtout(std_out,message,'COLL')
3747      if(degeneracy(imode)>=2)then
3748        if(degeneracy(imode)==2) write(message,'(a,i4)') &
3749 &       '        degenerate with vibration mode #',imode+1
3750        if(degeneracy(imode)>=3) write(message,'(a,i4,a,i4)') &
3751 &       '       degenerate with vibration modes #',imode+1,' to ',imode+degeneracy(imode)-1
3752        call wrtout(iout,message,'COLL')
3753        call wrtout(std_out,message,'COLL')
3754      end if
3755      do jj=1,(nsym-1)/16+1
3756        write(message,'(16f5.1)') (symph(isym,imode),isym=(jj-1)*16+1,min(nsym,jj*16))
3757        call wrtout(iout,message,'COLL')
3758        call wrtout(std_out,message,'COLL')
3759      end do
3760    end if
3761  end do !imode
3762 
3763  ABI_DEALLOCATE(degeneracy)
3764  ABI_DEALLOCATE(integer_characters)
3765  ABI_DEALLOCATE(eigvtr)
3766  ABI_DEALLOCATE(redvtr)
3767  ABI_DEALLOCATE(redvec)
3768  ABI_DEALLOCATE(symph)
3769  ABI_DEALLOCATE(symind)
3770 
3771 end subroutine dfpt_symph

m_phonons/freeze_displ_allmodes [ Functions ]

[ Top ] [ m_phonons ] [ Functions ]

NAME

 freeze_displ_allmodes

FUNCTION

  From a given set of phonon modes, generate and output supercells and
  displaced configurations of atoms.
  Typically useful to follow soft modes and see distorsions of crystal structures

INPUTS

 amu(ntypat) = mass of the atoms (atomic mass unit)
 displ(2,3*natom,3*natom) = phonon mode displacements (complex)
 freeze_displ = amplitude of the displacement to freeze into the supercell
 natom = number of atoms in the unit cell
 ntypat = number of atom types
 phfrq(3*natom) = phonon frequencies
 qphnrm = norm of phonon q vector (should be 1 or 0)
 qphon = phonon wavevector
 rprimd(3,3) = dimensionfull primitive translations in real space
 typat(natom) = integer label of each type of atom (1,2,...)
 xcart(3,natom) = cartesian coords of atoms in unit cell (bohr)

OUTPUT

 for the moment only prints to file, but could also return pointer to supercell object, with
 rprimd and atomic positions, for further use

NOTES

 freeze_displ could be determined automatically from a temperature and the phonon frequency,
 as the average displacement of the mode with a Bose distribution.

PARENTS

      m_phonons

CHILDREN

      destroy_supercell,freeze_displ_supercell,init_supercell_for_qpt
      prt_supercell_for_qpt

SOURCE

3815 subroutine freeze_displ_allmodes(displ, freeze_displ, natom, outfile_radix, phfreq,  &
3816 &         qphon, rprimd, typat, xcart, znucl)
3817 
3818 
3819 !Do not modify the following lines by hand.
3820 #undef ABI_FUNC
3821 #define ABI_FUNC 'freeze_displ_allmodes'
3822 !End of the abilint section
3823 
3824 !This section has been created automatically by the script Abilint (TD).
3825 !Do not modify the following lines by hand.
3826 #undef ABI_FUNC
3827 #define ABI_FUNC 'freeze_displ_allmodes'
3828 !End of the abilint section
3829 
3830  implicit none
3831 
3832 ! arguments
3833 ! scalar
3834  integer,intent(in) :: natom
3835  character(len=*),intent(in) :: outfile_radix
3836  real(dp), intent(in) :: freeze_displ
3837 
3838 !arrays
3839  integer,intent(in) :: typat(natom)
3840 
3841  real(dp),intent(in) :: displ(2,3*natom,3*natom)
3842  real(dp),intent(in) :: rprimd(3,3)
3843  real(dp),intent(in) :: phfreq(3*natom)
3844  real(dp),intent(in) :: qphon(3)
3845  real(dp),intent(in) :: xcart(3,natom)
3846  real(dp),intent(in) :: znucl(:)
3847 
3848 ! local vars
3849  integer :: jmode
3850  type(supercell_type) :: scell
3851 
3852 ! *************************************************************************
3853 
3854 !determine supercell needed to freeze phonon
3855  call init_supercell_for_qpt(natom, qphon, rprimd, typat, xcart, znucl, scell)
3856 
3857  do jmode = 1, 3*natom
3858 ! reset positions
3859    scell%xcart = scell%xcart_ref
3860 
3861 !  displace atoms according to phonon jmode
3862    call freeze_displ_supercell(displ(:,:,jmode), freeze_displ, scell)
3863 
3864 !  print out everything for this wavevector and mode
3865    call prt_supercell_for_qpt (phfreq(jmode), jmode, outfile_radix, scell)
3866  end do
3867 
3868  call destroy_supercell (scell)
3869 
3870 end subroutine freeze_displ_allmodes

m_phonons/ifc_mkphbs [ Functions ]

[ Top ] [ m_phonons ] [ Functions ]

NAME

 ifc_mkphbs

FUNCTION

 Compute the phonon band structure from the IFC and write data to file(s)

INPUTS

 ifc<ifc_type>=Interatomic force constants
 cryst<crystal_t> = Info on the crystalline structure.
 dtset=<datasets_type>: input: all input variables initialized from the input file.
 prefix=Prefix for output files.
 comm=MPI communicator

OUTPUT

  Only writing.

PARENTS

      eph

CHILDREN

      ifc_fourq,kpath_free,phonons_ncwrite,phonons_write_gnuplot
      phonons_write_phfrq,phonons_write_xmgrace,xmpi_sum_master

SOURCE

3448 subroutine ifc_mkphbs(ifc, cryst, dtset, prefix, comm)
3449 
3450 
3451 !This section has been created automatically by the script Abilint (TD).
3452 !Do not modify the following lines by hand.
3453 #undef ABI_FUNC
3454 #define ABI_FUNC 'ifc_mkphbs'
3455 !End of the abilint section
3456 
3457  implicit none
3458 
3459 !Arguments -------------------------------
3460 !scalars
3461  integer,intent(in) :: comm
3462  character(len=*),intent(in) :: prefix
3463  type(ifc_type),intent(in) :: ifc
3464  type(crystal_t),intent(in) :: cryst
3465  type(dataset_type),intent(in) :: dtset
3466 
3467 !Local variables -------------------------
3468 !scalars
3469  integer,parameter :: master=0
3470  integer :: iqpt,nqpts,natom,ncid,nprocs,my_rank,ierr
3471  !character(500) :: msg
3472  type(kpath_t) :: qpath
3473 !arrays
3474  real(dp),allocatable :: eigvec(:,:,:,:,:),phfrqs(:,:),phdispl_cart(:,:,:,:),weights(:)
3475 
3476 ! *********************************************************************
3477 
3478  if (dtset%prtphbands == 0) return
3479 
3480  if (dtset%ph_nqpath <= 0 .or. dtset%ph_ndivsm <= 0) then
3481    MSG_WARNING("ph_nqpath <=0 or ph_ndivsm <= 0, phonon bands won't be produced. returning")
3482    return
3483  end if
3484 
3485  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
3486  natom = cryst%natom
3487 
3488  qpath = kpath_new(dtset%ph_qpath(:,1:dtset%ph_nqpath), cryst%gprimd, dtset%ph_ndivsm)
3489  nqpts = qpath%npts
3490 
3491  ABI_CALLOC(phfrqs, (3*natom,nqpts))
3492  ABI_CALLOC(phdispl_cart, (2,3*natom,3*natom,nqpts))
3493  ABI_CALLOC(eigvec, (2,3,natom,3,natom))
3494 
3495  do iqpt=1,nqpts
3496    if (mod(iqpt, nprocs) /= my_rank) cycle ! mpi-parallelism
3497    ! Get phonon frequencies and displacements in cartesian coordinates for this q-point
3498    call ifc_fourq(ifc, cryst, qpath%points(:,iqpt), phfrqs(:,iqpt), phdispl_cart(:,:,:,iqpt), out_eigvec=eigvec)
3499  end do
3500 
3501  call xmpi_sum_master(phfrqs, master, comm, ierr)
3502  call xmpi_sum_master(phdispl_cart, master, comm, ierr)
3503 
3504  if (my_rank == master) then
3505    ABI_MALLOC(weights, (nqpts))
3506    weights = one
3507 
3508 #ifdef HAVE_NETCDF
3509    NCF_CHECK_MSG(nctk_open_create(ncid, strcat(prefix, "_PHBST.nc"), xmpi_comm_self), "Creating PHBST")
3510    NCF_CHECK(crystal_ncwrite(cryst, ncid))
3511    call phonons_ncwrite(ncid,natom,nqpts, qpath%points, weights, phfrqs, phdispl_cart)
3512    NCF_CHECK(nctk_def_arrays(ncid, [nctkarr_t('atomic_mass_units', "dp", "number_of_atom_species")],defmode=.True.))
3513    NCF_CHECK(nctk_set_datamode(ncid))
3514    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'atomic_mass_units'), ifc%amu))
3515    NCF_CHECK(nf90_close(ncid))
3516 #endif
3517 
3518    select case (dtset%prtphbands)
3519    case (1)
3520      call phonons_write_xmgrace(strcat(prefix, "_PHBANDS.agr"), natom, nqpts, qpath%points, phfrqs, qptbounds=qpath%bounds)
3521    case (2)
3522      call phonons_write_gnuplot(prefix, natom, nqpts, qpath%points, phfrqs, qptbounds=qpath%bounds)
3523    case (3)
3524      call phonons_write_phfrq(strcat(prefix, "_PHFRQ"), natom, nqpts, qpath%points, weights, phfrqs, phdispl_cart)
3525    case default
3526      MSG_WARNING(sjoin("Unsupported value for prtphbands:", itoa(dtset%prtphbands)))
3527    end select
3528 
3529    ABI_FREE(weights)
3530  end if ! master
3531 
3532  ABI_FREE(phfrqs)
3533  ABI_FREE(phdispl_cart)
3534  ABI_FREE(eigvec)
3535 
3536  call kpath_free(qpath)
3537 
3538 end subroutine ifc_mkphbs

m_phonons/mkphbs [ Functions ]

[ Top ] [ m_phonons ] [ Functions ]

NAME

 mkphbs

FUNCTION

 Function to calculate the phonon band structure, from the IFC

INPUTS

 Ifc<ifc_type>=Interatomic force constants
 crystal<type(crystal_t)> = Info on the crystalline structure.
 inp= (derived datatype) contains all the input variables
 ddb<type(ddb_type)>=Object storing the DDB results.
 asrq0<asrq0_t>=Object for the treatment of the ASR based on the q=0 block found in the DDB file.
 prefix=Prefix for output files.
 dielt(3,3)=dielectric tensor
 comm=MPI communicator

OUTPUT

  Only writing.

PARENTS

      anaddb

CHILDREN

      ifc_fourq,kpath_free,phonons_ncwrite,phonons_write_gnuplot
      phonons_write_phfrq,phonons_write_xmgrace,xmpi_sum_master

SOURCE

1889 subroutine mkphbs(Ifc,Crystal,inp,ddb,asrq0,prefix,comm)
1890 
1891 
1892 !This section has been created automatically by the script Abilint (TD).
1893 !Do not modify the following lines by hand.
1894 #undef ABI_FUNC
1895 #define ABI_FUNC 'mkphbs'
1896 !End of the abilint section
1897 
1898  implicit none
1899 
1900 !Arguments -------------------------------
1901 !scalars
1902  integer,intent(in) :: comm
1903  character(len=*),intent(in) :: prefix
1904  type(ifc_type),intent(in) :: Ifc
1905  type(crystal_t),intent(in) :: Crystal
1906  type(anaddb_dataset_type),target,intent(in) :: inp
1907  type(ddb_type),intent(in) :: ddb
1908  type(asrq0_t),intent(inout) :: asrq0
1909 
1910 !Local variables -------------------------
1911 !scalars
1912  integer,parameter :: master=0
1913  integer :: unt
1914  integer :: iphl1,iblok,rftyp, ii,nfineqpath,nsym,natom,ncid,nprocs,my_rank
1915  integer :: natprj_bs,eivec,enunit,ifcflag
1916  real(dp) :: freeze_displ
1917  real(dp) :: cfact
1918  character(500) :: msg
1919  character(len=8) :: unitname
1920 
1921 !arrays
1922  integer :: rfphon(4),rfelfd(4),rfstrs(4)
1923  integer :: nomega, imode, iomega
1924  integer,allocatable :: ndiv(:)
1925  real(dp) :: omega, omega_min, gaussmaxarg, gaussfactor, gaussprefactor, xx
1926  real(dp) :: speedofsound(3)
1927  real(dp) :: qphnrm(3), qphon(3), qphon_padded(3,3),res(3)
1928  real(dp) :: d2cart(2,ddb%msize),real_qphon(3)
1929  real(dp) :: displ(2*3*crystal%natom*3*crystal%natom),eigval(3,crystal%natom)
1930  real(dp),allocatable :: phfrq(:),eigvec(:,:,:,:,:)
1931  real(dp),allocatable :: save_phfrq(:,:),save_phdispl_cart(:,:,:,:),save_qpoints(:,:)
1932  real(dp),allocatable :: weights(:)
1933  real(dp),allocatable :: dos4bs(:)
1934  real(dp),allocatable,target :: alloc_path(:,:)
1935  real(dp),pointer :: fineqpath(:,:)
1936  type(atprj_type) :: atprj
1937 
1938 ! *********************************************************************
1939 
1940  ! Only master works for the time being
1941  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
1942  if (my_rank /= master) return
1943 
1944  nsym = Crystal%nsym; natom = Crystal%natom
1945 
1946  ! Copy parameters from inp (then I will try to remove inp from the API so that I can call mkphbs in eph)
1947  ifcflag = inp%ifcflag
1948  natprj_bs = inp%natprj_bs
1949  freeze_displ = inp%freeze_displ
1950  eivec = inp%eivec; enunit = inp%enunit
1951  rftyp=inp%rfmeth
1952 
1953  nullify(fineqpath)
1954  nfineqpath = inp%nph1l
1955  fineqpath => inp%qph1l
1956 
1957  if(inp%nph1l==0) then
1958    if (inp%nqpath==0) then
1959      return ! if there is nothing to do, return
1960    else
1961      ! allow override of nph1l with nqpath if the former is not set
1962      ! allocate and compute path here and make fineqpath points to it
1963      ABI_MALLOC(ndiv,(inp%nqpath-1))
1964      call make_path(inp%nqpath,inp%qpath,Crystal%gmet,'G',inp%ndivsm,ndiv,nfineqpath,alloc_path,std_out)
1965      ABI_FREE(ndiv)
1966      fineqpath => alloc_path
1967    end if
1968  end if
1969 
1970  write(msg, '(a,(80a),a,a,a,a)' ) ch10,('=',ii=1,80),ch10,ch10,' Treat the first list of vectors ',ch10
1971  call wrtout(std_out,msg,'COLL')
1972  call wrtout(ab_out,msg,'COLL')
1973 
1974  if (natprj_bs > 0) call atprj_init(atprj, natom, natprj_bs, inp%iatprj_bs, prefix)
1975 
1976  ABI_MALLOC(phfrq,(3*natom))
1977  ABI_MALLOC(eigvec,(2,3,natom,3,natom))
1978  ABI_MALLOC(save_qpoints,(3,nfineqpath))
1979  ABI_MALLOC(save_phfrq,(3*natom,nfineqpath))
1980  ABI_MALLOC(save_phdispl_cart,(2,3*natom,3*natom,nfineqpath))
1981  qphnrm = one
1982 
1983  do iphl1=1,nfineqpath
1984 
1985    ! Initialisation of the phonon wavevector
1986    qphon(:)=fineqpath(:,iphl1)
1987 
1988    if (inp%nph1l /= 0) qphnrm(1) = inp%qnrml1(iphl1)
1989 
1990    save_qpoints(:,iphl1) = qphon / qphnrm(1)
1991 
1992    ! Generation of the dynamical matrix in cartesian coordinates
1993    if (ifcflag == 1) then
1994 
1995      ! Get phonon frequencies and displacements in reduced coordinates for this q-point
1996      !call ifc_fourq(ifc, cryst, save_qpoints(:,iphl1), phfrq, displ, out_eigvec=eigvec)
1997 
1998      ! Get d2cart using the interatomic forces and the
1999      ! long-range coulomb interaction through Ewald summation
2000      call gtdyn9(ddb%acell,Ifc%atmfrc,Ifc%dielt,Ifc%dipdip,Ifc%dyewq0,d2cart,Crystal%gmet,ddb%gprim,ddb%mpert,natom,&
2001 &     Ifc%nrpt,qphnrm(1),qphon,Crystal%rmet,ddb%rprim,Ifc%rpt,Ifc%trans,Crystal%ucvol,Ifc%wghatm,Crystal%xred,ifc%zeff)
2002 
2003    else if (ifcflag == 0) then
2004 
2005      !call ddb_diagoq(ddb, crystal, save_qpoints(:,iphl1), asrq0, ifc%symdynmat, rftyp, phfrq, displ, &
2006      !                out_eigvec=eigvec)
2007 
2008      ! Look for the information in the DDB (no interpolation here!)
2009      rfphon(1:2)=1; rfelfd(1:2)=0; rfstrs(1:2)=0
2010      qphon_padded = zero; qphon_padded(:,1) = qphon
2011 
2012      call gtblk9(ddb,iblok,qphon_padded,qphnrm,rfphon,rfelfd,rfstrs,rftyp)
2013 
2014      ! Copy the dynamical matrix in d2cart
2015      d2cart(:,1:ddb%msize)=ddb%val(:,:,iblok)
2016 
2017      ! Eventually impose the acoustic sum rule based on previously calculated d2asr
2018      call asrq0_apply(asrq0, natom, ddb%mpert, ddb%msize, crystal%xcart, d2cart)
2019    end if
2020 
2021    ! Use inp%symdynmat instead of ifc because of ifcflag
2022    ! Calculation of the eigenvectors and eigenvalues of the dynamical matrix
2023    call dfpt_phfrq(ddb%amu,displ,d2cart,eigval,eigvec,Crystal%indsym,&
2024 &   ddb%mpert,Crystal%nsym,natom,nsym,Crystal%ntypat,phfrq,qphnrm(1),qphon,&
2025 &   crystal%rprimd,inp%symdynmat,Crystal%symrel,Crystal%symafm,Crystal%typat,Crystal%ucvol)
2026 
2027    if (abs(freeze_displ) > tol10) then
2028      real_qphon = zero
2029      if (abs(qphnrm(1)) > tol8) real_qphon = qphon / qphnrm(1)
2030      call freeze_displ_allmodes(displ, freeze_displ, natom, prefix, phfrq, &
2031 &     real_qphon, crystal%rprimd, Crystal%typat, crystal%xcart, crystal%znucl)
2032    end if
2033 
2034    ! If requested, output projection of each mode on given atoms
2035    if (natprj_bs > 0) call atprj_print(atprj, iphl1, phfrq, eigvec)
2036 
2037    ! In case eivec == 4, write output files for band2eps (visualization of phonon band structures)
2038    if (eivec == 4) then
2039      call sortph(eigvec,displ,strcat(prefix, "_B2EPS"),natom,phfrq)
2040    end if
2041 
2042    ! Write the phonon frequencies
2043    call dfpt_prtph(displ,eivec,enunit,ab_out,natom,phfrq,qphnrm(1),qphon)
2044 
2045    save_phfrq(:,iphl1) = phfrq
2046    save_phdispl_cart(:,:,:,iphl1) = RESHAPE(displ, [2, 3*natom, 3*natom])
2047 
2048    ! Determine the symmetries of the phonon mode at Gamma
2049    ! TODO: generalize for other q-point little groups.
2050    if (sum(abs(qphon)) < DDB_QTOL) then
2051      call dfpt_symph(ab_out,ddb%acell,eigvec,Crystal%indsym,natom,nsym,phfrq,ddb%rprim,Crystal%symrel)
2052    end if
2053 
2054    ! if we have an acoustic mode (small q and acoustic type displacements)
2055    ! extrapolate speed of sound in this direction, and Debye frequency
2056    call wrap2_pmhalf(qphon, real_qphon, res)
2057    if (sqrt(real_qphon(1)**2+real_qphon(2)**2+real_qphon(3)**2) < quarter .and. &
2058 &   sqrt(real_qphon(1)**2+real_qphon(2)**2+real_qphon(3)**2) > tol6) then
2059      call phdos_calc_vsound(eigvec, Crystal%gmet, natom, phfrq, real_qphon, speedofsound)
2060      if (my_rank == master) call phdos_print_vsound(ab_out, Crystal%ucvol, speedofsound)
2061    end if
2062 
2063  end do ! iphl1
2064 
2065 ! calculate dos for the specific q points along the BS calculated
2066 ! only Gaussians are possible - no interpolation
2067  omega_min = minval(save_phfrq(:,:))
2068  nomega=NINT( (maxval(save_phfrq(:,:))-omega_min) / inp%dosdeltae ) + 1
2069  nomega=MAX(6,nomega) ! Ensure Simpson integration will be ok
2070 
2071  ABI_MALLOC(dos4bs,(nomega))
2072  dos4bs = zero
2073  gaussmaxarg = sqrt(-log(1.d-90))
2074  gaussprefactor = one/(inp%dossmear*sqrt(two_pi))
2075  gaussfactor    = one/(sqrt2*inp%dossmear)
2076  do iphl1=1,nfineqpath
2077    do imode=1,3*natom
2078      do iomega=1, nomega
2079        omega = omega_min + (iomega-1) * inp%dosdeltae
2080        xx = (omega - save_phfrq(imode,iphl1)) * gaussfactor
2081        if(abs(xx) < gaussmaxarg) then
2082          dos4bs(iomega) = dos4bs(iomega) + gaussprefactor*exp(-xx*xx)
2083        end if
2084      end do
2085    end do
2086  end do
2087 
2088 
2089 !deallocate sortph array
2090  call end_sortph()
2091 
2092  if (natprj_bs > 0) call atprj_destroy(atprj)
2093 
2094 
2095 ! WRITE OUT FILES
2096  if (my_rank == master) then
2097    ABI_MALLOC(weights, (nfineqpath))
2098    weights = one
2099 
2100 #ifdef HAVE_NETCDF
2101    NCF_CHECK_MSG(nctk_open_create(ncid, strcat(prefix, "_PHBST.nc"), xmpi_comm_self), "Creating PHBST")
2102    NCF_CHECK(crystal_ncwrite(Crystal, ncid))
2103    call phonons_ncwrite(ncid,natom,nfineqpath,save_qpoints,weights,save_phfrq,save_phdispl_cart)
2104 
2105    ! Now treat the second list of vectors (only at the Gamma point, but can include non-analyticities)
2106    if (inp%nph2l /= 0 .and. inp%ifcflag == 1) then
2107      call ifc_calcnwrite_nana_terms(ifc, crystal, inp%nph2l, inp%qph2l, inp%qnrml2, ncid)
2108    end if
2109 
2110    NCF_CHECK(nf90_close(ncid))
2111 #endif
2112 
2113    call phonons_write_phfrq(strcat(prefix, "_PHFRQ"), natom,nfineqpath,save_qpoints,weights,save_phfrq,save_phdispl_cart)
2114 
2115    select case (inp%prtphbands)
2116    case (0)
2117      continue
2118 
2119    case (1)
2120      if (inp%nph1l == 0) then
2121        call phonons_write_xmgrace(strcat(prefix, "_PHBANDS.agr"), natom, nfineqpath, save_qpoints, save_phfrq, &
2122           qptbounds=inp%qpath)
2123      else
2124        call phonons_write_xmgrace(strcat(prefix, "_PHBANDS.agr"), natom, nfineqpath, save_qpoints, save_phfrq)
2125      end if
2126 
2127    case (2)
2128      if (inp%nph1l == 0) then
2129        call phonons_write_gnuplot(prefix, natom, nfineqpath, save_qpoints, save_phfrq, qptbounds=inp%qpath)
2130      else
2131        call phonons_write_gnuplot(prefix, natom, nfineqpath, save_qpoints, save_phfrq)
2132      end if
2133 
2134    !case (3)
2135      !call phonons_writeEPS(natom,nfineqpath,Crystal%ntypat,Crystal%typat, &
2136      !  save_phfrq,save_phdispl_cart)
2137 
2138    case default
2139      MSG_WARNING(sjoin("Don't know how to handle prtphbands:", itoa(inp%prtphbands)))
2140    end select
2141 
2142    ! write out DOS file for q along this path
2143    cfact=one
2144    unitname = 'Ha'
2145    if (open_file('PHBST_partial_DOS',msg,newunit=unt,form="formatted",action="write") /= 0) then
2146      MSG_ERROR(msg)
2147    end if
2148    write(msg,'(3a)')'# ',ch10,'# Partial phonon density of states for q along a band structure path'
2149    call wrtout(unt,msg,'COLL')
2150    write(msg,'(6a)')'# ',ch10,'# Energy in ',unitname,', DOS in states/',unitname
2151    call wrtout(unt,msg,'COLL')
2152    write(msg,'(a,E20.10,2a,i8)') '# Gaussian method with smearing = ',inp%dossmear*cfact,unitname, &
2153 &        ', nq =', nfineqpath
2154    call wrtout(unt,msg,'COLL')
2155    write(msg,'(5a)')'# ',ch10,'# omega     PHDOS ',ch10,'# '
2156    call wrtout(unt,msg,'COLL')
2157    do iomega=1,nomega
2158      omega = omega_min + (iomega-1) * inp%dosdeltae
2159      write(unt,'(2es17.8)',advance='NO')omega*cfact,dos4bs(iomega)/cfact
2160      write(unt,*)
2161    end do
2162    close(unt)
2163 
2164 
2165    ABI_FREE(weights)
2166  end if
2167 
2168  ABI_FREE(save_qpoints)
2169  ABI_FREE(save_phfrq)
2170  ABI_FREE(save_phdispl_cart)
2171  ABI_FREE(phfrq)
2172  ABI_FREE(eigvec)
2173  ABI_FREE(dos4bs)
2174 
2175  if (allocated(alloc_path)) then
2176    ABI_FREE(alloc_path)
2177  end if
2178 
2179 end subroutine mkphbs

m_phonons/mkphdos [ Functions ]

[ Top ] [ m_phonons ] [ Functions ]

NAME

 mkphdos

FUNCTION

 Calculate the phonon density of states as well as
 the contributions associated to the different types of atoms in the unit cell.
 Two methods are implemented: gaussian method and linear interpolation based on tetrahedrons.

INPUTS

 ifc<ifc_type>=Interatomic force constants
 crystal<crystal_t>=Info on the crystalline structure.
 prtdos=1 for gaussian method, 2 for tetrahedra.
 dosdeltae=Step of frequency mesh.
 dossmear=Gaussian broadening, used if prtdos==1.
 dos_ngqpt(3)=Divisions of the q-mesh used for computing the DOS
 nqshift=Number of shifts in Q-mesh
 dos_qshift(3, nqshift)=Shift of the q-mesh.
 comm=MPI communicator.

OUTPUT

 phdos<phonon_dos_type>=Container with phonon DOS, IDOS and atom-projected DOS.
 count_wminmax(2)=Number of (interpolated) phonon frequencies that are outside
   input range (see wminmax). Client code can use count_wminmax and wminmax to
   enlarge the mesh and call the routine again to recompute the DOS
   if all frequencies should be included.

SIDE EFFECTS

 wminmax(2)=
   In input: min and max value of frequency mesh. Used only if minmax(2) > minmax(1)
    else values are taken from ifc%omega_minmax (computed from ab-initio mesh + pad)
   In output: min and max frequency obtained after interpolating the IFCs on the dense q-mesh dos_ngqpt

PARENTS

      anaddb,eph,m_tdep_phdos

CHILDREN

      ifc_fourq,kpath_free,phonons_ncwrite,phonons_write_gnuplot
      phonons_write_phfrq,phonons_write_xmgrace,xmpi_sum_master

SOURCE

 693 subroutine mkphdos(phdos, crystal, ifc, prtdos, dosdeltae, dossmear, dos_ngqpt, nqshft, dos_qshift, &
 694                    wminmax, count_wminmax, comm)
 695 
 696 
 697 !This section has been created automatically by the script Abilint (TD).
 698 !Do not modify the following lines by hand.
 699 #undef ABI_FUNC
 700 #define ABI_FUNC 'mkphdos'
 701 !End of the abilint section
 702 
 703  implicit none
 704 
 705 !Arguments -------------------------------
 706 !scalars
 707  integer,intent(in) :: prtdos,nqshft,comm
 708  real(dp),intent(in) :: dosdeltae,dossmear
 709  type(crystal_t),intent(in) :: crystal
 710  type(ifc_type),intent(in) :: ifc
 711  type(phonon_dos_type),intent(out) :: phdos
 712 !arrays
 713  integer,intent(in) :: dos_ngqpt(3)
 714  integer,intent(out) :: count_wminmax(2)
 715  real(dp),intent(in) :: dos_qshift(3,nqshft)
 716  real(dp),intent(inout) :: wminmax(2)
 717 
 718 !Local variables -------------------------
 719 !scalars
 720  integer,parameter :: brav1=1,chksymbreak0=0,bcorr0=0,qptopt1=1,master=0
 721  integer :: iat,jat,idir,imode,io,iq_ibz,itype,nkpt_fullbz
 722  integer :: nqbz,ierr,natom,nomega,jdir, isym, nprocs, my_rank, ncid
 723  real(dp),parameter :: max_occ1=one, gaussmaxarg = sqrt(-log(1.d-90)), max_smallq = 0.0625_dp
 724  real(dp) :: nsmallq,gaussfactor,gaussprefactor,normq,debyefreq,rtmp
 725  real(dp) :: cpu, wall, gflops
 726  character(len=500) :: msg
 727  character(len=80) :: errstr
 728  character(len=80) ::  prefix = "freq_displ" ! FIXME
 729  type(t_tetrahedron) :: tetraq
 730 !arrays
 731  integer :: in_qptrlatt(3,3),new_qptrlatt(3,3)
 732  integer,allocatable :: bz2ibz(:)
 733  real(dp) :: speedofsound(3),speedofsound_(3)
 734  real(dp) :: displ(2*3*Crystal%natom*3*Crystal%natom)
 735  real(dp) :: eigvec(2,3,Crystal%natom,3*Crystal%natom),phfrq(3*Crystal%natom)
 736  real(dp) :: qlatt(3,3),rlatt(3,3)
 737  real(dp) :: msqd_atom_tmp(3,3),temp_33(3,3)
 738  real(dp) :: symcart(3,3,crystal%nsym)
 739  real(dp) :: syme2_xyza(3, crystal%natom)
 740  real(dp),allocatable :: dtweightde(:,:),full_eigvec(:,:,:,:,:),full_phfrq(:,:),new_shiftq(:,:)
 741  real(dp),allocatable :: kpt_fullbz(:,:),qbz(:,:),qibz(:,:),tmp_phfrq(:),tweight(:,:)
 742  real(dp),allocatable :: wtq_ibz(:),xvals(:), gvals_wtq(:), wdt(:,:)
 743 
 744 ! *********************************************************************
 745 
 746  DBG_ENTER("COLL")
 747 
 748  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
 749 
 750  ! Consistency check.
 751  if (all(prtdos /= [1, 2])) then
 752    MSG_BUG(sjoin('prtdos should be 1 or 2, but received', itoa(prtdos)))
 753  end if
 754  if (dosdeltae <= zero) then
 755    MSG_BUG(sjoin('dosdeltae should be positive, but received', ftoa(dosdeltae)))
 756  end if
 757  if (prtdos == 1 .and. dossmear <= zero) then
 758    MSG_BUG(sjoin('dossmear should be positive but received', ftoa(dossmear)))
 759  end if
 760 
 761  call cwtime(cpu, wall, gflops, "start")
 762 
 763  ! Get symmetries in cartesian coordinates
 764  do isym = 1, crystal%nsym
 765    call symredcart(crystal%rprimd,crystal%gprimd,symcart(:,:,isym),crystal%symrel(:,:,isym))
 766  end do
 767 
 768  natom = crystal%natom
 769  phdos%ntypat     = crystal%ntypat
 770  phdos%natom      = crystal%natom
 771  phdos%prtdos     = prtdos
 772  phdos%dossmear   = dossmear
 773  phdos%omega_step = dosdeltae
 774  ! Use values stored in ifc (obtained with ab-initio q-mesh + pad)
 775  if (wminmax(2) > wminmax(1)) then
 776    phdos%omega_min = wminmax(1)
 777    phdos%omega_max = wminmax(2)
 778  else
 779    phdos%omega_min = ifc%omega_minmax(1)
 780    phdos%omega_max = ifc%omega_minmax(2)
 781  end if
 782  ! Must be consistent with mesh computed in tetra routines!
 783  phdos%nomega = nint((phdos%omega_max - phdos%omega_min) / phdos%omega_step) + 1
 784  ! Ensure Simpson integration will be ok
 785  phdos%nomega = max(6, phdos%nomega)
 786  nomega = phdos%nomega
 787 
 788  ! Build frequency mesh.
 789  ABI_MALLOC(phdos%omega, (phdos%nomega))
 790  do io=1,phdos%nomega
 791    phdos%omega(io) = phdos%omega_min + phdos%omega_step * (io - 1)
 792  end do
 793  phdos%omega_min = phdos%omega(1)
 794  phdos%omega_max = phdos%omega(phdos%nomega)
 795 
 796  ! Allocate arrays that depend on nomega and set them to zero.
 797  ABI_CALLOC(phdos%phdos, (nomega))
 798  ABI_CALLOC(phdos%phdos_int, (nomega))
 799  ABI_CALLOC(phdos%pjdos, (nomega, 3, natom))
 800  ABI_CALLOC(phdos%pjdos_int, (nomega, 3, natom))
 801  ABI_CALLOC(phdos%msqd_dos_atom, (nomega, 3, 3, natom))
 802  ABI_MALLOC(phdos%atom_mass, (natom))
 803  phdos%atom_mass = crystal%amu(crystal%typat(:)) * amu_emass
 804 
 805  ABI_MALLOC(gvals_wtq, (nomega))
 806  ABI_MALLOC(xvals, (nomega))
 807 
 808  ! Parameters defining the gaussian approximant.
 809  if (prtdos == 1) then
 810    ! TODO: use dirac_delta and update reference files.
 811    gaussprefactor = one / (dossmear * sqrt(two_pi))
 812    gaussfactor = one / (sqrt2 * dossmear)
 813    write(msg, '(4a,f8.5,2a,f8.5)') ch10, &
 814     ' mkphdos: calculating phonon DOS using gaussian method:',ch10, &
 815     '    gaussian smearing [meV] = ',dossmear*Ha_meV,ch10, &
 816     '    frequency step    [meV] = ',phdos%omega_step*Ha_meV
 817  else if (prtdos == 2) then
 818    write(msg,'(2a)')ch10,' mkphdos: calculating phonon DOS using tetrahedron method'
 819  end if
 820  call wrtout(std_out, msg, 'COLL')
 821 
 822  ! TODO
 823  ! 1) fix bug in tetra if degenerate and update ref files
 824  ! 2) nshift > 1?
 825 
 826  ! This call will set %nqibz and IBZ and BZ arrays
 827  in_qptrlatt = 0; in_qptrlatt(1, 1) = dos_ngqpt(1); in_qptrlatt(2, 2) = dos_ngqpt(2); in_qptrlatt(3, 3) = dos_ngqpt(3)
 828 
 829  call kpts_ibz_from_kptrlatt(crystal, in_qptrlatt, qptopt1, nqshft, dos_qshift, &
 830    phdos%nqibz, qibz, wtq_ibz, nqbz, qbz, new_kptrlatt=new_qptrlatt, new_shiftk=new_shiftq)
 831 
 832  if (prtdos == 2) then
 833    ! Prepare tetrahedron method including workspace arrays.
 834    ! Convert kptrlatt to double and invert, qlatt here refer to the shortest qpt vectors
 835    rlatt = new_qptrlatt; call matr3inv(rlatt, qlatt)
 836 
 837    nkpt_fullbz = nqbz
 838    ABI_MALLOC(bz2ibz, (nkpt_fullbz))
 839    ABI_MALLOC(kpt_fullbz, (3, nkpt_fullbz))
 840 
 841    ! Make full kpoint grid and get equivalence to irred kpoints.
 842    call get_full_kgrid(bz2ibz, qibz, kpt_fullbz, new_qptrlatt, phdos%nqibz, &
 843        nkpt_fullbz, size(new_shiftq, dim=2), crystal%nsym, new_shiftq, crystal%symrel)
 844 
 845    ! Init tetrahedra, i.e. indexes of the full q-points at their summits
 846    call init_tetra(bz2ibz, crystal%gprimd, qlatt, kpt_fullbz, nqbz, tetraq, ierr, errstr)
 847    ABI_CHECK(ierr == 0, errstr)
 848 
 849    ABI_FREE(bz2ibz)
 850    ABI_FREE(kpt_fullbz)
 851 
 852    ! Allocate arrays used to store the entire spectrum, Required to calculate tetra weights.
 853    ! this may change in the future if Matteo refactorizes the tetra weights as sums over k instead of sums over bands
 854    ABI_CALLOC(full_phfrq, (3*natom, phdos%nqibz))
 855    ABI_STAT_MALLOC(full_eigvec, (2, 3, natom, 3*natom, phdos%nqibz), ierr)
 856    ABI_CHECK(ierr == 0, 'out-of-memory in full_eigvec')
 857    full_eigvec = zero
 858  end if ! tetra
 859  ABI_FREE(new_shiftq)
 860 
 861  ! MPI Sum over irreducible q-points then sync the following integrals:
 862  !   speedofsound, nsmallq
 863  !   wminmax and count_wminmax
 864  !   if gauss: %phdos, %msqd_dos_atom
 865  !   if tetra: full_phfrq, full_eigvec, %phdos_int
 866 
 867  nsmallq = zero; speedofsound = zero
 868  wminmax = [huge(one), -huge(one)]; count_wminmax = 0
 869  do iq_ibz=1,phdos%nqibz
 870    if (mod(iq_ibz, nprocs) /= my_rank) cycle ! mpi-parallelism
 871 
 872    ! Fourier interpolation (keep track of min/max to decide if initial mesh was large enough)
 873    call ifc_fourq(Ifc, crystal, qibz(:,iq_ibz), phfrq, displ, out_eigvec=eigvec)
 874    wminmax(1) = min(wminmax(1), minval(phfrq))
 875    if (wminmax(1) < phdos%omega(1)) count_wminmax(1) = count_wminmax(1) + 1
 876    wminmax(2) = max(wminmax(2), maxval(phfrq))
 877    if (wminmax(2) > phdos%omega(nomega)) count_wminmax(2) = count_wminmax(2) + 1
 878 
 879    normq = sum(qibz(:,iq_ibz) ** 2)
 880    if (normq < max_smallq .and. normq > tol6) then
 881      call phdos_calc_vsound(eigvec, crystal%gmet, natom, phfrq, qibz(:,iq_ibz), speedofsound_)
 882      speedofsound = speedofsound + speedofsound_ * wtq_ibz(iq_ibz)
 883      nsmallq = nsmallq + wtq_ibz(iq_ibz)
 884    end if
 885 
 886    select case (prtdos)
 887    case (1)
 888      do imode=1,3*natom
 889        ! Precompute \delta(w - w_{qnu}) * weight(q)
 890        xvals = (phdos%omega(:) - phfrq(imode)) * gaussfactor
 891        where (abs(xvals) < gaussmaxarg)
 892          gvals_wtq = gaussprefactor * exp(-xvals*xvals) * wtq_ibz(iq_ibz)
 893        elsewhere
 894          gvals_wtq = zero
 895        end where
 896 
 897        ! Accumulate PHDOS
 898        phdos%phdos(:) = phdos%phdos(:) + gvals_wtq
 899 
 900        ! Rotate e(q) to get e(Sq) to account for symmetrical q-points in BZ.
 901        ! eigenvectors indeed are not invariant under rotation. See e.g. Eq 39-40 of PhysRevB.76.165108 [[cite:Giustino2007]].
 902        ! In principle there's a phase due to nonsymmorphic translations
 903        ! but we here need |e(Sq)_iatom|**2
 904        syme2_xyza = zero
 905        do iat=1,natom
 906          do isym=1, crystal%nsym
 907            jat = crystal%indsym(4,isym,iat)
 908            syme2_xyza(:,jat) = syme2_xyza(:,jat) + &
 909              matmul(symcart(:,:,isym), eigvec(1,:,iat,imode)) ** 2 + &
 910              matmul(symcart(:,:,isym), eigvec(2,:,iat,imode)) ** 2
 911          end do
 912        end do
 913        !syme2_xyza = syme2_xyza / crystal%nsym
 914 
 915        ! Accumulate PJDOS
 916        do iat=1,natom
 917          do idir=1,3
 918            phdos%pjdos(:,idir,iat) = phdos%pjdos(:,idir,iat) + syme2_xyza(idir,iat) * gvals_wtq
 919          end do
 920        end do
 921 
 922        ! Accumulate outer product of displacement vectors
 923        ! NB: only accumulate real part. e(-q) = e(q)* thue full sum over the BZ guarantees Im=0
 924        ! this sum only does irreducible points: the matrix is symmetrized below
 925        ! msqd_atom_tmp has units of bohr^2 / Ha as gaussval ~ 1/smear ~ 1/Ha
 926        do iat=1,natom
 927          msqd_atom_tmp = zero
 928          do idir=1,3
 929            do jdir=1,3
 930              msqd_atom_tmp(jdir,idir) = msqd_atom_tmp(jdir,idir) + ( &
 931                    eigvec(1,idir,iat,imode)* eigvec(1,jdir,iat,imode) &
 932                 +  eigvec(2,idir,iat,imode)* eigvec(2,jdir,iat,imode) )
 933            end do
 934          end do
 935 
 936          ! Symmetrize matrices to get full sum of tensor over all BZ, not just IBZ.
 937          ! the atom is not necessarily invariant under symops, so these contributions should be added to each iat separately
 938          ! normalization by nsym is done at the end outside the iqpt loop and after the tetrahedron clause
 939          ! NB: looks consistent with the sym in harmonic thermo, just used in opposite
 940          ! direction for symops: symrel here instead of symrec and the inverse of
 941          ! indsym in harmonic_thermo
 942          do isym=1, crystal%nsym
 943            temp_33 = matmul( (symcart(:,:,isym)), matmul(msqd_atom_tmp, transpose(symcart(:,:,isym))) )
 944            jat = crystal%indsym(4,isym,iat)
 945            do idir=1,3
 946              do jdir=1,3
 947                phdos%msqd_dos_atom(:,idir,jdir,jat) = phdos%msqd_dos_atom(:,idir,jdir,jat) + &
 948                  temp_33(idir, jdir) * gvals_wtq
 949              end do
 950            end do
 951          end do
 952 
 953        end do ! iat
 954      end do ! imode
 955 
 956    case (2)
 957      ! Tetrahedra
 958      ! Save phonon frequencies and eigenvectors.
 959      ! Sum is done after the loops over the two meshes.
 960      full_phfrq(:,iq_ibz) = phfrq(:)
 961      full_eigvec(:,:,:,:,iq_ibz) = eigvec
 962 
 963    case default
 964      MSG_ERROR(sjoin("Wrong value for prtdos:", itoa(prtdos)))
 965    end select
 966  end do ! iq_ibz
 967 
 968  ABI_FREE(qbz)
 969  ABI_FREE(gvals_wtq)
 970  ABI_FREE(xvals)
 971 
 972  call xmpi_sum_master(nsmallq, master, comm, ierr)
 973  call xmpi_sum_master(speedofsound, master, comm, ierr)
 974 
 975  if (my_rank == master .and. nsmallq > tol10) then
 976    ! Write info about speed of sound
 977    speedofsound = speedofsound / nsmallq
 978    write (msg,'(a,E20.10,3a,F16.4,2a)') &
 979        ' Average speed of sound partial sums: ', third*sum(speedofsound), ' (at units)',ch10, &
 980        '-                                   = ', third*sum(speedofsound) * Bohr_Ang * 1.d-13 / Time_Sec, ' [km/s]',ch10
 981    call wrtout (ab_out,msg,"COLL")
 982    call wrtout (std_out,msg,"COLL")
 983 
 984    ! Debye frequency = vs * (6 pi^2 natom / ucvol)**1/3
 985    debyefreq = third*sum(speedofsound) * (six*pi**2/crystal%ucvol)**(1./3.)
 986    write (msg,'(a,E20.10,3a,E20.10,a)') &
 987       ' Debye frequency from partial sums: ', debyefreq, ' (Ha)',ch10, &
 988       '-                                 = ', debyefreq*Ha_THz, ' (THz)'
 989    call wrtout (ab_out,msg,"COLL")
 990    call wrtout (std_out,msg,"COLL")
 991 
 992    ! Debye temperature = hbar * Debye frequency / kb
 993    write (msg,'(a,E20.10,2a)') '-Debye temperature from partial sums: ', debyefreq*Ha_K, ' (K)', ch10
 994    call wrtout (ab_out,msg,"COLL")
 995    call wrtout (std_out,msg,"COLL")
 996  end if
 997 
 998  if (prtdos == 2) then
 999    ! Finalize integration with tetrahedra
1000    ! All the data are contained in full_phfrq and full_eigvec.
1001    call xmpi_sum(full_phfrq, comm, ierr)
1002    call xmpi_sum(full_eigvec, comm, ierr)
1003 
1004    ABI_MALLOC(tmp_phfrq, (phdos%nqibz))
1005    ABI_MALLOC(wdt, (nomega, 2))
1006    ABI_MALLOC(tweight, (nomega, phdos%nqibz))
1007    ABI_MALLOC(dtweightde, (nomega, phdos%nqibz))
1008 
1009    do imode=1,3*natom
1010      tmp_phfrq(:) = full_phfrq(imode,:)
1011 
1012      ! Parallelize weights computation inside comm, then distribute nqibz
1013      call tetra_blochl_weights(tetraq,tmp_phfrq,phdos%omega_min,phdos%omega_max,max_occ1,phdos%nomega,&
1014         phdos%nqibz,bcorr0,tweight,dtweightde,comm)
1015 
1016      !if (mod(imode, nprocs) /= my_rank) cycle ! mpi-parallelism
1017      !call tetra_blochl_weights(tetraq,tmp_phfrq,phdos%omega_min,phdos%omega_max,max_occ1,phdos%nomega,&
1018      !   phdos%nqibz,bcorr0,tweight,dtweightde,xmpi_comm_self)
1019 
1020      do iq_ibz=1,phdos%nqibz
1021        if (mod(iq_ibz, nprocs) /= my_rank) cycle ! mpi-parallelism
1022        wdt(:, 1) = dtweightde(:, iq_ibz)
1023        wdt(:, 2) = tweight(:, iq_ibz)
1024 
1025        !call tetra_get_onewk(tetraq, iq_ibz, bcorr0, phdos%nomega, phdos%nqibz, &
1026        !  tmp_phfrq, phdos%omega_min, phdos%omega_max, max_occ, wdt)
1027 
1028        ! Accumulate DOS/IDOS
1029        phdos%phdos(:) = phdos%phdos(:) + wdt(:, 1)
1030        phdos%phdos_int(:) = phdos%phdos_int(:) + wdt(:, 2)
1031 
1032        ! Rotate e(q) to get e(Sq) to account for other q-points in BZ. See notes in gaussian branch
1033        syme2_xyza = zero
1034        do iat=1,natom
1035          do isym=1, crystal%nsym
1036            jat = crystal%indsym(4,isym,iat)
1037            syme2_xyza(:,jat) = syme2_xyza(:,jat) + &
1038              matmul(symcart(:,:,isym), full_eigvec(1,:,iat,imode,iq_ibz)) ** 2 + &
1039              matmul(symcart(:,:,isym), full_eigvec(2,:,iat,imode,iq_ibz)) ** 2
1040          end do
1041        end do
1042        !syme2_xyza = syme2_xyza / crystal%nsym
1043 
1044        do iat=1,natom
1045          do idir=1,3
1046            phdos%pjdos(:,idir,iat) = phdos%pjdos(:,idir,iat) + syme2_xyza(idir,iat) * wdt(:,1)
1047            phdos%pjdos_int(:,idir,iat) = phdos%pjdos_int(:,idir,iat) + syme2_xyza(idir,iat) * wdt(:,2)
1048          end do
1049        end do
1050 
1051        do iat=1,natom
1052          ! Accumulate outer product of displacement vectors
1053          msqd_atom_tmp = zero
1054          do idir=1,3
1055            do jdir=1,3
1056              msqd_atom_tmp(jdir,idir) = msqd_atom_tmp(jdir,idir) + ( &
1057                    full_eigvec(1,idir,iat,imode,iq_ibz)* full_eigvec(1,jdir,iat,imode,iq_ibz) &
1058                 +  full_eigvec(2,idir,iat,imode,iq_ibz)* full_eigvec(2,jdir,iat,imode,iq_ibz) ) !* gvals_wtq
1059            end do ! jdie
1060          end do
1061 
1062          ! Symmetrize matrices to get full sum of tensor over all BZ, not just IBZ.
1063          ! the atom is not necessarily invariant under symops, so these contributions should be added to each iat separately
1064          ! normalization by nsym is done at the end outside the iqpt loop and after the tetrahedron clause
1065          ! from loops above only the eigvec are kept and not the displ, so we still have to divide by the masses
1066          ! TODO: need to check the direction of the symcart vs transpose or inverse, given that jat is the pre-image of iat...
1067          do isym=1, crystal%nsym
1068            temp_33 = matmul( (symcart(:,:,isym)), matmul(msqd_atom_tmp, transpose(symcart(:,:,isym))) )
1069            jat = crystal%indsym(4,isym,iat)
1070            do idir=1,3
1071              do jdir=1,3
1072                phdos%msqd_dos_atom(:,idir,jdir,jat) = phdos%msqd_dos_atom(:,idir,jdir,jat) + &
1073                  temp_33(idir, jdir) * wdt(:,1)
1074              end do
1075            end do
1076          end do
1077        end do ! iat
1078 
1079      end do ! imode
1080    end do ! iq_ibz
1081 
1082    ABI_FREE(wdt)
1083 
1084    ! Make eigvec into phonon displacements.
1085    do iat = 1, natom
1086      full_eigvec(:,:,iat,:,:) = full_eigvec(:,:,iat,:,:) / sqrt(phdos%atom_mass(iat))
1087    end do
1088 
1089    if (my_rank == master) then
1090 #ifdef HAVE_NETCDF
1091      ! TODO: should pass prefix as arg and make it optional
1092      NCF_CHECK_MSG(nctk_open_create(ncid, strcat(prefix, "_PHIBZ.nc"), xmpi_comm_self), "Creating PHIBZ")
1093      NCF_CHECK(crystal_ncwrite(crystal, ncid))
1094      call phonons_ncwrite(ncid, natom, phdos%nqibz, qibz, wtq_ibz, full_phfrq, full_eigvec)
1095      NCF_CHECK(nf90_close(ncid))
1096 #endif
1097    end if
1098 
1099    ! immediately free this - it contains displ and not eigvec at this stage
1100    ABI_FREE(full_eigvec)
1101    ABI_FREE(full_phfrq)
1102    ABI_FREE(tmp_phfrq)
1103    ABI_FREE(tweight)
1104    ABI_FREE(dtweightde)
1105    call destroy_tetra(tetraq)
1106  else
1107 #ifdef HAVE_NETCDF
1108    MSG_WARNING('The netcdf PHIBZ file is only output for tetrahedron integration and DOS calculations')
1109 #endif
1110  end if ! tetrahedra
1111 
1112  ! Test if the initial mesh was large enough
1113  call xmpi_sum(count_wminmax, comm, ierr)
1114  call xmpi_min(wminmax(1), rtmp, comm, ierr); wminmax(1) = rtmp
1115  call xmpi_max(wminmax(2), rtmp, comm, ierr); wminmax(2) = rtmp
1116 
1117  call xmpi_sum(phdos%phdos, comm, ierr)
1118  call xmpi_sum(phdos%msqd_dos_atom, comm, ierr)
1119  call xmpi_sum(phdos%pjdos, comm, ierr)
1120  if (prtdos == 2) then
1121    ! Reduce integrals if tetra, gauss will compute idos with simpson.
1122    call xmpi_sum(phdos%phdos_int, comm, ierr)
1123    call xmpi_sum(phdos%pjdos_int, comm, ierr)
1124  end if
1125 
1126  ! normalize by nsym: symmetrization is used in all prtdos cases
1127  phdos%msqd_dos_atom = phdos%msqd_dos_atom / crystal%nsym
1128  phdos%pjdos = phdos%pjdos / crystal%nsym
1129  if (prtdos == 2) phdos%pjdos_int = phdos%pjdos_int / crystal%nsym
1130  ! normalize by mass and factor of 2, now added in the printout to agree with harmonic_thermo
1131  ! do iat=1, natom
1132  !   phdos%msqd_dos_atom(:,:,:,iat) = phdos%msqd_dos_atom(:,:,:,iat) * invmass(iat) * half
1133  ! end do ! iat
1134 
1135  ! ===============================
1136  ! === Compute Integrated PDOS ===
1137  ! ===============================
1138  ABI_CALLOC(phdos%pjdos_rc_type, (nomega, 3, crystal%ntypat))
1139  ABI_CALLOC(phdos%pjdos_type, (nomega, crystal%ntypat))
1140  ABI_CALLOC(phdos%pjdos_type_int, (nomega, crystal%ntypat))
1141 
1142  do iat=1,natom
1143    itype = crystal%typat(iat)
1144    do io=1,phdos%nomega
1145      phdos%pjdos_rc_type(io,:,itype) = phdos%pjdos_rc_type(io,:,itype) + phdos%pjdos(io,:,iat)
1146      phdos%pjdos_type(io,itype) = phdos%pjdos_type(io,itype) + sum(phdos%pjdos(io,:,iat))
1147    end do
1148    if (prtdos == 2) then
1149      do io=1,phdos%nomega
1150        phdos%pjdos_type_int(io,itype) = phdos%pjdos_type_int(io,itype) + sum(phdos%pjdos_int(io,:,iat))
1151      end do
1152    end if
1153  end do
1154 
1155  ! Evaluate IDOS using simple simpson integration
1156  ! In principle one could use derf.F90, just to be consistent ...
1157  if (prtdos == 1) then
1158    call simpson_int(phdos%nomega,phdos%omega_step, phdos%phdos, phdos%phdos_int)
1159    do iat=1,natom
1160      do idir=1,3
1161        call simpson_int(phdos%nomega, phdos%omega_step, phdos%pjdos(:,idir,iat), phdos%pjdos_int(:,idir,iat))
1162      end do
1163    end do
1164    do itype=1,crystal%ntypat
1165      call simpson_int(phdos%nomega, phdos%omega_step, phdos%pjdos_type(:,itype), phdos%pjdos_type_int(:,itype))
1166    end do
1167  end if
1168 
1169  ABI_FREE(qibz)
1170  ABI_FREE(wtq_ibz)
1171 
1172  call cwtime(cpu, wall, gflops, "stop")
1173  write(msg,'(2(a,f8.2))')" mkphdos completed. cpu:", cpu, ", wall:", wall
1174  call wrtout(std_out, msg, do_flush=.True.)
1175 
1176  DBG_EXIT("COLL")
1177 
1178 end subroutine mkphdos

m_phonons/phdos_calc_vsound [ Functions ]

[ Top ] [ m_phonons ] [ Functions ]

NAME

 phdos_calc_vsound

FUNCTION

  From the frequencies for acoustic modes at small q, estimate speed of sound (which also gives Debye temperature)

INPUTS

 eigvec(2,3*natom,3*natom) = phonon eigenvectors at present q-point
 gmet(3,3) = metric tensor in reciprocal space.
 natom = number of atoms in the unit cell
 phfrq(3*natom) = phonon frequencies at present q-point
 qphon(3) = phonon q-point
 ucvol = unit cell volume

OUTPUT

PARENTS

      m_phonons

CHILDREN

      ifc_fourq,kpath_free,phonons_ncwrite,phonons_write_gnuplot
      phonons_write_phfrq,phonons_write_xmgrace,xmpi_sum_master

SOURCE

2209 subroutine phdos_calc_vsound(eigvec,gmet,natom,phfrq,qphon,speedofsound)
2210 
2211 
2212 !This section has been created automatically by the script Abilint (TD).
2213 !Do not modify the following lines by hand.
2214 #undef ABI_FUNC
2215 #define ABI_FUNC 'phdos_calc_vsound'
2216 !End of the abilint section
2217 
2218  implicit none
2219 
2220 !Arguments -------------------------------
2221 !scalras
2222  integer, intent(in) :: natom
2223 !arrays
2224  real(dp), intent(in) :: gmet(3,3),qphon(3)
2225  real(dp), intent(in) :: phfrq(3*natom),eigvec(2,3*natom,3*natom)
2226  real(dp), intent(out) :: speedofsound(3)
2227 
2228 !Local variables -------------------------
2229  integer :: iatref,imode, iatom, isacoustic, imode_acoustic
2230 ! character(len=500) :: msg
2231  real(dp) :: qnormcart
2232  real(dp) :: qtmp(3)
2233 
2234 ! *********************************************************************
2235 
2236  imode_acoustic = 0
2237  do imode = 1, 3*natom
2238 
2239 !  Check if this mode is acoustic like: scalar product of all displacement vectors are collinear
2240    isacoustic = 1
2241 !  Find reference atom with non-zero displacement
2242    do iatom=1,natom
2243      if(sum(eigvec(:,(iatom-1)*3+1:(iatom-1)*3+3,imode)**2) >tol16)iatref=iatom
2244    enddo
2245 !  Now compute scalar product, and check they are all positive
2246    do iatom = 1, natom
2247      if (sum(eigvec(:,(iatom-1)*3+1:(iatom-1)*3+3, imode)&
2248 &           *eigvec(:,(iatref-1)*3+1:(iatref-1)*3+3, imode)) < tol16 ) isacoustic = 0
2249    end do
2250    if (isacoustic == 0) cycle
2251    imode_acoustic = min(imode_acoustic + 1, 3)
2252 
2253 !   write (msg, '(a,I6,a,3F12.4)') ' Found acoustic mode ', imode, ' for |q| in red coord < 0.25 ; q = ', qphon
2254 !   call wrtout(std_out,msg,'COLL')
2255 
2256    qtmp = matmul(gmet, qphon)
2257    qnormcart = two * pi * sqrt(sum(qphon*qtmp))
2258    speedofsound(imode_acoustic) = phfrq(imode) / qnormcart
2259  end do
2260 
2261 end subroutine phdos_calc_vsound

m_phonons/phdos_free [ Functions ]

[ Top ] [ m_phonons ] [ Functions ]

NAME

 phdos_free

FUNCTION

 destructor function for phonon DOS object

INPUTS

 PHdos= container object for phonon DOS

OUTPUT

PARENTS

      anaddb,eph

CHILDREN

      ifc_fourq,kpath_free,phonons_ncwrite,phonons_write_gnuplot
      phonons_write_phfrq,phonons_write_xmgrace,xmpi_sum_master

SOURCE

597 subroutine phdos_free(PHdos)
598 
599 
600 !This section has been created automatically by the script Abilint (TD).
601 !Do not modify the following lines by hand.
602 #undef ABI_FUNC
603 #define ABI_FUNC 'phdos_free'
604 !End of the abilint section
605 
606  implicit none
607 
608 !Arguments -------------------------------
609  type(phonon_dos_type),intent(inout) ::PHdos
610 
611 ! *************************************************************************
612 
613  !@phonon_dos_type
614  if (allocated(PHdos%atom_mass)) then
615    ABI_FREE(PHdos%atom_mass)
616  end if
617  if (allocated(PHdos%omega)) then
618    ABI_FREE(PHdos%omega)
619  end if
620  if (allocated(PHdos%phdos)) then
621    ABI_FREE(PHdos%phdos)
622  end if
623  if (allocated(PHdos%phdos_int)) then
624    ABI_FREE(PHdos%phdos_int)
625  end if
626  if (allocated(PHdos%pjdos)) then
627    ABI_FREE(PHdos%pjdos)
628  end if
629  if (allocated(PHdos%pjdos_int)) then
630    ABI_FREE(PHdos%pjdos_int)
631  end if
632  if (allocated(PHdos%pjdos_type)) then
633    ABI_FREE(PHdos%pjdos_type)
634  end if
635  if (allocated(PHdos%pjdos_type_int)) then
636    ABI_FREE(PHdos%pjdos_type_int)
637  end if
638  if (allocated(PHdos%pjdos_rc_type)) then
639    ABI_FREE(PHdos%pjdos_rc_type)
640  end if
641  if (allocated(PHdos%msqd_dos_atom)) then
642    ABI_FREE(PHdos%msqd_dos_atom)
643  end if
644 
645 end subroutine phdos_free

m_phonons/phdos_ncwrite [ Functions ]

[ Top ] [ m_phonons ] [ Functions ]

NAME

 phdos_ncwrite

FUNCTION

  Save the content of the object in a netcdf file.

INPUTS

  ncid=NC file handle (open in the caller)
  phdos<phonon_dos_type>=Container object

OUTPUT

  Only writing

NOTES

  Frequencies are in eV, DOS are in states/eV.

PARENTS

      anaddb,eph

CHILDREN

      ifc_fourq,kpath_free,phonons_ncwrite,phonons_write_gnuplot
      phonons_write_phfrq,phonons_write_xmgrace,xmpi_sum_master

SOURCE

1780 subroutine phdos_ncwrite(phdos,ncid)
1781 
1782 
1783 !This section has been created automatically by the script Abilint (TD).
1784 !Do not modify the following lines by hand.
1785 #undef ABI_FUNC
1786 #define ABI_FUNC 'phdos_ncwrite'
1787 !End of the abilint section
1788 
1789  implicit none
1790 
1791 !Arguments ------------------------------------
1792 !scalars
1793  type(phonon_dos_type),intent(in) :: phdos
1794  integer,intent(in) :: ncid
1795 
1796 !Local variables-------------------------------
1797 !scalars
1798 #ifdef HAVE_NETCDF
1799  integer :: ncerr
1800 
1801 ! *************************************************************************
1802 
1803 ! Define dimensions
1804  NCF_CHECK(nctk_def_basedims(ncid, defmode=.True.))
1805 
1806  ncerr = nctk_def_dims(ncid, [nctkdim_t("three", 3), nctkdim_t("number_of_atoms", phdos%natom),&
1807    nctkdim_t("number_of_atom_species", phdos%ntypat), nctkdim_t("number_of_frequencies", phdos%nomega)])
1808  NCF_CHECK(ncerr)
1809 
1810 !scalars
1811  NCF_CHECK(nctk_def_iscalars(ncid, ["prtdos"]))
1812  NCF_CHECK(nctk_def_dpscalars(ncid, ["dossmear"]))
1813 
1814 !arrays
1815  ncerr = nctk_def_arrays(ncid, [&
1816    nctkarr_t('wmesh', "dp", 'number_of_frequencies'),&
1817    nctkarr_t('phdos', "dp", 'number_of_frequencies'),&
1818    nctkarr_t('pjdos', "dp", 'number_of_frequencies, three, number_of_atoms'),&
1819    nctkarr_t('pjdos_type', "dp", 'number_of_frequencies, number_of_atom_species'),&
1820    nctkarr_t('pjdos_rc_type', "dp", 'number_of_frequencies, three, number_of_atom_species'), &
1821    nctkarr_t('msqd_dos_atom', "dp", 'number_of_frequencies, three, three, number_of_atoms') &
1822  ])
1823  NCF_CHECK(ncerr)
1824 
1825  ! Write variables. Note unit conversion.
1826  NCF_CHECK(nctk_set_datamode(ncid))
1827  NCF_CHECK(nf90_put_var(ncid, vid("prtdos"), phdos%prtdos))
1828  NCF_CHECK(nf90_put_var(ncid, vid('dossmear'), phdos%dossmear*Ha_eV))
1829  NCF_CHECK(nf90_put_var(ncid, vid('wmesh'), phdos%omega*Ha_eV))
1830  NCF_CHECK(nf90_put_var(ncid, vid('phdos'), phdos%phdos/Ha_eV))
1831  NCF_CHECK(nf90_put_var(ncid, vid('pjdos'), phdos%pjdos/Ha_eV))
1832  NCF_CHECK(nf90_put_var(ncid, vid('pjdos_type'), phdos%pjdos_type/Ha_eV))
1833  NCF_CHECK(nf90_put_var(ncid, vid('pjdos_rc_type'), phdos%pjdos_rc_type/Ha_eV))
1834  NCF_CHECK(nf90_put_var(ncid, vid('msqd_dos_atom'), phdos%msqd_dos_atom/Ha_eV))
1835 
1836 #else
1837  MSG_ERROR("netcdf support not enabled")
1838  ABI_UNUSED((/ncid, phdos%nomega/))
1839 #endif
1840 
1841 contains
1842  integer function vid(vname)
1843 
1844 
1845 !This section has been created automatically by the script Abilint (TD).
1846 !Do not modify the following lines by hand.
1847 #undef ABI_FUNC
1848 #define ABI_FUNC 'vid'
1849 !End of the abilint section
1850 
1851    character(len=*),intent(in) :: vname
1852    vid = nctk_idname(ncid, vname)
1853  end function vid
1854 
1855 end subroutine phdos_ncwrite

m_phonons/phdos_print [ Functions ]

[ Top ] [ m_phonons ] [ Functions ]

NAME

 phdos_print

FUNCTION

 Print out phonon DOS (and partial DOS etc) in meV units

INPUTS

 PHdos= container object for phonon DOS
 fname=File name for output

OUTPUT

  Only writing.

PARENTS

      anaddb,eph,m_tdep_phdos

CHILDREN

      ifc_fourq,kpath_free,phonons_ncwrite,phonons_write_gnuplot
      phonons_write_phfrq,phonons_write_xmgrace,xmpi_sum_master

SOURCE

199 subroutine phdos_print(PHdos,fname)
200 
201 
202 !This section has been created automatically by the script Abilint (TD).
203 !Do not modify the following lines by hand.
204 #undef ABI_FUNC
205 #define ABI_FUNC 'phdos_print'
206 !End of the abilint section
207 
208  implicit none
209 
210 !Arguments ------------------------------------
211  character(len=*),intent(in) :: fname
212  type(phonon_dos_type),intent(in) :: PHdos
213 
214 !Local variables-------------------------------
215  integer :: io,itype,unt,unt_by_atom,unt_msqd,iatom
216  real(dp) :: tens(3,3)
217  character(len=500) :: msg
218  character(len=500) :: msg_method
219  character(len=fnlen) :: fname_by_atom
220  character(len=fnlen) :: fname_msqd
221  character(len=3) :: unitname
222 
223 ! *************************************************************************
224 
225 ! Use Ha units everywhere
226  unitname='Ha'
227 
228  select case (PHdos%prtdos)
229  case (1)
230    write(msg_method,'(a,es16.8,2a,i0)')&
231 &   '# Gaussian method with smearing = ',PHdos%dossmear,unitname,', nqibz =',PHdos%nqibz
232  case (2)
233    write(msg_method,'(a,i0)')'# Tetrahedron method, nqibz= ',PHdos%nqibz
234  case default
235    MSG_ERROR(sjoin(" Wrong prtdos: ",itoa(PHdos%prtdos)))
236  end select
237 
238 ! === Open external file and write results ===
239 ! TODO Here I have to rationalize how to write all this stuff!!
240 !
241  if (open_file(fname,msg,newunit=unt,form="formatted",action="write") /= 0) then
242    MSG_ERROR(msg)
243  end if
244  write(msg,'(3a)')'# ',ch10,'# Phonon density of states and atom type projected DOS'
245  call wrtout(unt,msg,'COLL')
246  write(msg,'(6a)')'# ',ch10,'# Energy in ',unitname,', DOS in states/',unitname
247  call wrtout(unt,msg,'COLL')
248  call wrtout(unt,msg_method,'COLL')
249  write(msg,'(5a)')'# ',ch10,'# omega     PHDOS    INT_PHDOS   PJDOS[atom_type=1]  INT_PJDOS[atom_type=1] ...  ',ch10,'# '
250  call wrtout(unt,msg,'COLL')
251  do io=1,PHdos%nomega
252    write(unt,'(3es17.8)',advance='NO')PHdos%omega(io),PHdos%phdos(io),PHdos%phdos_int(io)
253    do itype=1,PHdos%ntypat
254      write(unt,'(2es17.8,2x)',advance='NO')PHdos%pjdos_type(io,itype),PHdos%pjdos_type_int(io,itype)
255    end do
256    write(unt,*)
257  end do
258  close(unt)
259 
260  fname_by_atom = trim(fname) // "_by_atom"
261  if (open_file(fname_by_atom,msg,newunit=unt_by_atom,form="formatted",action="write") /= 0) then
262    MSG_ERROR(msg)
263  end if
264  write(msg,'(3a)')'# ',ch10,'# Phonon density of states and atom projected DOS'
265  call wrtout(unt_by_atom,msg,'COLL')
266  write(msg,'(6a)')'# ',ch10,'# Energy in ',unitname,', DOS in states/',unitname
267  call wrtout(unt_by_atom,msg,'COLL')
268  call wrtout(unt_by_atom,msg_method,'COLL')
269  write(msg,'(5a)')'# ',ch10,'# omega     PHDOS    PJDOS[atom=1]  PJDOS[atom=2] ...  ',ch10,'# '
270  call wrtout(unt_by_atom,msg,'COLL')
271  do io=1,PHdos%nomega
272    write(unt_by_atom,'(2es17.8)',advance='NO')PHdos%omega(io),PHdos%phdos(io)
273    do iatom=1,PHdos%natom
274      write(unt_by_atom,'(1es17.8,2x)',advance='NO') sum(PHdos%pjdos(io,1:3,iatom))
275    end do
276    write(unt_by_atom,*)
277  end do
278  close(unt_by_atom)
279 
280  fname_msqd = trim(fname) // "_msqd"
281  if (open_file(fname_msqd,msg,newunit=unt_msqd,form="formatted",action="write") /= 0) then
282    MSG_ERROR(msg)
283  end if
284  write(msg,'(3a)')'# ',ch10,'# Phonon density of states weighted msq displacement matrix (set to zero below 1e-12)'
285  call wrtout(unt_msqd,msg,'COLL')
286  write(msg,'(6a)')'# ',ch10,'# Energy in ',unitname,', DOS in bohr^2 states/',unitname
287  call wrtout(unt_msqd,msg,'COLL')
288  call wrtout(unt_msqd,msg_method,'COLL')
289  write(msg,'(5a)')'# ',ch10,'# omega     MSQDisp[atom=1, xx, yy, zz, yz, xz, xy]  MSQDisp[atom=2, xx, yy,...] ...  ',ch10,'# '
290  call wrtout(unt_msqd,msg,'COLL')
291  do io=1,PHdos%nomega
292    write(unt_msqd,'(2es17.8)',advance='NO')PHdos%omega(io)
293    do iatom=1,PHdos%natom
294      tens = PHdos%msqd_dos_atom(io,:,:,iatom)
295      where (abs(tens) < tol12)
296         tens = zero
297      end where
298      write(unt_msqd,'(6es17.8,2x)',advance='NO') &
299 &        tens(1,1), &
300 &        tens(2,2), &
301 &        tens(3,3), &
302 &        tens(2,3), &
303 &        tens(1,3), &
304 &        tens(1,2)
305    end do
306    write(unt_msqd,*)
307  end do
308  close(unt_msqd)
309 
310 end subroutine phdos_print

m_phonons/phdos_print_msqd [ Functions ]

[ Top ] [ m_phonons ] [ Functions ]

NAME

 phdos_print_msqd

FUNCTION

  Print out mean square displacement and velocity for each atom (trace and full matrix) as a function of T
  see for example https://atztogo.github.io/phonopy/thermal-displacement.html#thermal-displacement

INPUTS

   PHdos structure

OUTPUT

   to file only

PARENTS

      anaddb

CHILDREN

      ifc_fourq,kpath_free,phonons_ncwrite,phonons_write_gnuplot
      phonons_write_phfrq,phonons_write_xmgrace,xmpi_sum_master

SOURCE

2363 subroutine phdos_print_msqd(PHdos, fname, ntemper, tempermin, temperinc)
2364 
2365 
2366 !This section has been created automatically by the script Abilint (TD).
2367 !Do not modify the following lines by hand.
2368 #undef ABI_FUNC
2369 #define ABI_FUNC 'phdos_print_msqd'
2370 !End of the abilint section
2371 
2372  implicit none
2373 
2374 !Arguments -------------------------------
2375 !scalars
2376  integer, intent(in) :: ntemper
2377  type(phonon_dos_type),intent(in) :: PHdos
2378  character(len=*),intent(in) :: fname
2379  real(dp), intent(in) :: tempermin, temperinc
2380 
2381 !Local variables -------------------------
2382  integer :: io, iomin, itemp, iunit, junit, iatom
2383  real(dp) :: temper
2384  character(len=500) :: msg
2385  character(len=fnlen) :: fname_msqd
2386  character(len=fnlen) :: fname_veloc
2387 !arrays
2388  real(dp), allocatable :: bose_msqd(:,:), tmp_msqd(:,:), integ_msqd(:,:)
2389  real(dp), allocatable :: bose_msqv(:,:), tmp_msqv(:,:), integ_msqv(:,:)
2390 
2391 ! *********************************************************************
2392 
2393  fname_msqd = trim(fname) //"_MSQD_T"
2394  if (open_file(fname_msqd, msg, newunit=iunit, form="formatted", status="unknown", action="write") /= 0) then
2395    MSG_ERROR(msg)
2396  end if
2397  fname_veloc = trim(fname) // "_MSQV_T"
2398  if (open_file(fname_veloc, msg, newunit=junit, form="formatted", status="unknown", action="write") /= 0) then
2399    MSG_ERROR(msg)
2400  end if
2401 
2402 ! write a header
2403    write (msg, '(2a)') '# mean square displacement for each atom as a function of T (bohr^2)'
2404 
2405 ! NB: this call to wrtout does not seem to work from the eph executable, even in sequential, and whether within or outside a clause for me==master.
2406 !  Do not change this to wrtout without checking extensively.
2407    !call wrtout(iunit, msg, 'COLL')
2408    write (iunit, '(a)') trim(msg)
2409 
2410    write (msg, '(2a)') '# mean square velocity for each atom as a function of T (bohr^2/atomic time unit^2)'
2411    write (junit, '(a)') trim(msg)
2412 
2413    write (msg, '(a,F18.10,a,F18.10,a)') '#  T in Kelvin, from ', tempermin, ' to ', tempermin+(ntemper-1)*temperinc
2414    !call wrtout(iunit, msg, 'COLL')
2415    write (iunit, '(a)') trim(msg)
2416    write (junit, '(a)') trim(msg)
2417 
2418    write (msg, '(2a)') '#    T             |u^2|                u_xx                u_yy                u_zz',&
2419 &                                              '                u_yz                u_xz                u_xy in bohr^2'
2420    !call wrtout(iunit, msg, 'COLL')
2421    write (iunit, '(a)') trim(msg)
2422    write (msg, '(3a)') '#    T             |v^2|                v_xx                v_yy                v_zz',&
2423 &                                              '                v_yz                v_xz                v_xy',&
2424 &                                              ' in bohr^2/atomic time unit^2'
2425    !call wrtout(iunit, msg, 'COLL')
2426    write (junit, '(a)') trim(msg)
2427 
2428  ABI_ALLOCATE (tmp_msqd, (PHdos%nomega,9))
2429  ABI_ALLOCATE (tmp_msqv, (PHdos%nomega,9))
2430  ABI_ALLOCATE (integ_msqd, (9,ntemper))
2431  ABI_ALLOCATE (integ_msqv, (9,ntemper))
2432  ABI_ALLOCATE (bose_msqd, (PHdos%nomega, ntemper))
2433  ABI_ALLOCATE (bose_msqv, (PHdos%nomega, ntemper))
2434 
2435  do io=1, PHdos%nomega
2436      if ( PHdos%omega(io) >= 2._dp * 4.56d-6 ) exit ! 2 cm-1 TODO: make this an input parameter
2437  end do
2438  iomin = io
2439 
2440  ! calculate bose only once for each atom (instead of for each atom)
2441  bose_msqd = zero
2442  bose_msqv = zero
2443  do itemp = 1, ntemper
2444    temper = tempermin + (itemp-1) * temperinc
2445    if (temper < 1.e-3) cycle ! millikelvin at least to avoid exploding Bose factor(TM)
2446    do io = iomin, PHdos%nomega
2447 ! NB: factors follow convention in phonopy documentation
2448 !   the 1/sqrt(omega) factor in phonopy is contained in the displacement vector definition
2449 !   bose() is dimensionless
2450      !bose_msqd(io, itemp) =  (half + one  / ( exp(PHdos%omega(io)/(kb_HaK*temper)) - one )) / PHdos%omega(io)
2451      !bose_msqv(io, itemp) =  (half + one  / ( exp(PHdos%omega(io)/(kb_HaK*temper)) - one )) * PHdos%omega(io)
2452      bose_msqd(io, itemp) =  (half + bose_einstein(PHdos%omega(io),kb_HaK*temper)) / PHdos%omega(io)
2453      bose_msqv(io, itemp) =  (half + bose_einstein(PHdos%omega(io),kb_HaK*temper)) * PHdos%omega(io)
2454    end do
2455  end do
2456 
2457  do iatom=1,PHdos%natom
2458    write (msg, '(a,I8)') '# atom number ', iatom
2459    !call wrtout(iunit, msg, 'COLL')
2460    write (iunit, '(a)') trim(msg)
2461    write (junit, '(a)') trim(msg)
2462 
2463 ! for each T and each atom, integrate msqd matrix with Bose Einstein factor and output
2464    integ_msqd = zero
2465    tmp_msqd = reshape(PHdos%msqd_dos_atom(:,:,:,iatom), (/PHdos%nomega, 9/))
2466 ! perform all integrations as matrix multiplication: integ_msqd (idir, itemp) = [tmp_msqd(io,idir)]^T  * bose_msqd(io,itemp)
2467    call DGEMM('T','N', 9, ntemper, PHdos%nomega, one, tmp_msqd,PHdos%nomega,&
2468 &      bose_msqd, PHdos%nomega, zero, integ_msqd, 9)
2469 ! NB: this presumes an equidistant omega grid
2470    integ_msqd = integ_msqd * (PHdos%omega(2)-PHdos%omega(1)) / PHdos%atom_mass(iatom)
2471 
2472    integ_msqv = zero
2473    tmp_msqv = reshape(PHdos%msqd_dos_atom(:,:,:,iatom), (/PHdos%nomega, 9/))
2474 ! perform all integrations as matrix multiplication: integ_msqv (idir, itemp) = [tmp_msqv(io,idir)]^T  * bose_msqv(io,itemp)
2475    call DGEMM('T','N', 9, ntemper, PHdos%nomega, one, tmp_msqv,PHdos%nomega,&
2476 &      bose_msqv, PHdos%nomega, zero, integ_msqv, 9)
2477 ! NB: this presumes an equidistant omega grid
2478    integ_msqv = integ_msqv * (PHdos%omega(2)-PHdos%omega(1)) / PHdos%atom_mass(iatom)
2479 
2480 
2481 ! print out stuff
2482    do itemp = 1, ntemper
2483      temper = tempermin + (itemp-1) * temperinc
2484      write (msg, '(F10.2,4x,E22.10,2x,6E22.10)') temper, third*(integ_msqd(1,itemp)+integ_msqd(5,itemp)+integ_msqd(9,itemp)), &
2485 &                      integ_msqd(1,itemp),integ_msqd(5,itemp),integ_msqd(9,itemp), &
2486 &                      integ_msqd(6,itemp),integ_msqd(3,itemp),integ_msqd(2,itemp)
2487      !call wrtout(iunit, msg, 'COLL')
2488      write (iunit, '(a)') trim(msg)
2489      write (msg, '(F10.2,4x,E22.10,2x,6E22.10)') temper, third*(integ_msqv(1,itemp)+integ_msqv(5,itemp)+integ_msqv(9,itemp)), &
2490 &                      integ_msqv(1,itemp),integ_msqv(5,itemp),integ_msqv(9,itemp), &
2491 &                      integ_msqv(6,itemp),integ_msqv(3,itemp),integ_msqv(2,itemp)
2492      write (junit, '(a)') trim(msg)
2493    end do ! itemp
2494 
2495    !call wrtout(iunit, msg, 'COLL')
2496    write (iunit, '(a)') ''
2497    write (junit, '(a)') ''
2498  enddo ! iatom
2499 
2500  ABI_DEALLOCATE (tmp_msqd)
2501  ABI_DEALLOCATE (tmp_msqv)
2502  ABI_DEALLOCATE (bose_msqd)
2503  ABI_DEALLOCATE (bose_msqv)
2504  ABI_DEALLOCATE (integ_msqd)
2505  ABI_DEALLOCATE (integ_msqv)
2506  close(iunit)
2507  close(junit)
2508 
2509 end subroutine phdos_print_msqd

m_phonons/phdos_print_vsound [ Functions ]

[ Top ] [ m_phonons ] [ Functions ]

NAME

 phdos_print_vsound

FUNCTION

  Print out estimate speed of sound and Debye temperature at this (small) q
  should only be called by master proc for the hard unit number

INPUTS

 unit=Fortran unit number
 speedofsound(3)

OUTPUT

PARENTS

      m_phonons

CHILDREN

      ifc_fourq,kpath_free,phonons_ncwrite,phonons_write_gnuplot
      phonons_write_phfrq,phonons_write_xmgrace,xmpi_sum_master

SOURCE

2288 subroutine phdos_print_vsound(iunit,ucvol,speedofsound)
2289 
2290 
2291 !This section has been created automatically by the script Abilint (TD).
2292 !Do not modify the following lines by hand.
2293 #undef ABI_FUNC
2294 #define ABI_FUNC 'phdos_print_vsound'
2295 !End of the abilint section
2296 
2297  implicit none
2298 
2299 !Arguments -------------------------------
2300 !scalras
2301  integer, intent(in) :: iunit
2302  real(dp), intent(in) :: ucvol
2303 !arrays
2304  real(dp), intent(in) :: speedofsound(3)
2305 
2306 !Local variables -------------------------
2307  integer :: imode_acoustic
2308  character(len=500) :: msg
2309  real(dp) :: tdebye
2310 
2311 ! *********************************************************************
2312 
2313  do imode_acoustic = 1, 3
2314 !  from phonon frequency, estimate speed of sound by linear interpolation from Gamma
2315    write (msg, '(2a,a,E20.10,a,a,F20.5)') &
2316 &   ' Speed of sound for this q and mode:',ch10,&
2317 &   '   in atomic units: ', speedofsound(imode_acoustic), ch10,&
2318 &   '   in units km/s: ', speedofsound(imode_acoustic) * Bohr_Ang * 1.d-13 / Time_Sec
2319    call wrtout(iunit,msg,'COLL')
2320    call wrtout(std_out,msg,'COLL')
2321 
2322 !  also estimate partial Debye temperature, = energy if this band went to zone edge
2323    tdebye = speedofsound(imode_acoustic) * pi * (six / pi / ucvol)**(third)
2324    write (msg, '(2a,a,E20.10,a,a,F20.5)') &
2325 &   ' Partial Debye temperature for this q and mode:',ch10,&
2326 &   '   in atomic units: ', tdebye, ch10,&
2327 &   '   in SI units K  : ', tdebye * Ha_K
2328    call wrtout(iunit,msg,'COLL')
2329    call wrtout(iunit,"",'COLL')
2330    call wrtout(std_out,msg,'COLL')
2331    call wrtout(std_out,"",'COLL')
2332  end do
2333 
2334 end subroutine phdos_print_vsound

m_phonons/phonon_dos_type [ Types ]

[ Top ] [ m_phonons ] [ Types ]

NAME

 phonon_dos_type

FUNCTION

 Container for phonon DOS and atom projected contributions

SOURCE

 87  type,public :: phonon_dos_type
 88 
 89 ! Integer
 90   integer :: ntypat
 91   ! Number of type of atoms.
 92 
 93   integer :: natom
 94   ! Number of atoms is the unit cell.
 95 
 96   integer :: prtdos
 97   ! Option of DOS calculation (1 for Gaussian, 2 for tetrahedrons).
 98 
 99   integer :: nomega
100   ! Number of frequency points in DOS mesh.
101 
102   integer :: nqibz
103   ! Number of q-points in the IBZ.
104 
105 ! Reals
106   real(dp) :: omega_min
107   ! Min frequency for DOS calculation.
108 
109   real(dp) :: omega_max
110   ! Max frequency for DOS calculation.
111 
112   real(dp) :: omega_step
113   ! Frequency step.
114 
115   real(dp) :: dossmear
116   ! Gaussian broadening.
117 
118 ! Real pointers
119   real(dp),allocatable :: atom_mass(:)
120    ! atom_mass(natom)
121 
122   real(dp),allocatable :: omega(:)
123    ! omega(nomega)
124    ! Frequency grid.
125 
126   real(dp),allocatable :: phdos(:)
127    ! phdos(nomega)
128    ! phonon DOS.
129 
130   real(dp),allocatable :: phdos_int(:)
131    ! phdos_int(nomega)
132    ! integrated phonon DOS
133 
134   real(dp),allocatable :: pjdos(:,:,:)
135    ! pjdos(nomega,3,natom)
136    ! projected DOS (over atoms and cartesian directions)
137 
138   real(dp),allocatable :: pjdos_int(:,:,:)
139    ! pjdos_int(nomega,3,natom)
140    ! Integrated atomic PJDOS along the three cartesian directions.
141 
142   real(dp),allocatable :: pjdos_type(:,:)
143    ! pjdos_type(nomega,ntypat)
144    ! phonon DOS contribution arising from a particular atom-type.
145 
146   real(dp),allocatable :: pjdos_type_int(:,:)
147    ! pjdos_type_int(nomega,ntypat)
148    ! Integrate phonon DOS contribution arising from a particular atom-type.
149 
150   real(dp),allocatable :: pjdos_rc_type(:,:,:)
151    ! phdos(nomega,3,ntypat)
152    ! phonon DOS contribution arising from a particular atom-type
153    ! decomposed along the three cartesian directions.
154 
155   real(dp),allocatable :: msqd_dos_atom(:,:,:,:)
156    ! msqd_dos_atom(nomega,3,3,natom)
157    ! mean square displacement matrix, frequency dependent like a DOS, tensor in cartesian coords.
158    ! allows one to calculate Debye Waller factors by integration with 1/omega
159    ! and the Bose Einstein factor
160 
161  end type phonon_dos_type
162 
163  public :: mkphdos
164  public :: phdos_print
165  public :: phdos_print_debye
166  public :: phdos_print_msqd
167  public :: phdos_print_thermo
168  public :: phdos_free
169  public :: phdos_ncwrite
170 !!**
171 
172 CONTAINS  !===============================================================================

m_phonons/phonons_ncwrite [ Functions ]

[ Top ] [ m_phonons ] [ Functions ]

NAME

 phonons_ncwrite

FUNCTION

  Write phonon bandstructure in a netcdf file.

INPUTS

  ncid =NC file handle
  natom=Number of atoms
  nqpts=Number of q-points.
  qpoints=List of q-points in reduced coordinates
  weights(nqpts)= q-point weights
  phfreq=Phonon frequencies
  phdispl_cart=Phonon displacementent in Cartesian coordinates.

NOTES

  Input data is in a.u, whereas the netcdf files saves data in eV for frequencies
  and Angstrom for the displacements

OUTPUT

  Only writing

PARENTS

      m_phonons

CHILDREN

      ifc_fourq,kpath_free,phonons_ncwrite,phonons_write_gnuplot
      phonons_write_phfrq,phonons_write_xmgrace,xmpi_sum_master

SOURCE

2546 subroutine phonons_ncwrite(ncid,natom,nqpts,qpoints,weights,phfreq,phdispl_cart)
2547 
2548 
2549 !This section has been created automatically by the script Abilint (TD).
2550 !Do not modify the following lines by hand.
2551 #undef ABI_FUNC
2552 #define ABI_FUNC 'phonons_ncwrite'
2553 !End of the abilint section
2554 
2555  implicit none
2556 
2557 !Arguments ------------------------------------
2558 !scalars
2559  integer,intent(in) :: ncid,natom,nqpts
2560 !arrays
2561  real(dp),intent(in) :: qpoints(3,nqpts),weights(nqpts)
2562  real(dp),intent(in) :: phfreq(3*natom,nqpts),phdispl_cart(2,3*natom,3*natom,nqpts)
2563 
2564 !Local variables-------------------------------
2565 !scalars
2566 #ifdef HAVE_NETCDF
2567  integer :: nphmodes,ncerr
2568 
2569 ! *************************************************************************
2570 
2571  nphmodes = 3*natom
2572 
2573  NCF_CHECK(nctk_def_basedims(ncid, defmode=.True.))
2574 
2575  ncerr = nctk_def_dims(ncid, [&
2576    nctkdim_t("number_of_qpoints", nqpts), nctkdim_t('number_of_phonon_modes', nphmodes)])
2577  NCF_CHECK(ncerr)
2578 
2579 ! define arrays
2580  ncerr = nctk_def_arrays(ncid, [&
2581    nctkarr_t('qpoints', "dp" , 'number_of_reduced_dimensions, number_of_qpoints'),&
2582    nctkarr_t('qweights',"dp", 'number_of_qpoints'),&
2583    nctkarr_t('phfreqs',"dp", 'number_of_phonon_modes, number_of_qpoints'),&
2584    nctkarr_t('phdispl_cart',"dp", 'complex, number_of_phonon_modes, number_of_phonon_modes, number_of_qpoints')])
2585  NCF_CHECK(ncerr)
2586 
2587 !Write variables.
2588  NCF_CHECK(nctk_set_datamode(ncid))
2589  NCF_CHECK(nf90_put_var(ncid, vid('qpoints'), qpoints))
2590  NCF_CHECK(nf90_put_var(ncid, vid('qweights'), weights))
2591  NCF_CHECK(nf90_put_var(ncid, vid('phfreqs'), phfreq*Ha_eV))
2592  NCF_CHECK(nf90_put_var(ncid, vid('phdispl_cart'), phdispl_cart*Bohr_Ang))
2593 #endif
2594 
2595 contains
2596  integer function vid(vname)
2597 
2598 
2599 !This section has been created automatically by the script Abilint (TD).
2600 !Do not modify the following lines by hand.
2601 #undef ABI_FUNC
2602 #define ABI_FUNC 'vid'
2603 !End of the abilint section
2604 
2605    character(len=*),intent(in) :: vname
2606    vid = nctk_idname(ncid, vname)
2607  end function vid
2608 
2609 end subroutine phonons_ncwrite

m_phonons/phonons_write_gnuplot [ Functions ]

[ Top ] [ m_phonons ] [ Functions ]

NAME

 phonons_write_gnuplot

FUNCTION

  Write phonons bands in gnuplot format. This routine should be called by a single processor.

INPUTS

  prefix=prefix for files (.data, .gnuplot)
  natom=Number of atoms
  nqpts=Number of q-points
  qpts(3,nqpts)=Q-points
  phfreqs(3*natom,nqpts)=Phonon frequencies.
  [qptbounds(:,:)]=Optional argument giving the extrema of the q-path.

OUTPUT

  Only writing

PARENTS

      m_phonons

CHILDREN

      ifc_fourq,kpath_free,phonons_ncwrite,phonons_write_gnuplot
      phonons_write_phfrq,phonons_write_xmgrace,xmpi_sum_master

SOURCE

3309 subroutine phonons_write_gnuplot(prefix, natom, nqpts, qpts, phfreqs, qptbounds)
3310 
3311 
3312 !This section has been created automatically by the script Abilint (TD).
3313 !Do not modify the following lines by hand.
3314 #undef ABI_FUNC
3315 #define ABI_FUNC 'phonons_write_gnuplot'
3316 !End of the abilint section
3317 
3318  implicit none
3319 
3320 !Arguments ------------------------------------
3321 !scalars
3322  integer,intent(in) :: natom,nqpts
3323  real(dp),intent(in) :: qpts(3,nqpts),phfreqs(3*natom,nqpts)
3324  character(len=*),intent(in) :: prefix
3325 !arrays
3326  real(dp),optional,intent(in) :: qptbounds(:,:)
3327 
3328 !Local variables-------------------------------
3329 !scalars
3330  integer :: unt,iq,ii,start,nqbounds,gpl_unt
3331  character(len=500) :: msg,fmt
3332  character(len=fnlen) :: datafile,basefile
3333 !arrays
3334  integer :: g0(3)
3335  integer,allocatable :: bounds2qpt(:)
3336 
3337 ! *********************************************************************
3338 
3339  nqbounds = 0
3340  if (present(qptbounds)) then
3341    if (product(shape(qptbounds)) > 0 ) then
3342      ! Find correspondence between qptbounds and k-points in ebands.
3343      nqbounds = size(qptbounds, dim=2)
3344      ABI_MALLOC(bounds2qpt, (nqbounds))
3345      bounds2qpt = 1; start = 1
3346      do ii=1,nqbounds
3347         do iq=start,nqpts
3348           if (isamek(qpts(:, iq), qptbounds(:, ii), g0)) then
3349             bounds2qpt(ii) = iq; start = iq + 1; exit
3350           end if
3351         end do
3352      end do
3353    end if
3354  end if
3355 
3356  datafile = strcat(prefix, "_PHBANDS.data")
3357  if (open_file(datafile, msg, newunit=unt, form="formatted", action="write") /= 0) then
3358    MSG_ERROR(msg)
3359  end if
3360  if (open_file(strcat(prefix, "_PHBANDS.gnuplot"), msg, newunit=gpl_unt, form="formatted", action="write") /= 0) then
3361    MSG_ERROR(msg)
3362  end if
3363  basefile = basename(datafile)
3364 
3365  write(unt,'(a)') "# Phonon band structure data file"
3366  write(unt,'(a)') "# Generated by Abinit"
3367  write(unt,'(2(a,i0))') "# natom: ",natom,", nqpt: ",nqpts
3368  write(unt,'(a)') "# Frequencies are in meV"
3369  write(unt,'(a)')"# List of q-points and their index (C notation i.e. count from 0)"
3370  do iq=1,nqpts
3371    write(unt, "(a)")sjoin("#", itoa(iq-1), ktoa(qpts(:,iq)))
3372  end do
3373 
3374  fmt = sjoin("(i0,1x,", itoa(3*natom), "(es16.8,1x))")
3375  write(unt,'(a)')"# [kpt-index, mode_1, mode_2 ...]"
3376  do iq=1,nqpts
3377    write(unt, fmt) iq-1, phfreqs(:, iq) * Ha_meV
3378  end do
3379 
3380  ! gnuplot script file
3381   write(gpl_unt,'(a)') '# File to plot electron bandstructure with gnuplot'
3382   !write(gpl_unt,'(a)') "#set terminal postscript eps enhanced color font 'Times-Roman,26' lw 2"
3383   write(gpl_unt,'(a)') '#use the next lines to make a nice figure for a paper'
3384   write(gpl_unt,'(a)') '#set term postscript enhanced eps color lw 0.5 dl 0.5'
3385   write(gpl_unt,'(a)') '#set pointsize 0.275'
3386   write(gpl_unt,'(a)') 'set palette defined ( 0 "blue", 3 "green", 6 "yellow", 10 "red" )'
3387   write(gpl_unt,'(a)') 'unset key'
3388   write(gpl_unt,'(a)') '# can make pointsize smaller (~0.5). Too small and nothing is printed'
3389   write(gpl_unt,'(a)') 'set pointsize 0.8'
3390   write(gpl_unt,'(a)') 'set view 0,0'
3391   write(gpl_unt,'(a,i0,a)') 'set xrange [0:',nqpts-1,']'
3392   write(gpl_unt,'(2(a,es16.8),a)')&
3393     'set yrange [',minval(phfreqs * Ha_meV),':',maxval(phfreqs * Ha_meV),']'
3394   write(gpl_unt,'(a)') 'set xlabel "Momentum"'
3395   write(gpl_unt,'(a)') 'set ylabel "Energy [meV]"'
3396   write(gpl_unt,'(a)') strcat('set title "', replace(basefile, "_", "\\_"),'"')
3397   if (nqbounds == 0) then
3398      write(gpl_unt,'(a)') 'set grid xtics'
3399   else
3400     write(gpl_unt,"(a)")"# Add vertical lines in correspondence of high-symmetry points."
3401     write(gpl_unt,'(a)') 'unset xtics'
3402     do ii=1,nqbounds
3403       write(gpl_unt,"(a,2(i0,a))") &
3404         "set arrow from ",bounds2qpt(ii)-1,",graph(0,0) to ",bounds2qpt(ii)-1,",graph(1,1) nohead"
3405       !write(gpl_unt,"(a)")sjoin("set xtics add ('kname'", itoa(bounds2kpt(ii)-1), ")")
3406     end do
3407   end if
3408   write(gpl_unt,"(a)")sjoin("nbranch =", itoa(3*natom))
3409   write(gpl_unt,"(a)")strcat('plot for [i=2:nbranch] "', basefile, '" u 1:i every :1 with lines linetype -1')
3410   write(gpl_unt,"(a)")"pause -1"
3411 
3412  close(unt)
3413  close(gpl_unt)
3414 
3415  if (allocated(bounds2qpt)) then
3416    ABI_FREE(bounds2qpt)
3417  end if
3418 
3419 end subroutine phonons_write_gnuplot

m_phonons/phonons_write_phfrq [ Functions ]

[ Top ] [ m_phonons ] [ Functions ]

NAME

 phonons_write_phfrq

FUNCTION

  Write phonon bandstructure in a text file. Fixed file name for the moment

INPUTS

  natom=Number of atoms
  nqpts=Number of q-points.
  qpoints=List of q-points in reduced coordinates
  weights(nqpts)= q-point weights
  phfreq=Phonon frequencies
  phdispl_cart=Phonon displacementent in Cartesian coordinates.

NOTES

  Input data is in a.u, output too

OUTPUT

  Only writing

PARENTS

      m_phonons

CHILDREN

      ifc_fourq,kpath_free,phonons_ncwrite,phonons_write_gnuplot
      phonons_write_phfrq,phonons_write_xmgrace,xmpi_sum_master

SOURCE

2644  subroutine phonons_write_phfrq(path,natom,nqpts,qpoints,weights,phfreq,phdispl_cart)
2645 
2646 
2647 !This section has been created automatically by the script Abilint (TD).
2648 !Do not modify the following lines by hand.
2649 #undef ABI_FUNC
2650 #define ABI_FUNC 'phonons_write_phfrq'
2651 !End of the abilint section
2652 
2653  implicit none
2654 
2655 !Arguments ------------------------------------
2656 !scalars
2657  integer,intent(in) :: natom,nqpts
2658  character(len=*),intent(in) :: path
2659 !arrays
2660  real(dp),intent(in) :: qpoints(3,nqpts),weights(nqpts)
2661  real(dp),intent(in) :: phfreq(3*natom,nqpts)
2662  real(dp),intent(in) :: phdispl_cart(2,3*natom,3*natom,nqpts)
2663 
2664 !Local variables-------------------------------
2665 !scalars
2666  integer :: nphmodes, iq, iunit, imod, icomp
2667  real(dp) :: dummy
2668  character(len=300) :: formt
2669  character(len=500) :: msg
2670 
2671 ! *************************************************************************
2672 
2673  nphmodes = 3*natom
2674 
2675  dummy = qpoints(1,1); dummy = weights(1)
2676 
2677  if (open_file(path, msg, newunit=iunit, form="formatted", status="unknown", action="write") /= 0) then
2678    MSG_ERROR(msg)
2679  end if
2680 
2681  write (iunit, '(a)')  '# ABINIT generated phonon band structure file. All in Ha atomic units'
2682  write (iunit, '(a)')  '# '
2683  write (iunit, '(a,i0)')  '# number_of_qpoints ', nqpts
2684  write (iunit, '(a,i0)')  '# number_of_phonon_modes ', nphmodes
2685  write (iunit, '(a)')  '# '
2686 
2687  write (formt,'(a,i0,a)') "(I5, ", nphmodes, "E20.10)"
2688  do iq= 1, nqpts
2689    write (iunit, formt)  iq, phfreq(:,iq)
2690  end do
2691 
2692  close(iunit)
2693 
2694  if (.False.) then
2695    if (open_file(strcat(path, "_PHDISPL"), msg, unit=iunit, form="formatted", status="unknown", action="write") /= 0) then
2696      MSG_ERROR(msg)
2697    end if
2698 
2699    write (iunit, '(a)')     '# ABINIT generated phonon displacements, along points in PHFRQ file. All in Ha atomic units'
2700    write (iunit, '(a)')     '# '
2701    write (iunit, '(a)')     '# displacements in cartesian coordinates, Re and Im parts '
2702    write (iunit, '(a,i0)')  '# number_of_qpoints ', nqpts
2703    write (iunit, '(a,i0)')  '# number_of_phonon_modes ', nphmodes
2704    write (iunit, '(a)')     '# '
2705 
2706    !write (formt,'(a,I3,a)') "( ", nphmodes, "(2E20.10,2x))"
2707    formt = "(2E20.10,2x)"
2708 
2709    do iq = 1, nqpts
2710      write (iunit, '(a, i0)') '# iq ', iq
2711      do imod = 1, nphmodes
2712        write (iunit, '(a, i0)') '# imode ', imod
2713        do icomp = 1, nphmodes
2714          write (iunit, formt, ADVANCE='NO') phdispl_cart(:,icomp,imod,iq)
2715        end do
2716        write (iunit, '(a)') ' '
2717      end do
2718    end do
2719 
2720    close(iunit)
2721  end if
2722 
2723 end subroutine phonons_write_phfrq

m_phonons/phonons_write_xmgrace [ Functions ]

[ Top ] [ m_phonons ] [ Functions ]

NAME

 phonons_write_xmgrace

FUNCTION

  Write phonons bands in Xmgrace format. This routine should be called by a single processor.

INPUTS

  filename=Filename
  natom=Number of atoms
  nqpts=Number of q-points
  qpts(3,nqpts)=Q-points
  phfreqs(3*natom,nqpts)=Phonon frequencies.
  [qptbounds(:,:)]=Optional argument giving the extrema of the q-path.

OUTPUT

  Only writing

PARENTS

      m_phonons

CHILDREN

      ifc_fourq,kpath_free,phonons_ncwrite,phonons_write_gnuplot
      phonons_write_phfrq,phonons_write_xmgrace,xmpi_sum_master

SOURCE

3171 subroutine phonons_write_xmgrace(filename, natom, nqpts, qpts, phfreqs, qptbounds)
3172 
3173 
3174 !This section has been created automatically by the script Abilint (TD).
3175 !Do not modify the following lines by hand.
3176 #undef ABI_FUNC
3177 #define ABI_FUNC 'phonons_write_xmgrace'
3178 !End of the abilint section
3179 
3180  implicit none
3181 
3182 !Arguments ------------------------------------
3183 !scalars
3184  integer,intent(in) :: natom,nqpts
3185  real(dp),intent(in) :: qpts(3,nqpts),phfreqs(3*natom,nqpts)
3186  character(len=*),intent(in) :: filename
3187 !arrays
3188  real(dp),optional,intent(in) :: qptbounds(:,:)
3189 
3190 !Local variables-------------------------------
3191 !scalars
3192  integer :: unt,iq,nu,ii,start,nqbounds
3193  character(len=500) :: msg
3194 !arrays
3195  integer :: g0(3)
3196  integer,allocatable :: bounds2qpt(:)
3197 
3198 ! *********************************************************************
3199 
3200  nqbounds = 0
3201  if (present(qptbounds)) then
3202    if (product(shape(qptbounds)) > 0 ) then
3203      ! Find correspondence between qptbounds and k-points in ebands.
3204      nqbounds = size(qptbounds, dim=2)
3205      ABI_MALLOC(bounds2qpt, (nqbounds))
3206      bounds2qpt = 1; start = 1
3207      do ii=1,nqbounds
3208         do iq=start,nqpts
3209           if (isamek(qpts(:, iq), qptbounds(:, ii), g0)) then
3210             bounds2qpt(ii) = iq; start = iq + 1; exit
3211           end if
3212         end do
3213      end do
3214    end if
3215  end if
3216 
3217  if (open_file(filename, msg, newunit=unt, form="formatted", action="write") /= 0) then
3218    MSG_ERROR(msg)
3219  end if
3220 
3221  write(unt,'(a)') "# Grace project file"
3222  write(unt,'(a)') "# Generated by Abinit"
3223  write(unt,'(2(a,i0))') "# natom: ",natom,", nqpt: ",nqpts
3224  write(unt,'(a)') "# Frequencies are in meV"
3225  write(unt,'(a)')"# List of q-points and their index (C notation i.e. count from 0)"
3226  do iq=1,nqpts
3227    write(unt, "(a)")sjoin("#", itoa(iq-1), ktoa(qpts(:,iq)))
3228  end do
3229 
3230  write(unt,'(a)') "@page size 792, 612"
3231  write(unt,'(a)') "@page scroll 5%"
3232  write(unt,'(a)') "@page inout 5%"
3233  write(unt,'(a)') "@link page off"
3234  write(unt,'(a)') "@with g0"
3235  write(unt,'(a)') "@world xmin 0.00"
3236  write(unt,'(a,i0)') '@world xmax ',nqpts
3237  write(unt,'(a,e16.8)') '@world ymin ',minval(phfreqs * Ha_meV)
3238  write(unt,'(a,e16.8)') '@world ymax ',maxval(phfreqs * Ha_meV)
3239  write(unt,'(a)') '@default linewidth 1.5'
3240  write(unt,'(a)') '@xaxis  tick on'
3241  write(unt,'(a)') '@xaxis  tick major 1'
3242  write(unt,'(a)') '@xaxis  tick major color 1'
3243  write(unt,'(a)') '@xaxis  tick major linestyle 3'
3244  write(unt,'(a)') '@xaxis  tick major grid on'
3245  write(unt,'(a)') '@xaxis  tick spec type both'
3246  write(unt,'(a)') '@xaxis  tick major 0, 0'
3247  if (nqbounds /= 0) then
3248    write(unt,'(a,i0)') '@xaxis  tick spec ',nqbounds
3249    do iq=1,nqbounds
3250      !write(unt,'(a,i0,a,a)') '@xaxis  ticklabel ',iq-1,',', "foo"
3251      write(unt,'(a,i0,a,i0)') '@xaxis  tick major ',iq-1,' , ',bounds2qpt(iq) - 1
3252    end do
3253  end if
3254  write(unt,'(a)') '@xaxis  ticklabel char size 1.500000'
3255  write(unt,'(a)') '@yaxis  tick major 10'
3256  write(unt,'(a)') '@yaxis  label "Phonon Energy [meV]"'
3257  write(unt,'(a)') '@yaxis  label char size 1.500000'
3258  write(unt,'(a)') '@yaxis  ticklabel char size 1.500000'
3259  do nu=1,3*natom
3260    write(unt,'(a,i0,a)') '@    s',nu-1,' line color 1'
3261  end do
3262  do nu=1,3*natom
3263    write(unt,'(a,i0)') '@target G0.S',nu-1
3264    write(unt,'(a)') '@type xy'
3265    do iq=1,nqpts
3266       write(unt,'(i0,1x,e16.8)') iq-1, phfreqs(nu, iq) * Ha_meV
3267    end do
3268    write(unt,'(a)') '&'
3269  end do
3270 
3271  close(unt)
3272 
3273  if (allocated(bounds2qpt)) then
3274    ABI_FREE(bounds2qpt)
3275  end if
3276 
3277 end subroutine phonons_write_xmgrace

m_phonons/phonons_writeEPS [ Functions ]

[ Top ] [ m_phonons ] [ Functions ]

NAME

 phonons_writeEPS

FUNCTION

  Write phonons bands in EPS format. This routine should be called by a single processor.

INPUTS

OUTPUT

PARENTS

CHILDREN

      ifc_fourq,kpath_free,phonons_ncwrite,phonons_write_gnuplot
      phonons_write_phfrq,phonons_write_xmgrace,xmpi_sum_master

SOURCE

2747 subroutine phonons_writeEPS(natom,nqpts,ntypat,typat,phfreq,phdispl_cart)
2748 
2749 
2750 !This section has been created automatically by the script Abilint (TD).
2751 !Do not modify the following lines by hand.
2752 #undef ABI_FUNC
2753 #define ABI_FUNC 'phonons_writeEPS'
2754 !End of the abilint section
2755 
2756  implicit none
2757 
2758 !Arguments ------------------------------------
2759 !scalars
2760  integer,intent(in) :: natom,nqpts,ntypat
2761 !arrays
2762  integer,intent(in) :: typat(natom)
2763  real(dp),intent(in) :: phfreq(3*natom,nqpts)
2764  real(dp),intent(in) :: phdispl_cart(2,3*natom,3*natom,nqpts)
2765 
2766 !Local variables-------------------------------
2767 !scalars
2768  integer :: cunits,EmaxN,EminN,gradRes,kmaxN,kminN,lastPos,pos,posk
2769  integer :: iatom,ii,imode,iqpt,jj,nqpt
2770  integer :: option,unt
2771  real(dp) :: E,Emax,Emin,deltaE
2772  real(dp) :: facUnit,norm,renorm
2773  character(len=500) :: msg
2774  logical :: set_color = .true.
2775 !array
2776  complex(dpc) :: displcpx(3*natom,3*natom,nqpts)
2777  integer,allocatable :: nqptl(:)
2778  real(dp),allocatable :: phfrq(:),phfrqqm1(:),scale(:)
2779  real(dp),allocatable :: colorAtom(:,:),color(:,:)
2780  real(dp),allocatable :: displ(:,:)
2781  character(len=6),allocatable :: qname(:)
2782 
2783 ! *********************************************************************
2784 
2785 
2786  if (open_file("PHFRQ.eps", msg, unit=unt, form="formatted", status="unknown", action="write") /= 0) then
2787    MSG_ERROR(msg)
2788  end if
2789 
2790 !Multiplication factor for units (from Hartree to cm-1 or THz)
2791  if(cunits==1) then
2792    facUnit=Ha_cmm1
2793  elseif(cunits==2) then
2794    facUnit=Ha_THz
2795  else
2796  end if
2797 
2798 !Boundings of the plot (only the plot and not what is around)
2799  EminN=6900
2800  EmaxN=2400
2801  kminN=2400
2802  kmaxN=9600
2803 
2804 !convert phdispl_cart in cpx array
2805  displcpx = dcmplx(phdispl_cart(1,:,:,:),phdispl_cart(2,:,:,:))
2806 
2807 !Read the input file, and store the information in a long string of characters
2808 !strlen from defs_basis module
2809  option = 1
2810 
2811 !Allocate dynamique variables
2812  ABI_ALLOCATE(phfrqqm1,(3*natom))
2813  ABI_ALLOCATE(phfrq,(3*natom))
2814  ABI_ALLOCATE(color,(3,3*natom))
2815  ABI_ALLOCATE(qname,(nqpts+1))
2816  ABI_ALLOCATE(scale,(nqpts))
2817  ABI_ALLOCATE(nqptl,(nqpts))
2818  ABI_ALLOCATE(colorAtom,(3,natom))
2819 !colorAtom(1,1:5) : atoms contributing to red (ex : [1 0 0 0 0])
2820 !colorAtom(2,1:5) : atoms contributing to green (ex : [0 1 0 0 0])
2821 !colorAtom(3,1:5) : atoms contributing to blue (ex : [0 0 1 1 1])
2822  ABI_ALLOCATE(displ,(natom,3*natom))
2823 
2824 
2825 
2826 !TEST_AM TO DO
2827 !Set Values
2828  if(ntypat /= 3) then
2829    set_color = .false.
2830  else
2831    color = zero
2832    do ii=1,natom
2833      if(typat(ii)==1) colorAtom(1,ii) = one
2834      if(typat(ii)==2) colorAtom(2,ii) = one
2835      if(typat(ii)==3) colorAtom(3,ii) = one
2836    end do
2837  end if
2838 
2839  Emin = -300.0
2840  Emax =   800.0
2841  gradRes = 8
2842  cunits = 1
2843  qname(:) = "T"
2844 !Read end of input file
2845  ! read(21,*)
2846  ! read(21,*) (qname(ii),ii=1,nqpts+1)
2847  ! read(21,*)
2848  ! read(21,*) (nqptl(ii),ii=1,nqpts)
2849  ! read(21,*)
2850  ! read(21,*) (scale(ii),ii=1,nqpts)
2851  ! read(21,*)
2852  ! read(21,*)
2853  ! read(21,*)
2854  ! read(21,*) (colorAtom(1,ii),ii=1,natom)
2855  ! read(21,*)
2856  ! read(21,*) (colorAtom(2,ii),ii=1,natom)
2857  ! read(21,*)
2858  ! read(21,*) (colorAtom(3,ii),ii=1,natom)
2859 !calculate nqpt
2860  nqpt=0
2861  do ii=1,nqpts
2862    nqpt=nqpt+nqptl(ii)
2863  end do
2864 !compute normalisation factor
2865  renorm=0
2866  do ii=1,nqpts
2867    renorm=renorm+nqptl(ii)*scale(ii)
2868  end do
2869  renorm=renorm/nqpt
2870 !Calculate Emin and Emax
2871  Emin=Emin/FacUnit
2872  Emax=Emax/FacUnit
2873 
2874 !*******************************************************
2875 !Begin to write some comments in the eps file
2876 !This is based to 'xfig'
2877 
2878  write(unt,'(a)') '% !PS-Adobe-2.0 EPSF-2.0'
2879  write(unt,'(a)') '%%Title: band.ps'
2880  write(unt,'(a)') '%%BoundingBox: 0 0 581 310'
2881  write(unt,'(a)') '%%Magnification: 1.0000'
2882 
2883  write(unt,'(a)') '/$F2psDict 200 dict def'
2884  write(unt,'(a)') '$F2psDict begin'
2885  write(unt,'(a)') '$F2psDict /mtrx matrix put'
2886  write(unt,'(a)') '/col-1 {0 setgray} bind def'
2887  write(unt,'(a)') '/col0 {0.000 0.000 0.000 srgb} bind def'
2888  write(unt,'(a)') 'end'
2889  write(unt,'(a)') 'save'
2890  write(unt,'(a)') 'newpath 0 310 moveto 0 0 lineto 581 0 lineto 581 310 lineto closepath clip newpath'
2891  write(unt,'(a)') '-36.0 446.0 translate'
2892  write(unt,'(a)') '1 -1 scale'
2893 
2894  write(unt,'(a)') '/cp {closepath} bind def'
2895  write(unt,'(a)') '/ef {eofill} bind def'
2896  write(unt,'(a)') '/gr {grestore} bind def'
2897  write(unt,'(a)') '/gs {gsave} bind def'
2898  write(unt,'(a)') '/sa {save} bind def'
2899  write(unt,'(a)') '/rs {restore} bind def'
2900  write(unt,'(a)') '/l {lineto} bind def'
2901  write(unt,'(a)') '/m {moveto} bind def'
2902  write(unt,'(a)') '/rm {rmoveto} bind def'
2903  write(unt,'(a)') '/n {newpath} bind def'
2904  write(unt,'(a)') '/s {stroke} bind def'
2905  write(unt,'(a)') '/sh {show} bind def'
2906  write(unt,'(a)') '/slc {setlinecap} bind def'
2907  write(unt,'(a)') '/slj {setlinejoin} bind def'
2908  write(unt,'(a)') '/slw {setlinewidth} bind def'
2909  write(unt,'(a)') '/srgb {setrgbcolor} bind def'
2910  write(unt,'(a)') '/rot {rotate} bind def'
2911  write(unt,'(a)') '/sc {scale} bind def'
2912  write(unt,'(a)') '/sd {setdash} bind def'
2913  write(unt,'(a)') '/ff {findfont} bind def'
2914  write(unt,'(a)') '/sf {setfont} bind def'
2915  write(unt,'(a)') '/scf {scalefont} bind def'
2916  write(unt,'(a)') '/sw {stringwidth} bind def'
2917  write(unt,'(a)') '/tr {translate} bind def'
2918  write(unt,'(a)') '/tnt {dup dup currentrgbcolor'
2919 
2920  write(unt,'(a)') '4 -2 roll dup 1 exch sub 3 -1 roll mul add'
2921  write(unt,'(a)') '4 -2 roll dup 1 exch sub 3 -1 roll mul add'
2922  write(unt,'(a)') '4 -2 roll dup 1 exch sub 3 -1 roll mul add srgb}'
2923  write(unt,'(a)') 'bind def'
2924  write(unt,'(a)') '/shd {dup dup currentrgbcolor 4 -2 roll mul 4 -2 roll mul'
2925  write(unt,'(a)') ' 4 -2 roll mul srgb} bind def'
2926  write(unt,'(a)') '/$F2psBegin {$F2psDict begin /$F2psEnteredState save def} def'
2927  write(unt,'(a)') '/$F2psEnd {$F2psEnteredState restore end} def'
2928  write(unt,'(a)') '$F2psBegin'
2929  write(unt,'(a)') '%%Page: 1 1'
2930  write(unt,'(a)') '10 setmiterlimit'
2931  write(unt,'(a)') '0.06000 0.06000 sc'
2932 
2933 !****************************************************************
2934 !Begin of the intelligible part of the postcript document
2935 
2936  write(unt,'(a)') '%**************************************'
2937 !****************************************************************
2938 !Draw the box containing the plot
2939  write(unt,'(a)') '%****Big Box****'
2940  write(unt,'(a)') '16 slw'
2941  write(unt,'(a,i4,a,i4,a,i4,a,i4,a,i4,a,i4,a,i4,a,i4,a)') 'n ', kminN,' ', EmaxN,&
2942 & ' m ', kmaxN,' ', EmaxN, ' l ', &
2943 & kmaxN,' ', EminN, ' l ', kminN,' ', EminN, ' l'
2944  write(unt,'(a)') 'cp gs col0 s gr'
2945 
2946 !****************************************************************
2947 !Write unit on the middle left of the vertical axe
2948  write(unt,'(a)') '%****Units****'
2949  if(cunits==1) then
2950 !  1/lambda
2951    write(unt,'(a)') '/Times-Roman ff 270.00 scf sf'
2952    write(unt,'(a)') '1425 5650 m'
2953    write(unt,'(3a)') 'gs 1 -1 sc  90.0 rot (Frequency ',achar(92),'(cm) col0 sh gr'
2954 !  cm-1
2955    write(unt,'(a)') '/Times-Roman ff 200.00 scf sf'
2956    write(unt,'(a)') '1325 4030 m'
2957    write(unt,'(a)') 'gs 1 -1 sc 90.0 rot  (-1) col0 sh gr'
2958    write(unt,'(a)') '/Times-Roman ff 270.00 scf sf'
2959    write(unt,'(a)') '1425 3850 m'
2960    write(unt,'(3a)') 'gs 1 -1 sc  90.0 rot (',achar(92),')) col0 sh gr'
2961  else
2962 !  Freq
2963    write(unt,'(a)') '/Times-Roman ff 270.00 scf sf'
2964    write(unt,'(a)') '825 4850 m'
2965    write(unt,'(a)') 'gs 1 -1 sc  90.0 rot (Freq) col0 sh gr'
2966 !  THz
2967    write(unt,'(a)') '/Times-Roman ff 270.00 scf sf'
2968    write(unt,'(a)') '825 4350 m'
2969    write(unt,'(a)') 'gs 1 -1 sc 90.0 rot  (THz) col0 sh gr'
2970  end if
2971 !*****************************************************************
2972 !Write graduation on the vertical axe
2973  write(unt,'(a)') '%****Vertical graduation****'
2974  deltaE=(Emax-Emin)/gradRes
2975 
2976 !Replacing do loop with real variables with standard g95 do loop
2977  E=Emin
2978  do
2979 !  do E=Emin,(Emax-deltaE/2),deltaE
2980    if (E >= (Emax-deltaE/2)-tol6) exit
2981    pos=int(((EminN-EmaxN)*E &
2982 &   +EmaxN*Emin -EminN*Emax)/(Emin-Emax))
2983 
2984 !  write the value of energy(or frequence)
2985    write(unt,'(a)') '/Times-Roman ff 270.00 scf sf'
2986    write(unt,'(i4,a,i4,a)') kminN-800,' ',pos+60,' m'        !-1300 must be CHANGED
2987 !  as a function of the width of E
2988    write(unt,'(a,i6,a)') 'gs 1 -1 sc (', nint(E*facUnit),') col0 sh gr'
2989 
2990 !  write a little bar
2991    write(unt,'(a,i4,a,i4,a,i4,a,i4,a)') 'n ', kminN,' ',pos ,' m ', kminN+100,' ', pos, ' l'
2992    write(unt,'(a)') 'gs col0 s gr '
2993 
2994    E = E+deltaE
2995  end do
2996 
2997 !do the same thing for E=Emax (floating point error)
2998  write(unt,'(a)') '/Times-Roman ff 270.00 scf sf'
2999  write(unt,'(i4,a,i4,a)') kminN-800,' ',EmaxN+60,' m'        !-1300 must be changed as E
3000  write(unt,'(a,i6,a)') 'gs 1 -1 sc (', nint(Emax*facUnit),') col0 sh gr'
3001 
3002 
3003 !draw zero line
3004  E=0
3005  pos=int(((EminN-EmaxN)*E &
3006 & +EmaxN*Emin -EminN*Emax)/(Emin-Emax))
3007  write(unt,'(a,i4,a,i4,a,i4,a,i4,a)') 'n ', kminN,' ',pos ,' m ', kmaxN,' ', pos, ' l'
3008  write(unt,'(a)') 'gs col0 s gr '
3009 
3010 
3011 !******************************************************
3012 !draw legend of horizontal axe
3013 !+vertical line
3014 
3015  write(unt,'(a)') '%****Horizontal graduation****'
3016 
3017  lastPos=kminN
3018 
3019  do ii=0,nqpts
3020 
3021    if(ii/=0) then
3022      posk=int(((kminN-kmaxN)*(nqptl(ii))) &
3023 &     *scale(ii)/renorm/(-nqpt))
3024    else
3025      posk=0
3026    end if
3027 
3028    posk=posk+lastPos
3029    lastPos=posk
3030 
3031    if(qname(ii+1)=='gamma') then             !GAMMA
3032      write(unt,'(a)') '/Symbol ff 270.00 scf sf'
3033      write(unt,'(i4,a,i4,a)') posk-100,' ', 7150, ' m'
3034      write(unt,'(a)') 'gs 1 -1 sc (G) col0 sh gr'
3035    elseif(qname(ii+1)=='lambda') then              !LAMBDA
3036      write(unt,'(a)') '/Symbol ff 270.00 scf sf'
3037      write(unt,'(i4,a,i4,a)') posk-100,' ', 7150, ' m'
3038      write(unt,'(a)') 'gs 1 -1 sc (L) col0 sh gr'
3039    else                                     !autre
3040      write(unt,'(a)') '/Times-Roman ff 270.00 scf sf'
3041      write(unt,'(i4,a,i4,a)') posk-100,' ', 7150, ' m'
3042      write(unt,'(a,a1,a)') 'gs 1 -1 sc (',qname(ii+1),') col0 sh gr'
3043    end if
3044 
3045 
3046 !  draw vertical line
3047    write(unt,'(a,i4,a,i4,a,i4,a,i4,a)') 'n ', posk,' ',EminN ,' m ', posk,' ', EmaxN, ' l'
3048    write(unt,'(a)') 'gs col0 s gr '
3049 
3050 
3051  end do
3052 
3053 
3054 
3055 
3056 !***********************************************************
3057 !Write the bands (the most important part actually)
3058 
3059  write(unt,'(a)') '%****Write Bands****'
3060 
3061 
3062 ! read(19,*) (phfrqqm1(ii),ii=1,3*natom)
3063  jj = 1
3064  lastPos=kminN
3065  do iqpt=1,nqpts
3066 !  Copy frequency of the qpoint
3067    phfrqqm1(:) = phfreq(:,iqpt)
3068 !  Set displacement
3069    do ii=1,3*natom
3070      do iatom=1,natom
3071        displ(iatom,ii) =  real(sqrt(displcpx(3*(iatom-1)+1,ii,iqpt)*   &
3072            conjg(displcpx(3*(iatom-1)+1,ii,iqpt)) + &
3073 &                displcpx(3*(iatom-1)+2,ii,iqpt)*   &
3074 &          conjg(displcpx(3*(iatom-1)+2,ii,iqpt)) + &
3075 &               displcpx(3*(iatom-1)+3,ii,iqpt)*   &
3076 &          conjg(displcpx(3*(iatom-1)+3,ii,iqpt)) ))
3077      end do
3078    end do
3079 
3080 
3081    do imode=1,3*natom
3082 !    normalize displ
3083      norm=0
3084      do iatom=1,natom
3085        norm=norm+displ(iatom,imode)
3086      end do
3087 
3088      do iatom=1,natom
3089        displ(iatom,imode)=displ(iatom,imode)/norm
3090      end do
3091 
3092 !    Treat color
3093      color(:,imode)=0
3094      if(set_color)then
3095        do ii=1,natom
3096 !        Red
3097          color(1,imode)=color(1,imode)+displ(ii,imode)*colorAtom(1,ii)
3098 !        Green
3099          color(2,imode)=color(2,imode)+displ(ii,imode)*colorAtom(2,ii)
3100 !        Blue
3101          color(3,imode)=color(3,imode)+displ(ii,imode)*colorAtom(3,ii)
3102        end do
3103      end if
3104 
3105      pos=int(((EminN-EmaxN)*phfrqqm1(imode) &
3106 &     +EmaxN*Emin -EminN*Emax)/(Emin-Emax))
3107 
3108      posk=int(((kminN-kmaxN)*(iqpt-1) &
3109 &        *scale(jj)/renorm/(-nqpts)))
3110      posk=posk+lastPos
3111      write(unt,'(a,i4,a,i4,a)') 'n ',posk,' ',pos,' m'
3112      pos=int(((EminN-EmaxN)*phfrq(imode) &
3113 &       +EmaxN*Emin -EminN*Emax)/(Emin-Emax))
3114      posk=int(((kminN-kmaxN)*(iqpt) &
3115 &       *scale(jj)/renorm/(-nqpts)))
3116      posk=posk+lastPos
3117      write(unt,'(i4,a,i4,a)') posk,' ',pos,' l gs'
3118 
3119      if(set_color) then     !(in color)
3120        write(unt,'(f6.3,a,f6.3,a,f6.3,a)') color(1,imode),' ', &
3121 &        color(2,imode),' ',color(3,imode), ' srgb s gr'
3122      else
3123        write(unt,'(f6.3,a,f6.3,a,f6.3,a)') 0.0,' ', &
3124 &        0.0,' ',0.0, ' srgb s gr'
3125      end if
3126    end do
3127    lastPos=posk
3128  end do
3129 
3130 
3131 !**********************************************************
3132 !Ending the poscript document
3133  write(unt,'(a)') '$F2psEnd'
3134  write(unt,'(a)') 'rs'
3135 
3136 ! *************************************************************************
3137  close(unt)
3138 
3139 end subroutine phonons_writeEPS

m_phonons/thermal_supercell_free [ Functions ]

[ Top ] [ m_phonons ] [ Functions ]

NAME

 thermal_supercell_free

FUNCTION

  deallocate thermal array of supercells

INPUTS

OUTPUT

NOTES

PARENTS

      anaddb

CHILDREN

      ifc_fourq,kpath_free,phonons_ncwrite,phonons_write_gnuplot
      phonons_write_phfrq,phonons_write_xmgrace,xmpi_sum_master

SOURCE

1606 subroutine thermal_supercell_free(nscells, thm_scells)
1607 
1608 
1609 !This section has been created automatically by the script Abilint (TD).
1610 !Do not modify the following lines by hand.
1611 #undef ABI_FUNC
1612 #define ABI_FUNC 'thermal_supercell_free'
1613 !End of the abilint section
1614 
1615  implicit none
1616 
1617 !Arguments ------------------------------------
1618 !scalars
1619  integer, intent(in) :: nscells
1620  type(supercell_type), allocatable, intent(inout) :: thm_scells(:)
1621 
1622 ! local
1623  integer :: icell
1624 
1625  if(allocated(thm_scells)) then
1626    do icell = 1, nscells
1627      call destroy_supercell(thm_scells(icell))
1628    end do
1629  end if
1630 end subroutine thermal_supercell_free

m_phonons/thermal_supercell_make [ Functions ]

[ Top ] [ m_phonons ] [ Functions ]

NAME

 thermal_supercell_make

FUNCTION

  Construct an random thermalized supercell configuration, as in TDEP
  main function is for training set generation in multibinit

INPUTS

   Crystal = crystal object with rprim etc...
   Ifc = interatomic force constants object from anaddb
   option = option to deal with negative frequency -> Bose factor explodes (eg acoustic at Gamma)
      several philosophies to be implemented for the unstable modes:
      option == 1 =>  ignore
      option == 2 =>  populate them according to a default amplitude
      option == 3 =>  populate according to their modulus squared
      option == 4 =>  USER defined value(s), require namplitude and amplitude
   nconfig = numer of requested configurations
   rlatt = matrix of conversion for supercell (3 0 0   0 3 0   0 0 3 for example)
   temperature_K =  temperature in Kelvin
   nqpt = number of q-point 
   namplitude = number of amplitude provided by the user
   amplitudes(namplitude) = list of the amplitudes of the unstable phonons
                            amplitudes(1:3,iamplitude) = qpt
                            amplitudes(4,iamplitude)   = mode
                            amplitudes(5,iamplitude)   = amplitude

OUTPUT

   thm_scells = array of configurations with thermalized supercells

NOTES

PARENTS

CHILDREN

      ifc_fourq,kpath_free,phonons_ncwrite,phonons_write_gnuplot
      phonons_write_phfrq,phonons_write_xmgrace,xmpi_sum_master

SOURCE

1396 subroutine thermal_supercell_make(amplitudes,Crystal, Ifc,namplitude, nconfig,option,&
1397 &                                 rlatt, temperature_K, thm_scells)
1398 
1399 
1400 !This section has been created automatically by the script Abilint (TD).
1401 !Do not modify the following lines by hand.
1402 #undef ABI_FUNC
1403 #define ABI_FUNC 'thermal_supercell_make'
1404 !End of the abilint section
1405 
1406  implicit none
1407 
1408 !Arguments ------------------------------------
1409 !scalars
1410  integer, intent(in) :: option,nconfig
1411  integer, intent(in) :: rlatt(3,3)
1412  real(dp), intent(in) :: temperature_K
1413  type(crystal_t),intent(in) :: Crystal
1414  type(ifc_type),intent(in) :: Ifc
1415  type(supercell_type), intent(out) :: thm_scells(nconfig)
1416  integer,intent(in) :: namplitude
1417 !Local variables-------------------------------
1418 !scalars
1419  integer :: iq, nqibz, nqbz, qptopt1, iampl ,imode, ierr, iconfig
1420  real(dp) :: temperature, sigma, freeze_displ
1421  real(dp) :: rand !, rand1, rand2
1422  real(dp),intent(in):: amplitudes(5,namplitude)
1423  !arrays
1424  real(dp), allocatable :: qshft(:,:) ! dummy with 2 dimensions for call to kpts_ibz_from_kptrlatt
1425  real(dp), allocatable :: qbz(:,:), qibz(:,:), wtqibz(:)
1426  real(dp), allocatable :: phfrq_allq(:,:), phdispl_allq(:,:,:,:,:)
1427  real(dp), allocatable :: phfrq(:), phdispl(:,:,:,:),pheigvec(:,:,:,:)
1428  real(dp), allocatable :: phdispl1(:,:,:)
1429  character (len=500) :: msg
1430 
1431 ! *************************************************************************
1432 ! check inputs
1433 ! TODO: add check that all rlatt are the same on input
1434  if (rlatt(1,2)/=0 .or.  rlatt(1,3)/=0 .or.  rlatt(2,3)/=0 .or. &
1435 &    rlatt(2,1)/=0 .or.  rlatt(3,1)/=0 .or.  rlatt(3,2)/=0) then
1436    write (msg, '(4a, 9I6, a)') ' for the moment I have not implemented ', &
1437 &    ' non diagonal supercells.',ch10,' rlatt for temp 1 = ', rlatt, ' Returning '
1438    MSG_WARNING(msg)
1439    return
1440  end if
1441 
1442  temperature = temperature_K /  Ha_K
1443 
1444  ! build qpoint grid used for the Fourier interpolation.
1445  !(use no syms for the moment!)
1446  qptopt1 = 3
1447 
1448  ! for the moment do not allow shifted q grids.
1449  ! We are interpolating anyway, so it will always work
1450  ABI_MALLOC(qshft,(3,1))
1451  qshft(:,1)=zero
1452 
1453  ! This call will set nqibz, IBZ and BZ arrays
1454  call kpts_ibz_from_kptrlatt(crystal, rlatt, qptopt1, 1, qshft, &
1455 &   nqibz, qibz, wtqibz, nqbz, qbz) ! new_kptrlatt, new_shiftk)  ! Optional
1456  ABI_FREE(qshft)
1457 
1458  ! allocate arrays wzith all of the q, omega, and displacement vectors
1459  ABI_STAT_MALLOC(phfrq_allq, (3*Crystal%natom, nqibz), ierr)
1460  ABI_CHECK(ierr==0, 'out-of-memory in phfrq_allq')
1461  ABI_STAT_MALLOC(phdispl_allq, (2, 3, Crystal%natom, 3*Crystal%natom, nqibz), ierr)
1462  ABI_CHECK(ierr==0, 'out-of-memory in phdispl_allq')
1463 
1464  ABI_STAT_MALLOC(phfrq, (3*Crystal%natom), ierr)
1465  ABI_CHECK(ierr==0, 'out-of-memory in phfrq_allq')
1466  ABI_STAT_MALLOC(phdispl, (2, 3, Crystal%natom, 3*Crystal%natom), ierr)
1467  ABI_CHECK(ierr==0, 'out-of-memory in phdispl_allq')
1468  ABI_STAT_MALLOC(pheigvec, (2, 3, Crystal%natom, 3*Crystal%natom), ierr)
1469  ABI_CHECK(ierr==0, 'out-of-memory in phdispl_allq')
1470 
1471  ! loop over q to get all frequencies and displacement vectors
1472  imode = 0
1473  do iq = 1, nqibz
1474    ! Fourier interpolation.
1475    call ifc_fourq(Ifc, Crystal, qibz(:,iq), phfrq, phdispl, out_eigvec=pheigvec)
1476    phfrq_allq(1:3*Crystal%natom, iq) = phfrq
1477    phdispl_allq(1:2, 1:3, 1:Crystal%natom, 1:3*Crystal%natom, iq) = phdispl
1478  end do
1479  ABI_FREE(phfrq)
1480  ABI_FREE(pheigvec)
1481  ABI_FREE(phdispl)
1482 
1483  ! only diagonal supercell case for the moment
1484  do iconfig = 1, nconfig
1485    call init_supercell(Crystal%natom, rlatt, Crystal%rprimd, Crystal%typat, Crystal%xcart, Crystal%znucl, thm_scells(iconfig))
1486  end do
1487 
1488  ! precalculate phase factors???
1489 
1490  ABI_STAT_MALLOC(phdispl1, (2, 3, Crystal%natom), ierr)
1491 
1492  ! for all modes at all q in whole list, sorted
1493  do iq = 1, nqibz
1494    do imode = 1, 3*Crystal%natom
1495 
1496      ! skip modes with too low or negative frequency -> Bose factor explodes (eg acoustic at Gamma)
1497      ! TODO: check the convergence wrt the tolerance
1498      ! several philosophies to be implemented for the unstable modes:
1499      ! 1) ignore
1500      ! 2) populate them according to a default amplitude
1501      ! 3) populate according to their modulus squared
1502      if (abs(phfrq_allq(imode, iq))<tol6) cycle
1503 
1504      phdispl1 = phdispl_allq(:,:,:,imode,iq)
1505 
1506      ! loop over configurations
1507      do iconfig = 1, nconfig
1508 
1509        ! trick supercell object into using present q point
1510        thm_scells(iconfig)%qphon(:) = qibz(:,iq)
1511 
1512        ! find thermal displacement amplitude eq 4 of Zacharias
1513        !   combined with l_nu,q expression in paragraph before
1514        if (phfrq_allq(imode, iq) > tol6) then
1515          sigma = sqrt( (bose_einstein(phfrq_allq(imode,iq), temperature) + half)/phfrq_allq(imode,iq))
1516        else
1517          !Treat negative frequencies
1518          select case (option)
1519          case(1)
1520            !Do not populate
1521            sigma = 0._dp
1522          case(2)
1523            !Default amplitude for all the frequencies
1524            sigma = 100._dp
1525          case(3)
1526            !Absolute value of the frequencie 
1527            sigma=sqrt((bose_einstein(abs(phfrq_allq(imode,iq)),temperature)+half)/&
1528 &                abs(phfrq_allq(imode,iq)))
1529          case(4)
1530            sigma = 0._dp
1531            !Search if the amplitude of this unstable phonon is in the input argument amplitudes
1532            do iampl=1,namplitude
1533              if(abs(thm_scells(iconfig)%qphon(1) - amplitudes(1,iampl)) < tol8.and.&
1534 &               abs(thm_scells(iconfig)%qphon(2) - amplitudes(2,iampl)) < tol8.and.&
1535 &               abs(thm_scells(iconfig)%qphon(3) - amplitudes(3,iampl)) < tol8.and.&
1536 &               abs(imode - amplitudes(4,iampl)) < tol8) then
1537                sigma = amplitudes(5,iampl)
1538              end if
1539            end do
1540            !If not, the amplitude is zero
1541            if(abs(sigma) < tol8)then
1542              write (msg, '(a,I0,a,3es12.5,2a,I0)') ' The amplitude of the unstable mode ',&
1543 &                int(imode),' of the qpt ',thm_scells(iconfig)%qphon(:), ch10,&
1544 &                'is set to zero for the configuration ',iconfig
1545              MSG_WARNING(msg)
1546            end if
1547          end select
1548        end if
1549 
1550        ! add displacement for this mode to supercell positions eq 5 of Zacharias
1551        call RANDOM_NUMBER(rand)
1552        rand = two * rand - one
1553 
1554        ! from TDEP documentation for gaussian distribution of displacements
1555        !call RANDOM_NUMBER(rand1)
1556        !call RANDOM_NUMBER(rand2)
1557        ! rand = sqrt(-two*rand1) * sin(twopi*rand2)
1558 
1559        ! if (rand > half) then
1560        !   rand = one
1561        ! else
1562        !   rand = -one
1563        ! end if
1564 
1565        freeze_displ =  rand * sigma
1566 
1567        call freeze_displ_supercell (phdispl1(:,:,:), freeze_displ, thm_scells(iconfig))       
1568      end do !iconfig
1569    end do !imode
1570  end do !iq
1571 
1572  ABI_FREE(phfrq_allq)
1573  ABI_FREE(phdispl_allq)
1574  ABI_FREE(phdispl1)
1575  ABI_FREE(qibz)
1576  ABI_FREE(qbz)
1577  ABI_FREE(wtqibz)
1578 
1579 end subroutine thermal_supercell_make

m_phonons/thermal_supercell_print [ Functions ]

[ Top ] [ m_phonons ] [ Functions ]

NAME

 thermal_supercell_print

FUNCTION

  print files with thermalized array of random supercell configurations

INPUTS

OUTPUT

NOTES

PARENTS

      anaddb

CHILDREN

      ifc_fourq,kpath_free,phonons_ncwrite,phonons_write_gnuplot
      phonons_write_phfrq,phonons_write_xmgrace,xmpi_sum_master

SOURCE

1717 subroutine thermal_supercell_print(fname, nconfig, temperature_K, thm_scells)
1718 
1719 
1720 !This section has been created automatically by the script Abilint (TD).
1721 !Do not modify the following lines by hand.
1722 #undef ABI_FUNC
1723 #define ABI_FUNC 'thermal_supercell_print'
1724 !End of the abilint section
1725 
1726  implicit none
1727 
1728 !Arguments ------------------------------------
1729 !scalars
1730  integer, intent(in) :: nconfig
1731  type(supercell_type), intent(in) :: thm_scells(nconfig)
1732  character(len=fnlen), intent(in) :: fname
1733  real(dp), intent(in) :: temperature_K
1734 
1735 ! local
1736  integer :: iconfig,itemp
1737  character(len=80) :: title1, title2
1738  character(len=fnlen) :: filename
1739  character(len=10) :: config_str
1740 
1741  do iconfig = 1, nconfig
1742    write (config_str,'(I8)') iconfig
1743    write (filename, '(3a)') trim(fname), "_cf_", trim(adjustl(config_str))
1744    write (title1, '(a,I6,a)') "#  thermalized supercell at temperature T= ", temperature_K, " Kelvin"
1745    title2 = "#  generated with random thermal displacements of all phonons"
1746    call prt_supercell (filename, thm_scells(itemp), title1, title2)
1747  end do
1748 
1749 end subroutine thermal_supercell_print

m_phonons/zacharias_supercell_make [ Functions ]

[ Top ] [ m_phonons ] [ Functions ]

NAME

 zacharias_supercell_make

FUNCTION

  Construct an optimally thermalized supercell following Zacharias and Giustino
 PRB 94 075125 (2016) [[cite:Zacharias2016]]

INPUTS

OUTPUT

NOTES

PARENTS

      anaddb

CHILDREN

      ifc_fourq,kpath_free,phonons_ncwrite,phonons_write_gnuplot
      phonons_write_phfrq,phonons_write_xmgrace,xmpi_sum_master

SOURCE

1206 subroutine zacharias_supercell_make(Crystal, Ifc, ntemper, &
1207 &    rlatt, tempermin, temperinc, thm_scells)
1208 
1209 
1210 !This section has been created automatically by the script Abilint (TD).
1211 !Do not modify the following lines by hand.
1212 #undef ABI_FUNC
1213 #define ABI_FUNC 'zacharias_supercell_make'
1214 !End of the abilint section
1215 
1216  implicit none
1217 
1218 !Arguments ------------------------------------
1219 !scalars
1220  integer, intent(in) :: ntemper
1221  integer, intent(in) :: rlatt(3,3)
1222  real(dp), intent(in) :: tempermin, temperinc
1223  type(crystal_t),intent(in) :: Crystal
1224  type(ifc_type),intent(in) :: Ifc
1225  type(supercell_type), intent(out) :: thm_scells(ntemper)
1226 
1227 !Local variables-------------------------------
1228 !scalars
1229  integer :: iq, nqibz, nqbz, qptopt1, imode, itemper, ierr, jmode
1230  real(dp) :: temperature_K, temperature, modesign, sigma, freeze_displ
1231 
1232 !arrays
1233  integer, allocatable :: modeindex(:)
1234  real(dp), allocatable :: qshft(:,:) ! dummy with 2 dimensions for call to kpts_ibz_from_kptrlatt
1235  real(dp), allocatable :: qbz(:,:), qibz(:,:), wtq_ibz(:)
1236  real(dp), allocatable :: phfrq_allq(:), phdispl_allq(:,:,:,:,:)
1237  real(dp), allocatable :: phfrq(:), phdispl(:,:,:,:),pheigvec(:,:,:,:)
1238  real(dp), allocatable :: phdispl1(:,:,:)
1239  character (len=500) :: msg
1240 
1241 ! *************************************************************************
1242 
1243  ! check inputs
1244 ! TODO: add check that all rlatt are the same on input
1245 
1246  if (rlatt(1,2)/=0 .or.  rlatt(1,3)/=0 .or.  rlatt(2,3)/=0 .or. &
1247 &    rlatt(2,1)/=0 .or.  rlatt(3,1)/=0 .or.  rlatt(3,2)/=0) then
1248    write (msg, '(4a, 9I6, a)') ' for the moment I have not implemented ', &
1249 &    ' non diagonal supercells.',ch10,' rlatt for temp 1 = ', rlatt, ' Returning '
1250    MSG_WARNING(msg)
1251    return
1252  end if
1253 
1254  ! build qpoint grid used for the Fourier interpolation.
1255  !(use no syms for the moment!)
1256  qptopt1 = 3
1257 
1258  ! for the moment do not allow shifted q grids.
1259  ! We are interpolating anyway, so it will always work
1260  ABI_MALLOC(qshft,(3,1))
1261  qshft(:,1)=zero
1262 
1263  ! This call will set nqibz, IBZ and BZ arrays
1264  call kpts_ibz_from_kptrlatt(crystal, rlatt, qptopt1, 1, qshft, &
1265 &   nqibz, qibz, wtq_ibz, nqbz, qbz) ! new_kptrlatt, new_shiftk)  ! Optional
1266  ABI_FREE(qshft)
1267 
1268  ! allocate arrays with all of the q, omega, and displacement vectors
1269  ABI_STAT_MALLOC(phfrq_allq, (3*Crystal%natom*nqibz), ierr)
1270  ABI_CHECK(ierr==0, 'out-of-memory in phfrq_allq')
1271  ABI_STAT_MALLOC(phdispl_allq, (2, 3, Crystal%natom, 3*Crystal%natom, nqibz), ierr)
1272  ABI_CHECK(ierr==0, 'out-of-memory in phdispl_allq')
1273 
1274  ABI_STAT_MALLOC(phfrq, (3*Crystal%natom), ierr)
1275  ABI_CHECK(ierr==0, 'out-of-memory in phfrq_allq')
1276  ABI_STAT_MALLOC(phdispl, (2, 3, Crystal%natom, 3*Crystal%natom), ierr)
1277  ABI_CHECK(ierr==0, 'out-of-memory in phdispl_allq')
1278  ABI_STAT_MALLOC(pheigvec, (2, 3, Crystal%natom, 3*Crystal%natom), ierr)
1279  ABI_CHECK(ierr==0, 'out-of-memory in phdispl_allq')
1280 
1281  ! loop over q to get all frequencies and displacement vectors
1282  ABI_ALLOCATE(modeindex, (nqibz*3*Crystal%natom))
1283  imode = 0
1284  do iq = 1, nqibz
1285    ! Fourier interpolation.
1286    call ifc_fourq(Ifc, Crystal, qibz(:,iq), phfrq, phdispl, out_eigvec=pheigvec)
1287    phfrq_allq((iq-1)*3*Crystal%natom+1 : iq*3*Crystal%natom) = phfrq
1288    phdispl_allq(1:2, 1:3, 1:Crystal%natom, 1:3*Crystal%natom, iq) = phdispl
1289    do jmode = 1, 3*Crystal%natom
1290      imode = imode + 1
1291      modeindex(imode) = imode
1292    end do
1293  end do
1294  ABI_FREE(phfrq)
1295  ABI_FREE(pheigvec)
1296  ABI_FREE(phdispl)
1297 
1298  ! sort modes in whole list: get indirect indexing for qbz and displ
1299  call sort_dp(nqibz*3*Crystal%natom, phfrq_allq, modeindex, tol10)
1300  ! NB: phfrq is sorted now, but displ and qibz will have to be indexed indirectly with modeindex
1301 
1302  ! only diagonal supercell case for the moment
1303  do itemper = 1, ntemper
1304    call init_supercell(Crystal%natom, rlatt, Crystal%rprimd, Crystal%typat, Crystal%xcart, Crystal%znucl, thm_scells(itemper))
1305  end do
1306 
1307  ! precalculate phase factors???
1308 
1309  ABI_STAT_MALLOC(phdispl1, (2, 3, Crystal%natom), ierr)
1310  ! for all modes at all q in whole list, sorted
1311  modesign=one
1312  do imode = 1, 3*Crystal%natom*nqibz
1313    ! skip modes with too low or negative frequency -> Bose factor explodes (eg acoustic at Gamma)
1314    if (phfrq_allq(imode) < tol10) cycle
1315 
1316    iq = ceiling(dble(modeindex(imode))/dble(3*Crystal%natom))
1317    jmode = modeindex(imode) - (iq-1)*3*Crystal%natom
1318    phdispl1 = phdispl_allq(:,:,:,jmode,iq)
1319 
1320    ! loop over temperatures
1321    do itemper = 1, ntemper
1322      temperature_K = tempermin + dble(itemper-1)*temperinc  ! this is in Kelvin
1323      temperature = temperature_K / Ha_K !=315774.65_dp
1324 
1325      ! trick supercell object into using present q point
1326      thm_scells(itemper)%qphon(:) = qibz(:,iq)
1327 
1328      ! find thermal displacement amplitude eq 4 of Zacharias
1329      !   combined with l_nu,q expression in paragraph before
1330      sigma = sqrt( (bose_einstein(phfrq_allq(imode), temperature) + half)/phfrq_allq(imode) )
1331 
1332      ! add displacement for this mode to supercell positions eq 5 of Zacharias
1333      freeze_displ = modesign * sigma
1334      call freeze_displ_supercell (phdispl1(:,:,:), freeze_displ, thm_scells(itemper))
1335    end do !itemper
1336 
1337    ! this is the prescription: flip sign for each successive mode in full
1338    ! spectrum, to cancel electron phonon coupling to 1st order
1339    ! (hopeflly 3rd order as well)
1340    modesign=-modesign
1341 
1342  end do !imode
1343 
1344  ABI_FREE(modeindex)
1345  ABI_FREE(phfrq_allq)
1346  ABI_FREE(phdispl_allq)
1347  ABI_FREE(phdispl1)
1348  ABI_FREE(qibz)
1349  ABI_FREE(qbz)
1350  ABI_FREE(wtq_ibz)
1351 
1352 end subroutine zacharias_supercell_make

m_phonons/zacharias_supercell_print [ Functions ]

[ Top ] [ m_phonons ] [ Functions ]

NAME

 zacharias_supercell_print

FUNCTION

  print files with thermal array of supercells

INPUTS

OUTPUT

NOTES

PARENTS

      anaddb

CHILDREN

      ifc_fourq,kpath_free,phonons_ncwrite,phonons_write_gnuplot
      phonons_write_phfrq,phonons_write_xmgrace,xmpi_sum_master

SOURCE

1657 subroutine zacharias_supercell_print(fname, ntemper, tempermin, temperinc, thm_scells)
1658 
1659 
1660 !This section has been created automatically by the script Abilint (TD).
1661 !Do not modify the following lines by hand.
1662 #undef ABI_FUNC
1663 #define ABI_FUNC 'zacharias_supercell_print'
1664 !End of the abilint section
1665 
1666  implicit none
1667 
1668 !Arguments ------------------------------------
1669 !scalars
1670  integer, intent(in) :: ntemper
1671  real(dp), intent(in) :: tempermin
1672  real(dp), intent(in) :: temperinc
1673  type(supercell_type), intent(in) :: thm_scells(ntemper)
1674  character(len=fnlen), intent(in) :: fname
1675 
1676 ! local
1677  integer :: itemp
1678  character(len=80) :: title1, title2
1679  character(len=fnlen) :: filename
1680  real(dp) :: temper
1681  character(len=10) :: temper_str
1682 
1683  do itemp = 1, ntemper
1684    temper = dble(itemp-1)*temperinc+tempermin
1685    write (temper_str,'(I8)') int(temper)
1686    write (filename, '(3a)') trim(fname), "_T_", trim(adjustl(temper_str))
1687    write (title1, '(3a)') "#  Zacharias thermalized supercell at temperature T= ", trim(temper_str), " Kelvin"
1688    title2 = "#  generated with alternating thermal displacements of all phonons"
1689    call prt_supercell (filename, thm_scells(itemp), title1, title2)
1690  end do
1691 
1692 end subroutine zacharias_supercell_print