TABLE OF CONTENTS
- ABINIT/m_gruneisen
- m_gruneisen/gruns_anaddb
- m_gruneisen/gruns_fourq
- m_gruneisen/gruns_free
- m_gruneisen/gruns_new
- m_gruneisen/gruns_qmesh
- m_gruneisen/gruns_qpath
- m_gruneisen/gruns_t
ABINIT/m_gruneisen [ 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-2024 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 .
SOURCE
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 MODULE m_gruneisen 24 25 use defs_basis 26 use m_errors 27 use m_abicore 28 use m_xmpi 29 use m_crystal 30 use m_htetra 31 use m_ddb 32 use m_ddb_hdr 33 use m_ifc 34 use m_cgtools 35 use m_nctk 36 #ifdef HAVE_NETCDF 37 use netcdf 38 #endif 39 40 use m_io_tools, only : get_unit, open_file 41 use m_time, only : cwtime, cwtime_report 42 use m_fstrings, only : sjoin, itoa, ltoa, ftoa, strcat 43 use m_numeric_tools, only : central_finite_diff, arth 44 use m_kpts, only : kpts_ibz_from_kptrlatt, tetra_from_kptrlatt 45 use m_bz_mesh, only : kpath_t, kpath_new 46 use m_anaddb_dataset, only : anaddb_dataset_type 47 use m_dynmat, only : massmult_and_breaksym, dfpt_phfrq, gtdyn9 48 49 implicit none 50 51 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.
SOURCE
717 subroutine gruns_anaddb(inp, prefix, comm) 718 719 !Arguments ------------------------------------ 720 integer,intent(in) :: comm 721 character(len=*),intent(in) :: prefix 722 type(anaddb_dataset_type) :: inp 723 724 !Local variables------------------------------- 725 !scalars 726 integer,parameter :: master=0 727 integer :: ii,nprocs,my_rank,ncid,iv0 728 #ifdef HAVE_NETCDF 729 integer :: ncerr 730 #endif 731 real(dp) :: cpu,wall,gflops 732 type(gruns_t),target :: gruns 733 type(kpath_t) :: qpath 734 735 ! ************************************************************************ 736 737 nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 738 739 ABI_CHECK(inp%ifcflag == 1, "Gruneisen requires ifcflag == 1") 740 741 call cwtime(cpu, wall, gflops, "start") 742 743 gruns = gruns_new(inp%gruns_ddbs, inp, comm) 744 iv0 = gruns%iv0 745 746 ncid = nctk_noid 747 #ifdef HAVE_NETCDF 748 if (my_rank == master) then 749 NCF_CHECK_MSG(nctk_open_create(ncid, strcat(prefix, "_GRUNS.nc"), xmpi_comm_self), "Creating _GRUNS.nc") 750 751 ! Write structure corresponding to iv0 752 NCF_CHECK(gruns%cryst_vol(iv0)%ncwrite(ncid)) 753 754 ! Add important dimensions and additional metadata. 755 ncerr = nctk_def_dims(ncid, [ & 756 nctkdim_t("gruns_nvols", gruns%nvols), & 757 nctkdim_t('number_of_phonon_modes', gruns%natom3)], defmode=.True.) 758 NCF_CHECK(ncerr) 759 ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "gruns_iv0"]) 760 NCF_CHECK(ncerr) 761 762 ! Add lattice parameters and positions of the `gruns_nvols` structures so 763 ! that we can easily reconstruct structure objects in AbiPy. 764 ncerr = nctk_def_arrays(ncid, [ & 765 nctkarr_t("gruns_rprimd", "dp", "three, three, gruns_nvols"), & 766 nctkarr_t("gruns_xred", "dp", "three, number_of_atoms, gruns_nvols") & 767 ]) 768 NCF_CHECK(ncerr) 769 770 ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: "gruns_iv0"], [iv0], datamode=.True.) 771 NCF_CHECK(ncerr) 772 773 do ii=1,gruns%nvols 774 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_rprimd"), gruns%cryst_vol(ii)%rprimd, start=[1,1,ii])) 775 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_xred"), gruns%cryst_vol(ii)%xred, start=[1,1,ii])) 776 end do 777 778 !call phonons_ncwrite(ncid,natom,nfineqpath,save_qpoints,weights,save_phfrq,save_phdispl_cart) 779 !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'atomic_mass_units'), ddb%amu)) 780 end if 781 #endif 782 783 ! Compute gruneisen parameters on the q-mesh. 784 if (all(inp%ng2qpt /= 0)) then 785 call gruns_qmesh(gruns, prefix, inp%dosdeltae, inp%ng2qpt, 1, inp%q2shft, ncid, comm) 786 else 787 ABI_WARNING("Cannot compute Gruneisen parameters on q-mesh because ng2qpt == 0") 788 end if 789 790 ! Compute gruneisen on the q-path. 791 if (inp%nqpath /= 0) then 792 qpath = kpath_new(inp%qpath, gruns%cryst_vol(iv0)%gprimd, inp%ndivsm) 793 call gruns_qpath(gruns, prefix, qpath, ncid, comm) 794 call qpath%free() 795 else 796 ABI_WARNING("Cannot compute Gruneisen parameters on q-path because nqpath == 0") 797 end if 798 799 ! Compute speed of sound for V0. 800 if (inp%vs_qrad_tolkms(1) > zero) then 801 call gruns%ifc_vol(iv0)%speedofsound(gruns%cryst_vol(iv0), inp%vs_qrad_tolkms, ncid, comm) 802 end if 803 804 ! Now treat the second list of vectors (only at the Gamma point, but can include non-analyticities) 805 if (my_rank == master .and. inp%nph2l /= 0 .and. inp%ifcflag == 1) then 806 call gruns%ifc_vol(iv0)%calcnwrite_nana_terms(gruns%cryst_vol(iv0), inp%nph2l, inp%qph2l, inp%qnrml2, ncid) 807 end if 808 809 #ifdef HAVE_NETCDF 810 if (my_rank == master) then 811 NCF_CHECK(nf90_close(ncid)) 812 end if 813 #endif 814 815 call gruns_free(gruns) 816 call cwtime_report("gruns_anaddb", cpu, wall, gflops) 817 818 end subroutine gruns_anaddb
m_gruneisen/gruns_fourq [ Functions ]
[ Top ] [ m_gruneisen ] [ Functions ]
NAME
gruns_fourq
FUNCTION
Compute gruneisen 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 Gruneisen 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.
SOURCE
238 subroutine gruns_fourq(gruns, qpt, wvols, gvals, dwdq, phdispl_cart) 239 240 !Arguments ------------------------------------ 241 !scalars 242 type(gruns_t),intent(in) :: gruns 243 !arrays 244 real(dp),intent(in) :: qpt(3) 245 real(dp),intent(out) :: wvols(gruns%natom3,gruns%nvols),gvals(gruns%natom3),dwdq(3,gruns%natom3) 246 real(dp),intent(out) :: phdispl_cart(2, gruns%natom3, gruns%natom3, gruns%nvols) 247 248 !Local variables------------------------------- 249 !scalars 250 integer :: ivol,natom3,nu 251 real(dp) :: fact 252 !arrays 253 real(dp) :: dot(2) 254 real(dp) :: eigvec(2,gruns%natom3,gruns%natom3,gruns%nvols),d2cart(2,gruns%natom3,gruns%natom3,gruns%nvols) 255 real(dp) :: dddv(2,gruns%natom3,gruns%natom3) 256 real(dp) :: omat(2,gruns%natom3, gruns%natom3) 257 258 ! ************************************************************************ 259 260 natom3 = gruns%natom3 261 do ivol=1,gruns%nvols 262 if (ivol == gruns%iv0) then 263 ! Compute group velocities for V=V0 264 call gruns%ifc_vol(ivol)%fourq(gruns%cryst_vol(ivol), qpt, wvols(:,ivol), phdispl_cart(:,:,:,ivol), & 265 out_d2cart=d2cart(:,:,:,ivol), out_eigvec=eigvec(:,:,:,ivol), dwdq=dwdq) 266 else 267 call gruns%ifc_vol(ivol)%fourq(gruns%cryst_vol(ivol), qpt, wvols(:,ivol), phdispl_cart(:,:,:,ivol), & 268 out_d2cart=d2cart(:,:,:,ivol), out_eigvec=eigvec(:,:,:,ivol)) 269 end if 270 271 call massmult_and_breaksym(gruns%cryst_vol(ivol)%natom, gruns%cryst_vol(ivol)%ntypat, & 272 gruns%cryst_vol(ivol)%typat, gruns%ifc_vol(ivol)%amu, d2cart(:,:,:,ivol)) 273 274 !call zgemm('N','N',natom3,natom3,natom3,cone,d2cart(:,:,:,ivol),natom3,eigvec(:,:,:,ivol),natom3,czero,omat,natom3) 275 !do nu=1,natom3 276 ! write(std_out,*)"H|psi> - w**2 |psi>",maxval(abs(omat(:,:,nu) - wvols(nu,ivol) ** 2 * eigvec(:,:,nu,ivol))) 277 !end do 278 end do 279 280 ! Compute dD(q)/dV with central finite difference. 281 dddv = zero 282 do ivol=1,gruns%nvols 283 fact = central_finite_diff(1, ivol, gruns%nvols) 284 if (fact /= zero) dddv = dddv + fact * d2cart(:,:,:,ivol) 285 end do 286 dddv = dddv / gruns%delta_vol 287 288 ! Compute -V0/(2w(q)**2) <u(q)|dD(q)/dq|u(q)> 289 call zgemm('N','N',natom3,natom3,natom3,cone,dddv,natom3,eigvec(:,:,:,gruns%iv0),natom3,czero,omat,natom3) 290 do nu=1,natom3 291 if (abs(wvols(nu, gruns%iv0)) > tol12) then 292 dot = cg_zdotc(natom3, eigvec(1,1,nu, gruns%iv0), omat(1,1,nu)) 293 ! Must change sign if we have a purely imaginary solution i.e. z = iw 294 gvals(nu) = -sign(one, wvols(nu, gruns%iv0)) * gruns%v0 * dot(1) / (two * wvols(nu, gruns%iv0)**2) 295 else 296 gvals(nu) = zero 297 end if 298 end do 299 300 end subroutine gruns_fourq
m_gruneisen/gruns_free [ Functions ]
[ Top ] [ m_gruneisen ] [ Functions ]
NAME
gruns_free
FUNCTION
Free dynamic memory.
SOURCE
662 subroutine gruns_free(gruns) 663 664 !Arguments ------------------------------------ 665 !array 666 type(gruns_t),intent(inout) :: gruns 667 668 !Local variables------------------------------- 669 !scalars 670 integer :: ii 671 672 ! ************************************************************************ 673 674 if (allocated(gruns%ifc_vol)) then 675 do ii=1,size(gruns%cryst_vol) 676 call gruns%cryst_vol(ii)%free() 677 end do 678 ABI_FREE(gruns%cryst_vol) 679 end if 680 681 if (allocated(gruns%ddb_vol)) then 682 do ii=1,size(gruns%ddb_vol) 683 call gruns%ddb_vol(ii)%free() 684 end do 685 ABI_FREE(gruns%ddb_vol) 686 end if 687 688 if (allocated(gruns%ifc_vol)) then 689 do ii=1,size(gruns%ifc_vol) 690 call gruns%ifc_vol(ii)%free() 691 end do 692 ABI_FREE(gruns%ifc_vol) 693 end if 694 695 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_filepaths(:)=Paths of the DDB files (must be ordered by volume) inp<anaddb_dataset_type>=Anaddb dataset with input variables comm=MPI communicator
SOURCE
121 type(gruns_t) function gruns_new(ddb_filepaths, inp, comm) result(new) 122 123 !Arguments ------------------------------------ 124 integer,intent(in) :: comm 125 type(anaddb_dataset_type),intent(in) :: inp 126 !arrays 127 character(len=*),intent(in) :: ddb_filepaths(:) 128 129 !Local variables------------------------------- 130 integer,parameter :: master=0 131 integer :: ivol,iblock,natom,ddbun, nprocs,my_rank,ierr 132 character(len=500) :: msg 133 type(ddb_hdr_type) :: ddb_hdr 134 !arrays 135 real(dp) :: dielt(3,3) 136 real(dp),allocatable :: zeff(:,:,:), qdrp_cart(:,:,:,:) 137 138 ! ************************************************************************ 139 140 nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 141 142 new%nvols = size(ddb_filepaths) 143 ABI_MALLOC(new%cryst_vol, (new%nvols)) 144 ABI_MALLOC(new%ddb_vol, (new%nvols)) 145 ABI_MALLOC(new%ifc_vol, (new%nvols)) 146 147 call wrtout(ab_out, "Computation of Gruneisen parameter with central finite difference:") 148 149 ddbun = get_unit() 150 do ivol=1,new%nvols 151 call wrtout(ab_out, sjoin(" Reading DDB file:", ddb_filepaths(ivol))) 152 153 call new%ddb_vol(ivol)%from_file(ddb_filepaths(ivol), ddb_hdr, new%cryst_vol(ivol), comm) 154 call new%ddb_vol(ivol)%set_brav(inp%brav) 155 natom = ddb_hdr%natom 156 call ddb_hdr%free() 157 158 if (my_rank == master) then 159 call new%cryst_vol(ivol)%print(header=sjoin("Structure for ivol:", itoa(ivol)), unit=ab_out, prtvol=-1) 160 end if 161 162 ! Get Dielectric Tensor and Effective Charges 163 ! (initialized to one_3D and zero if the derivatives are not available in the DDB file) 164 ABI_MALLOC(zeff, (3,3,natom)) 165 ABI_CALLOC(qdrp_cart, (3,3,3,natom)) 166 iblock = new%ddb_vol(ivol)%get_dielt_zeff(new%cryst_vol(ivol), inp%rfmeth, inp%chneut, inp%selectz, dielt, zeff) 167 if (iblock == 0) then 168 call wrtout(ab_out, sjoin("- Cannot find dielectric tensor and Born effective charges in DDB file:", ddb_filepaths(ivol))) 169 call wrtout(ab_out, "Values initialized with zeros") 170 else 171 call wrtout(ab_out, sjoin("- Found dielectric tensor and Born effective charges in DDB file:", ddb_filepaths(ivol))) 172 end if 173 174 call new%ifc_vol(ivol)%init(new%cryst_vol(ivol), new%ddb_vol(ivol),& 175 inp%brav,inp%asr,inp%symdynmat,inp%dipdip,inp%rfmeth,inp%ngqpt(1:3),inp%nqshft,inp%q1shft,dielt,zeff,& 176 qdrp_cart,inp%nsphere,inp%rifcsph,inp%prtsrlr,inp%enunit,comm) 177 ABI_FREE(zeff) 178 ABI_FREE(qdrp_cart) 179 end do 180 181 ! Consistency check 182 ! TODO: Add more tests. 183 ABI_CHECK(any(new%nvols == [3, 5, 7, 9]), "Central finite difference requires [3,5,7,9] DDB files") 184 185 new%natom3 = 3 * new%cryst_vol(1)%natom 186 new%iv0 = 1 + new%nvols / 2 187 new%v0 = new%cryst_vol(new%iv0)%ucvol 188 new%delta_vol = new%cryst_vol(new%iv0+1)%ucvol - new%cryst_vol(new%iv0)%ucvol 189 190 ierr = 0 191 do ivol=1,new%nvols 192 if (abs(new%cryst_vol(1)%ucvol + new%delta_vol * (ivol-1) - new%cryst_vol(ivol)%ucvol) > tol4) then 193 write(std_out,*)"ucvol, delta_vol, diff_vol", new%cryst_vol(ivol)%ucvol, new%delta_vol, & 194 abs(new%cryst_vol(1)%ucvol + new%delta_vol * (ivol-1) - new%cryst_vol(ivol)%ucvol) 195 ierr = ierr + 1 196 end if 197 end do 198 if (ierr /= 0) then 199 msg = ltoa([(new%cryst_vol(ivol)%ucvol, ivol=1,new%nvols)]) 200 ABI_ERROR(sjoin("Gruneisen calculations requires linear mesh of volumes but received:", msg)) 201 end if 202 203 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
SOURCE
448 subroutine gruns_qmesh(gruns, prefix, dosdeltae, ngqpt, nshiftq, shiftq, ncid, comm) 449 450 !Arguments ------------------------------------ 451 !scalars 452 integer,intent(in) :: nshiftq,ncid,comm 453 real(dp),intent(in) :: dosdeltae !,dossmear 454 type(gruns_t),intent(in) :: gruns 455 character(len=*),intent(in) :: prefix 456 !arrays 457 integer,intent(in) :: ngqpt(3) 458 real(dp),intent(in) :: shiftq(3,nshiftq) 459 460 !Local variables------------------------------- 461 !scalars 462 integer,parameter :: master=0,qptopt1=1,bcorr0=0 463 integer :: nprocs,my_rank,iqibz,nqbz,nqibz,ierr,ii,nu,ncerr,nomega,cnt,unt,io 464 real(dp) :: gavg,omega_min,omega_max,v2 465 type(htetra_t) :: tetra 466 character(len=500) :: msg 467 !arrays 468 integer :: qptrlatt(3,3) 469 real(dp),allocatable :: gvals_qibz(:,:),wvols_qibz(:,:,:),dwdq_qibz(:,:,:) 470 real(dp),allocatable :: qibz(:,:),qbz(:,:),wtq(:) 471 real(dp),allocatable :: wdt(:,:),wdos(:,:),grdos(:,:),gr2dos(:,:),wibz(:),omega_mesh(:) 472 real(dp),allocatable :: vdos(:,:),v2dos(:,:) 473 real(dp),allocatable :: phdispl_cart_qibz(:,:,:,:,:) 474 475 ! ************************************************************************ 476 477 nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 478 479 write(msg,'(a,(80a),4a)')ch10,('=',ii=1,80),ch10,ch10,' Calculation of Gruneisen DOSes ',ch10 480 call wrtout(std_out, msg) 481 482 ! Generate the q-mesh by finding the IBZ and the corresponding weights. 483 ABI_CHECK(all(ngqpt > 0), sjoin("invalid ngqpt:", ltoa(ngqpt))) 484 qptrlatt = 0 485 do ii=1,3 486 qptrlatt(ii,ii) = ngqpt(ii) 487 end do 488 489 ! Get IBZ and BZ. 490 call kpts_ibz_from_kptrlatt(gruns%cryst_vol(gruns%iv0), qptrlatt, qptopt1, nshiftq, shiftq, & 491 nqibz, qibz, wtq, nqbz, qbz) 492 493 ! Build tetrahedra 494 tetra = tetra_from_kptrlatt(gruns%cryst_vol(gruns%iv0), qptopt1, qptrlatt, nshiftq, shiftq, nqibz, qibz, comm, msg, ierr) 495 if (ierr /= 0) ABI_ERROR(msg) 496 497 ABI_CALLOC(wvols_qibz, (gruns%natom3, gruns%nvols, nqibz)) 498 ABI_CALLOC(gvals_qibz, (gruns%natom3, nqibz)) 499 ABI_CALLOC(dwdq_qibz, (3, gruns%natom3, nqibz)) 500 ABI_CALLOC(phdispl_cart_qibz, (2, gruns%natom3, gruns%natom3, gruns%nvols, nqibz)) 501 502 gavg = zero 503 do iqibz=1,nqibz 504 if (mod(iqibz, nprocs) /= my_rank) cycle ! mpi-parallelism 505 call gruns_fourq(gruns, qibz(:,iqibz), wvols_qibz(:,:,iqibz), gvals_qibz(:,iqibz), & 506 dwdq_qibz(:,:,iqibz), phdispl_cart_qibz(:,:,:,:,iqibz)) 507 gavg = gavg + wtq(iqibz) * sum(gvals_qibz(:,iqibz)) 508 end do 509 gavg = gavg / gruns%natom3 510 511 call xmpi_sum(gavg, comm, ierr) 512 call xmpi_sum(wvols_qibz, comm, ierr) 513 call xmpi_sum(gvals_qibz, comm, ierr) 514 call xmpi_sum(dwdq_qibz, comm, ierr) 515 call xmpi_sum(phdispl_cart_qibz, comm, ierr) 516 517 omega_min = gruns%ifc_vol(gruns%iv0)%omega_minmax(1) 518 omega_max = gruns%ifc_vol(gruns%iv0)%omega_minmax(2) 519 nomega = nint((omega_max - omega_min) / dosdeltae) + 1 520 nomega = max(6, nomega) ! Ensure Simpson integration will be ok 521 522 ABI_MALLOC(omega_mesh, (nomega)) 523 omega_mesh = arth(omega_min, dosdeltae, nomega) 524 omega_max = omega_mesh(nomega) 525 !write(std_out,*)"hello",omega_min,omega_max,dosdeltae,(omega_max-omega_min) / (nomega-1) 526 ABI_MALLOC(wibz, (nqibz)) 527 ABI_MALLOC(wdt, (nomega, 2)) 528 ABI_CALLOC(wdos, (nomega, 2)) 529 ABI_CALLOC(grdos, (nomega, 2)) 530 ABI_CALLOC(gr2dos, (nomega, 2)) 531 ABI_CALLOC(vdos, (nomega, 2)) 532 ABI_CALLOC(v2dos, (nomega, 2)) 533 534 ! Compute DOSes. 535 cnt = 0 536 do iqibz=1,nqibz 537 do nu=1,gruns%natom3 538 cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! mpi-parallelism 539 wibz = wvols_qibz(nu, gruns%iv0, :) 540 call tetra%get_onewk(iqibz,bcorr0,nomega,nqibz,wibz,omega_min,omega_max,one,wdt) 541 wdt = wdt*wtq(iqibz) 542 wdos = wdos + wdt 543 grdos = grdos + wdt * gvals_qibz(nu,iqibz) 544 gr2dos = gr2dos + wdt * gvals_qibz(nu,iqibz) ** 2 545 v2 = sum(dwdq_qibz(1:3,nu,iqibz) ** 2) 546 vdos = vdos + wdt * sqrt(v2) 547 v2dos = v2dos + wdt * v2 548 end do 549 end do 550 551 call xmpi_sum(wdos, comm, ierr) 552 call xmpi_sum(grdos, comm, ierr) 553 call xmpi_sum(gr2dos, comm, ierr) 554 call xmpi_sum(vdos, comm, ierr) 555 call xmpi_sum(v2dos, comm, ierr) 556 557 if (my_rank == master) then 558 call wrtout(ab_out, sjoin(" Average Gruneisen parameter:", ftoa(gavg, fmt="f8.5"))) 559 560 ! Write text files with Gruneisen and DOSes. 561 if (open_file(strcat(prefix, "_GRUNS_DOS"), msg, newunit=unt, form="formatted", action="write") /= 0) then 562 ABI_ERROR(msg) 563 end if 564 write(unt,'(a)')'# Phonon density of states, Gruneisen DOS and phonon group velocity DOS' 565 write(unt,'(a)')"# Energy in Hartree, DOS in states/Hartree" 566 write(unt,'(a,i0)')'# Tetrahedron method with nqibz= ',nqibz 567 write(unt,"(a,f8.5)")"# Average Gruneisen parameter:", gavg 568 write(unt,'(5a)') & 569 "# 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" 570 do io=1,nomega 571 write(unt, "(11es17.8)")omega_mesh(io), & 572 wdos(io,1), grdos(io,1), gr2dos(io,1), vdos(io,1), v2dos(io,1), & 573 wdos(io,2), grdos(io,2), gr2dos(io,2), vdos(io,2), v2dos(io,2) 574 end do 575 close(unt) 576 end if 577 578 #ifdef HAVE_NETCDF 579 ! Write netcdf files. 580 if (my_rank == master .and. ncid /= nctk_noid) then 581 ncerr = nctk_def_dims(ncid, [ & 582 nctkdim_t("gruns_nqibz", nqibz), nctkdim_t('gruns_nshiftq', nshiftq), & 583 nctkdim_t('gruns_nomega', nomega)], defmode=.True.) 584 NCF_CHECK(ncerr) 585 586 ncerr = nctk_def_arrays(ncid, [ & 587 ! q-point sampling in IBZ, 588 nctkarr_t("gruns_qptrlatt", "int", "three, three"), & 589 nctkarr_t("gruns_shiftq", "dp", "three, gruns_nshiftq"), & 590 nctkarr_t("gruns_qibz", "dp", "three, gruns_nqibz"), & 591 nctkarr_t("gruns_wtq", "dp", "gruns_nqibz"), & 592 ! gruneisen parameters in IBZ 593 ! phonon frequencies at the different volumes, 594 ! group velocities at V0 in Cartesian coordinates. 595 nctkarr_t("gruns_gvals_qibz", "dp", "number_of_phonon_modes, gruns_nqibz"), & 596 nctkarr_t("gruns_wvols_qibz", "dp", "number_of_phonon_modes, gruns_nvols, gruns_nqibz"), & 597 nctkarr_t("gruns_dwdq_qibz", "dp", "three, number_of_phonon_modes, gruns_nqibz"), & 598 ! displacements for the different volumes. 599 nctkarr_t("gruns_phdispl_cart_qibz", "dp", & 600 "two, number_of_phonon_modes, number_of_phonon_modes, gruns_nvols, gruns_nqibz"), & 601 ! DOSes and IDOSes 602 nctkarr_t("gruns_omega_mesh", "dp", "gruns_nomega"), & 603 nctkarr_t("gruns_wdos", "dp", "gruns_nomega, two"), & 604 nctkarr_t("gruns_grdos", "dp", "gruns_nomega, two"), & 605 nctkarr_t("gruns_gr2dos", "dp", "gruns_nomega, two"), & 606 nctkarr_t("gruns_v2dos", "dp", "gruns_nomega, two"), & 607 nctkarr_t("gruns_vdos", "dp", "gruns_nomega, two") & 608 ]) 609 NCF_CHECK(ncerr) 610 611 ! Write data. 612 NCF_CHECK(nctk_set_datamode(ncid)) 613 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_qptrlatt"), qptrlatt)) 614 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_shiftq"), shiftq)) 615 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_qibz"), qibz)) 616 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_wtq"), wtq)) 617 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_gvals_qibz"), gvals_qibz)) 618 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_wvols_qibz"), wvols_qibz)) 619 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_dwdq_qibz"), dwdq_qibz)) 620 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_phdispl_cart_qibz"), phdispl_cart_qibz)) 621 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_omega_mesh"), omega_mesh)) 622 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_wdos"), wdos)) 623 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_grdos"), grdos)) 624 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_gr2dos"), gr2dos)) 625 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_v2dos"), v2dos)) 626 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_vdos"), vdos)) 627 end if 628 #endif 629 630 ABI_FREE(qibz) 631 ABI_FREE(wtq) 632 ABI_FREE(qbz) 633 ABI_FREE(wvols_qibz) 634 ABI_FREE(gvals_qibz) 635 ABI_FREE(dwdq_qibz) 636 ABI_FREE(phdispl_cart_qibz) 637 ABI_FREE(omega_mesh) 638 ABI_FREE(wibz) 639 ABI_FREE(wdt) 640 ABI_FREE(wdos) 641 ABI_FREE(grdos) 642 ABI_FREE(gr2dos) 643 ABI_FREE(v2dos) 644 ABI_FREE(vdos) 645 646 call tetra%free() 647 648 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
SOURCE
324 subroutine gruns_qpath(gruns, prefix, qpath, ncid, comm) 325 326 !Arguments ------------------------------------ 327 !scalars 328 integer,intent(in) :: ncid,comm 329 type(gruns_t),intent(in) :: gruns 330 type(kpath_t),intent(in) :: qpath 331 character(len=*),intent(in) :: prefix 332 333 !Local variables------------------------------- 334 !scalars 335 integer,parameter :: master=0 336 integer :: nprocs,my_rank,iqpt,ierr,ncerr,unt,iv0,nu,ii 337 character(len=500) :: msg 338 !arrays 339 real(dp),allocatable :: gvals_qpath(:,:),wvols_qpath(:,:,:),dwdq_qpath(:,:,:) 340 real(dp),allocatable :: phdispl_cart_qpath(:,:,:,:,:) 341 342 ! ************************************************************************ 343 344 nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 345 346 write(msg,'(a,(80a),4a)')ch10,('=',ii=1,80),ch10,ch10,' Calculation of Gruneisen parameters along q-path ',ch10 347 call wrtout(std_out, msg) 348 !call wrtout(ab_out, msg) 349 350 iv0 = gruns%iv0 351 ABI_CALLOC(wvols_qpath, (gruns%natom3, gruns%nvols, qpath%npts)) 352 ABI_CALLOC(gvals_qpath, (gruns%natom3, qpath%npts)) 353 ABI_CALLOC(dwdq_qpath, (3, gruns%natom3, qpath%npts)) 354 ABI_CALLOC(phdispl_cart_qpath, (2, gruns%natom3, gruns%natom3, gruns%nvols, qpath%npts)) 355 356 do iqpt=1,qpath%npts 357 if (mod(iqpt, nprocs) /= my_rank) cycle ! mpi-parallelism 358 call gruns_fourq(gruns, qpath%points(:,iqpt), wvols_qpath(:,:,iqpt), gvals_qpath(:,iqpt), & 359 dwdq_qpath(:,:,iqpt), phdispl_cart_qpath(:,:,:,:,iqpt)) 360 end do 361 362 call xmpi_sum(wvols_qpath, comm, ierr) 363 call xmpi_sum(gvals_qpath, comm, ierr) 364 call xmpi_sum(dwdq_qpath, comm, ierr) 365 call xmpi_sum(phdispl_cart_qpath, comm, ierr) 366 367 ! Write text files with phonon frequencies and gruneisen on the path. 368 if (my_rank == master) then 369 if (open_file(strcat(prefix, "_GRUNS_QPATH"), msg, newunit=unt, form="formatted", action="write") /= 0) then 370 ABI_ERROR(msg) 371 end if 372 write(unt,'(a)')'# Phonon band structure, Gruneisen parameters and group velocity' 373 write(unt,'(a)')"# Energy in Hartree, DOS in states/Hartree" 374 call qpath%print(unit=unt, pre="#") 375 write(unt,'(5a)')& 376 "# phfreq(mode=1) gruneisen(mode=1) velocity(mode=1) phfreq(mode=2) gruneisen(mode=2) velocity(mode=2) ..." 377 do iqpt=1,qpath%npts 378 do nu=1,gruns%natom3 379 write(unt, "(3es17.8)", advance="no") & 380 wvols_qpath(nu, iv0, iqpt), gvals_qpath(nu, iqpt), sum(dwdq_qpath(1:3, nu, iqpt) ** 2) 381 end do 382 write(unt, "(a)")" " 383 end do 384 close(unt) 385 end if 386 387 #ifdef HAVE_NETCDF 388 if (my_rank == master .and. ncid /= nctk_noid) then 389 ncerr = nctk_def_dims(ncid, [nctkdim_t("gruns_nqpath", qpath%npts)], defmode=.True.) 390 NCF_CHECK(ncerr) 391 392 ncerr = nctk_def_arrays(ncid, [ & 393 ! q-points of the path 394 nctkarr_t("gruns_qpath", "dp", "three, gruns_nqpath"), & 395 ! gruneisen parameters on the path 396 nctkarr_t("gruns_gvals_qpath", "dp", "number_of_phonon_modes, gruns_nqpath"), & 397 ! phonon frequencies at the different volumes 398 nctkarr_t("gruns_wvols_qpath", "dp", "number_of_phonon_modes, gruns_nvols, gruns_nqpath"), & 399 ! group velocities at V0 in Cartesian coordinates. 400 nctkarr_t("gruns_dwdq_qpath", "dp", "three, number_of_phonon_modes, gruns_nqpath"), & 401 ! displacements for the different volumes. 402 nctkarr_t("gruns_phdispl_cart_qpath", "dp", & 403 "two, number_of_phonon_modes, number_of_phonon_modes, gruns_nvols, gruns_nqpath") & 404 ]) 405 NCF_CHECK(ncerr) 406 407 ! Write data. 408 NCF_CHECK(nctk_set_datamode(ncid)) 409 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_qpath"), qpath%points)) 410 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_gvals_qpath"), gvals_qpath)) 411 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_wvols_qpath"), wvols_qpath)) 412 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_dwdq_qpath"), dwdq_qpath)) 413 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_phdispl_cart_qpath"), phdispl_cart_qpath)) 414 end if 415 #endif 416 417 ABI_FREE(wvols_qpath) 418 ABI_FREE(gvals_qpath) 419 ABI_FREE(dwdq_qpath) 420 ABI_FREE(phdispl_cart_qpath) 421 422 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
64 type,public :: gruns_t 65 66 integer :: natom3 67 ! 3 * natom 68 69 integer :: nvols 70 ! Number of volumes. 71 72 integer :: iv0 73 ! Index of the DDB file corresponding to the equilibrium volume V0. 74 75 real(dp) :: v0 76 ! Equilibrium volume. 77 78 real(dp) :: delta_vol 79 ! Uniform grid spacing for finite difference. 80 81 type(crystal_t),allocatable :: cryst_vol(:) 82 ! cryst_vol(nvols) 83 ! crystalline structure for the different volumes. 84 85 type(ddb_type),allocatable :: ddb_vol(:) 86 ! dbb_vol(nvols) 87 ! DDB objects for the different volumes. 88 89 type(ifc_type),allocatable :: ifc_vol(:) 90 ! ifc_vol(nvols) 91 ! interatomic force constants for the different volumes. 92 93 end type gruns_t 94 95 public :: gruns_new ! Constructor. 96 public :: gruns_qpath ! Compute Gruneisen parameters on a q-path. 97 public :: gruns_qmesh ! Compute Gruneisen parameters on a q-mesh. 98 public :: gruns_free ! Release memory. 99 public :: gruns_anaddb ! Driver routine called in anaddb.