TABLE OF CONTENTS


ABINIT/m_gruneisen [ Modules ]

[ Top ] [ Modules ]

NAME

  m_gruneisen

FUNCTION

  Objects and methods to compute Gruneisen parameters with central finite differences
  of dynamical matrices obtained with different unit cell volumes.

COPYRIGHT

 Copyright (C) 2011-2018 ABINIT group (MG)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

22 #if defined HAVE_CONFIG_H
23 #include "config.h"
24 #endif
25 
26 #include "abi_common.h"
27 
28 MODULE m_gruneisen
29 
30  use defs_basis
31  use m_errors
32  use m_abicore
33  use m_xmpi
34  use m_crystal
35  use m_crystal_io
36  use m_tetrahedron
37  use m_ddb
38  use m_ddb_hdr
39  use m_ifc
40  use m_cgtools
41  use m_nctk
42 #ifdef HAVE_NETCDF
43  use netcdf
44 #endif
45 
46  use m_io_tools,            only : get_unit, open_file
47  use m_time,                only : cwtime
48  use m_fstrings,            only : sjoin, itoa, ltoa, ftoa, strcat
49  use m_numeric_tools,       only : central_finite_diff, arth
50  use m_kpts,                only : kpts_ibz_from_kptrlatt, tetra_from_kptrlatt
51  use m_bz_mesh,             only : kpath_t, kpath_new, kpath_free, kpath_print
52  use m_anaddb_dataset,      only : anaddb_dataset_type
53  use m_dynmat,              only : massmult_and_breaksym, dfpt_phfrq, gtdyn9
54 
55  implicit none
56 
57  private

m_gruneisen/gruns_anaddb [ Functions ]

[ Top ] [ m_gruneisen ] [ Functions ]

NAME

  gruns_anaddb

FUNCTION

  Driver routine called in anaddb to compute Gruneisen parameters.

INPUTS

  inp<anaddb_dataset_type: :: anaddb input variables.
  prefix=Prefix for output files
  comm=MPI communicator

OUTPUT

  Only writing.

PARENTS

      anaddb

CHILDREN

      cwtime,gruns_free,gruns_qmesh,gruns_qpath,ifc_calcnwrite_nana_terms
      ifc_speedofsound,kpath_free,wrtout

SOURCE

811 subroutine gruns_anaddb(inp, prefix, comm)
812 
813 
814 !This section has been created automatically by the script Abilint (TD).
815 !Do not modify the following lines by hand.
816 #undef ABI_FUNC
817 #define ABI_FUNC 'gruns_anaddb'
818 !End of the abilint section
819 
820  implicit none
821 
822 !Arguments ------------------------------------
823  integer,intent(in) :: comm
824  character(len=*),intent(in) :: prefix
825  type(anaddb_dataset_type) :: inp
826 
827 !Local variables-------------------------------
828 !scalars
829  integer,parameter :: master=0
830  integer :: ii,nprocs,my_rank,ncid,iv0
831 #ifdef HAVE_NETCDF
832  integer :: ncerr
833 #endif
834  real(dp) :: cpu,wall,gflops
835  type(gruns_t),target :: gruns
836  type(kpath_t) :: qpath
837  character(len=500) :: msg
838 !arrays
839 
840 ! ************************************************************************
841 
842  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
843 
844  ABI_CHECK(inp%ifcflag == 1, "Gruneisen requires ifcflag == 1")
845 
846  call cwtime(cpu, wall, gflops, "start")
847 
848  gruns = gruns_new(inp%gruns_ddbs, inp, comm)
849  iv0 = gruns%iv0
850 
851  ncid = nctk_noid
852 #ifdef HAVE_NETCDF
853  if (my_rank == master) then
854    NCF_CHECK_MSG(nctk_open_create(ncid, strcat(prefix, "_GRUNS.nc"), xmpi_comm_self), "Creating _GRUNS.nc")
855 
856    ! Write structure corresponding to iv0
857    NCF_CHECK(crystal_ncwrite(gruns%cryst_vol(iv0), ncid))
858 
859    ! Add important dimensions and additional metadata.
860    ncerr = nctk_def_dims(ncid, [ &
861      nctkdim_t("gruns_nvols", gruns%nvols), &
862      nctkdim_t('number_of_phonon_modes', gruns%natom3)], defmode=.True.)
863    NCF_CHECK(ncerr)
864    ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "gruns_iv0"])
865    NCF_CHECK(ncerr)
866 
867    ! Add lattice parameters and positions of the `gruns_nvols` structures so
868    ! that we can easily reconstruct structure objects in AbiPy.
869    ncerr = nctk_def_arrays(ncid, [ &
870      nctkarr_t("gruns_rprimd", "dp", "three, three, gruns_nvols"), &
871      nctkarr_t("gruns_xred", "dp", "three, number_of_atoms, gruns_nvols") &
872    ])
873    NCF_CHECK(ncerr)
874 
875    ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: "gruns_iv0"], [iv0], datamode=.True.)
876    NCF_CHECK(ncerr)
877 
878    do ii=1,gruns%nvols
879      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_rprimd"), gruns%cryst_vol(ii)%rprimd, start=[1,1,ii]))
880      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_xred"), gruns%cryst_vol(ii)%xred, start=[1,1,ii]))
881    end do
882 
883    !call phonons_ncwrite(ncid,natom,nfineqpath,save_qpoints,weights,save_phfrq,save_phdispl_cart)
884    !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'atomic_mass_units'), ddb%amu))
885  end if
886 #endif
887 
888  ! Compute grunesein parameters on the q-mesh.
889  if (all(inp%ng2qpt /= 0)) then
890    call gruns_qmesh(gruns, prefix, inp%dosdeltae, inp%ng2qpt, 1, inp%q2shft, ncid, comm)
891  else
892    MSG_WARNING("Cannot compute Grunesein parameters on q-mesh because ng2qpt == 0")
893  end if
894 
895  ! Compute grunesein on the q-path.
896  if (inp%nqpath /= 0) then
897    qpath = kpath_new(inp%qpath, gruns%cryst_vol(iv0)%gprimd, inp%ndivsm)
898    call gruns_qpath(gruns, prefix, qpath, ncid, comm)
899    call kpath_free(qpath)
900  else
901    MSG_WARNING("Cannot compute Grunesein parameters on q-path because nqpath == 0")
902  end if
903 
904  ! Compute speed of sound for V0.
905  if (inp%vs_qrad_tolkms(1) > zero) then
906    call ifc_speedofsound(gruns%ifc_vol(iv0), gruns%cryst_vol(iv0), inp%vs_qrad_tolkms, ncid, comm)
907  end if
908 
909  ! Now treat the second list of vectors (only at the Gamma point, but can include non-analyticities)
910  if (my_rank == master .and. inp%nph2l /= 0 .and. inp%ifcflag == 1) then
911    call ifc_calcnwrite_nana_terms(gruns%ifc_vol(iv0), gruns%cryst_vol(iv0), inp%nph2l, inp%qph2l, inp%qnrml2, ncid)
912  end if
913 
914 #ifdef HAVE_NETCDF
915  if (my_rank == master) then
916    NCF_CHECK(nf90_close(ncid))
917  end if
918 #endif
919 
920  call gruns_free(gruns)
921 
922  call cwtime(cpu,wall,gflops,"stop")
923  write(msg,'(2(a,f8.2))')" gruns_anaddb completed. cpu:",cpu,", wall:",wall
924  call wrtout(std_out, msg, do_flush=.True.)
925 
926 end subroutine gruns_anaddb

