TABLE OF CONTENTS
- ABINIT/m_phonons
- m_phonons/dfpt_symph
- m_phonons/freeze_displ_allmodes
- m_phonons/ifc_mkphbs
- m_phonons/mkphbs
- m_phonons/phdos_calc_vsound
- m_phonons/phdos_free
- m_phonons/phdos_init
- m_phonons/phdos_malloc
- m_phonons/phdos_ncwrite
- m_phonons/phdos_print
- m_phonons/phdos_print_msqd
- m_phonons/phdos_print_vsound
- m_phonons/phdos_t
- m_phonons/pheigvec_rotate
- m_phonons/phonons_ncwrite
- m_phonons/phonons_write_gnuplot
- m_phonons/phonons_write_phfrq
- m_phonons/phonons_write_xmgrace
- m_phonons/phstore_async_rotate
- m_phonons/phstore_free
- m_phonons/phstore_new
- m_phonons/phstore_t
- m_phonons/phstore_wait
- m_phonons/test_phrotation
- m_phonons/thermal_supercell_free
- m_phonons/thermal_supercell_make
- m_phonons/thermal_supercell_print
- m_phonons/zacharias_supercell_make
- m_phonons/zacharias_supercell_print
ABINIT/m_phonons [ 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 phdos_init
COPYRIGHT
Copyright (C) 1999-2024 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
19 #if defined HAVE_CONFIG_H 20 #include "config.h" 21 #endif 22 23 #include "abi_common.h" 24 25 module m_phonons 26 27 use defs_basis 28 use m_errors 29 use m_xmpi 30 use m_abicore 31 use m_htetra 32 use m_numeric_tools 33 use m_cgtools 34 use m_crystal 35 use m_nctk 36 use, intrinsic :: iso_c_binding 37 use m_atprj 38 use m_sortph 39 use m_ddb 40 use netcdf 41 use m_supercell 42 use m_dtset 43 use m_krank 44 45 use m_fstrings, only : itoa, ftoa, sjoin, ltoa, ktoa, strcat, basename, replace 46 use m_symtk, only : matr3inv 47 use m_time, only : cwtime, cwtime_report 48 use m_io_tools, only : open_file 49 use m_geometry, only : mkrdim, symredcart, normv, phdispl_cart2red 50 use m_dynmat, only : gtdyn9, dfpt_phfrq, dfpt_prtph, & 51 pheigvec_normalize, massmult_and_breaksym, & 52 phdispl_from_eigvec, phangmom_from_eigvec 53 use m_bz_mesh, only : isamek, make_path, kpath_t, kpath_new 54 use m_ifc, only : ifc_type 55 use m_anaddb_dataset, only : anaddb_dataset_type 56 use m_kpts, only : kpts_ibz_from_kptrlatt, get_full_kgrid, kpts_map, kpts_timrev_from_kptopt 57 use m_special_funcs, only : bose_einstein 58 use m_sort, only : sort_dp 59 use m_symfind, only : symanal 60 61 implicit none 62 63 private 64 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
SOURCE
2963 subroutine dfpt_symph(iout,acell,eigvec,indsym,natom,nsym,phfrq,rprim,symrel) 2964 2965 !Arguments ------------------------------------ 2966 !scalars 2967 integer,intent(in) :: iout,natom,nsym 2968 !arrays 2969 integer,intent(in) :: indsym(4,nsym,natom),symrel(3,3,nsym) 2970 real(dp),intent(in) :: acell(3),eigvec(2*3*natom*3*natom),phfrq(3*natom) 2971 real(dp),intent(in) :: rprim(3,3) 2972 2973 !Local variables ------------------------- 2974 !scalars 2975 integer :: iad1,iad2,iad3,iatom,idir,ii1,ii2,ii3,imode,isym,itol,jad,jatom,jj 2976 integer :: jmode,kk,ntol 2977 character(len=500) :: message 2978 !arrays 2979 integer,allocatable :: degeneracy(:),integer_characters(:),symind(:,:) 2980 real(dp) :: gprimd(3,3),rprimd(3,3) 2981 real(dp),allocatable :: eigvtr(:),redvec(:),redvtr(:),symph(:,:) 2982 2983 !****************************************************************** 2984 2985 ! Compute dimensional primitive translations rprimd and its inverse gprimd 2986 call mkrdim(acell,rprim,rprimd) 2987 call matr3inv(rprimd,gprimd) 2988 2989 ! Build the symmetry index (inverse of indsym(4,:,:)) 2990 ABI_MALLOC(symind, (nsym,natom)) 2991 do isym=1,nsym 2992 do iatom=1,natom 2993 symind(isym,indsym(4,isym,iatom))=iatom 2994 end do 2995 end do 2996 2997 ABI_MALLOC(symph,(nsym,3*natom)) 2998 ABI_MALLOC(redvec,(2*3*natom)) 2999 ABI_MALLOC(redvtr,(2*3*natom)) 3000 ABI_MALLOC(eigvtr,(2*3*natom)) 3001 3002 ! Loop over the vibration modes 3003 do imode=1,3*natom 3004 3005 ! Compute eigvec for this mode in reduced coordinates redvec 3006 do iatom=1,natom 3007 iad1=3*(iatom-1)+1 3008 ii1=2*3*natom*(imode-1)+2*(iad1-1)+1 3009 iad2=3*(iatom-1)+2 3010 ii2=2*3*natom*(imode-1)+2*(iad2-1)+1 3011 iad3=3*(iatom-1)+3 3012 ii3=2*3*natom*(imode-1)+2*(iad3-1)+1 3013 do idir=1,3 3014 jad=3*(iatom-1)+idir 3015 jj=2*(jad-1)+1 3016 redvec(jj)=gprimd(1,idir)*eigvec(ii1)+& 3017 gprimd(2,idir)*eigvec(ii2)+& 3018 gprimd(3,idir)*eigvec(ii3) 3019 redvec(jj+1)=gprimd(1,idir)*eigvec(ii1+1)+& 3020 gprimd(2,idir)*eigvec(ii2+1)+& 3021 gprimd(3,idir)*eigvec(ii3+1) 3022 end do !idir 3023 end do !iatom 3024 3025 ! Apply each transformation to redvec and store at the correct location in redvtr (iatom -> jatom) 3026 do isym=1,nsym 3027 do iatom=1,natom 3028 jatom=symind(isym,iatom) 3029 iad1=3*(iatom-1)+1 3030 ii1=2*(iad1-1)+1 3031 iad2=3*(iatom-1)+2 3032 ii2=2*(iad2-1)+1 3033 iad3=3*(iatom-1)+3 3034 ii3=2*(iad3-1)+1 3035 do idir=1,3 3036 jad=3*(jatom-1)+idir 3037 jj=2*(jad-1)+1 3038 redvtr(jj)=dble(symrel(idir,1,isym))*redvec(ii1)+& 3039 dble(symrel(idir,2,isym))*redvec(ii2)+& 3040 dble(symrel(idir,3,isym))*redvec(ii3) 3041 redvtr(jj+1)=dble(symrel(idir,1,isym))*redvec(ii1+1)+& 3042 dble(symrel(idir,2,isym))*redvec(ii2+1)+& 3043 dble(symrel(idir,3,isym))*redvec(ii3+1) 3044 end do !idir 3045 end do !iatom 3046 3047 ! Compute redvtr in cartesian coordinates eigvtr 3048 do iatom=1,natom 3049 iad1=3*(iatom-1)+1 3050 ii1=2*(iad1-1)+1 3051 iad2=3*(iatom-1)+2 3052 ii2=2*(iad2-1)+1 3053 iad3=3*(iatom-1)+3 3054 ii3=2*(iad3-1)+1 3055 do idir=1,3 3056 jad=3*(iatom-1)+idir 3057 jj=2*(jad-1)+1 3058 eigvtr(jj)=rprimd(idir,1)*redvtr(ii1)+& 3059 rprimd(idir,2)*redvtr(ii2)+& 3060 rprimd(idir,3)*redvtr(ii3) 3061 eigvtr(jj+1)=rprimd(idir,1)*redvtr(ii1+1)+& 3062 rprimd(idir,2)*redvtr(ii2+1)+& 3063 rprimd(idir,3)*redvtr(ii3+1) 3064 end do !idir 3065 end do !iatom 3066 3067 ! Compute scalar product... 3068 symph(isym,imode)=0.0_dp 3069 do jad=1,3*natom 3070 jj=2*(jad-1)+1 3071 kk=2*3*natom*(imode-1)+2*(jad-1)+1 3072 symph(isym,imode)=symph(isym,imode)+eigvtr(jj)*eigvec(kk)+eigvtr(jj+1)*eigvec(kk+1) 3073 end do 3074 3075 end do !isym 3076 end do !imode 3077 3078 !Treat degeneracies (different tolerances will be tried) 3079 !Compute the order of the degeneracy, and 3080 !attribute it to the lowest of the degenerate modes 3081 !Also attribute the characters to the lowest mode 3082 !When all the characters are integers, consider that the 3083 !mode is non-degenerate. The maximum difference in frequency 3084 !that is tolerated is on the order of 4cm-1 (which is large...) 3085 ABI_MALLOC(degeneracy, (3*natom)) 3086 ABI_MALLOC(integer_characters, (3*natom)) 3087 degeneracy(:)=1 3088 integer_characters(:)=0 3089 do itol=1,20 3090 ntol=itol 3091 do imode=3*natom,2,-1 3092 if(integer_characters(imode)==0)then 3093 do jmode=imode-1,1,-1 3094 if(integer_characters(jmode)==0)then 3095 if(abs(phfrq(imode)-phfrq(jmode))<itol*tol6)then 3096 degeneracy(jmode)=degeneracy(jmode)+degeneracy(imode) 3097 degeneracy(imode)=0 3098 symph(:,jmode)=symph(:,jmode)+symph(:,imode) 3099 symph(:,imode)=0.0_dp 3100 end if 3101 end if !integer_characters(jmode)==0 3102 end do !jmode 3103 end if !integer_characters(imode)==0 3104 end do !imode 3105 do imode=1,3*natom 3106 if(maxval(abs( symph(:,imode)-nint(symph(:,imode)) ))<0.05_dp)then 3107 integer_characters(imode)=1 3108 end if 3109 end do 3110 if(sum(integer_characters(:))==3*natom)exit 3111 end do !itol 3112 3113 !write(std_out,*)' dfpt_symph : degeneracy=',degeneracy(:) 3114 3115 write(message,'(a,a,es8.2,5a)')ch10,' Analysis of degeneracies and characters (maximum tolerance=',ntol*tol6,' a.u.)',ch10,& 3116 ' For each vibration mode, or group of modes if degenerate,',ch10,& 3117 ' the characters are given for each symmetry operation (see the list in the log file).' 3118 call wrtout([std_out, iout], message) 3119 3120 do imode=1,3*natom 3121 if(degeneracy(imode)/=0)then 3122 write(message,'(a,i4)') ' Symmetry characters of vibration mode #',imode 3123 call wrtout([std_out, iout], message) 3124 if(degeneracy(imode)>=2)then 3125 if(degeneracy(imode)==2) write(message,'(a,i4)') ' degenerate with vibration mode #',imode+1 3126 if(degeneracy(imode)>=3) write(message,'(a,i4,a,i4)') & 3127 ' degenerate with vibration modes #',imode+1,' to ',imode+degeneracy(imode)-1 3128 call wrtout([std_out, iout], message) 3129 end if 3130 do jj=1,(nsym-1)/16+1 3131 write(message,'(16f5.1)') (symph(isym,imode),isym=(jj-1)*16+1,min(nsym,jj*16)) 3132 call wrtout([std_out, iout], message) 3133 end do 3134 end if 3135 end do !imode 3136 3137 ABI_FREE(degeneracy) 3138 ABI_FREE(integer_characters) 3139 ABI_FREE(eigvtr) 3140 ABI_FREE(redvtr) 3141 ABI_FREE(redvec) 3142 ABI_FREE(symph) 3143 ABI_FREE(symind) 3144 3145 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.
SOURCE
3182 subroutine freeze_displ_allmodes(displ, freeze_displ, natom, outfile_radix, phfreq, & 3183 qphon, rprimd, typat, xcart, znucl) 3184 3185 !Arguments ------------------------------------ 3186 !scalars 3187 integer,intent(in) :: natom 3188 character(len=*),intent(in) :: outfile_radix 3189 real(dp), intent(in) :: freeze_displ 3190 !arrays 3191 integer,intent(in) :: typat(natom) 3192 real(dp),intent(in) :: displ(2,3*natom,3*natom) 3193 real(dp),intent(in) :: rprimd(3,3) 3194 real(dp),intent(in) :: phfreq(3*natom) 3195 real(dp),intent(in) :: qphon(3) 3196 real(dp),intent(in) :: xcart(3,natom) 3197 real(dp),intent(in) :: znucl(:) 3198 3199 ! local vars 3200 integer :: jmode 3201 type(supercell_type) :: scell 3202 3203 ! ************************************************************************* 3204 3205 !determine supercell needed to freeze phonon 3206 call init_supercell_for_qpt(natom, qphon, rprimd, typat, xcart, znucl, scell) 3207 3208 do jmode = 1, 3*natom 3209 ! reset positions 3210 scell%xcart = scell%xcart_ref 3211 3212 ! displace atoms according to phonon jmode 3213 call freeze_displ_supercell(displ(:,:,jmode), freeze_displ, scell) 3214 3215 ! print out everything for this wavevector and mode 3216 call prt_supercell_for_qpt (phfreq(jmode), jmode, outfile_radix, scell) 3217 end do 3218 3219 call destroy_supercell(scell) 3220 3221 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.
SOURCE
2812 subroutine ifc_mkphbs(ifc, cryst, dtset, prefix, comm) 2813 2814 !Arguments ------------------------------- 2815 !scalars 2816 integer,intent(in) :: comm 2817 character(len=*),intent(in) :: prefix 2818 type(ifc_type),intent(in) :: ifc 2819 type(crystal_t),intent(in) :: cryst 2820 type(dataset_type),intent(in) :: dtset 2821 2822 !Local variables ------------------------- 2823 !scalars 2824 integer,parameter :: master = 0 2825 integer :: iqpt, nqpts, natom, ncid, nprocs, my_rank, ierr, nph2l, ncerr 2826 type(kpath_t) :: qpath 2827 !arrays 2828 real(dp),allocatable :: qph2l(:,:), qnrml2(:), eigvec(:,:,:,:,:),phfrqs(:,:),phdispl_cart(:,:,:,:),phangmom(:,:,:),weights(:) 2829 2830 ! ********************************************************************* 2831 2832 if (dtset%prtphbands == 0) return 2833 2834 if (dtset%ph_nqpath <= 0 .or. dtset%ph_ndivsm <= 0) then 2835 ABI_COMMENT("ph_nqpath <= 0 or ph_ndivsm <= 0. Phonon bands won't be produced. Returning") 2836 return 2837 end if 2838 2839 call wrtout(std_out, " Writing phonon bands, use prtphbands 0 to disable this part") 2840 2841 nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 2842 natom = cryst%natom 2843 2844 qpath = kpath_new(dtset%ph_qpath(:,1:dtset%ph_nqpath), cryst%gprimd, dtset%ph_ndivsm) 2845 nqpts = qpath%npts 2846 2847 ABI_CALLOC(phfrqs, (3*natom,nqpts)) 2848 ABI_CALLOC(phdispl_cart, (2,3*natom,3*natom,nqpts)) 2849 ABI_CALLOC(phangmom, (3,3*natom,nqpts)) 2850 ABI_CALLOC(eigvec, (2,3,natom,3,natom)) 2851 2852 do iqpt=1,nqpts 2853 if (mod(iqpt, nprocs) /= my_rank) cycle ! MPI-parallelism 2854 ! Get phonon frequencies and displacements in cartesian coordinates for this q-point 2855 call ifc%fourq(cryst, qpath%points(:,iqpt), phfrqs(:,iqpt), phdispl_cart(:,:,:,iqpt), out_eigvec=eigvec) 2856 call phangmom_from_eigvec(natom, eigvec, phangmom(:,:,iqpt)) 2857 end do 2858 2859 call xmpi_sum_master(phfrqs, master, comm, ierr) 2860 call xmpi_sum_master(phdispl_cart, master, comm, ierr) 2861 call xmpi_sum_master(phangmom, master, comm, ierr) 2862 2863 if (my_rank == master) then 2864 ABI_MALLOC(weights, (nqpts)) 2865 weights = one 2866 2867 ! Compute directions for non-analytical behaviour. 2868 ! TODO: The same approach should be used in anaddb at the level of the parser. 2869 ABI_MALLOC(qph2l, (3, 2*dtset%ph_nqpath)) 2870 ABI_MALLOC(qnrml2, (2*dtset%ph_nqpath)) 2871 2872 nph2l = 0 2873 if (any(ifc%zeff /= zero)) then 2874 do iqpt=1,dtset%ph_nqpath 2875 if (sum(dtset%ph_qpath(:, iqpt)**2) < tol14) then 2876 nph2l = nph2l + 1 2877 if (iqpt == 1) then 2878 qph2l(:, nph2l) = dtset%ph_qpath(:, 2) - dtset%ph_qpath(:, 1) 2879 else if (iqpt == dtset%ph_nqpath) then 2880 qph2l(:, nph2l) = dtset%ph_qpath(:, dtset%ph_nqpath - 1) - dtset%ph_qpath(:, dtset%ph_nqpath) 2881 else 2882 qph2l(:, nph2l) = dtset%ph_qpath(:, iqpt - 1) - dtset%ph_qpath(:, iqpt) 2883 nph2l = nph2l + 1 2884 qph2l(:, nph2l) = dtset%ph_qpath(:, iqpt + 1) - dtset%ph_qpath(:, iqpt) 2885 end if 2886 end if 2887 end do 2888 2889 ! Convert to Cartesian coordinates. 2890 do iqpt=1,nph2l 2891 qph2l(:, iqpt) = matmul(cryst%gprimd, qph2l(:, iqpt)) 2892 end do 2893 qnrml2 = zero 2894 end if 2895 2896 ! TODO: A similar piece of code is used in anaddb (mkpbs + ifc_calcnwrite_nana_terms). 2897 ! Should centralize everything in a single routine 2898 NCF_CHECK_MSG(nctk_open_create(ncid, strcat(prefix, "_PHBST.nc"), xmpi_comm_self), "Creating PHBST") 2899 NCF_CHECK(cryst%ncwrite(ncid)) 2900 call phonons_ncwrite(ncid, natom, nqpts, qpath%points, weights, phfrqs, phdispl_cart, phangmom) 2901 ! This flag tells AbiPy that all the non-analytic directions have been computed. 2902 NCF_CHECK(nctk_defnwrite_ivars(ncid, ["has_abipy_non_anal_ph"], [1])) 2903 ncerr = nctk_def_arrays(ncid, [nctkarr_t("atomic_mass_units", "dp", "number_of_atom_species")], defmode=.True.) 2904 NCF_CHECK(ncerr) 2905 NCF_CHECK(nctk_set_datamode(ncid)) 2906 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "atomic_mass_units"), ifc%amu)) 2907 if (nph2l /= 0) call ifc%calcnwrite_nana_terms(cryst, nph2l, qph2l, qnrml2, ncid=ncid) 2908 NCF_CHECK(nf90_close(ncid)) 2909 2910 ABI_FREE(qph2l) 2911 ABI_FREE(qnrml2) 2912 2913 select case (dtset%prtphbands) 2914 case (1) 2915 call phonons_write_xmgrace(strcat(prefix, "_PHBANDS.agr"), natom, nqpts, qpath%points, phfrqs, qptbounds=qpath%bounds) 2916 case (2) 2917 call phonons_write_gnuplot(prefix, natom, nqpts, qpath%points, phfrqs, qptbounds=qpath%bounds) 2918 case (3) 2919 call phonons_write_phfrq(prefix, natom, nqpts, qpath%points, weights, phfrqs, phdispl_cart, phangmom) 2920 case default 2921 ABI_WARNING(sjoin("Unsupported value for prtphbands:", itoa(dtset%prtphbands))) 2922 end select 2923 2924 ABI_FREE(weights) 2925 end if ! master 2926 2927 ABI_FREE(phfrqs) 2928 ABI_FREE(phdispl_cart) 2929 ABI_FREE(phangmom) 2930 ABI_FREE(eigvec) 2931 2932 call qpath%free() 2933 2934 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.
SOURCE
1796 subroutine mkphbs(Ifc,Crystal,inp,ddb,asrq0,prefix,comm) 1797 1798 !Arguments ------------------------------- 1799 !scalars 1800 integer,intent(in) :: comm 1801 character(len=*),intent(in) :: prefix 1802 type(ifc_type),intent(in) :: Ifc 1803 type(crystal_t),intent(in) :: Crystal 1804 type(anaddb_dataset_type),target,intent(in) :: inp 1805 type(ddb_type),intent(in) :: ddb 1806 type(asrq0_t),intent(inout) :: asrq0 1807 1808 !Local variables ------------------------- 1809 !scalars 1810 integer,parameter :: master=0 1811 integer :: unt 1812 integer :: iphl1,iblok,rftyp, ii,nfineqpath,nsym,natom,ncid,nprocs,my_rank 1813 integer :: natprj_bs,eivec,enunit,ifcflag,ptgroupma,spgroup 1814 real(dp) :: freeze_displ 1815 real(dp) :: cfact 1816 character(500) :: msg 1817 character(len=8) :: unitname 1818 !arrays 1819 integer :: bravais(11),rfphon(4),rfelfd(4),rfstrs(4) 1820 integer :: nomega, imode, iomega 1821 integer,allocatable :: ndiv(:) 1822 real(dp) :: omega, omega_min, gaussmaxarg, gaussfactor, gaussprefactor, xx 1823 real(dp) :: speedofsound(3),genafm(3) 1824 real(dp) :: qphnrm(3), qphon(3), qphon_padded(3,3),res(3) 1825 real(dp) :: d2cart(2,ddb%msize),real_qphon(3) 1826 real(dp) :: displ(2*3*Crystal%natom*3*Crystal%natom),eigval(3,Crystal%natom),phangmom(3,3*Crystal%natom) 1827 real(dp),allocatable :: phfrq(:),eigvec(:,:,:,:,:) 1828 real(dp),allocatable :: save_phfrq(:,:),save_phdispl_cart(:,:,:,:),save_qpoints(:,:),save_phangmom(:,:,:) 1829 real(dp),allocatable :: weights(:) 1830 real(dp),allocatable :: dos4bs(:) 1831 real(dp),allocatable,target :: alloc_path(:,:) 1832 real(dp),pointer :: fineqpath(:,:) 1833 type(atprj_type) :: atprj 1834 1835 ! ********************************************************************* 1836 1837 ! Only master works for the time being 1838 nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 1839 if (my_rank /= master) return 1840 1841 nsym = Crystal%nsym; natom = Crystal%natom 1842 1843 ! Copy parameters from inp (then I will try to remove inp from the API so that I can call mkphbs in eph) 1844 ifcflag = inp%ifcflag 1845 natprj_bs = inp%natprj_bs 1846 freeze_displ = inp%freeze_displ 1847 eivec = inp%eivec; enunit = inp%enunit 1848 rftyp=inp%rfmeth 1849 1850 nullify(fineqpath) 1851 nfineqpath = inp%nph1l 1852 fineqpath => inp%qph1l 1853 1854 if(inp%nph1l==0) then 1855 if (inp%nqpath==0) then 1856 return ! if there is nothing to do, return 1857 else 1858 ! allow override of nph1l with nqpath if the former is not set 1859 ! allocate and compute path here and make fineqpath points to it 1860 ABI_MALLOC(ndiv,(inp%nqpath-1)) 1861 call make_path(inp%nqpath,inp%qpath,Crystal%gmet,'G',inp%ndivsm,ndiv,nfineqpath,alloc_path,std_out) 1862 ABI_FREE(ndiv) 1863 fineqpath => alloc_path 1864 end if 1865 end if 1866 1867 write(msg, '(a,(80a),a,a,a,a)' ) ch10,('=',ii=1,80),ch10,ch10,' Treat the first list of vectors ',ch10 1868 call wrtout([std_out, ab_out], msg) 1869 1870 if (natprj_bs > 0) call atprj_init(atprj, natom, natprj_bs, inp%iatprj_bs, prefix) 1871 1872 ABI_MALLOC(phfrq, (3*natom)) 1873 ABI_MALLOC(eigvec, (2,3,natom,3,natom)) 1874 ABI_MALLOC(save_qpoints, (3,nfineqpath)) 1875 ABI_MALLOC(save_phfrq, (3*natom,nfineqpath)) 1876 ABI_MALLOC(save_phdispl_cart, (2,3*natom,3*natom,nfineqpath)) 1877 ABI_MALLOC(save_phangmom, (3,3*natom,nfineqpath)) 1878 qphnrm = one 1879 1880 do iphl1=1,nfineqpath 1881 1882 ! Initialisation of the phonon wavevector 1883 qphon(:)=fineqpath(:,iphl1) 1884 1885 if (inp%nph1l /= 0) qphnrm(1) = inp%qnrml1(iphl1) 1886 1887 save_qpoints(:,iphl1) = qphon / qphnrm(1) 1888 1889 ! Generation of the dynamical matrix in cartesian coordinates 1890 if (ifcflag == 1) then 1891 1892 ! Get phonon frequencies and displacements in reduced coordinates for this q-point 1893 !call ifc%fourq(cryst, save_qpoints(:,iphl1), phfrq, displ, out_eigvec=eigvec) 1894 1895 ! Get d2cart using the interatomic forces and the 1896 ! long-range coulomb interaction through Ewald summation 1897 call gtdyn9(ddb%acell,Ifc%atmfrc,Ifc%dielt,Ifc%dipdip,Ifc%dyewq0,d2cart,Crystal%gmet,ddb%gprim,ddb%mpert,natom, & 1898 Ifc%nrpt,qphnrm(1),qphon,Crystal%rmet,ddb%rprim,Ifc%rpt,Ifc%trans,Crystal%ucvol,Ifc%wghatm,Crystal%xred,ifc%zeff,& 1899 ifc%qdrp_cart,ifc%ewald_option,xmpi_comm_self,dipquad=Ifc%dipquad,quadquad=Ifc%quadquad) 1900 1901 else if (ifcflag == 0) then 1902 1903 !call ddb_diagoq(ddb, crystal, save_qpoints(:,iphl1), asrq0, ifc%symdynmat, rftyp, phfrq, displ, & 1904 ! out_eigvec=eigvec) 1905 1906 ! Look for the information in the DDB (no interpolation here!) 1907 rfphon(1:2)=1; rfelfd(1:2)=0; rfstrs(1:2)=0 1908 qphon_padded = zero; qphon_padded(:,1) = qphon 1909 1910 call ddb%get_block(iblok,qphon_padded,qphnrm,rfphon,rfelfd,rfstrs,rftyp) 1911 1912 ! Copy the dynamical matrix in d2cart 1913 d2cart(:,1:ddb%msize)=ddb%val(:,:,iblok) 1914 1915 ! Eventually impose the acoustic sum rule based on previously calculated d2asr 1916 call asrq0%apply(natom, ddb%mpert, ddb%msize, crystal%xcart, d2cart) 1917 end if 1918 1919 ! Use inp%symdynmat instead of ifc because of ifcflag 1920 ! Calculation of the eigenvectors and eigenvalues of the dynamical matrix 1921 call dfpt_phfrq(ddb%amu,displ,d2cart,eigval,eigvec,Crystal%indsym,& 1922 ddb%mpert,Crystal%nsym,natom,nsym,Crystal%ntypat,phfrq,qphnrm(1),qphon,& 1923 crystal%rprimd,inp%symdynmat,Crystal%symrel,Crystal%symafm,Crystal%typat,Crystal%ucvol) 1924 ! Calculation of the phonon angular momentum 1925 ! maybe add it in dpft_phfrq directly ? 1926 call phangmom_from_eigvec(natom, eigvec, phangmom) 1927 1928 if (abs(freeze_displ) > tol10) then 1929 real_qphon = zero 1930 if (abs(qphnrm(1)) > tol8) real_qphon = qphon / qphnrm(1) 1931 call freeze_displ_allmodes(displ, freeze_displ, natom, prefix, phfrq, & 1932 real_qphon, crystal%rprimd, Crystal%typat, crystal%xcart, crystal%znucl) 1933 end if 1934 1935 ! If requested, output projection of each mode on given atoms 1936 if (natprj_bs > 0) call atprj_print(atprj, iphl1, phfrq, eigvec) 1937 1938 ! In case eivec == 4, write output files for band2eps (visualization of phonon band structures) 1939 if (eivec == 4) then 1940 call sortph(eigvec,displ,strcat(prefix, "_B2EPS"),natom,phfrq) 1941 end if 1942 1943 ! Write the phonon frequencies 1944 call dfpt_prtph(displ,eivec,enunit,ab_out,natom,phfrq,qphnrm(1),qphon) 1945 1946 save_phfrq(:,iphl1) = phfrq 1947 save_phdispl_cart(:,:,:,iphl1) = RESHAPE(displ, [2, 3*natom, 3*natom]) 1948 save_phangmom(:,:,iphl1) = RESHAPE(phangmom, [3, 3*natom]) 1949 1950 ! Determine the symmetries of the phonon mode at Gamma 1951 ! TODO: generalize for other q-point little groups. 1952 if (sum(abs(qphon)) < DDB_QTOL) then 1953 call symanal(bravais,0,genafm,nsym,nsym,ptgroupma,Crystal%rprimd,spgroup, & 1954 Crystal%symafm,Crystal%symrel,Crystal%tnons,tol5,verbose=.TRUE.) 1955 call dfpt_symph(ab_out,ddb%acell,eigvec,Crystal%indsym,natom,nsym,phfrq,ddb%rprim,Crystal%symrel) 1956 end if 1957 1958 ! if we have an acoustic mode (small q and acoustic type displacements) 1959 ! extrapolate speed of sound in this direction, and Debye frequency 1960 call wrap2_pmhalf(qphon, real_qphon, res) 1961 if (sqrt(real_qphon(1)**2+real_qphon(2)**2+real_qphon(3)**2) < quarter .and. & 1962 sqrt(real_qphon(1)**2+real_qphon(2)**2+real_qphon(3)**2) > tol6) then 1963 call phdos_calc_vsound(eigvec, Crystal%gmet, natom, phfrq, real_qphon, speedofsound) 1964 if (my_rank == master) call phdos_print_vsound(ab_out, Crystal%ucvol, speedofsound) 1965 end if 1966 1967 end do ! iphl1 1968 1969 ! calculate dos for the specific q points along the BS calculated 1970 ! only Gaussians are possible - no interpolation 1971 omega_min = minval(save_phfrq(:,:)) 1972 nomega=NINT( (maxval(save_phfrq(:,:))-omega_min) / inp%dosdeltae ) + 1 1973 nomega=MAX(6,nomega) ! Ensure Simpson integration will be ok 1974 1975 ABI_MALLOC(dos4bs,(nomega)) 1976 dos4bs = zero 1977 gaussmaxarg = sqrt(-log(1.d-90)) 1978 gaussprefactor = one/(inp%dossmear*sqrt(two_pi)) 1979 gaussfactor = one/(sqrt2*inp%dossmear) 1980 do iphl1=1,nfineqpath 1981 do imode=1,3*natom 1982 do iomega=1, nomega 1983 omega = omega_min + (iomega-1) * inp%dosdeltae 1984 xx = (omega - save_phfrq(imode,iphl1)) * gaussfactor 1985 if(abs(xx) < gaussmaxarg) then 1986 dos4bs(iomega) = dos4bs(iomega) + gaussprefactor*exp(-xx*xx) 1987 end if 1988 end do 1989 end do 1990 end do 1991 1992 1993 !deallocate sortph array 1994 call end_sortph() 1995 1996 if (natprj_bs > 0) call atprj_destroy(atprj) 1997 1998 ! WRITE OUT FILES 1999 if (my_rank == master) then 2000 ABI_MALLOC(weights, (nfineqpath)) 2001 weights = one 2002 2003 NCF_CHECK_MSG(nctk_open_create(ncid, strcat(prefix, "_PHBST.nc"), xmpi_comm_self), "Creating PHBST") 2004 NCF_CHECK(crystal%ncwrite(ncid)) 2005 call phonons_ncwrite(ncid,natom,nfineqpath,save_qpoints,weights,save_phfrq,save_phdispl_cart,save_phangmom) 2006 2007 ! Now treat the second list of vectors (only at the Gamma point, but can include non-analyticities) 2008 if (inp%nph2l /= 0 .and. inp%ifcflag == 1) then 2009 call ifc%calcnwrite_nana_terms(crystal, inp%nph2l, inp%qph2l, inp%qnrml2, ncid) 2010 end if 2011 NCF_CHECK(nf90_close(ncid)) 2012 2013 call phonons_write_phfrq(prefix, natom,nfineqpath,save_qpoints,weights,save_phfrq,save_phdispl_cart, save_phangmom) 2014 2015 select case (inp%prtphbands) 2016 case (0) 2017 continue 2018 2019 case (1) 2020 if (inp%nph1l == 0) then 2021 call phonons_write_xmgrace(strcat(prefix, "_PHBANDS.agr"), natom, nfineqpath, save_qpoints, save_phfrq, & 2022 qptbounds=inp%qpath) 2023 else 2024 call phonons_write_xmgrace(strcat(prefix, "_PHBANDS.agr"), natom, nfineqpath, save_qpoints, save_phfrq) 2025 end if 2026 2027 case (2) 2028 if (inp%nph1l == 0) then 2029 call phonons_write_gnuplot(prefix, natom, nfineqpath, save_qpoints, save_phfrq, qptbounds=inp%qpath) 2030 else 2031 call phonons_write_gnuplot(prefix, natom, nfineqpath, save_qpoints, save_phfrq) 2032 end if 2033 2034 case default 2035 ABI_WARNING(sjoin("Don't know how to handle prtphbands:", itoa(inp%prtphbands))) 2036 end select 2037 2038 ! write out DOS file for q along this path 2039 cfact=one 2040 unitname = 'Ha' 2041 if (open_file('PHBST_partial_DOS',msg,newunit=unt,form="formatted",action="write") /= 0) then 2042 ABI_ERROR(msg) 2043 end if 2044 write(msg,'(3a)')'# ',ch10,'# Partial phonon density of states for q along a band structure path' 2045 call wrtout(unt, msg) 2046 write(msg,'(6a)')'# ',ch10,'# Energy in ',unitname,', DOS in states/',unitname 2047 call wrtout(unt, msg) 2048 write(msg,'(a,E20.10,2a,i8)') '# Gaussian method with smearing = ',inp%dossmear*cfact,unitname, ', nq =', nfineqpath 2049 call wrtout(unt, msg) 2050 write(msg,'(5a)')'# ',ch10,'# omega PHDOS ',ch10,'# ' 2051 call wrtout(unt, msg) 2052 do iomega=1,nomega 2053 omega = omega_min + (iomega-1) * inp%dosdeltae 2054 write(unt,'(2es17.8)',advance='NO')omega*cfact,dos4bs(iomega)/cfact 2055 write(unt,*) 2056 end do 2057 close(unt) 2058 2059 ABI_FREE(weights) 2060 end if 2061 2062 ABI_FREE(save_qpoints) 2063 ABI_FREE(save_phfrq) 2064 ABI_FREE(save_phdispl_cart) 2065 ABI_FREE(save_phangmom) 2066 ABI_FREE(phfrq) 2067 ABI_FREE(eigvec) 2068 ABI_FREE(dos4bs) 2069 ABI_SFREE(alloc_path) 2070 2071 end subroutine mkphbs
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
SOURCE
2094 subroutine phdos_calc_vsound(eigvec,gmet,natom,phfrq,qphon,speedofsound) 2095 2096 !Arguments ------------------------------- 2097 !scalras 2098 integer, intent(in) :: natom 2099 !arrays 2100 real(dp), intent(in) :: gmet(3,3),qphon(3) 2101 real(dp), intent(in) :: phfrq(3*natom),eigvec(2,3*natom,3*natom) 2102 real(dp), intent(out) :: speedofsound(3) 2103 2104 !Local variables ------------------------- 2105 integer :: iatref,imode, iatom, isacoustic, imode_acoustic 2106 ! character(len=500) :: msg 2107 real(dp) :: qnormcart 2108 real(dp) :: qtmp(3) 2109 2110 ! ********************************************************************* 2111 2112 imode_acoustic = 0 2113 do imode = 1, 3*natom 2114 2115 ! Check if this mode is acoustic like: scalar product of all displacement vectors are collinear 2116 isacoustic = 1 2117 ! Find reference atom with non-zero displacement 2118 do iatom=1,natom 2119 if(sum(eigvec(:,(iatom-1)*3+1:(iatom-1)*3+3,imode)**2) >tol16)iatref=iatom 2120 enddo 2121 ! Now compute scalar product, and check they are all positive 2122 do iatom = 1, natom 2123 if (sum(eigvec(:,(iatom-1)*3+1:(iatom-1)*3+3, imode)& 2124 *eigvec(:,(iatref-1)*3+1:(iatref-1)*3+3, imode)) < tol16 ) isacoustic = 0 2125 end do 2126 if (isacoustic == 0) cycle 2127 imode_acoustic = min(imode_acoustic + 1, 3) 2128 2129 ! write (msg, '(a,I6,a,3F12.4)') ' Found acoustic mode ', imode, ' for |q| in red coord < 0.25 ; q = ', qphon 2130 ! call wrtout(std_out, msg) 2131 qtmp = matmul(gmet, qphon) 2132 qnormcart = two * pi * sqrt(sum(qphon*qtmp)) 2133 speedofsound(imode_acoustic) = phfrq(imode) / qnormcart 2134 end do 2135 2136 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
SOURCE
595 subroutine phdos_free(PHdos) 596 597 !Arguments ------------------------------- 598 class(phdos_t),intent(inout) ::PHdos 599 600 ! ************************************************************************* 601 602 !@phdos_t 603 ABI_SFREE(PHdos%atom_mass) 604 ABI_SFREE(PHdos%omega) 605 ABI_SFREE(PHdos%phdos) 606 ABI_SFREE(PHdos%phdos_int) 607 ABI_SFREE(PHdos%pjdos) 608 ABI_SFREE(PHdos%pjdos_int) 609 ABI_SFREE(PHdos%pjdos_type) 610 ABI_SFREE(PHdos%pjdos_type_int) 611 ABI_SFREE(PHdos%pjdos_rc_type) 612 ABI_SFREE(PHdos%msqd_dos_atom) 613 614 end subroutine phdos_free
m_phonons/phdos_init [ Functions ]
[ Top ] [ m_phonons ] [ Functions ]
NAME
phdos_init
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 tetrahedra.
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. prefix=Prefix for PHBIZ output file. Empty string to deactivate output. comm=MPI communicator.
OUTPUT
phdos<phdos_t>=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
SOURCE
721 subroutine phdos_init(phdos, crystal, ifc, prtdos, dosdeltae_in, dossmear, dos_ngqpt, nqshft, dos_qshift, prefix, & 722 wminmax, count_wminmax, comm, dos_maxmode) 723 724 !Arguments ------------------------------- 725 !scalars 726 class(phdos_t),intent(out) :: phdos 727 integer,intent(in) :: prtdos,nqshft,comm 728 real(dp),intent(in) :: dosdeltae_in,dossmear 729 character(len=*),intent(in) :: prefix 730 type(crystal_t),intent(in) :: crystal 731 type(ifc_type),intent(in) :: ifc 732 integer, optional, intent(in) :: dos_maxmode 733 !arrays 734 integer,intent(in) :: dos_ngqpt(3) 735 integer,intent(out) :: count_wminmax(2) 736 real(dp),intent(in) :: dos_qshift(3,nqshft) 737 real(dp),intent(inout) :: wminmax(2) 738 739 !Local variables ------------------------- 740 !scalars 741 integer,parameter :: bcorr0 = 0, master = 0 742 integer :: iat,jat,idir,imode,io,iq_ibz,itype, my_qptopt, my_nsym 743 integer :: nqbz,ierr,natom,nomega,jdir, isym, nprocs, my_rank, ncid 744 logical :: refine_dosdeltae 745 real(dp),parameter :: max_occ1=one, gaussmaxarg = sqrt(-log(1.d-90)), max_smallq = 0.0625_dp 746 real(dp) :: nsmallq,gaussfactor,gaussprefactor,normq,debyefreq,rtmp 747 real(dp) :: cpu, wall, gflops, cpu_all, wall_all, gflops_all 748 real(dp) :: dosdeltae, phdos_int 749 character(len=500) :: msg 750 character(len=80) :: errstr 751 type(htetra_t) :: htetraq 752 !arrays 753 integer :: in_qptrlatt(3,3),new_qptrlatt(3,3) 754 integer :: dos_maxmode_ 755 integer,allocatable :: bz2ibz_smap(:,:), bz2ibz(:) 756 real(dp) :: speedofsound(3),speedofsound_(3) 757 real(dp) :: displ(2*3*Crystal%natom*3*Crystal%natom) 758 real(dp) :: eigvec(2,3,Crystal%natom,3*Crystal%natom),phfrq(3*Crystal%natom),phangmom(3,3*Crystal%natom) 759 real(dp) :: qlatt(3,3),rlatt(3,3), msqd_atom_tmp(3,3),temp_33(3,3) 760 real(dp) :: symcart(3,3,crystal%nsym), syme2_xyza(3, crystal%natom) 761 real(dp),allocatable :: full_eigvec(:,:,:,:,:),full_phfrq(:,:),full_phangmom(:,:,:),new_shiftq(:,:) 762 real(dp),allocatable :: qbz(:,:),qibz(:,:),tmp_phfrq(:) !, work_msqd(:,:,:,:) 763 real(dp),allocatable :: wtq_ibz(:),xvals(:), gvals_wtq(:), wdt(:,:), energies(:) 764 765 ! ********************************************************************* 766 767 DBG_ENTER("COLL") 768 769 nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 770 771 ! Consistency check. 772 if (all(prtdos /= [1, 2])) then 773 ABI_BUG(sjoin('prtdos should be 1 or 2, but received', itoa(prtdos))) 774 end if 775 dosdeltae = dosdeltae_in; refine_dosdeltae = .false. 776 if (dosdeltae <= zero) then 777 dosdeltae = -dosdeltae; refine_dosdeltae = .true. 778 ABI_CHECK(nprocs == 1, "refine_dosdeltae cannot be used with nprocs > 1") 779 end if 780 if (prtdos == 1 .and. dossmear <= zero) then 781 ABI_BUG(sjoin('dossmear should be positive but received', ftoa(dossmear))) 782 end if 783 784 call cwtime(cpu_all, wall_all, gflops_all, "start") 785 786 ! Get symmetries in cartesian coordinates 787 do isym=1,crystal%nsym 788 call symredcart(crystal%rprimd,crystal%gprimd,symcart(:,:,isym),crystal%symrel(:,:,isym)) 789 end do 790 791 natom = crystal%natom 792 call phdos_malloc(phdos, crystal, ifc, dosdeltae, dossmear, wminmax, prtdos) 793 nomega = phdos%nomega 794 795 ABI_MALLOC(gvals_wtq, (nomega)) 796 ABI_MALLOC(xvals, (nomega)) 797 798 ! Parameters defining the gaussian approximant. 799 if (prtdos == 1) then 800 ! TODO: use gaussian and update reference files. 801 gaussprefactor = one / (dossmear * sqrt(two_pi)) 802 gaussfactor = one / (sqrt2 * dossmear) 803 write(msg, '(4a,f8.5,2a,f8.5,a,i0)') ch10, & 804 ' phdos_init: calculating phonon DOS using gaussian method:', ch10, & 805 ' gaussian smearing [meV] = ', dossmear * Ha_meV, ch10, & 806 ' frequency step [meV] = ', phdos%omega_step * Ha_meV, ", nomega = ",phdos%nomega 807 else if (prtdos == 2) then 808 write(msg, '(4a,f8.5,a,i0)') ch10, & 809 ' phdos_init: calculating phonon DOS using tetrahedron method:', ch10, & 810 ' frequency step [meV] = ',phdos%omega_step * Ha_meV, ", nomega = ",phdos%nomega 811 end if 812 call wrtout(std_out, msg) 813 814 ! This call will set %nqibz and IBZ and BZ arrays 815 in_qptrlatt = 0; in_qptrlatt(1, 1) = dos_ngqpt(1); in_qptrlatt(2, 2) = dos_ngqpt(2); in_qptrlatt(3, 3) = dos_ngqpt(3) 816 817 my_qptopt = 1 818 !my_qptopt = 3 ! This to deactivate the use of symmetries for debugging purposes. 819 call kpts_ibz_from_kptrlatt(crystal, in_qptrlatt, my_qptopt, nqshft, dos_qshift, & 820 phdos%nqibz, qibz, wtq_ibz, nqbz, qbz, new_kptrlatt=new_qptrlatt, new_shiftk=new_shiftq, bz2ibz=bz2ibz_smap) 821 822 phdos%qptrlatt = new_qptrlatt 823 phdos%shiftq(:) = new_shiftq(:, 1) ! only one shift in output 824 825 if (my_rank == master) then 826 write(msg, "(3a, i0)")" DOS ngqpt: ", trim(ltoa(dos_ngqpt)), ", qptopt: ", my_qptopt 827 call wrtout(std_out, msg) 828 write(msg, "(2(a, i0))")" Number of q-points in the IBZ: ", phdos%nqibz, ", number of MPI processes: ", nprocs 829 call wrtout(std_out, msg) 830 end if 831 !call cwtime_report(" kpts_ibz_from_kptrlatt", cpu, wall, gflops) 832 833 if (prtdos == 2) then 834 ! Prepare tetrahedron method including workspace arrays. 835 ! Convert kptrlatt to double and invert, qlatt here refer to the shortest qpt vectors 836 rlatt = new_qptrlatt; call matr3inv(rlatt, qlatt) 837 838 ABI_MALLOC(bz2ibz, (nqbz)) 839 bz2ibz = bz2ibz_smap(1,:) 840 841 call htetra_init(htetraq, bz2ibz, crystal%gprimd, qlatt, qbz, nqbz, qibz, phdos%nqibz, ierr, errstr, comm) 842 !call cwtime_report(" init_tetra", cpu, wall, gflops) 843 ABI_CHECK(ierr == 0, errstr) 844 ABI_FREE(bz2ibz) 845 846 ! Allocate arrays used to store the entire spectrum, Required to calculate tetra weights. 847 ! this may change in the future if Matteo refactorizes the tetra weights as sums over k instead of sums over bands 848 ABI_CALLOC(full_phfrq, (3*natom, phdos%nqibz)) 849 ABI_MALLOC(full_phangmom, (3, 3*natom, phdos%nqibz)) 850 ABI_MALLOC_OR_DIE(full_eigvec, (2, 3, natom, 3*natom, phdos%nqibz), ierr) 851 full_eigvec = zero 852 end if ! tetra 853 854 ABI_SFREE(bz2ibz_smap) 855 ABI_FREE(new_shiftq) 856 857 ! MPI Sum over irreducible q-points then sync the following integrals: 858 ! speedofsound, nsmallq 859 ! wminmax and count_wminmax 860 ! if gauss: %phdos, %msqd_dos_atom 861 ! if tetra: full_phfrq, full_eigvec, full_phangmom, %phdos_int 862 863 nsmallq = zero; speedofsound = zero 864 wminmax = [huge(one), -huge(one)]; count_wminmax = 0 865 866 do iq_ibz=1,phdos%nqibz 867 if (mod(iq_ibz, nprocs) /= my_rank) cycle ! mpi-parallelism 868 869 ! Fourier interpolation (keep track of min/max to decide if initial mesh was large enough) 870 call ifc%fourq(crystal, qibz(:,iq_ibz), phfrq, displ, out_eigvec=eigvec) 871 wminmax(1) = min(wminmax(1), minval(phfrq)) 872 if (wminmax(1) < phdos%omega(1)) count_wminmax(1) = count_wminmax(1) + 1 873 wminmax(2) = max(wminmax(2), maxval(phfrq)) 874 if (wminmax(2) > phdos%omega(nomega)) count_wminmax(2) = count_wminmax(2) + 1 875 876 normq = sum(qibz(:,iq_ibz) ** 2) 877 if (normq < max_smallq .and. normq > tol6) then 878 call phdos_calc_vsound(eigvec, crystal%gmet, natom, phfrq, qibz(:,iq_ibz), speedofsound_) 879 speedofsound = speedofsound + speedofsound_ * wtq_ibz(iq_ibz) 880 nsmallq = nsmallq + wtq_ibz(iq_ibz) 881 end if 882 883 dos_maxmode_ = 3*natom 884 if (present(dos_maxmode)) then 885 if (dos_maxmode > 0 .and. dos_maxmode < 3*natom) then 886 dos_maxmode_ = dos_maxmode 887 end if 888 end if 889 890 select case (prtdos) 891 case (1) 892 do imode=1,dos_maxmode_ 893 ! Precompute \delta(w - w_{qnu}) * weight(q) 894 xvals = (phdos%omega(:) - phfrq(imode)) * gaussfactor 895 where (abs(xvals) < gaussmaxarg) 896 gvals_wtq = gaussprefactor * exp(-xvals*xvals) * wtq_ibz(iq_ibz) 897 elsewhere 898 gvals_wtq = zero 899 end where 900 901 ! Accumulate PHDOS 902 phdos%phdos(:) = phdos%phdos(:) + gvals_wtq 903 904 ! Rotate e(q) to get e(Sq) to account for symmetrical q-points in BZ. 905 ! eigenvectors indeed are not invariant under rotation. See e.g. Eq 39-40 of PhysRevB.76.165108 [[cite:Giustino2007]]. 906 ! In principle there's a phase due to nonsymmorphic translations but we here need |e(Sq)_iatom|**2 907 syme2_xyza = zero 908 my_nsym = crystal%nsym; if (my_qptopt == 3) my_nsym = 1 909 do iat=1,natom 910 do isym=1,my_nsym 911 jat = crystal%indsym(4,isym,iat) 912 syme2_xyza(:,jat) = syme2_xyza(:,jat) + & 913 matmul(symcart(:,:,isym), eigvec(1,:,iat,imode)) ** 2 + & 914 matmul(symcart(:,:,isym), eigvec(2,:,iat,imode)) ** 2 915 end do 916 end do 917 syme2_xyza = syme2_xyza / my_nsym 918 919 ! Accumulate PJDOS 920 do iat=1,natom 921 do idir=1,3 922 phdos%pjdos(:,idir,iat) = phdos%pjdos(:,idir,iat) + syme2_xyza(idir,iat) * gvals_wtq 923 end do 924 end do 925 926 ! Accumulate outer product of displacement vectors 927 ! NB: only accumulate real part. e(-q) = e(q)* the full sum over the BZ guarantees Im=0 928 ! this sum only does irreducible points: the matrix is symmetrized below 929 ! msqd_atom_tmp has units of bohr^2 / Ha as gaussval ~ 1/smear ~ 1/Ha 930 do iat=1,natom 931 msqd_atom_tmp = zero 932 do idir=1,3 933 do jdir=1,3 934 msqd_atom_tmp(jdir,idir) = msqd_atom_tmp(jdir,idir) + ( & 935 eigvec(1,idir,iat,imode) * eigvec(1,jdir,iat,imode) & 936 + eigvec(2,idir,iat,imode) * eigvec(2,jdir,iat,imode) ) 937 end do 938 end do 939 940 ! Symmetrize matrices to get full sum of tensor over all BZ, not just IBZ. 941 ! the atom is not necessarily invariant under symops, so these contributions should be added to each iat separately 942 ! normalization by nsym is done at the end outside the iqpt loop and after the tetrahedron clause 943 ! NB: looks consistent with the sym in harmonic thermo, just used in opposite 944 ! direction for symops: symrel here instead of symrec and the inverse of indsym in harmonic_thermo 945 my_nsym = crystal%nsym; if (my_qptopt == 3) my_nsym = 1 946 do isym=1, my_nsym 947 !temp_33 = matmul( (symcart(:,:,isym)), matmul(msqd_atom_tmp, transpose(symcart(:,:,isym))) ) 948 ! MG Version 949 temp_33 = matmul( (transpose(symcart(:,:,isym))), matmul(msqd_atom_tmp, symcart(:,:,isym)) ) 950 temp_33 = temp_33 / my_nsym 951 jat = crystal%indsym(4,isym,iat) 952 do idir=1,3 953 do jdir=1,3 954 phdos%msqd_dos_atom(:,idir,jdir,jat) = phdos%msqd_dos_atom(:,idir,jdir,jat) + & 955 temp_33(idir, jdir) * gvals_wtq 956 end do 957 end do 958 end do 959 960 end do ! iat 961 end do ! imode 962 963 case (2) 964 ! Tetrahedra; Save phonon frequencies, eigenvectors and angular momentum. 965 ! Sum is done after the loops over the two meshes. 966 full_phfrq(:,iq_ibz) = phfrq(:) 967 full_eigvec(:,:,:,:,iq_ibz) = eigvec 968 full_phangmom(:,:,iq_ibz) = phangmom 969 970 case default 971 ABI_ERROR(sjoin("Wrong value for prtdos:", itoa(prtdos))) 972 end select 973 end do ! iq_ibz 974 975 ABI_FREE(qbz) 976 ABI_FREE(gvals_wtq) 977 ABI_FREE(xvals) 978 979 call xmpi_sum_master(nsmallq, master, comm, ierr) 980 call xmpi_sum_master(speedofsound, master, comm, ierr) 981 982 !call cwtime_report(" phdos", cpu, wall, gflops) 983 984 if (my_rank == master .and. nsmallq > tol10) then 985 ! Write info about speed of sound 986 speedofsound = speedofsound / nsmallq 987 write (msg,'(a,E20.10,3a,F16.4,2a)') & 988 ' Average speed of sound partial sums: ', third*sum(speedofsound), ' (at units)',ch10, & 989 '- = ', third*sum(speedofsound) * Bohr_Ang * 1.d-13 / Time_Sec, ' [km/s]',ch10 990 call wrtout([std_out, ab_out], msg) 991 992 ! Debye frequency = vs * (6 pi^2 natom / ucvol)**1/3 993 debyefreq = third*sum(speedofsound) * (six*pi**2/crystal%ucvol)**(1./3.) 994 write (msg,'(a,E20.10,3a,E20.10,a)') & 995 ' Debye frequency from partial sums: ', debyefreq, ' (Ha)',ch10, & 996 '- = ', debyefreq*Ha_THz, ' (THz)' 997 call wrtout([std_out, ab_out], msg) 998 999 ! Debye temperature = hbar * Debye frequency / kb 1000 write (msg,'(a,E20.10,2a)') '-Debye temperature from partial sums: ', debyefreq*Ha_K, ' (K)', ch10 1001 call wrtout([std_out, ab_out], msg) 1002 end if 1003 1004 if (prtdos == 2) then 1005 call cwtime(cpu, wall, gflops, "start") 1006 ! Finalize integration with tetrahedra 1007 ! All the data are contained in full_phfrq, full_eigvec and full_phangmom. 1008 call xmpi_sum(full_phfrq, comm, ierr) 1009 call xmpi_sum(full_eigvec, comm, ierr) 1010 call xmpi_sum(full_phangmom, comm, ierr) 1011 1012 ABI_MALLOC(tmp_phfrq, (phdos%nqibz)) 1013 1014 do 1015 ABI_MALLOC(wdt, (phdos%nomega, 2)) 1016 ABI_MALLOC(energies, (phdos%nomega)) 1017 energies = linspace(phdos%omega_min, phdos%omega_max, phdos%nomega) 1018 1019 do iq_ibz=1,phdos%nqibz 1020 if (mod(iq_ibz, nprocs) /= my_rank) cycle ! mpi-parallelism 1021 1022 do imode=1,dos_maxmode_ 1023 ! Compute the weights for this q-point using tetrahedron 1024 tmp_phfrq(:) = full_phfrq(imode,:) 1025 call htetraq%get_onewk_wvals(iq_ibz,bcorr0,phdos%nomega,energies,max_occ1,phdos%nqibz,tmp_phfrq,wdt) 1026 wdt = wdt * wtq_ibz(iq_ibz) 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 my_nsym = crystal%nsym; if (my_qptopt == 3) my_nsym = 1 1035 do iat=1,natom 1036 do isym=1,my_nsym 1037 jat = crystal%indsym(4,isym,iat) 1038 syme2_xyza(:,jat) = syme2_xyza(:,jat) + & 1039 matmul(symcart(:,:,isym), full_eigvec(1,:,iat,imode,iq_ibz)) ** 2 + & 1040 matmul(symcart(:,:,isym), full_eigvec(2,:,iat,imode,iq_ibz)) ** 2 1041 end do 1042 end do 1043 syme2_xyza = syme2_xyza / my_nsym 1044 1045 do iat=1,natom 1046 do idir=1,3 1047 phdos%pjdos(:,idir,iat) = phdos%pjdos(:,idir,iat) + syme2_xyza(idir,iat) * wdt(:,1) 1048 phdos%pjdos_int(:,idir,iat) = phdos%pjdos_int(:,idir,iat) + syme2_xyza(idir,iat) * wdt(:,2) 1049 end do 1050 end do 1051 1052 do iat=1,natom 1053 1054 ! Accumulate outer product of displacement vectors 1055 msqd_atom_tmp = zero 1056 do idir=1,3 1057 do jdir=1,3 1058 msqd_atom_tmp(jdir,idir) = msqd_atom_tmp(jdir,idir) + ( & 1059 full_eigvec(1,idir,iat,imode,iq_ibz) * full_eigvec(1,jdir,iat,imode,iq_ibz) & 1060 + full_eigvec(2,idir,iat,imode,iq_ibz) * full_eigvec(2,jdir,iat,imode,iq_ibz) ) 1061 end do 1062 end do 1063 1064 ! Symmetrize matrices to get full sum of tensor over all BZ, not just IBZ. 1065 ! the atom is not necessarily invariant under symops, so these contributions should be added to each iat separately 1066 ! normalization by nsym is done at the end outside the iqpt loop and after the tetrahedron clause 1067 ! from loops above only the eigvec are kept and not the displ, so we still have to divide by the masses 1068 ! TODO: need to check the direction of the symcart vs transpose or inverse, given that jat is the pre-image of iat... 1069 my_nsym = crystal%nsym; if (my_qptopt == 3) my_nsym = 1 1070 !my_nsym = 1 1071 do isym=1,my_nsym 1072 !temp_33 = matmul((symcart(:,:,isym)), matmul(msqd_atom_tmp, transpose(symcart(:,:,isym)))) 1073 ! MG Version 1074 temp_33 = matmul( (transpose(symcart(:,:,isym))), matmul(msqd_atom_tmp, symcart(:,:,isym)) ) 1075 temp_33 = temp_33 / my_nsym 1076 jat = crystal%indsym(4,isym,iat) 1077 do idir=1,3 1078 do jdir=1,3 1079 phdos%msqd_dos_atom(:,idir,jdir,jat) = phdos%msqd_dos_atom(:,idir,jdir,jat) + & 1080 temp_33(idir, jdir) * wdt(:,1) 1081 end do 1082 end do 1083 end do 1084 end do ! iat 1085 1086 end do ! imode 1087 end do ! iq_ibz 1088 1089 if (refine_dosdeltae) then 1090 ! HM: Check if the integration of the DOS is correct, otherwise half dos%deltae and re-run 1091 ! MG FIXME: This won't work in parallel because we still have to call xmpi_sum 1092 call ctrap(phdos%nomega, phdos%phdos, phdos%omega_step, phdos_int) 1093 if (abs(phdos_int - crystal%natom*3) > tol2) then 1094 write(msg,'(a,f6.2,a,i4,2a,e10.3,a,e10.3)') "The value of the integral is", phdos_int, & 1095 " but it should be", crystal%natom*3, ch10,& 1096 "I will decrease dosdeltae from: ", dosdeltae, " to: ", dosdeltae/two 1097 ABI_WARNING(msg) 1098 dosdeltae = dosdeltae / two 1099 call phdos%free() 1100 call phdos_malloc(phdos, crystal, ifc, dosdeltae, dossmear, wminmax, prtdos) 1101 nomega = phdos%nomega 1102 ABI_FREE(wdt) 1103 ABI_FREE(energies) 1104 cycle 1105 endif 1106 endif 1107 1108 exit 1109 end do 1110 1111 ABI_FREE(energies) 1112 ABI_FREE(wdt) 1113 1114 ! Make eigvec into phonon displacements. 1115 do iat = 1, natom 1116 full_eigvec(:,:,iat,:,:) = full_eigvec(:,:,iat,:,:) / sqrt(phdos%atom_mass(iat)) 1117 end do 1118 1119 call cwtime_report(" tetra accumulate", cpu, wall, gflops) 1120 1121 if (my_rank == master .and. len_trim(prefix) > 0) then 1122 NCF_CHECK_MSG(nctk_open_create(ncid, strcat(prefix, "_PHIBZ.nc"), xmpi_comm_self), "Creating PHIBZ") 1123 NCF_CHECK(crystal%ncwrite(ncid)) 1124 call phonons_ncwrite(ncid, natom, phdos%nqibz, qibz, wtq_ibz, full_phfrq, full_eigvec, full_phangmom) 1125 NCF_CHECK(nf90_close(ncid)) 1126 end if 1127 1128 ! Immediately free this - it contains displ and not eigvec at this stage 1129 ABI_FREE(full_eigvec) 1130 ABI_FREE(full_phfrq) 1131 ABI_FREE(full_phangmom) 1132 ABI_FREE(tmp_phfrq) 1133 call htetraq%free() 1134 else 1135 ABI_WARNING('The netcdf PHIBZ file is only output for tetrahedron integration and DOS calculations') 1136 end if ! tetrahedra 1137 1138 ! Test if the initial mesh was large enough 1139 call xmpi_sum(count_wminmax, comm, ierr) 1140 call xmpi_min(wminmax(1), rtmp, comm, ierr); wminmax(1) = rtmp 1141 call xmpi_max(wminmax(2), rtmp, comm, ierr); wminmax(2) = rtmp 1142 1143 call xmpi_sum(phdos%phdos, comm, ierr) 1144 call xmpi_sum(phdos%msqd_dos_atom, comm, ierr) 1145 call xmpi_sum(phdos%pjdos, comm, ierr) 1146 if (prtdos == 2) then 1147 ! Reduce integrals if tetra, gauss will compute idos with simpson. 1148 call xmpi_sum(phdos%phdos_int, comm, ierr) 1149 call xmpi_sum(phdos%pjdos_int, comm, ierr) 1150 end if 1151 1152 #if 0 1153 !my_nsym = crystal%nsym; if (my_qptopt == 3) my_nsym = 1 1154 !ABI_MALLOC(work_msqd, (phdos%nomega, 3, 3, crystal%natom)) 1155 !work_msqd = phdos%msqd_dos_atom 1156 !phdos%msqd_dos_atom = zero 1157 !do iat=1,natom 1158 ! do isym=1,my_nsym 1159 ! jat = crystal%indsym(4,isym,iat) 1160 ! do io=1,phdos%nomega 1161 ! phdos%msqd_dos_atom(io,:,:,jat) = phdos%msqd_dos_atom(io,:,:,jat) + & 1162 ! matmul(transpose(symcart(:,:,isym)), matmul(work_msqd(io,:,:,iat), symcart(:,:,isym))) 1163 ! end do 1164 ! end do 1165 !end do 1166 !ABI_FREE(work_msqd) 1167 !phdos%msqd_dos_atom = phdos%msqd_dos_atom / my_nsym 1168 #endif 1169 1170 ! normalize by mass and factor of 2, now added in the printout to agree with harmonic_thermo 1171 ! do iat=1, natom 1172 ! phdos%msqd_dos_atom(:,:,:,iat) = phdos%msqd_dos_atom(:,:,:,iat) * invmass(iat) * half 1173 ! end do ! iat 1174 1175 ! =============================== 1176 ! === Compute Integrated PDOS === 1177 ! =============================== 1178 ABI_CALLOC(phdos%pjdos_rc_type, (nomega, 3, crystal%ntypat)) 1179 ABI_CALLOC(phdos%pjdos_type, (nomega, crystal%ntypat)) 1180 ABI_CALLOC(phdos%pjdos_type_int, (nomega, crystal%ntypat)) 1181 1182 do iat=1,natom 1183 itype = crystal%typat(iat) 1184 do io=1,phdos%nomega 1185 phdos%pjdos_rc_type(io,:,itype) = phdos%pjdos_rc_type(io,:,itype) + phdos%pjdos(io,:,iat) 1186 phdos%pjdos_type(io,itype) = phdos%pjdos_type(io,itype) + sum(phdos%pjdos(io,:,iat)) 1187 end do 1188 if (prtdos == 2) then 1189 do io=1,phdos%nomega 1190 phdos%pjdos_type_int(io,itype) = phdos%pjdos_type_int(io,itype) + sum(phdos%pjdos_int(io,:,iat)) 1191 end do 1192 end if 1193 end do 1194 1195 ! Evaluate IDOS using simple simpson integration 1196 ! In principle one could use derf.F90, just to be consistent ... 1197 if (prtdos == 1) then 1198 call simpson_int(phdos%nomega,phdos%omega_step, phdos%phdos, phdos%phdos_int) 1199 do iat=1,natom 1200 do idir=1,3 1201 call simpson_int(phdos%nomega, phdos%omega_step, phdos%pjdos(:,idir,iat), phdos%pjdos_int(:,idir,iat)) 1202 end do 1203 end do 1204 do itype=1,crystal%ntypat 1205 call simpson_int(phdos%nomega, phdos%omega_step, phdos%pjdos_type(:,itype), phdos%pjdos_type_int(:,itype)) 1206 end do 1207 end if 1208 1209 ABI_FREE(qibz) 1210 ABI_FREE(wtq_ibz) 1211 1212 call cwtime_report(" phdos_init", cpu_all, wall_all, gflops_all) 1213 DBG_EXIT("COLL") 1214 1215 end subroutine phdos_init
m_phonons/phdos_malloc [ Functions ]
[ Top ] [ m_phonons ] [ Functions ]
NAME
phdos_malloc
FUNCTION
INPUTS
OUTPUT
SOURCE
631 subroutine phdos_malloc(phdos, crystal, ifc, dosdeltae, dossmear, wminmax, prtdos) 632 633 ! Arguments ------------------------------------------------------ 634 class(phdos_t),intent(out) :: phdos 635 type(crystal_t),intent(in) :: crystal 636 type(ifc_type),intent(in) :: ifc 637 integer,intent(in) :: prtdos 638 real(dp),intent(in) :: dosdeltae,dossmear 639 real(dp),intent(in) :: wminmax(2) 640 641 !Local variables ------------------------- 642 integer :: io 643 ! ********************************************************************* 644 645 phdos%ntypat = crystal%ntypat 646 phdos%natom = crystal%natom 647 phdos%prtdos = prtdos 648 phdos%dossmear = dossmear 649 phdos%omega_step = dosdeltae 650 ! Use values stored in ifc (obtained with ab-initio q-mesh + pad) 651 if (wminmax(2) > wminmax(1)) then 652 phdos%omega_min = wminmax(1) 653 phdos%omega_max = wminmax(2) 654 else 655 phdos%omega_min = ifc%omega_minmax(1) 656 phdos%omega_max = ifc%omega_minmax(2) 657 end if 658 ! Must be consistent with mesh computed in tetra routines! 659 phdos%nomega = nint((phdos%omega_max - phdos%omega_min) / phdos%omega_step) + 1 660 ! Ensure Simpson integration will be ok 661 phdos%nomega = max(6, phdos%nomega) 662 663 ! Build frequency mesh. 664 ABI_MALLOC(phdos%omega, (phdos%nomega)) 665 do io=1,phdos%nomega 666 phdos%omega(io) = phdos%omega_min + phdos%omega_step * (io - 1) 667 end do 668 phdos%omega_min = phdos%omega(1) 669 phdos%omega_max = phdos%omega(phdos%nomega) 670 671 ! Allocate arrays that depend on nomega and set them to zero. 672 ABI_CALLOC(phdos%phdos, (phdos%nomega)) 673 ABI_CALLOC(phdos%phdos_int, (phdos%nomega)) 674 ABI_CALLOC(phdos%pjdos, (phdos%nomega, 3, crystal%natom)) 675 ABI_CALLOC(phdos%pjdos_int, (phdos%nomega, 3, crystal%natom)) 676 ABI_CALLOC(phdos%msqd_dos_atom, (phdos%nomega, 3, 3, crystal%natom)) 677 ABI_MALLOC(phdos%atom_mass, (crystal%natom)) 678 phdos%atom_mass = crystal%amu(crystal%typat(:)) * amu_emass 679 680 end subroutine phdos_malloc
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<phdos_t>=Container object
OUTPUT
Only writing
NOTES
Frequencies are in eV, DOS are in states/eV.
SOURCE
1718 subroutine phdos_ncwrite(phdos, ncid) 1719 1720 !Arguments ------------------------------------ 1721 !scalars 1722 class(phdos_t),intent(in) :: phdos 1723 integer,intent(in) :: ncid 1724 1725 !Local variables------------------------------- 1726 !scalars 1727 integer :: ncerr 1728 1729 ! ************************************************************************* 1730 1731 ! Define dimensions 1732 NCF_CHECK(nctk_def_basedims(ncid, defmode=.True.)) 1733 1734 ncerr = nctk_def_dims(ncid, [nctkdim_t("number_of_atoms", phdos%natom),& 1735 nctkdim_t("number_of_atom_species", phdos%ntypat), nctkdim_t("number_of_frequencies", phdos%nomega), & 1736 nctkdim_t("nqibz", phdos%nqibz)]) 1737 NCF_CHECK(ncerr) 1738 1739 !scalars 1740 NCF_CHECK(nctk_def_iscalars(ncid, ["prtdos"])) 1741 NCF_CHECK(nctk_def_dpscalars(ncid, ["dossmear"])) 1742 1743 !arrays 1744 ncerr = nctk_def_arrays(ncid, [& 1745 nctkarr_t('wmesh', "dp", 'number_of_frequencies'),& 1746 nctkarr_t('phdos', "dp", 'number_of_frequencies'),& 1747 nctkarr_t('pjdos', "dp", 'number_of_frequencies, three, number_of_atoms'),& 1748 nctkarr_t('pjdos_type', "dp", 'number_of_frequencies, number_of_atom_species'),& 1749 nctkarr_t('pjdos_rc_type', "dp", 'number_of_frequencies, three, number_of_atom_species'), & 1750 nctkarr_t('msqd_dos_atom', "dp", 'number_of_frequencies, three, three, number_of_atoms'), & 1751 nctkarr_t('qptrlatt', "int", 'three, three'), & 1752 nctkarr_t('shiftq', "dp", 'three') & 1753 ]) 1754 NCF_CHECK(ncerr) 1755 1756 ! Write variables. Note unit conversion. 1757 NCF_CHECK(nctk_set_datamode(ncid)) 1758 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "prtdos"), phdos%prtdos)) 1759 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'dossmear'), phdos%dossmear*Ha_eV)) 1760 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'wmesh'), phdos%omega*Ha_eV)) 1761 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'phdos'), phdos%phdos/Ha_eV)) 1762 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'pjdos'), phdos%pjdos/Ha_eV)) 1763 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'pjdos_type'), phdos%pjdos_type/Ha_eV)) 1764 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'pjdos_rc_type'), phdos%pjdos_rc_type/Ha_eV)) 1765 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'msqd_dos_atom'), phdos%msqd_dos_atom/Ha_eV)) 1766 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'qptrlatt'), phdos%qptrlatt)) 1767 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'shiftq'), phdos%shiftq)) 1768 1769 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.
SOURCE
271 subroutine phdos_print(PHdos, fname) 272 273 !Arguments ------------------------------------ 274 class(phdos_t),intent(in) :: PHdos 275 character(len=*),intent(in) :: fname 276 277 !Local variables------------------------------- 278 integer :: io,itype,unt,unt_by_atom,unt_msqd,iatom 279 real(dp) :: tens(3,3) 280 character(len=500) :: msg, msg_method 281 character(len=fnlen) :: fname_by_atom, fname_msqd 282 character(len=3) :: unitname 283 284 ! ************************************************************************* 285 286 ! Use Ha units everywhere 287 unitname='Ha' 288 289 select case (PHdos%prtdos) 290 case (1) 291 write(msg_method,'(a,es16.8,2a,i0)')& 292 '# Gaussian method with smearing = ',PHdos%dossmear,unitname,', nqibz =',PHdos%nqibz 293 case (2) 294 write(msg_method,'(a,i0)')'# Tetrahedron method, nqibz= ',PHdos%nqibz 295 case default 296 ABI_ERROR(sjoin(" Wrong prtdos: ",itoa(PHdos%prtdos))) 297 end select 298 299 ! Open external file and write results 300 if (open_file(fname,msg,newunit=unt,form="formatted",action="write") /= 0) then 301 ABI_ERROR(msg) 302 end if 303 write(msg,'(3a)')'# ',ch10,'# Phonon density of states and atom type projected DOS' 304 call wrtout(unt,msg) 305 write(msg,'(6a)')'# ',ch10,'# Energy in ',unitname,', DOS in states/',unitname 306 call wrtout(unt,msg) 307 call wrtout(unt,msg_method,'COLL') 308 write(msg,'(5a)')'# ',ch10,'# omega PHDOS INT_PHDOS PJDOS[atom_type=1] INT_PJDOS[atom_type=1] ... ',ch10,'# ' 309 call wrtout(unt,msg) 310 do io=1,PHdos%nomega 311 write(unt,'(3es17.8)',advance='NO')PHdos%omega(io),PHdos%phdos(io),PHdos%phdos_int(io) 312 do itype=1,PHdos%ntypat 313 write(unt,'(2es17.8,2x)',advance='NO')PHdos%pjdos_type(io,itype),PHdos%pjdos_type_int(io,itype) 314 end do 315 write(unt,*) 316 end do 317 close(unt) 318 319 fname_by_atom = trim(fname) // "_by_atom" 320 if (open_file(fname_by_atom,msg,newunit=unt_by_atom,form="formatted",action="write") /= 0) then 321 ABI_ERROR(msg) 322 end if 323 write(msg,'(3a)')'# ',ch10,'# Phonon density of states and atom projected DOS' 324 call wrtout(unt_by_atom,msg) 325 write(msg,'(6a)')'# ',ch10,'# Energy in ',unitname,', DOS in states/',unitname 326 call wrtout(unt_by_atom,msg) 327 call wrtout(unt_by_atom,msg_method) 328 write(msg,'(5a)')'# ',ch10,'# omega PHDOS PJDOS[atom=1] PJDOS[atom=2] ... ',ch10,'# ' 329 call wrtout(unt_by_atom,msg) 330 do io=1,PHdos%nomega 331 write(unt_by_atom,'(2es17.8)',advance='NO')PHdos%omega(io),PHdos%phdos(io) 332 do iatom=1,PHdos%natom 333 write(unt_by_atom,'(1es17.8,2x)',advance='NO') sum(PHdos%pjdos(io,1:3,iatom)) 334 end do 335 write(unt_by_atom,*) 336 end do 337 close(unt_by_atom) 338 339 fname_msqd = trim(fname) // "_msqd" 340 if (open_file(fname_msqd,msg,newunit=unt_msqd,form="formatted",action="write") /= 0) then 341 ABI_ERROR(msg) 342 end if 343 write(msg,'(3a)')'# ',ch10,'# Phonon density of states weighted msq displacement matrix (set to zero below 1e-12)' 344 call wrtout(unt_msqd,msg) 345 write(msg,'(6a)')'# ',ch10,'# Energy in ',unitname,', DOS in bohr^2 states/',unitname 346 call wrtout(unt_msqd,msg) 347 call wrtout(unt_msqd,msg_method) 348 write(msg,'(5a)')'# ',ch10,'# omega MSQDisp[atom=1, xx, yy, zz, yz, xz, xy] MSQDisp[atom=2, xx, yy,...] ... ',ch10,'# ' 349 call wrtout(unt_msqd,msg) 350 do io=1,PHdos%nomega 351 write(unt_msqd,'(2es17.8)',advance='NO')PHdos%omega(io) 352 do iatom=1,PHdos%natom 353 tens = PHdos%msqd_dos_atom(io,:,:,iatom) 354 where (abs(tens) < tol12) 355 tens = zero 356 end where 357 write(unt_msqd,'(6es17.8,2x)',advance='NO') & 358 tens(1,1), & 359 tens(2,2), & 360 tens(3,3), & 361 tens(2,3), & 362 tens(1,3), & 363 tens(1,2) 364 end do 365 write(unt_msqd,*) 366 end do 367 close(unt_msqd) 368 369 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 Only master node should call this routine.
INPUTS
PHdos structure
OUTPUT
to file only
SOURCE
2216 subroutine phdos_print_msqd(PHdos, fname, ntemper, tempermin, temperinc) 2217 2218 !Arguments ------------------------------- 2219 !scalars 2220 integer, intent(in) :: ntemper 2221 class(phdos_t),intent(in) :: PHdos 2222 character(len=*),intent(in) :: fname 2223 real(dp), intent(in) :: tempermin, temperinc 2224 2225 !Local variables ------------------------- 2226 integer :: io, iomin, itemp, iunit, junit, iatom 2227 real(dp) :: temper 2228 character(len=500) :: msg 2229 character(len=fnlen) :: fname_msqd 2230 character(len=fnlen) :: fname_veloc 2231 !arrays 2232 real(dp), allocatable :: bose_msqd(:,:), tmp_msqd(:,:), integ_msqd(:,:) 2233 real(dp), allocatable :: bose_msqv(:,:), tmp_msqv(:,:), integ_msqv(:,:) 2234 2235 ! ********************************************************************* 2236 2237 fname_msqd = trim(fname) //"_MSQD_T" 2238 if (open_file(fname_msqd, msg, newunit=iunit, form="formatted", status="unknown", action="write") /= 0) then 2239 ABI_ERROR(msg) 2240 end if 2241 fname_veloc = trim(fname) // "_MSQV_T" 2242 if (open_file(fname_veloc, msg, newunit=junit, form="formatted", status="unknown", action="write") /= 0) then 2243 ABI_ERROR(msg) 2244 end if 2245 2246 ! write the header 2247 write (iunit, '(a)') '# mean square displacement for each atom as a function of T (bohr^2)' 2248 write (junit, '(a)') "# mean square velocity for each atom as a function of T (bohr^2/atomic time unit^2)" 2249 2250 write (msg, '(a,F18.10,a,F18.10,a)') '# T in Kelvin, from ', tempermin, ' to ', tempermin+(ntemper-1)*temperinc 2251 write (iunit, '(a)') trim(msg) 2252 write (junit, '(a)') trim(msg) 2253 2254 write (msg, '(2a)') '# T |u^2| u_xx u_yy u_zz',& 2255 ' u_yz u_xz u_xy in bohr^2' 2256 write (iunit, '(a)') trim(msg) 2257 write (msg, '(3a)') '# T |v^2| v_xx v_yy v_zz',& 2258 ' v_yz v_xz v_xy',& 2259 ' in bohr^2/atomic time unit^2' 2260 write (junit, '(a)') trim(msg) 2261 2262 ABI_MALLOC(tmp_msqd, (PHdos%nomega,9)) 2263 ABI_MALLOC(tmp_msqv, (PHdos%nomega,9)) 2264 ABI_MALLOC(integ_msqd, (9,ntemper)) 2265 ABI_MALLOC(integ_msqv, (9,ntemper)) 2266 ABI_MALLOC(bose_msqd, (PHdos%nomega, ntemper)) 2267 ABI_MALLOC(bose_msqv, (PHdos%nomega, ntemper)) 2268 2269 do io=1, PHdos%nomega 2270 if ( PHdos%omega(io) >= 2._dp * 4.56d-6 ) exit ! 2 cm-1 TODO: make this an input parameter 2271 end do 2272 iomin = io 2273 2274 ! calculate bose only once for each atom (instead of for each atom) 2275 bose_msqd = zero 2276 bose_msqv = zero 2277 do itemp = 1, ntemper 2278 temper = tempermin + (itemp-1) * temperinc 2279 if (temper < 1.e-3) cycle ! millikelvin at least to avoid exploding Bose factor(TM) 2280 do io = iomin, PHdos%nomega 2281 ! NB: factors follow convention in phonopy documentation 2282 ! the 1/sqrt(omega) factor in phonopy is contained in the displacement vector definition 2283 ! bose() is dimensionless 2284 !bose_msqd(io, itemp) = (half + one / ( exp(PHdos%omega(io)/(kb_HaK*temper)) - one )) / PHdos%omega(io) 2285 !bose_msqv(io, itemp) = (half + one / ( exp(PHdos%omega(io)/(kb_HaK*temper)) - one )) * PHdos%omega(io) 2286 bose_msqd(io, itemp) = (half + bose_einstein(PHdos%omega(io),kb_HaK*temper)) / PHdos%omega(io) 2287 bose_msqv(io, itemp) = (half + bose_einstein(PHdos%omega(io),kb_HaK*temper)) * PHdos%omega(io) 2288 end do 2289 end do 2290 2291 do iatom=1,PHdos%natom 2292 write (msg, '(a,I8)') '# atom number ', iatom 2293 write (iunit, '(a)') trim(msg) 2294 write (junit, '(a)') trim(msg) 2295 2296 ! for each T and each atom, integrate msqd matrix with Bose Einstein factor and output 2297 integ_msqd = zero 2298 tmp_msqd = reshape(PHdos%msqd_dos_atom(:,:,:,iatom), (/PHdos%nomega, 9/)) 2299 2300 ! perform all integrations as matrix multiplication: integ_msqd (idir, itemp) = [tmp_msqd(io,idir)]^T * bose_msqd(io,itemp) 2301 call DGEMM('T','N', 9, ntemper, PHdos%nomega, one, tmp_msqd,PHdos%nomega, bose_msqd, PHdos%nomega, zero, integ_msqd, 9) 2302 ! NB: this presumes an equidistant omega grid 2303 integ_msqd = integ_msqd * (PHdos%omega(2)-PHdos%omega(1)) / PHdos%atom_mass(iatom) 2304 2305 integ_msqv = zero 2306 tmp_msqv = reshape(PHdos%msqd_dos_atom(:,:,:,iatom), (/PHdos%nomega, 9/)) 2307 2308 ! perform all integrations as matrix multiplication: integ_msqv (idir, itemp) = [tmp_msqv(io,idir)]^T * bose_msqv(io,itemp) 2309 call DGEMM('T','N', 9, ntemper, PHdos%nomega, one, tmp_msqv,PHdos%nomega, bose_msqv, PHdos%nomega, zero, integ_msqv, 9) 2310 ! NB: this presumes an equidistant omega grid 2311 integ_msqv = integ_msqv * (PHdos%omega(2)-PHdos%omega(1)) / PHdos%atom_mass(iatom) 2312 2313 ! print out stuff 2314 do itemp = 1, ntemper 2315 temper = tempermin + (itemp-1) * temperinc 2316 write (msg, '(F10.2,4x,E22.10,2x,6E22.10)') & 2317 temper, third*(integ_msqd(1,itemp)+integ_msqd(5,itemp)+integ_msqd(9,itemp)), & 2318 integ_msqd(1,itemp),integ_msqd(5,itemp),integ_msqd(9,itemp), & 2319 integ_msqd(6,itemp),integ_msqd(3,itemp),integ_msqd(2,itemp) 2320 write (iunit, '(a)') trim(msg) 2321 write (msg, '(F10.2,4x,E22.10,2x,6E22.10)') & 2322 temper, third*(integ_msqv(1,itemp)+integ_msqv(5,itemp)+integ_msqv(9,itemp)), & 2323 integ_msqv(1,itemp),integ_msqv(5,itemp),integ_msqv(9,itemp), & 2324 integ_msqv(6,itemp),integ_msqv(3,itemp),integ_msqv(2,itemp) 2325 write (junit, '(a)') trim(msg) 2326 end do ! itemp 2327 2328 write (iunit, '(a)') '' 2329 write (junit, '(a)') '' 2330 enddo ! iatom 2331 2332 ABI_FREE(tmp_msqd) 2333 ABI_FREE(tmp_msqv) 2334 ABI_FREE(bose_msqd) 2335 ABI_FREE(bose_msqv) 2336 ABI_FREE(integ_msqd) 2337 ABI_FREE(integ_msqv) 2338 2339 close(iunit) 2340 close(junit) 2341 2342 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
SOURCE
2156 subroutine phdos_print_vsound(iunit,ucvol,speedofsound) 2157 2158 !Arguments ------------------------------- 2159 !scalras 2160 integer, intent(in) :: iunit 2161 real(dp), intent(in) :: ucvol 2162 !arrays 2163 real(dp), intent(in) :: speedofsound(3) 2164 2165 !Local variables ------------------------- 2166 integer :: imode_acoustic 2167 character(len=500) :: msg 2168 real(dp) :: tdebye 2169 2170 ! ********************************************************************* 2171 2172 do imode_acoustic = 1, 3 2173 ! from phonon frequency, estimate speed of sound by linear interpolation from Gamma 2174 write (msg, '(2a,a,E20.10,a,a,F20.5)') & 2175 ' Speed of sound for this q and mode:',ch10,& 2176 ' in atomic units: ', speedofsound(imode_acoustic), ch10,& 2177 ' in units km/s: ', speedofsound(imode_acoustic) * Bohr_Ang * 1.d-13 / Time_Sec 2178 call wrtout(iunit, msg) 2179 call wrtout(std_out, msg) 2180 2181 ! also estimate partial Debye temperature, = energy if this band went to zone edge 2182 tdebye = speedofsound(imode_acoustic) * pi * (six / pi / ucvol)**(third) 2183 write (msg, '(2a,a,E20.10,a,a,F20.5)') & 2184 ' Partial Debye temperature for this q and mode:',ch10,& 2185 ' in atomic units: ', tdebye, ch10,& 2186 ' in SI units K : ', tdebye * Ha_K 2187 call wrtout(iunit,msg) 2188 call wrtout(iunit,"") 2189 call wrtout(std_out,msg) 2190 call wrtout(std_out,"") 2191 end do 2192 2193 end subroutine phdos_print_vsound
m_phonons/phdos_t [ Types ]
[ Top ] [ m_phonons ] [ Types ]
NAME
phdos_t
FUNCTION
Container for phonon DOS and atom projected contributions
SOURCE
87 type,public :: phdos_t 88 89 integer :: ntypat 90 ! Number of type of atoms. 91 92 integer :: natom 93 ! Number of atoms is the unit cell. 94 95 integer :: prtdos 96 ! Option of DOS calculation (1 for Gaussian, 2 for tetrahedrons). 97 98 integer :: nomega 99 ! Number of frequency points in DOS mesh. 100 101 integer :: nqibz 102 ! Number of q-points in the IBZ. 103 104 real(dp) :: omega_min 105 ! Min frequency for DOS calculation. 106 107 real(dp) :: omega_max 108 ! Max frequency for DOS calculation. 109 110 real(dp) :: omega_step 111 ! Frequency step. 112 113 real(dp) :: dossmear 114 ! Gaussian broadening. 115 116 integer :: qptrlatt(3,3) = 0 117 ! q-mesh as computed in getkgrid_low 118 119 real(dp) :: shiftq(3) 120 ! Shigt of Q-mesh computed by getkgrid_low (1 shift is enough) 121 122 real(dp),allocatable :: atom_mass(:) 123 ! atom_mass(natom) 124 125 real(dp),allocatable :: omega(:) 126 ! omega(nomega) 127 ! Frequency grid. 128 129 real(dp),allocatable :: phdos(:) 130 ! phdos(nomega) 131 ! phonon DOS. 132 133 real(dp),allocatable :: phdos_int(:) 134 ! phdos_int(nomega) 135 ! integrated phonon DOS 136 137 real(dp),allocatable :: pjdos(:,:,:) 138 ! pjdos(nomega,3,natom) 139 ! projected DOS (over atoms and cartesian directions) 140 141 real(dp),allocatable :: pjdos_int(:,:,:) 142 ! pjdos_int(nomega,3,natom) 143 ! Integrated atomic PJDOS along the three cartesian directions. 144 145 real(dp),allocatable :: pjdos_type(:,:) 146 ! pjdos_type(nomega,ntypat) 147 ! phonon DOS contribution arising from a particular atom-type. 148 149 real(dp),allocatable :: pjdos_type_int(:,:) 150 ! pjdos_type_int(nomega,ntypat) 151 ! Integrate phonon DOS contribution arising from a particular atom-type. 152 153 real(dp),allocatable :: pjdos_rc_type(:,:,:) 154 ! phdos(nomega,3,ntypat) 155 ! phonon DOS contribution arising from a particular atom-type 156 ! decomposed along the three cartesian directions. 157 158 real(dp),allocatable :: msqd_dos_atom(:,:,:,:) 159 ! msqd_dos_atom(nomega,3,3,natom) 160 ! mean square displacement matrix, frequency dependent like a DOS, tensor in cartesian coords. 161 ! allows one to calculate Debye Waller factors by integration with 1/omega 162 ! and the Bose Einstein factor 163 164 contains 165 166 procedure :: print => phdos_print 167 procedure :: print_debye => phdos_print_debye 168 procedure :: print_msqd => phdos_print_msqd 169 procedure :: print_thermo => phdos_print_thermo 170 procedure :: free => phdos_free 171 procedure :: ncwrite => phdos_ncwrite 172 procedure :: init => phdos_init ! Constructor 173 end type phdos_t
m_phonons/pheigvec_rotate [ Functions ]
[ Top ] [ m_phonons ] [ Functions ]
NAME
pheigvec_rotate
FUNCTION
Return phonon eigenvectors for q in the BZ from the symmetrical image in the IBZ.
INPUTS
cryst: crystal structure qq_ibz: q-point in the IBZ isym itimrev eigvec_ibz: Input phonon eigenvectors at qq_ibz.
OUTPUT
eigvec_bz: phonon eigenvectors at q_bz. displ_cart_qbz: phonon displacement at q_bz in Cartesian coordinates. [displ_red_qbz]: phonon displacement at q_bz in reduced coordinates.
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. phangmom= Phonon angular momentum 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 The angular momentum is output in units of hbar
OUTPUT
Only writing
SOURCE
2374 subroutine phonons_ncwrite(ncid,natom,nqpts,qpoints,weights,phfreq,phdispl_cart,phangmom) 2375 2376 !Arguments ------------------------------------ 2377 !scalars 2378 integer,intent(in) :: ncid,natom,nqpts 2379 !arrays 2380 real(dp),intent(in) :: qpoints(3,nqpts),weights(nqpts) 2381 real(dp),intent(in) :: phfreq(3*natom,nqpts),phdispl_cart(2,3*natom,3*natom,nqpts),phangmom(3,3*natom,nqpts) 2382 2383 !Local variables------------------------------- 2384 !scalars 2385 integer :: nphmodes,ncerr 2386 2387 ! ************************************************************************* 2388 2389 nphmodes = 3*natom 2390 2391 NCF_CHECK(nctk_def_basedims(ncid, defmode=.True.)) 2392 2393 ncerr = nctk_def_dims(ncid, [& 2394 nctkdim_t("number_of_qpoints", nqpts), nctkdim_t('number_of_phonon_modes', nphmodes), nctkdim_t('three', 3)]) 2395 NCF_CHECK(ncerr) 2396 2397 ! Define arrays 2398 ncerr = nctk_def_arrays(ncid, [& 2399 nctkarr_t('qpoints', "dp" , 'number_of_reduced_dimensions, number_of_qpoints'),& 2400 nctkarr_t('qweights',"dp", 'number_of_qpoints'),& 2401 nctkarr_t('phfreqs',"dp", 'number_of_phonon_modes, number_of_qpoints'),& 2402 nctkarr_t('phdispl_cart',"dp", 'complex, number_of_phonon_modes, number_of_phonon_modes, number_of_qpoints'),& 2403 nctkarr_t('phangmom',"dp", 'three, number_of_phonon_modes, number_of_qpoints')]) 2404 NCF_CHECK(ncerr) 2405 2406 ! Write variables. 2407 NCF_CHECK(nctk_set_datamode(ncid)) 2408 NCF_CHECK(nf90_put_var(ncid, vid('qpoints'), qpoints)) 2409 NCF_CHECK(nf90_put_var(ncid, vid('qweights'), weights)) 2410 NCF_CHECK(nf90_put_var(ncid, vid('phfreqs'), phfreq*Ha_eV)) 2411 NCF_CHECK(nf90_put_var(ncid, vid('phdispl_cart'), phdispl_cart*Bohr_Ang)) 2412 NCF_CHECK(nf90_put_var(ncid, vid('phangmom'), phangmom)) 2413 2414 2415 contains 2416 integer function vid(vname) 2417 character(len=*),intent(in) :: vname 2418 vid = nctk_idname(ncid, vname) 2419 end function vid 2420 2421 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
SOURCE
2691 subroutine phonons_write_gnuplot(prefix, natom, nqpts, qpts, phfreqs, qptbounds) 2692 2693 !Arguments ------------------------------------ 2694 !scalars 2695 integer,intent(in) :: natom,nqpts 2696 real(dp),intent(in) :: qpts(3,nqpts),phfreqs(3*natom,nqpts) 2697 character(len=*),intent(in) :: prefix 2698 !arrays 2699 real(dp),optional,intent(in) :: qptbounds(:,:) 2700 2701 !Local variables------------------------------- 2702 !scalars 2703 integer :: unt,iq,ii,start,nqbounds,gpl_unt 2704 character(len=500) :: msg,fmt 2705 character(len=fnlen) :: datafile,basefile 2706 !arrays 2707 integer :: g0(3) 2708 integer,allocatable :: bounds2qpt(:) 2709 2710 ! ********************************************************************* 2711 2712 nqbounds = 0 2713 if (present(qptbounds)) then 2714 if (product(shape(qptbounds)) > 0 ) then 2715 ! Find correspondence between qptbounds and k-points in ebands. 2716 nqbounds = size(qptbounds, dim=2) 2717 ABI_MALLOC(bounds2qpt, (nqbounds)) 2718 bounds2qpt = 1; start = 1 2719 do ii=1,nqbounds 2720 do iq=start,nqpts 2721 if (isamek(qpts(:, iq), qptbounds(:, ii), g0)) then 2722 bounds2qpt(ii) = iq; start = iq + 1; exit 2723 end if 2724 end do 2725 end do 2726 end if 2727 end if 2728 2729 datafile = strcat(prefix, "_PHBANDS.data") 2730 if (open_file(datafile, msg, newunit=unt, form="formatted", action="write") /= 0) then 2731 ABI_ERROR(msg) 2732 end if 2733 if (open_file(strcat(prefix, "_PHBANDS.gnuplot"), msg, newunit=gpl_unt, form="formatted", action="write") /= 0) then 2734 ABI_ERROR(msg) 2735 end if 2736 basefile = basename(datafile) 2737 2738 write(unt,'(a)') "# Phonon band structure data file" 2739 write(unt,'(a)') "# Generated by Abinit" 2740 write(unt,'(2(a,i0))') "# natom: ",natom,", nqpt: ",nqpts 2741 write(unt,'(a)') "# Frequencies are in meV" 2742 write(unt,'(a)')"# List of q-points and their index (C notation i.e. count from 0)" 2743 do iq=1,nqpts 2744 write(unt, "(a)")sjoin("#", itoa(iq-1), ktoa(qpts(:,iq))) 2745 end do 2746 2747 fmt = sjoin("(i0,1x,", itoa(3*natom), "(es16.8,1x))") 2748 write(unt,'(a)')"# [kpt-index, mode_1, mode_2 ...]" 2749 do iq=1,nqpts 2750 write(unt, fmt) iq-1, phfreqs(:, iq) * Ha_meV 2751 end do 2752 2753 ! gnuplot script file 2754 write(gpl_unt,'(a)') '# File to plot electron bandstructure with gnuplot' 2755 !write(gpl_unt,'(a)') "#set terminal postscript eps enhanced color font 'Times-Roman,26' lw 2" 2756 write(gpl_unt,'(a)') '#use the next lines to make a nice figure for a paper' 2757 write(gpl_unt,'(a)') '#set term postscript enhanced eps color lw 0.5 dl 0.5' 2758 write(gpl_unt,'(a)') '#set pointsize 0.275' 2759 write(gpl_unt,'(a)') 'set palette defined ( 0 "blue", 3 "green", 6 "yellow", 10 "red" )' 2760 write(gpl_unt,'(a)') 'unset key' 2761 write(gpl_unt,'(a)') '# can make pointsize smaller (~0.5). Too small and nothing is printed' 2762 write(gpl_unt,'(a)') 'set pointsize 0.8' 2763 write(gpl_unt,'(a)') 'set view 0,0' 2764 write(gpl_unt,'(a,i0,a)') 'set xrange [0:',nqpts-1,']' 2765 write(gpl_unt,'(2(a,es16.8),a)')& 2766 'set yrange [',minval(phfreqs * Ha_meV),':',maxval(phfreqs * Ha_meV),']' 2767 write(gpl_unt,'(a)') 'set xlabel "Momentum"' 2768 write(gpl_unt,'(a)') 'set ylabel "Energy [meV]"' 2769 write(gpl_unt,'(a)') strcat('set title "', replace(basefile, "_", "\\_"),'"') 2770 if (nqbounds == 0) then 2771 write(gpl_unt,'(a)') 'set grid xtics' 2772 else 2773 write(gpl_unt,"(a)")"# Add vertical lines in correspondence of high-symmetry points." 2774 write(gpl_unt,'(a)') 'unset xtics' 2775 do ii=1,nqbounds 2776 write(gpl_unt,"(a,2(i0,a))") & 2777 "set arrow from ",bounds2qpt(ii)-1,",graph(0,0) to ",bounds2qpt(ii)-1,",graph(1,1) nohead" 2778 !write(gpl_unt,"(a)")sjoin("set xtics add ('kname'", itoa(bounds2kpt(ii)-1), ")") 2779 end do 2780 end if 2781 write(gpl_unt,"(a)")sjoin("nbranch =", itoa(3*natom)) 2782 write(gpl_unt,"(a)")strcat('plot for [i=2:nbranch] "', basefile, '" u 1:i every :1 with lines linetype -1') 2783 write(gpl_unt,"(a)")"pause -1" 2784 2785 close(unt) 2786 close(gpl_unt) 2787 2788 ABI_SFREE(bounds2qpt) 2789 2790 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
SOURCE
2449 subroutine phonons_write_phfrq(path,natom,nqpts,qpoints,weights,phfreq,phdispl_cart,phangmom) 2450 2451 !Arguments ------------------------------------ 2452 !scalars 2453 integer,intent(in) :: natom,nqpts 2454 character(len=*),intent(in) :: path 2455 !arrays 2456 real(dp),intent(in) :: qpoints(3,nqpts),weights(nqpts) 2457 real(dp),intent(in) :: phfreq(3*natom,nqpts) 2458 real(dp),intent(in) :: phdispl_cart(2,3*natom,3*natom,nqpts) 2459 real(dp),intent(in) :: phangmom(3,3*natom,nqpts) 2460 2461 !Local variables------------------------------- 2462 !scalars 2463 integer :: nphmodes, iq, iunit, imod, icomp 2464 real(dp) :: dummy 2465 character(len=300) :: formt 2466 character(len=500) :: msg 2467 2468 ! ************************************************************************* 2469 2470 nphmodes = 3*natom 2471 2472 dummy = qpoints(1,1); dummy = weights(1) 2473 2474 ! Write phonon frequencies 2475 if (open_file(strcat(path, "_PHFRQ"), msg, newunit=iunit, form="formatted", status="unknown", action="write") /= 0) then 2476 ABI_ERROR(msg) 2477 end if 2478 2479 write (iunit, '(a)') '# ABINIT generated phonon band structure file. All in Ha atomic units' 2480 write (iunit, '(a)') '# ' 2481 write (iunit, '(a,i0)') '# number_of_qpoints ', nqpts 2482 write (iunit, '(a,i0)') '# number_of_phonon_modes ', nphmodes 2483 write (iunit, '(a)') '# ' 2484 2485 write (formt,'(a,i0,a)') "(I5, ", nphmodes, "E20.10)" 2486 do iq= 1, nqpts 2487 write (iunit, formt) iq, phfreq(:,iq) 2488 end do 2489 2490 close(iunit) 2491 2492 ! Does not Write phonon displacement ? 2493 if (.False.) then 2494 if (open_file(strcat(path, "_PHDISPL"), msg, unit=iunit, form="formatted", status="unknown", action="write") /= 0) then 2495 ABI_ERROR(msg) 2496 end if 2497 2498 write (iunit, '(a)') '# ABINIT generated phonon displacements, along points in PHFRQ file. All in Ha atomic units' 2499 write (iunit, '(a)') '# ' 2500 write (iunit, '(a)') '# displacements in cartesian coordinates, Re and Im parts ' 2501 write (iunit, '(a,i0)') '# number_of_qpoints ', nqpts 2502 write (iunit, '(a,i0)') '# number_of_phonon_modes ', nphmodes 2503 write (iunit, '(a)') '# ' 2504 2505 !write (formt,'(a,I3,a)') "( ", nphmodes, "(2E20.10,2x))" 2506 formt = "(2E20.10,2x)" 2507 2508 do iq = 1, nqpts 2509 write (iunit, '(a, i0)') '# iq ', iq 2510 do imod = 1, nphmodes 2511 write (iunit, '(a, i0)') '# imode ', imod 2512 do icomp = 1, nphmodes 2513 write (iunit, formt, ADVANCE='NO') phdispl_cart(:,icomp,imod,iq) 2514 end do 2515 write (iunit, '(a)') ' ' 2516 end do 2517 end do 2518 2519 close(iunit) 2520 end if 2521 2522 ! Write phonon angular momentum 2523 if (open_file(strcat(path, "_PHANGMOM"), msg, unit=iunit, form="formatted", status="unknown", action="write") /= 0) then 2524 ABI_ERROR(msg) 2525 end if 2526 2527 write (iunit, '(a)') '# ABINIT generated phonon angular momentum, along points in PHFRQ file. All in Ha atomic units' 2528 write (iunit, '(a)') '# ' 2529 write (iunit, '(a)') '# angular momentum in cartesian coordinates ' 2530 write (iunit, '(a,i0)') '# number_of_qpoints ', nqpts 2531 write (iunit, '(a,i0)') '# number_of_phonon_modes ', nphmodes 2532 write (iunit, '(a)') '# ' 2533 2534 write (formt,'(a,i0,a)') "(I5, ", nphmodes, "E20.10)" 2535 do icomp = 1, 3 2536 do iq= 1, nqpts 2537 write (iunit, formt) iq, phangmom(icomp,:,iq) 2538 end do 2539 if (icomp /= 3) then 2540 write (iunit, '(a,a)') '' 2541 end if 2542 end do 2543 2544 close(iunit) 2545 2546 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
SOURCE
2571 subroutine phonons_write_xmgrace(filename, natom, nqpts, qpts, phfreqs, qptbounds) 2572 2573 !Arguments ------------------------------------ 2574 !scalars 2575 integer,intent(in) :: natom,nqpts 2576 real(dp),intent(in) :: qpts(3,nqpts),phfreqs(3*natom,nqpts) 2577 character(len=*),intent(in) :: filename 2578 !arrays 2579 real(dp),optional,intent(in) :: qptbounds(:,:) 2580 2581 !Local variables------------------------------- 2582 !scalars 2583 integer :: unt,iq,nu,ii,start,nqbounds 2584 character(len=500) :: msg 2585 !arrays 2586 integer :: g0(3) 2587 integer,allocatable :: bounds2qpt(:) 2588 2589 ! ********************************************************************* 2590 2591 nqbounds = 0 2592 if (present(qptbounds)) then 2593 if (product(shape(qptbounds)) > 0 ) then 2594 ! Find correspondence between qptbounds and k-points in ebands. 2595 nqbounds = size(qptbounds, dim=2) 2596 ABI_MALLOC(bounds2qpt, (nqbounds)) 2597 bounds2qpt = 1; start = 1 2598 do ii=1,nqbounds 2599 do iq=start,nqpts 2600 if (isamek(qpts(:, iq), qptbounds(:, ii), g0)) then 2601 bounds2qpt(ii) = iq; start = iq + 1; exit 2602 end if 2603 end do 2604 end do 2605 end if 2606 end if 2607 2608 if (open_file(filename, msg, newunit=unt, form="formatted", action="write") /= 0) then 2609 ABI_ERROR(msg) 2610 end if 2611 2612 write(unt,'(a)') "# Grace project file" 2613 write(unt,'(a)') "# Generated by Abinit" 2614 write(unt,'(2(a,i0))') "# natom: ",natom,", nqpt: ",nqpts 2615 write(unt,'(a)') "# Frequencies are in meV" 2616 write(unt,'(a)')"# List of q-points and their index (C notation i.e. count from 0)" 2617 do iq=1,nqpts 2618 write(unt, "(a)")sjoin("#", itoa(iq-1), ktoa(qpts(:,iq))) 2619 end do 2620 2621 write(unt,'(a)') "@page size 792, 612" 2622 write(unt,'(a)') "@page scroll 5%" 2623 write(unt,'(a)') "@page inout 5%" 2624 write(unt,'(a)') "@link page off" 2625 write(unt,'(a)') "@with g0" 2626 write(unt,'(a)') "@world xmin 0.00" 2627 write(unt,'(a,i0)') '@world xmax ',nqpts 2628 write(unt,'(a,e16.8)') '@world ymin ',minval(phfreqs * Ha_meV) 2629 write(unt,'(a,e16.8)') '@world ymax ',maxval(phfreqs * Ha_meV) 2630 write(unt,'(a)') '@default linewidth 1.5' 2631 write(unt,'(a)') '@xaxis tick on' 2632 write(unt,'(a)') '@xaxis tick major 1' 2633 write(unt,'(a)') '@xaxis tick major color 1' 2634 write(unt,'(a)') '@xaxis tick major linestyle 3' 2635 write(unt,'(a)') '@xaxis tick major grid on' 2636 write(unt,'(a)') '@xaxis tick spec type both' 2637 write(unt,'(a)') '@xaxis tick major 0, 0' 2638 if (nqbounds /= 0) then 2639 write(unt,'(a,i0)') '@xaxis tick spec ',nqbounds 2640 do iq=1,nqbounds 2641 !write(unt,'(a,i0,a,a)') '@xaxis ticklabel ',iq-1,',', "foo" 2642 write(unt,'(a,i0,a,i0)') '@xaxis tick major ',iq-1,' , ',bounds2qpt(iq) - 1 2643 end do 2644 end if 2645 write(unt,'(a)') '@xaxis ticklabel char size 1.500000' 2646 write(unt,'(a)') '@yaxis tick major 10' 2647 write(unt,'(a)') '@yaxis label "Phonon Energy [meV]"' 2648 write(unt,'(a)') '@yaxis label char size 1.500000' 2649 write(unt,'(a)') '@yaxis ticklabel char size 1.500000' 2650 do nu=1,3*natom 2651 write(unt,'(a,i0,a)') '@ s',nu-1,' line color 1' 2652 end do 2653 do nu=1,3*natom 2654 write(unt,'(a,i0)') '@target G0.S',nu-1 2655 write(unt,'(a)') '@type xy' 2656 do iq=1,nqpts 2657 write(unt,'(i0,1x,e16.8)') iq-1, phfreqs(nu, iq) * Ha_meV 2658 end do 2659 write(unt,'(a)') '&' 2660 end do 2661 2662 close(unt) 2663 2664 ABI_SFREE(bounds2qpt) 2665 2666 end subroutine phonons_write_xmgrace
m_phonons/phstore_async_rotate [ Functions ]
[ Top ] [ m_phonons ] [ Functions ]
NAME
phstore_async_rotate
FUNCTION
Begin non-blocking collective MPI communication inside self%comm to obtain phonon frequencies and eigenvectors in the BZ from data in the IBZ.
INPUTS
m_phonons/phstore_free [ Functions ]
[ Top ] [ m_phonons ] [ Functions ]
NAME
phstore_free
FUNCTION
Free dynamic memory.
INPUTS
m_phonons/phstore_new [ Functions ]
[ Top ] [ m_phonons ] [ Functions ]
NAME
phstore_new
FUNCTION
Create new object with phonon quantities in the IBZ.
INPUTS
cryst: Crystal structure ifc: Interatomic force constants. nqibz: Number of q-points in the IBZ. qibz: q-points in the IBZ use_ifc_fourq: True to replace symmetrization with call to ifc_fourq (debugging option) comm: MPI communicator in which phonon arrays in the IBZ will be MPI distributed.
m_phonons/phstore_t [ Types ]
[ Top ] [ m_phonons ] [ Types ]
NAME
phstore_t
FUNCTION
This object stores ph eigenvalues and eigenvectors in the IBZ and provides methods to compute the corresponding quantites in the full BZ using symmetries. Useful for very intensive loops of q-points in the BZ in which the call to ifc_fourq may become a significant bottleneck. Note that IBZ quantities are memory-distributed inside the MPI communicator comm. and the symmetrization is performed in a non-blocking fashion so that it's possible to overlap the symmetrization with computations. See sigmaph for usage.
SOURCE
192 type,public :: phstore_t 193 194 integer :: nqibz 195 ! Number of q-points in the IBZ 196 197 integer :: comm 198 ! MPI communicator used to distribute memory. 199 200 integer :: nprocs 201 ! Number of MPI procs in comm. 202 203 integer :: my_rank 204 ! Rank of this MPI proc inside comm 205 206 integer :: natom, natom3 207 208 logical :: use_ifc_fourq = .False. 209 ! Debuggin flag. If True, replace symmetrization with call to ifc_fourq. 210 211 integer :: requests(2) 212 ! MPI requests 213 214 integer,allocatable :: qibz_start(:), qibz_stop(:) 215 ! (0:%nprocs-1)) 216 ! Initial and final index of the IBZ qpoint treated by this MPI proc inside comm. 217 218 real(dp), ABI_CONTIGUOUS pointer :: qibz(:,:) 219 ! q-points in the IBZ. 220 221 real(dp),allocatable :: phfreqs_qibz(:,:) 222 ! (natom3, %nqibz)) 223 ! Ph frequencies in the IBZ 224 225 real(dp),allocatable :: pheigvec_qibz(:,:,:,:) 226 ! (2, natom3, natom3, %nqibz)) 227 ! Ph eigenvectors in the IBZ 228 229 real(dp),allocatable :: phfrq(:) 230 ! (natom3) 231 ! Ph frequencies for q in the full BZ 232 233 real(dp),allocatable :: displ_cart(:,:,:,:) 234 ! displ_cart(2, 3, natom, %natom3) 235 ! Ph displacement for q in the full BZ 236 237 contains 238 239 procedure :: async_rotate => phstore_async_rotate ! Begin non-blocking collective MPI communication to symmetrize stuff 240 procedure :: wait => phstore_wait ! Wait from non-blocking MPI BCAST started in phstore_async_rotate, 241 ! return ph frequencies and displacements. 242 procedure :: free => phstore_free ! Free dynamic memory 243 244 end type phstore_t
m_phonons/phstore_wait [ Functions ]
[ Top ] [ m_phonons ] [ Functions ]
NAME
phstore_wait
FUNCTION
Wait from non-blocking MPI BCAST started in phstore_async_rotate, returns phonon frequencies and displacements in Cartesian and reduced coordinates.
INPUTS
m_phonons/test_phrotation [ Functions ]
[ Top ] [ m_phonons ] [ Functions ]
NAME
test_phrotation
FUNCTION
Test the symmetrization of the phonon eigenvalues and eigenvectors. INPUT cryst=Crystalline structure ifc<ifc_type>=interatomic force constants and corresponding real space grid info. ngqpt(3)=Divisions of the ab-initio q-mesh. qptopt=option for the generation of q points (defines whether spatial symmetries and/or time-reversal can be used) comm= MPI communicator
m_phonons/thermal_supercell_free [ Functions ]
[ Top ] [ m_phonons ] [ Functions ]
NAME
thermal_supercell_free
FUNCTION
deallocate thermal array of supercells
INPUTS
OUTPUT
NOTES
SOURCE
1593 subroutine thermal_supercell_free(nscells, thm_scells) 1594 1595 !Arguments ------------------------------------ 1596 !scalars 1597 integer, intent(in) :: nscells 1598 type(supercell_type), allocatable, intent(inout) :: thm_scells(:) 1599 1600 ! local 1601 integer :: icell 1602 1603 if (allocated(thm_scells)) then 1604 do icell = 1, nscells 1605 call destroy_supercell(thm_scells(icell)) 1606 end do 1607 end if 1608 1609 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
SOURCE
1404 subroutine thermal_supercell_make(amplitudes,Crystal, Ifc,namplitude, nconfig,option,& 1405 & rlatt, temperature_K, thm_scells) 1406 1407 !Arguments ------------------------------------ 1408 !scalars 1409 integer, intent(in) :: option,nconfig 1410 integer, intent(in) :: rlatt(3,3) 1411 real(dp), intent(in) :: temperature_K 1412 type(crystal_t),intent(in) :: Crystal 1413 type(ifc_type),intent(in) :: Ifc 1414 type(supercell_type), intent(out) :: thm_scells(nconfig) 1415 integer,intent(in) :: namplitude 1416 !Local variables------------------------------- 1417 !scalars 1418 integer :: iq, nqibz, nqbz, qptopt1, iampl ,imode, ierr, iconfig 1419 real(dp) :: temperature, sigma, freeze_displ 1420 real(dp) :: rand !, rand1, rand2 1421 real(dp),intent(in):: amplitudes(5,namplitude) 1422 !arrays 1423 real(dp), allocatable :: qshft(:,:) ! dummy with 2 dimensions for call to kpts_ibz_from_kptrlatt 1424 real(dp), allocatable :: qbz(:,:), qibz(:,:), wtqibz(:) 1425 real(dp), allocatable :: phfrq_allq(:,:), phdispl_allq(:,:,:,:,:) 1426 real(dp), allocatable :: phfrq(:), phdispl(:,:,:,:),pheigvec(:,:,:,:) 1427 real(dp), allocatable :: phdispl1(:,:,:) 1428 character (len=500) :: msg 1429 1430 ! ************************************************************************* 1431 ! check inputs 1432 ! TODO: add check that all rlatt are the same on input 1433 if (rlatt(1,2)/=0 .or. rlatt(1,3)/=0 .or. rlatt(2,3)/=0 .or. & 1434 & rlatt(2,1)/=0 .or. rlatt(3,1)/=0 .or. rlatt(3,2)/=0) then 1435 write (msg, '(4a, 9I6, a)') ' for the moment I have not implemented ', & 1436 & ' non diagonal supercells.',ch10,' rlatt for temp 1 = ', rlatt, ' Returning ' 1437 ABI_WARNING(msg) 1438 return 1439 end if 1440 1441 temperature = temperature_K / Ha_K 1442 1443 ! build qpoint grid used for the Fourier interpolation. 1444 !(use no syms for the moment!) 1445 qptopt1 = 3 1446 1447 ! for the moment do not allow shifted q grids. 1448 ! We are interpolating anyway, so it will always work 1449 ABI_MALLOC(qshft,(3,1)) 1450 qshft(:,1)=zero 1451 1452 ! This call will set nqibz, IBZ and BZ arrays 1453 call kpts_ibz_from_kptrlatt(crystal, rlatt, qptopt1, 1, qshft, & 1454 nqibz, qibz, wtqibz, nqbz, qbz) ! new_kptrlatt, new_shiftk) ! Optional 1455 ABI_FREE(qshft) 1456 1457 ! allocate arrays wzith all of the q, omega, and displacement vectors 1458 ABI_MALLOC_OR_DIE(phfrq_allq, (3*Crystal%natom, nqibz), ierr) 1459 ABI_MALLOC_OR_DIE(phdispl_allq, (2, 3, Crystal%natom, 3*Crystal%natom, nqibz), ierr) 1460 1461 ABI_MALLOC(phfrq, (3*Crystal%natom)) 1462 ABI_MALLOC(phdispl, (2, 3, Crystal%natom, 3*Crystal%natom)) 1463 ABI_MALLOC(pheigvec, (2, 3, Crystal%natom, 3*Crystal%natom)) 1464 1465 ! loop over q to get all frequencies and displacement vectors 1466 imode = 0 1467 do iq = 1, nqibz 1468 ! Fourier interpolation. 1469 call ifc%fourq(Crystal, qibz(:,iq), phfrq, phdispl, out_eigvec=pheigvec) 1470 phfrq_allq(1:3*Crystal%natom, iq) = phfrq 1471 phdispl_allq(1:2, 1:3, 1:Crystal%natom, 1:3*Crystal%natom, iq) = phdispl 1472 end do 1473 ABI_FREE(phfrq) 1474 ABI_FREE(pheigvec) 1475 ABI_FREE(phdispl) 1476 1477 ! only diagonal supercell case for the moment 1478 do iconfig = 1, nconfig 1479 call init_supercell(Crystal%natom, rlatt, Crystal%rprimd, Crystal%typat, Crystal%xcart, Crystal%znucl, thm_scells(iconfig)) 1480 end do 1481 1482 ! precalculate phase factors??? 1483 1484 ABI_MALLOC(phdispl1, (2, 3, Crystal%natom)) 1485 1486 ! for all modes at all q in whole list, sorted 1487 do iq = 1, nqibz 1488 do imode = 1, 3*Crystal%natom 1489 1490 ! skip modes with too low or negative frequency -> Bose factor explodes (eg acoustic at Gamma) 1491 ! TODO: check the convergence wrt the tolerance 1492 ! several philosophies to be implemented for the unstable modes: 1493 ! 1) ignore 1494 ! 2) populate them according to a default amplitude 1495 ! 3) populate according to their modulus squared 1496 if (abs(phfrq_allq(imode, iq))<tol6) cycle 1497 1498 phdispl1 = phdispl_allq(:,:,:,imode,iq) 1499 1500 ! loop over configurations 1501 do iconfig = 1, nconfig 1502 1503 ! trick supercell object into using present q point 1504 thm_scells(iconfig)%qphon(:) = qibz(:,iq) 1505 1506 ! find thermal displacement amplitude eq 4 of Zacharias 1507 ! combined with l_nu,q expression in paragraph before 1508 if (phfrq_allq(imode, iq) > tol6) then 1509 sigma = sqrt( (bose_einstein(phfrq_allq(imode,iq), temperature) + half)/phfrq_allq(imode,iq)) 1510 else 1511 !Treat negative frequencies 1512 select case (option) 1513 case(1) 1514 !Do not populate 1515 sigma = 0._dp 1516 case(2) 1517 !Default amplitude for all the frequencies 1518 sigma = 100._dp 1519 case(3) 1520 !Absolute value of the frequencies 1521 sigma=sqrt((bose_einstein(abs(phfrq_allq(imode,iq)),temperature)+half)/& 1522 & abs(phfrq_allq(imode,iq))) 1523 case(4) 1524 sigma = 0._dp 1525 !Search if the amplitude of this unstable phonon is in the input argument amplitudes 1526 do iampl=1,namplitude 1527 if(abs(thm_scells(iconfig)%qphon(1) - amplitudes(1,iampl)) < tol8.and.& 1528 & abs(thm_scells(iconfig)%qphon(2) - amplitudes(2,iampl)) < tol8.and.& 1529 & abs(thm_scells(iconfig)%qphon(3) - amplitudes(3,iampl)) < tol8.and.& 1530 & abs(imode - amplitudes(4,iampl)) < tol8) then 1531 sigma = amplitudes(5,iampl) 1532 end if 1533 end do 1534 !If not, the amplitude is zero 1535 if(abs(sigma) < tol8)then 1536 write (msg, '(a,I0,a,3es12.5,2a,I0)') ' The amplitude of the unstable mode ',& 1537 & int(imode),' of the qpt ',thm_scells(iconfig)%qphon(:), ch10,& 1538 & 'is set to zero for the configuration ',iconfig 1539 ABI_WARNING(msg) 1540 end if 1541 end select 1542 end if 1543 1544 ! add displacement for this mode to supercell positions eq 5 of Zacharias 1545 call RANDOM_NUMBER(rand) 1546 rand = two * rand - one 1547 1548 ! from TDEP documentation for gaussian distribution of displacements 1549 !call RANDOM_NUMBER(rand1) 1550 !call RANDOM_NUMBER(rand2) 1551 ! rand = sqrt(-two*rand1) * sin(twopi*rand2) 1552 1553 ! if (rand > half) then 1554 ! rand = one 1555 ! else 1556 ! rand = -one 1557 ! end if 1558 1559 freeze_displ = rand * sigma 1560 1561 call freeze_displ_supercell (phdispl1(:,:,:), freeze_displ, thm_scells(iconfig)) 1562 end do !iconfig 1563 end do !imode 1564 end do !iq 1565 1566 ABI_FREE(phfrq_allq) 1567 ABI_FREE(phdispl_allq) 1568 ABI_FREE(phdispl1) 1569 ABI_FREE(qibz) 1570 ABI_FREE(qbz) 1571 ABI_FREE(wtqibz) 1572 1573 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
SOURCE
1671 subroutine thermal_supercell_print(fname, nconfig, temperature_K, thm_scells) 1672 1673 !Arguments ------------------------------------ 1674 !scalars 1675 integer, intent(in) :: nconfig 1676 type(supercell_type), intent(in) :: thm_scells(nconfig) 1677 character(len=fnlen), intent(in) :: fname 1678 real(dp), intent(in) :: temperature_K 1679 1680 ! local 1681 integer :: iconfig,itemp 1682 character(len=80) :: title1, title2 1683 character(len=fnlen) :: filename 1684 character(len=10) :: config_str 1685 1686 do iconfig = 1, nconfig 1687 write (config_str,'(I8)') iconfig 1688 write (filename, '(3a)') trim(fname), "_cf_", trim(adjustl(config_str)) 1689 write (title1, '(a,I6,a)') "# thermalized supercell at temperature T= ", temperature_K, " Kelvin" 1690 title2 = "# generated with random thermal displacements of all phonons" 1691 call prt_supercell (filename, thm_scells(itemp), title1, title2) 1692 end do 1693 1694 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
SOURCE
1234 subroutine zacharias_supercell_make(Crystal, Ifc, ntemper, rlatt, tempermin, temperinc, thm_scells) 1235 1236 !Arguments ------------------------------------ 1237 !scalars 1238 integer, intent(in) :: ntemper 1239 integer, intent(in) :: rlatt(3,3) 1240 real(dp), intent(in) :: tempermin, temperinc 1241 type(crystal_t),intent(in) :: Crystal 1242 type(ifc_type),intent(in) :: Ifc 1243 type(supercell_type), intent(out) :: thm_scells(ntemper) 1244 1245 !Local variables------------------------------- 1246 !scalars 1247 integer :: iq, nqibz, nqbz, qptopt1, imode, itemper, ierr, jmode 1248 real(dp) :: temperature_K, temperature, modesign, sigma, freeze_displ 1249 !arrays 1250 integer, allocatable :: modeindex(:) 1251 real(dp), allocatable :: qshft(:,:) ! dummy with 2 dimensions for call to kpts_ibz_from_kptrlatt 1252 real(dp), allocatable :: qbz(:,:), qibz(:,:), wtq_ibz(:) 1253 real(dp), allocatable :: phfrq_allq(:), phdispl_allq(:,:,:,:,:) 1254 real(dp), allocatable :: phfrq(:), phdispl(:,:,:,:),pheigvec(:,:,:,:) 1255 real(dp), allocatable :: phdispl1(:,:,:) 1256 character (len=500) :: msg 1257 1258 ! ************************************************************************* 1259 1260 ! check inputs 1261 ! TODO: add check that all rlatt are the same on input 1262 1263 if (rlatt(1,2)/=0 .or. rlatt(1,3)/=0 .or. rlatt(2,3)/=0 .or. & 1264 rlatt(2,1)/=0 .or. rlatt(3,1)/=0 .or. rlatt(3,2)/=0) then 1265 write (msg, '(4a, 9I6, a)') ' for the moment I have not implemented ', & 1266 ' non diagonal supercells.',ch10,' rlatt for temp 1 = ', rlatt, ' Returning ' 1267 ABI_WARNING(msg) 1268 return 1269 end if 1270 1271 ! build qpoint grid used for the Fourier interpolation. 1272 !(use no syms for the moment!) 1273 qptopt1 = 3 1274 1275 ! for the moment do not allow shifted q grids. 1276 ! We are interpolating anyway, so it will always work 1277 ABI_MALLOC(qshft,(3,1)) 1278 qshft(:,1)=zero 1279 1280 ! This call will set nqibz, IBZ and BZ arrays 1281 call kpts_ibz_from_kptrlatt(crystal, rlatt, qptopt1, 1, qshft, & 1282 nqibz, qibz, wtq_ibz, nqbz, qbz) ! new_kptrlatt, new_shiftk) ! Optional 1283 ABI_FREE(qshft) 1284 1285 ! allocate arrays with all of the q, omega, and displacement vectors 1286 ABI_MALLOC_OR_DIE(phfrq_allq, (3*Crystal%natom*nqibz), ierr) 1287 ABI_MALLOC_OR_DIE(phdispl_allq, (2, 3, Crystal%natom, 3*Crystal%natom, nqibz), ierr) 1288 1289 ABI_MALLOC(phfrq, (3*Crystal%natom)) 1290 ABI_MALLOC(phdispl, (2, 3, Crystal%natom, 3*Crystal%natom)) 1291 ABI_MALLOC(pheigvec, (2, 3, Crystal%natom, 3*Crystal%natom)) 1292 1293 ! loop over q to get all frequencies and displacement vectors 1294 ABI_MALLOC(modeindex, (nqibz*3*Crystal%natom)) 1295 imode = 0 1296 do iq = 1, nqibz 1297 ! Fourier interpolation. 1298 call ifc%fourq(Crystal, qibz(:,iq), phfrq, phdispl, out_eigvec=pheigvec) 1299 phfrq_allq((iq-1)*3*Crystal%natom+1 : iq*3*Crystal%natom) = phfrq 1300 phdispl_allq(1:2, 1:3, 1:Crystal%natom, 1:3*Crystal%natom, iq) = phdispl 1301 do jmode = 1, 3*Crystal%natom 1302 imode = imode + 1 1303 modeindex(imode) = imode 1304 end do 1305 end do 1306 ABI_FREE(phfrq) 1307 ABI_FREE(pheigvec) 1308 ABI_FREE(phdispl) 1309 1310 ! sort modes in whole list: get indirect indexing for qbz and displ 1311 call sort_dp(nqibz*3*Crystal%natom, phfrq_allq, modeindex, tol10) 1312 ! NB: phfrq is sorted now, but displ and qibz will have to be indexed indirectly with modeindex 1313 1314 ! only diagonal supercell case for the moment 1315 do itemper = 1, ntemper 1316 call init_supercell(Crystal%natom, rlatt, Crystal%rprimd, Crystal%typat, Crystal%xcart, Crystal%znucl, thm_scells(itemper)) 1317 end do 1318 1319 ! precalculate phase factors??? 1320 1321 ABI_MALLOC(phdispl1, (2, 3, Crystal%natom)) 1322 ! for all modes at all q in whole list, sorted 1323 modesign=one 1324 do imode = 1, 3*Crystal%natom*nqibz 1325 ! skip modes with too low or negative frequency -> Bose factor explodes (eg acoustic at Gamma) 1326 if (phfrq_allq(imode) < tol10) cycle 1327 1328 iq = ceiling(dble(modeindex(imode))/dble(3*Crystal%natom)) 1329 jmode = modeindex(imode) - (iq-1)*3*Crystal%natom 1330 phdispl1 = phdispl_allq(:,:,:,jmode,iq) 1331 1332 ! loop over temperatures 1333 do itemper = 1, ntemper 1334 temperature_K = tempermin + dble(itemper-1)*temperinc ! this is in Kelvin 1335 temperature = temperature_K / Ha_K !=315774.65_dp 1336 1337 ! trick supercell object into using present q point 1338 thm_scells(itemper)%qphon(:) = qibz(:,iq) 1339 1340 ! find thermal displacement amplitude eq 4 of Zacharias 1341 ! combined with l_nu,q expression in paragraph before 1342 sigma = sqrt( (bose_einstein(phfrq_allq(imode), temperature) + half)/phfrq_allq(imode) ) 1343 1344 ! add displacement for this mode to supercell positions eq 5 of Zacharias 1345 freeze_displ = modesign * sigma 1346 call freeze_displ_supercell (phdispl1(:,:,:), freeze_displ, thm_scells(itemper)) 1347 1348 end do !itemper 1349 1350 ! this is the prescription: flip sign for each successive mode in full 1351 ! spectrum, to cancel electron phonon coupling to 1st order 1352 ! (hopeflly 3rd order as well) 1353 modesign=-modesign 1354 1355 end do !imode 1356 1357 ABI_FREE(modeindex) 1358 ABI_FREE(phfrq_allq) 1359 ABI_FREE(phdispl_allq) 1360 ABI_FREE(phdispl1) 1361 ABI_FREE(qibz) 1362 ABI_FREE(qbz) 1363 ABI_FREE(wtq_ibz) 1364 1365 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
SOURCE
1627 subroutine zacharias_supercell_print(fname, ntemper, tempermin, temperinc, thm_scells) 1628 1629 !Arguments ------------------------------------ 1630 !scalars 1631 integer, intent(in) :: ntemper 1632 real(dp), intent(in) :: tempermin 1633 real(dp), intent(in) :: temperinc 1634 type(supercell_type), intent(in) :: thm_scells(ntemper) 1635 character(len=fnlen), intent(in) :: fname 1636 1637 ! local 1638 integer :: itemp 1639 character(len=80) :: title1, title2 1640 character(len=fnlen) :: filename 1641 real(dp) :: temper 1642 character(len=10) :: temper_str 1643 1644 do itemp = 1, ntemper 1645 temper = dble(itemp-1)*temperinc+tempermin 1646 write (temper_str,'(I8)') int(temper) 1647 write (filename, '(3a)') trim(fname), "_T_", trim(adjustl(temper_str)) 1648 write (title1, '(3a)') "# Zacharias thermalized supercell at temperature T= ", trim(temper_str), " Kelvin" 1649 title2 = "# generated with alternating thermal displacements of all phonons" 1650 call prt_supercell (filename, thm_scells(itemp), title1, title2) 1651 end do 1652 1653 end subroutine zacharias_supercell_print