m_gruneisen/gruns_fourq [ Functions ]

[ Top ] [ m_gruneisen ] [ Functions ]

NAME

  gruns_fourq

FUNCTION

  Compute grunesein parameters at an arbitrary q-point.

INPUTS

  qpt(3)=q-point in reduced coordinates.

OUTPUT

  wvols(3*natom, nvols) = Phonon frequencies for the different volumen.
  gvals(3*natom)=Gruneisen parameters evaluated at V0.
  dwdq(3,3*natom)=Group velocities at V0 in Cartesian coordinates.
  phdispl_cart(2, natom3, natom3, nvols)=Phonon displacement in Cartesian coordinates for the different volumes

NOTES

  The Grunesein parameter is given by:

     gamma(q,nu) = - (V / w(q,nu)) dw(q,nu)/dV

  Using w*2 = <u|D|u> and the Hellmann-Feynmann theorem, one obtains:

     gamma(q,nu) = - (V / 2 w(q,nu)**2) <u(q,nu)|dD(q)/dV|u(q,nu)>

  The derivative dD/dV is computed via central finite difference.

PARENTS

      m_gruneisen

CHILDREN

      cwtime,gruns_free,gruns_qmesh,gruns_qpath,ifc_calcnwrite_nana_terms
      ifc_speedofsound,kpath_free,wrtout

SOURCE

268 subroutine gruns_fourq(gruns, qpt, wvols, gvals, dwdq, phdispl_cart)
269 
270 
271 !This section has been created automatically by the script Abilint (TD).
272 !Do not modify the following lines by hand.
273 #undef ABI_FUNC
274 #define ABI_FUNC 'gruns_fourq'
275 !End of the abilint section
276 
277  implicit none
278 
279 !Arguments ------------------------------------
280 !scalars
281  type(gruns_t),intent(in) :: gruns
282 !arrays
283  real(dp),intent(in) :: qpt(3)
284  real(dp),intent(out) :: wvols(gruns%natom3,gruns%nvols),gvals(gruns%natom3),dwdq(3,gruns%natom3)
285  real(dp),intent(out) :: phdispl_cart(2, gruns%natom3, gruns%natom3, gruns%nvols)
286 
287 !Local variables-------------------------------
288 !scalars
289  integer :: ivol,natom3,nu
290  real(dp) :: fact
291 !arrays
292  real(dp) :: dot(2)
293  real(dp) :: eigvec(2,gruns%natom3,gruns%natom3,gruns%nvols),d2cart(2,gruns%natom3,gruns%natom3,gruns%nvols)
294  real(dp) :: dddv(2,gruns%natom3,gruns%natom3)
295  real(dp) :: omat(2,gruns%natom3, gruns%natom3)
296 
297 ! ************************************************************************
298 
299  natom3 = gruns%natom3
300  do ivol=1,gruns%nvols
301    if (ivol == gruns%iv0) then
302      ! Compute group velocities for V=V0
303      call ifc_fourq(gruns%ifc_vol(ivol), gruns%cryst_vol(ivol), qpt, wvols(:,ivol), phdispl_cart(:,:,:,ivol), &
304                     out_d2cart=d2cart(:,:,:,ivol), out_eigvec=eigvec(:,:,:,ivol), dwdq=dwdq)
305    else
306      call ifc_fourq(gruns%ifc_vol(ivol), gruns%cryst_vol(ivol), qpt, wvols(:,ivol), phdispl_cart(:,:,:,ivol), &
307                     out_d2cart=d2cart(:,:,:,ivol), out_eigvec=eigvec(:,:,:,ivol))
308    end if
309 
310    call massmult_and_breaksym(gruns%cryst_vol(ivol)%natom, gruns%cryst_vol(ivol)%ntypat, &
311      gruns%cryst_vol(ivol)%typat, gruns%ifc_vol(ivol)%amu, d2cart(:,:,:,ivol))
312 
313    !call zgemm('N','N',natom3,natom3,natom3,cone,d2cart(:,:,:,ivol),natom3,eigvec(:,:,:,ivol),natom3,czero,omat,natom3)
314    !do nu=1,natom3
315    !  write(std_out,*)"H|psi> - w**2 |psi>",maxval(abs(omat(:,:,nu) - wvols(nu,ivol) ** 2 * eigvec(:,:,nu,ivol)))
316    !end do
317  end do
318 
319  ! Compute dD(q)/dV with central finite difference.
320  dddv = zero
321  do ivol=1,gruns%nvols
322    fact = central_finite_diff(1, ivol, gruns%nvols)
323    if (fact /= zero) dddv = dddv + fact * d2cart(:,:,:,ivol)
324  end do
325  dddv = dddv / gruns%delta_vol
326 
327  ! Compute -V0/(2w(q)**2) <u(q)|dD(q)/dq|u(q)>
328  call zgemm('N','N',natom3,natom3,natom3,cone,dddv,natom3,eigvec(:,:,:,gruns%iv0),natom3,czero,omat,natom3)
329  do nu=1,natom3
330    if (abs(wvols(nu, gruns%iv0)) > tol12) then
331      dot = cg_zdotc(natom3, eigvec(1,1,nu, gruns%iv0), omat(1,1,nu))
332      ! Must change sign if we have a purely imaginary solution i.e. z = iw
333      gvals(nu) = -sign(one, wvols(nu, gruns%iv0)) * gruns%v0 * dot(1) / (two * wvols(nu, gruns%iv0)**2)
334    else
335      gvals(nu) = zero
336    end if
337  end do
338 
339 end subroutine gruns_fourq

m_gruneisen/gruns_free [ Functions ]

[ Top ] [ m_gruneisen ] [ Functions ]

NAME

  gruns_free

FUNCTION

  Free dynamic memory.

PARENTS

      m_gruneisen

CHILDREN

      cwtime,gruns_free,gruns_qmesh,gruns_qpath,ifc_calcnwrite_nana_terms
      ifc_speedofsound,kpath_free,wrtout

SOURCE

740 subroutine gruns_free(gruns)
741 
742 
743 !This section has been created automatically by the script Abilint (TD).
744 !Do not modify the following lines by hand.
745 #undef ABI_FUNC
746 #define ABI_FUNC 'gruns_free'
747 !End of the abilint section
748 
749  implicit none
750 
751 !Arguments ------------------------------------
752 !array
753  type(gruns_t),intent(inout) :: gruns
754 
755 !Local variables-------------------------------
756 !scalars
757  integer :: ii
758 
759 ! ************************************************************************
760 
761  if (allocated(gruns%ifc_vol)) then
762    do ii=1,size(gruns%cryst_vol)
763      call crystal_free(gruns%cryst_vol(ii))
764    end do
765    ABI_FREE(gruns%cryst_vol)
766  end if
767 
768  if (allocated(gruns%ddb_vol)) then
769    do ii=1,size(gruns%ddb_vol)
770      call ddb_free(gruns%ddb_vol(ii))
771    end do
772    ABI_FREE(gruns%ddb_vol)
773  end if
774 
775  if (allocated(gruns%ifc_vol)) then
776    do ii=1,size(gruns%ifc_vol)
777      call ifc_free(gruns%ifc_vol(ii))
778    end do
779    ABI_FREE(gruns%ifc_vol)
780  end if
781 
782 end subroutine gruns_free

m_gruneisen/gruns_new [ Functions ]

[ Top ] [ m_gruneisen ] [ Functions ]

NAME

  gruns_new

FUNCTION

  Construct new object from a list of DDB files.

INPUTS

  ddb_paths(:)=Paths of the DDB files (must be ordered by volume)
  inp<anaddb_dataset_type>=Anaddb dataset with input variables
  comm=MPI communicator

PARENTS

CHILDREN

SOURCE

131 type(gruns_t) function gruns_new(ddb_paths, inp, comm) result(new)
132 
133 
134 !This section has been created automatically by the script Abilint (TD).
135 !Do not modify the following lines by hand.
136 #undef ABI_FUNC
137 #define ABI_FUNC 'gruns_new'
138 !End of the abilint section
139 
140  implicit none
141 
142 !Arguments ------------------------------------
143  integer,intent(in) :: comm
144  type(anaddb_dataset_type),intent(in) :: inp
145 !arrays
146  character(len=*),intent(in) :: ddb_paths(:)
147 
148 !Local variables-------------------------------
149  integer,parameter :: natifc0=0,master=0
150  integer :: ivol,iblock,natom,ddbun
151  integer :: nprocs,my_rank,ierr
152  character(len=500) :: msg
153  type(ddb_hdr_type) :: ddb_hdr
154 !arrays
155  integer,allocatable :: atifc0(:)
156  real(dp) :: dielt(3,3)
157  real(dp),allocatable :: zeff(:,:,:)
158 
159 ! ************************************************************************
160 
161  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
162 
163  new%nvols = size(ddb_paths)
164  ABI_MALLOC(new%cryst_vol, (new%nvols))
165  ABI_MALLOC(new%ddb_vol, (new%nvols))
166  ABI_MALLOC(new%ifc_vol, (new%nvols))
167 
168  call wrtout(ab_out, "Computation of Gruneisen parameter with central finite difference:")
169 
170  ddbun = get_unit()
171  do ivol=1,new%nvols
172    call wrtout(ab_out, sjoin(" Reading DDB file:", ddb_paths(ivol)))
173    call ddb_hdr_open_read(ddb_hdr,ddb_paths(ivol),ddbun,DDB_VERSION,&
174 &                         dimonly=1)
175    natom = ddb_hdr%natom
176 
177    call ddb_hdr_free(ddb_hdr)
178 
179    ABI_MALLOC(atifc0, (natom))
180    atifc0 = 0
181    call ddb_from_file(new%ddb_vol(ivol), ddb_paths(ivol), inp%brav, natom, natifc0, atifc0, new%cryst_vol(ivol), comm)
182    ABI_FREE(atifc0)
183    if (my_rank == master) then
184      call crystal_print(new%cryst_vol(ivol), header=sjoin("Structure for ivol:", itoa(ivol)), unit=ab_out, prtvol=-1)
185    end if
186 
187    ! Get Dielectric Tensor and Effective Charges
188    ! (initialized to one_3D and zero if the derivatives are not available in the DDB file)
189    ABI_MALLOC(zeff, (3,3,natom))
190    iblock = ddb_get_dielt_zeff(new%ddb_vol(ivol), new%cryst_vol(ivol), inp%rfmeth, inp%chneut, inp%selectz, dielt, zeff)
191    if (iblock == 0) then
192      call wrtout(ab_out, sjoin("- Cannot find dielectric tensor and Born effective charges in DDB file:", ddb_paths(ivol)))
193      call wrtout(ab_out, "Values initialized with zeros")
194    else
195      call wrtout(ab_out, sjoin("- Found dielectric tensor and Born effective charges in DDB file:", ddb_paths(ivol)))
196    end if
197 
198    call ifc_init(new%ifc_vol(ivol), new%cryst_vol(ivol), new%ddb_vol(ivol),&
199      inp%brav,inp%asr,inp%symdynmat,inp%dipdip,inp%rfmeth,inp%ngqpt(1:3),inp%nqshft,inp%q1shft,dielt,zeff,&
200      inp%nsphere,inp%rifcsph,inp%prtsrlr,inp%enunit,comm)
201    ABI_FREE(zeff)
202  end do
203 
204  ! Consistency check
205  ! TODO: Add more tests.
206  ABI_CHECK(any(new%nvols == [3, 5, 7, 9]), "Central finite difference requires [3,5,7,9] DDB files")
207 
208  new%natom3 = 3 * new%cryst_vol(1)%natom
209  new%iv0 = 1 + new%nvols / 2
210  new%v0 = new%cryst_vol(new%iv0)%ucvol
211  new%delta_vol = new%cryst_vol(new%iv0+1)%ucvol - new%cryst_vol(new%iv0)%ucvol
212 
213  ierr = 0
214  do ivol=1,new%nvols
215    if (abs(new%cryst_vol(1)%ucvol + new%delta_vol * (ivol-1) - new%cryst_vol(ivol)%ucvol) > tol4) then
216       write(std_out,*)"ucvol, delta_vol, diff_vol", new%cryst_vol(ivol)%ucvol, new%delta_vol, &
217         abs(new%cryst_vol(1)%ucvol + new%delta_vol * (ivol-1) - new%cryst_vol(ivol)%ucvol)
218       ierr = ierr + 1
219    end if
220  end do
221  if (ierr /= 0) then
222    msg = ltoa([(new%cryst_vol(ivol)%ucvol, ivol=1,new%nvols)])
223    MSG_ERROR(sjoin("Gruneisen calculations requires linear mesh of volumes but received:", msg))
224  end if
225 
226 end function gruns_new

m_gruneisen/gruns_qmesh [ Functions ]

[ Top ] [ m_gruneisen ] [ Functions ]

NAME

  gruns_qmesh

FUNCTION

  Compute gruneisen parameters and group velocities on a q-mesh.
  Save results to file.

INPUTS

  prefix=Prefix for output files.
  dosdeltae=Step for the frequency mesh.
  ngqpt(3)=q-mesh divisions
  nshiftq=Number of shifts used to generated the ab-initio q-mesh.
  shiftq(3,nshiftq)=The shifts of the ab-initio q-mesh.
  ncid=netcdf file id.
  comm=MPI communicator

OUTPUT

PARENTS

      m_gruneisen

CHILDREN

      cwtime,gruns_free,gruns_qmesh,gruns_qpath,ifc_calcnwrite_nana_terms
      ifc_speedofsound,kpath_free,wrtout

SOURCE

510 subroutine gruns_qmesh(gruns, prefix, dosdeltae, ngqpt, nshiftq, shiftq, ncid, comm)
511 
512 
513 !This section has been created automatically by the script Abilint (TD).
514 !Do not modify the following lines by hand.
515 #undef ABI_FUNC
516 #define ABI_FUNC 'gruns_qmesh'
517 !End of the abilint section
518 
519  implicit none
520 
521 !Arguments ------------------------------------
522 !scalars
523  integer,intent(in) :: nshiftq,ncid,comm
524  real(dp),intent(in) :: dosdeltae !,dossmear
525  type(gruns_t),intent(in) :: gruns
526  character(len=*),intent(in) :: prefix
527 !arrays
528  integer,intent(in) :: ngqpt(3)
529  real(dp),intent(in) :: shiftq(3,nshiftq)
530 
531 !Local variables-------------------------------
532 !scalars
533  integer,parameter :: master=0,qptopt1=1,bcorr0=0
534  integer :: nprocs,my_rank,iqibz,nqbz,nqibz,ierr,ii,nu,ncerr,nomega,cnt,unt,io
535  real(dp) :: gavg,omega_min,omega_max,v2
536  type(t_tetrahedron) :: tetra
537  character(len=500) :: msg
538 !arrays
539  integer :: qptrlatt(3,3)
540  real(dp),allocatable :: gvals_qibz(:,:),wvols_qibz(:,:,:),dwdq_qibz(:,:,:)
541  real(dp),allocatable :: qibz(:,:),qbz(:,:),wtq(:)
542  real(dp),allocatable :: wdt(:,:),wdos(:,:),grdos(:,:),gr2dos(:,:),wibz(:),omega_mesh(:)
543  real(dp),allocatable :: vdos(:,:),v2dos(:,:)
544  real(dp),allocatable :: phdispl_cart_qibz(:,:,:,:,:)
545 
546 ! ************************************************************************
547 
548  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
549 
550  write(msg,'(a,(80a),4a)')ch10,('=',ii=1,80),ch10,ch10,' Calculation of Grunesein DOSes ',ch10
551  call wrtout(std_out, msg)
552  !call wrtout(ab_out, msg)
553 
554  ! Generate the q-mesh by finding the IBZ and the corresponding weights.
555  ABI_CHECK(all(ngqpt > 0), sjoin("invalid ngqpt:", ltoa(ngqpt)))
556  qptrlatt = 0
557  do ii=1,3
558    qptrlatt(ii,ii) = ngqpt(ii)
559  end do
560 
561  ! Get IBZ and BZ.
562  call kpts_ibz_from_kptrlatt(gruns%cryst_vol(gruns%iv0), qptrlatt, qptopt1, nshiftq, shiftq, &
563    nqibz, qibz, wtq, nqbz, qbz)
564 
565  ! Build tetrahedra
566  tetra = tetra_from_kptrlatt(gruns%cryst_vol(gruns%iv0), qptopt1, qptrlatt, nshiftq, shiftq, nqibz, qibz, msg, ierr)
567  if (ierr /= 0) MSG_ERROR(msg)
568 
569  ABI_CALLOC(wvols_qibz, (gruns%natom3, gruns%nvols, nqibz))
570  ABI_CALLOC(gvals_qibz, (gruns%natom3, nqibz))
571  ABI_CALLOC(dwdq_qibz, (3, gruns%natom3, nqibz))
572  ABI_CALLOC(phdispl_cart_qibz, (2, gruns%natom3, gruns%natom3, gruns%nvols, nqibz))
573 
574  gavg = zero
575  do iqibz=1,nqibz
576    if (mod(iqibz, nprocs) /= my_rank) cycle ! mpi-parallelism
577    call gruns_fourq(gruns, qibz(:,iqibz), wvols_qibz(:,:,iqibz), gvals_qibz(:,iqibz), &
578                     dwdq_qibz(:,:,iqibz), phdispl_cart_qibz(:,:,:,:,iqibz))
579    gavg = gavg + wtq(iqibz) * sum(gvals_qibz(:,iqibz))
580  end do
581  gavg = gavg / gruns%natom3
582 
583  call xmpi_sum(gavg, comm, ierr)
584  call xmpi_sum(wvols_qibz, comm, ierr)
585  call xmpi_sum(gvals_qibz, comm, ierr)
586  call xmpi_sum(dwdq_qibz, comm, ierr)
587  call xmpi_sum(phdispl_cart_qibz, comm, ierr)
588 
589  omega_min = gruns%ifc_vol(gruns%iv0)%omega_minmax(1)
590  omega_max = gruns%ifc_vol(gruns%iv0)%omega_minmax(2)
591  nomega = nint((omega_max - omega_min) / dosdeltae) + 1
592  nomega = max(6, nomega) ! Ensure Simpson integration will be ok
593 
594  ABI_MALLOC(omega_mesh, (nomega))
595  omega_mesh = arth(omega_min, dosdeltae, nomega)
596  omega_max = omega_mesh(nomega)
597  !write(std_out,*)"hello",omega_min,omega_max,dosdeltae,(omega_max-omega_min) / (nomega-1)
598  ABI_MALLOC(wibz, (nqibz))
599  ABI_MALLOC(wdt, (nomega, 2))
600  ABI_CALLOC(wdos, (nomega, 2))
601  ABI_CALLOC(grdos, (nomega, 2))
602  ABI_CALLOC(gr2dos, (nomega, 2))
603  ABI_CALLOC(vdos, (nomega, 2))
604  ABI_CALLOC(v2dos, (nomega, 2))
605 
606  ! Compute DOSes.
607  cnt = 0
608  do iqibz=1,nqibz
609    do nu=1,gruns%natom3
610      cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! mpi-parallelism
611      wibz = wvols_qibz(nu, gruns%iv0, :)
612      call tetra_get_onewk(tetra,iqibz,bcorr0,nomega,nqibz,wibz,omega_min,omega_max,one,wdt)
613      wdos = wdos + wdt
614      grdos = grdos + wdt * gvals_qibz(nu,iqibz)
615      gr2dos = gr2dos + wdt * gvals_qibz(nu,iqibz) ** 2
616      v2 = sum(dwdq_qibz(1:3,nu,iqibz) ** 2)
617      vdos = vdos + wdt * sqrt(v2)
618      v2dos = v2dos + wdt * v2
619    end do
620  end do
621 
622  call xmpi_sum(wdos, comm, ierr)
623  call xmpi_sum(grdos, comm, ierr)
624  call xmpi_sum(gr2dos, comm, ierr)
625  call xmpi_sum(vdos, comm, ierr)
626  call xmpi_sum(v2dos, comm, ierr)
627 
628  if (my_rank == master) then
629    call wrtout(ab_out, sjoin(" Average Gruneisen parameter:", ftoa(gavg, fmt="f8.5")))
630 
631    ! Write text files with Grunesein and DOSes.
632    if (open_file(strcat(prefix, "_GRUNS_DOS"), msg, newunit=unt, form="formatted", action="write") /= 0) then
633      MSG_ERROR(msg)
634    end if
635    write(unt,'(a)')'# Phonon density of states, Gruneisen DOS and phonon group velocity DOS'
636    write(unt,'(a)')"# Energy in Hartree, DOS in states/Hartree"
637    write(unt,'(a,i0)')'# Tetrahedron method with nqibz= ',nqibz
638    write(unt,"(a,f8.5)")"# Average Gruneisen parameter:", gavg
639    write(unt,'(5a)') &
640      "# omega PH_DOS Gruns_DOS Gruns**2_DOS Vel_DOS  Vel**2_DOS  PH_IDOS Gruns_IDOS Gruns**2_IDOS Vel_IDOS Vel**2_IDOS"
641    do io=1,nomega
642      write(unt, "(11es17.8)")omega_mesh, &
643        wdos(io,1), grdos(io,1), gr2dos(io,1), vdos(io,1), v2dos(io,1), &
644        wdos(io,2), grdos(io,2), gr2dos(io,2), vdos(io,2), v2dos(io,2)
645    end do
646    close(unt)
647  end if
648 
649 #ifdef HAVE_NETCDF
650  ! Write netcdf files.
651  if (my_rank == master .and. ncid /= nctk_noid) then
652    ncerr = nctk_def_dims(ncid, [ &
653      nctkdim_t("gruns_nqibz", nqibz), nctkdim_t('gruns_nshiftq', nshiftq), &
654      nctkdim_t('gruns_nomega', nomega)], defmode=.True.)
655    NCF_CHECK(ncerr)
656 
657    ncerr = nctk_def_arrays(ncid, [ &
658     ! q-point sampling in IBZ,
659     nctkarr_t("gruns_qptrlatt", "int", "three, three"), &
660     nctkarr_t("gruns_shiftq", "dp", "three, gruns_nshiftq"), &
661     nctkarr_t("gruns_qibz", "dp", "three, gruns_nqibz"), &
662     nctkarr_t("gruns_wtq", "dp", "gruns_nqibz"), &
663     ! grunesein parameters in IBZ
664     ! phonon frequencies at the different volumes,
665     ! group velocities at V0 in Cartesian coordinates.
666     nctkarr_t("gruns_gvals_qibz", "dp", "number_of_phonon_modes, gruns_nqibz"), &
667     nctkarr_t("gruns_wvols_qibz", "dp", "number_of_phonon_modes, gruns_nvols, gruns_nqibz"), &
668     nctkarr_t("gruns_dwdq_qibz", "dp", "three, number_of_phonon_modes, gruns_nqibz"), &
669     ! displacements for the different volumes.
670     nctkarr_t("gruns_phdispl_cart_qibz", "dp", &
671         "two, number_of_phonon_modes, number_of_phonon_modes, gruns_nvols, gruns_nqibz"), &
672     ! DOSes and IDOSes
673     nctkarr_t("gruns_omega_mesh", "dp", "gruns_nomega"), &
674     nctkarr_t("gruns_wdos", "dp", "gruns_nomega, two"), &
675     nctkarr_t("gruns_grdos", "dp", "gruns_nomega, two"), &
676     nctkarr_t("gruns_gr2dos", "dp", "gruns_nomega, two"), &
677     nctkarr_t("gruns_v2dos", "dp", "gruns_nomega, two"), &
678     nctkarr_t("gruns_vdos", "dp", "gruns_nomega, two") &
679    ])
680    NCF_CHECK(ncerr)
681 
682    ! Write data.
683    NCF_CHECK(nctk_set_datamode(ncid))
684    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_qptrlatt"), qptrlatt))
685    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_shiftq"), shiftq))
686    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_qibz"), qibz))
687    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_wtq"), wtq))
688    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_gvals_qibz"), gvals_qibz))
689    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_wvols_qibz"), wvols_qibz))
690    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_dwdq_qibz"), dwdq_qibz))
691    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_phdispl_cart_qibz"), phdispl_cart_qibz))
692    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_omega_mesh"), omega_mesh))
693    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_wdos"), wdos))
694    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_grdos"), grdos))
695    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_gr2dos"), gr2dos))
696    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_v2dos"), v2dos))
697    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_vdos"), vdos))
698  end if
699 #endif
700 
701  ABI_FREE(qibz)
702  ABI_FREE(wtq)
703  ABI_FREE(qbz)
704  ABI_FREE(wvols_qibz)
705  ABI_FREE(gvals_qibz)
706  ABI_FREE(dwdq_qibz)
707  ABI_FREE(phdispl_cart_qibz)
708  ABI_FREE(omega_mesh)
709  ABI_FREE(wibz)
710  ABI_FREE(wdt)
711  ABI_FREE(wdos)
712  ABI_FREE(grdos)
713  ABI_FREE(gr2dos)
714  ABI_FREE(v2dos)
715  ABI_FREE(vdos)
716 
717  call destroy_tetra(tetra)
718 
719 end subroutine gruns_qmesh

m_gruneisen/gruns_qpath [ Functions ]

[ Top ] [ m_gruneisen ] [ Functions ]

NAME

  gruns_qpath

FUNCTION

  Compute gruneisen parameters and group velocities on a q-path
  Write results to file.

INPUTS

  prefix=Prefix for output files.
  qpath<kpath_t>=Object describing the q-path.
  ncid=netcdf file id.
  comm=MPI communicator.

OUTPUT

  Only writing

PARENTS

      m_gruneisen

CHILDREN

      cwtime,gruns_free,gruns_qmesh,gruns_qpath,ifc_calcnwrite_nana_terms
      ifc_speedofsound,kpath_free,wrtout

SOURCE

370 subroutine gruns_qpath(gruns, prefix, qpath, ncid, comm)
371 
372 
373 !This section has been created automatically by the script Abilint (TD).
374 !Do not modify the following lines by hand.
375 #undef ABI_FUNC
376 #define ABI_FUNC 'gruns_qpath'
377 !End of the abilint section
378 
379  implicit none
380 
381 !Arguments ------------------------------------
382 !scalars
383  integer,intent(in) :: ncid,comm
384  type(gruns_t),intent(in) :: gruns
385  type(kpath_t),intent(in) :: qpath
386  character(len=*),intent(in) :: prefix
387 
388 !Local variables-------------------------------
389 !scalars
390  integer,parameter :: master=0
391  integer :: nprocs,my_rank,iqpt,ierr,ncerr,unt,iv0,nu,ii
392  character(len=500) :: msg
393 !arrays
394  real(dp),allocatable :: gvals_qpath(:,:),wvols_qpath(:,:,:),dwdq_qpath(:,:,:)
395  real(dp),allocatable :: phdispl_cart_qpath(:,:,:,:,:)
396 
397 ! ************************************************************************
398 
399  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
400 
401  write(msg,'(a,(80a),4a)')ch10,('=',ii=1,80),ch10,ch10,' Calculation of Grunesein parameters along q-path ',ch10
402  call wrtout(std_out, msg)
403  !call wrtout(ab_out, msg)
404 
405  iv0 = gruns%iv0
406  ABI_CALLOC(wvols_qpath, (gruns%natom3, gruns%nvols, qpath%npts))
407  ABI_CALLOC(gvals_qpath, (gruns%natom3, qpath%npts))
408  ABI_CALLOC(dwdq_qpath, (3, gruns%natom3, qpath%npts))
409  ABI_CALLOC(phdispl_cart_qpath, (2, gruns%natom3, gruns%natom3, gruns%nvols, qpath%npts))
410 
411  do iqpt=1,qpath%npts
412    if (mod(iqpt, nprocs) /= my_rank) cycle ! mpi-parallelism
413    call gruns_fourq(gruns, qpath%points(:,iqpt), wvols_qpath(:,:,iqpt), gvals_qpath(:,iqpt), &
414                     dwdq_qpath(:,:,iqpt), phdispl_cart_qpath(:,:,:,:,iqpt))
415  end do
416 
417  call xmpi_sum(wvols_qpath, comm, ierr)
418  call xmpi_sum(gvals_qpath, comm, ierr)
419  call xmpi_sum(dwdq_qpath, comm, ierr)
420  call xmpi_sum(phdispl_cart_qpath, comm, ierr)
421 
422  ! Write text files with phonon frequencies and grunesein on the path.
423  if (my_rank == master) then
424    if (open_file(strcat(prefix, "_GRUNS_QPATH"), msg, newunit=unt, form="formatted", action="write") /= 0) then
425      MSG_ERROR(msg)
426    end if
427    write(unt,'(a)')'# Phonon band structure, Gruneseinen parameters and group velocity'
428    write(unt,'(a)')"# Energy in Hartree, DOS in states/Hartree"
429    call kpath_print(qpath, unit=unt, pre="#")
430    write(unt,'(5a)')&
431      "# phfreq(mode=1) grunesein(mode=1) velocity(mode=1) phfreq(mode=2) grunesein(mode=2) velocity(mode=2)"
432    do iqpt=1,qpath%npts
433      do nu=1,gruns%natom3
434        write(unt, "(3es17.8)", advance="no") &
435          wvols_qpath(nu, iv0, iqpt), gvals_qpath(nu, iqpt), sum(dwdq_qpath(1:3, nu, iqpt) ** 2)
436      end do
437      write(unt, "(a)")" "
438    end do
439    close(unt)
440  end if
441 
442 #ifdef HAVE_NETCDF
443  if (my_rank == master .and. ncid /= nctk_noid) then
444    ncerr = nctk_def_dims(ncid, [nctkdim_t("gruns_nqpath", qpath%npts)], defmode=.True.)
445    NCF_CHECK(ncerr)
446 
447    ncerr = nctk_def_arrays(ncid, [ &
448      ! q-points of the path
449      nctkarr_t("gruns_qpath", "dp", "three, gruns_nqpath"), &
450      ! grunesein parameters on the path
451      nctkarr_t("gruns_gvals_qpath", "dp", "number_of_phonon_modes, gruns_nqpath"), &
452      ! phonon frequencies at the different volumes
453      nctkarr_t("gruns_wvols_qpath", "dp", "number_of_phonon_modes, gruns_nvols, gruns_nqpath"), &
454      ! group velocities at V0 in Cartesian coordinates.
455      nctkarr_t("gruns_dwdq_qpath", "dp", "three, number_of_phonon_modes, gruns_nqpath"), &
456      ! displacements for the different volumes.
457      nctkarr_t("gruns_phdispl_cart_qpath", "dp", &
458          "two, number_of_phonon_modes, number_of_phonon_modes, gruns_nvols, gruns_nqpath") &
459    ])
460    NCF_CHECK(ncerr)
461 
462    ! Write data.
463    NCF_CHECK(nctk_set_datamode(ncid))
464    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_qpath"), qpath%points))
465    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_gvals_qpath"), gvals_qpath))
466    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_wvols_qpath"), wvols_qpath))
467    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_dwdq_qpath"), dwdq_qpath))
468    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_phdispl_cart_qpath"), phdispl_cart_qpath))
469  end if
470 #endif
471 
472  ABI_FREE(wvols_qpath)
473  ABI_FREE(gvals_qpath)
474  ABI_FREE(dwdq_qpath)
475  ABI_FREE(phdispl_cart_qpath)
476 
477 end subroutine gruns_qpath

m_gruneisen/gruns_t [ Types ]

[ Top ] [ m_gruneisen ] [ Types ]

NAME

 gruns_t

FUNCTION

  Contains the interatomic force constants for different volumes.
  Provides methods to compute phonon bandstructures and gruneisen parameters.

SOURCE

 70  type,public :: gruns_t
 71 
 72    integer :: natom3
 73     ! 3 * natom
 74 
 75    integer :: nvols
 76     ! Number of volumes.
 77 
 78    integer :: iv0
 79     ! Index of the DDB file corresponding to the equilibrium volume V0.
 80 
 81    real(dp) :: v0
 82     ! Equilibrium volume.
 83 
 84    real(dp) :: delta_vol
 85     ! Uniform grid spacing for finite difference.
 86 
 87    type(crystal_t),allocatable :: cryst_vol(:)
 88     ! cryst_vol(nvols)
 89     ! crystalline structure for the different volumes.
 90 
 91    type(ddb_type),allocatable :: ddb_vol(:)
 92     ! dbb_vol(nvols)
 93     ! DDB objects for the different volumes.
 94 
 95    type(ifc_type),allocatable :: ifc_vol(:)
 96     ! ifc_vol(nvols)
 97     ! interatomic force constants for the different volumes.
 98 
 99  end type gruns_t
100 
101  public :: gruns_new        ! Constructor.
102  public :: gruns_qpath      ! Compute Grunesein parameters on a q-path.
103  public :: gruns_qmesh      ! Compute Grunesein parameters on a q-mesh.
104  public :: gruns_free       ! Release memory.
105  public :: gruns_anaddb     ! Driver routine called in anaddb.