TABLE OF CONTENTS


ABINIT/m_ifc [ Modules ]

[ Top ] [ Modules ]

NAME

  m_ifc

FUNCTION

  This module contains the declaration of data types and methods
  used to handle interatomic force constant sets

COPYRIGHT

 Copyright (C) 2011-2018 ABINIT group (XG,MJV,EB,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_ifc
29 
30  use defs_basis
31  use m_errors
32  use m_abicore
33  use m_xmpi
34  use m_sort
35  use m_cgtools
36  use m_bspline
37  use m_skw
38  use m_lebedev
39  use m_nctk
40 #ifdef HAVE_NETCDF
41  use netcdf
42 #endif
43 
44  use m_io_tools,    only : open_file
45  use m_fstrings,    only : ktoa, int2char4, sjoin, itoa, ltoa, ftoa
46  use m_symtk,       only : matr3inv
47  use m_special_funcs,  only : abi_derfc
48  use m_time,        only : cwtime
49  use m_numeric_tools, only : wrap2_zero_one, wrap2_pmhalf
50  use m_copy,        only : alloc_copy
51  use m_pptools,     only : printbxsf
52  use m_ewald,       only : ewald9
53  use m_crystal,     only : crystal_t,crystal_free
54  use m_geometry,    only : phdispl_cart2red, normv, mkrdim
55  use m_kpts,        only : kpts_ibz_from_kptrlatt, listkk, smpbz
56  use m_dynmat,      only : canct9, dist9 , ifclo9, axial9, q0dy3_apply, q0dy3_calc, asrif9, dynmat_dq, &
57 &                          make_bigbox, canat9, chkrp9, ftifc_q2r, wght9, symdm9, nanal9, gtdyn9, dymfz9, &
58 &                          massmult_and_breaksym, dfpt_phfrq, dfpt_prtph
59  use m_ddb
60  use m_ddb_hdr
61  use m_symkpt
62 
63  implicit none
64 
65  private

m_ifc/corsifc9 [ Functions ]

[ Top ] [ m_ifc ] [ Functions ]

NAME

 corsifc9

FUNCTION

 Applies a cutoff on the ifc in real space

INPUTS

 acell(3)=length scales by which rprim is to be multiplied
 gprim(3,3)=dimensionless primitive translations in reciprocal space
 natom=number of atoms in unit cell
 nrpt= Number of R points in the Big Box
 rcan(3,natom)=canonical coordinates of atoms
 rprim(3,3)=dimensionless primitive translations in real space
 rpt(3,nrpt)=canonical coordinates of the points in the BigBox.
 nsphere=number of atoms to be included in the cut-off sphere for interatomic
  force constant; if = 0 : maximum extent allowed by the grid.
 rifcsph=radius for cutoff of IFC
 wghatm(natom,natom,nrpt) = Weights associated to a pair of atoms and to a R vector.

OUTPUT

 wghatm(natom,natom,nrpt) = Weights associated to a pair of atoms and to a R vector
  with the required cutoff applied.
 rcut_min=Effective cutoff. Defined by the minimum cutoff radius over the natom sites.

PARENTS

      m_ifc

CHILDREN

      dfpt_phfrq,gtdyn9,nctk_defwrite_nonana_terms

SOURCE

1603 subroutine corsifc9(acell,gprim,natom,nrpt,nsphere,rifcsph,rcan,rprim,rpt,rcut_min,wghatm)
1604 
1605 
1606 !This section has been created automatically by the script Abilint (TD).
1607 !Do not modify the following lines by hand.
1608 #undef ABI_FUNC
1609 #define ABI_FUNC 'corsifc9'
1610 !End of the abilint section
1611 
1612  implicit none
1613 
1614 !Arguments -------------------------------
1615 !scalars
1616  integer,intent(in) :: natom,nrpt,nsphere
1617  real(dp),intent(in) :: rifcsph
1618  real(dp),intent(out) :: rcut_min
1619 !arrays
1620  real(dp),intent(in) :: acell(3)
1621  real(dp),intent(in) :: gprim(3,3),rcan(3,natom)
1622  real(dp),intent(in) :: rprim(3,3),rpt(3,nrpt)
1623  real(dp),intent(inout) :: wghatm(natom,natom,nrpt)
1624 
1625 !Local variables -------------------------
1626 !scalars
1627  integer :: ia,ib,ii,index,irpt
1628  real(dp) :: rmax,rsigma,r0
1629 !arrays
1630  integer,allocatable :: list(:)
1631  real(dp),allocatable :: dist(:,:,:),wkdist(:)
1632 
1633 ! *********************************************************************
1634 
1635  ! Compute the distances between atoms
1636  ! dist(ia,ib,irpt) contains the distance from atom ia to atom ib in unit cell irpt.
1637  ABI_MALLOC(dist,(natom,natom,nrpt))
1638  call dist9(acell,dist,gprim,natom,nrpt,rcan,rprim,rpt)
1639 
1640  ABI_MALLOC(list,(natom*nrpt))
1641  ABI_MALLOC(wkdist,(natom*nrpt))
1642 
1643  ! loop on all generic atoms.
1644  rcut_min = huge(one)
1645  do ia=1,natom
1646 
1647    wkdist = reshape(dist(ia,:,:), [natom*nrpt])
1648    do ii=1,natom*nrpt
1649      list(ii)=ii
1650    end do
1651    ! This sorting algorithm is slow ...
1652    call sort_dp(natom*nrpt,wkdist,list,tol14)
1653    rmax = wkdist(natom*nrpt)
1654 
1655    ! zero the outside IFCs: act on wghatm
1656 
1657    ! fix number of spheres
1658    if(nsphere/=0.and.nsphere<natom*nrpt)then
1659      rcut_min = min(rcut_min, wkdist(nsphere+1))
1660      do ii=nsphere+1,natom*nrpt
1661        index=list(ii)
1662        irpt=(index-1)/natom+1
1663        ib=index-natom*(irpt-1)
1664        wghatm(ia,ib,irpt)=zero
1665      end do
1666    end if
1667 
1668    ! or fix radius of maximum ifc
1669    if(rifcsph>tol10)then
1670      do ii=nsphere+1,natom*nrpt
1671        index=list(ii)
1672        ! preserve weights for atoms inside sphere of radius rifcsph
1673        if (wkdist(ii) < rifcsph) cycle
1674        rcut_min = min(rcut_min, wkdist(ii))
1675        irpt=(index-1)/natom+1
1676        ib=index-natom*(irpt-1)
1677        wghatm(ia,ib,irpt)=zero
1678      end do
1679    end if
1680 
1681    ! filter smoothly to 0 at edge of WS cell
1682    if (rifcsph < -tol10) then
1683      ! Use different filter
1684      r0 = abs(rifcsph) * rmax; rsigma = half*(rmax-r0) !one
1685      rcut_min = r0 ! Set it to r0
1686      do ii=nsphere+1,natom*nrpt
1687        index=list(ii)
1688        irpt=(index-1)/natom+1
1689        ib=index-natom*(irpt-1)
1690        wghatm(ia,ib,irpt) = wghatm(ia,ib,irpt) * half * abi_derfc((wkdist(ii) - r0) / rsigma)
1691      end do
1692    end if
1693 
1694  end do
1695 
1696  ABI_FREE(dist)
1697  ABI_FREE(list)
1698  ABI_FREE(wkdist)
1699 
1700 end subroutine corsifc9

m_ifc/ifc_autocutoff [ Functions ]

[ Top ] [ m_ifc ] [ Functions ]

NAME

  ifc_autocutoff

FUNCTION

 Find the value of nsphere that gives non-negative frequencies around Gamma
 in a small sphere of radius qrad.
 Use bisection to reduce the number of attemps although there's no guarantee
 that the number of negative frequencies is monotonic.

INPUTS

  crystal<crystal_t> = Information on the crystalline structure.
  comm=MPI communicator

SIDE EFFECTS

  ifc%wghatm(natom,natom,nrpt) = Weights associated to a pair of atoms and to a R vector
    with the last cutoff found by the bisection algorithm applied.
  ifc%atmfrc(2,3,natom,3,natom,nrpt)= ASR-imposed Interatomic Forces

PARENTS

      m_ifc

CHILDREN

      dfpt_phfrq,gtdyn9,nctk_defwrite_nonana_terms

SOURCE

1441 subroutine ifc_autocutoff(ifc, crystal, comm)
1442 
1443 
1444 !This section has been created automatically by the script Abilint (TD).
1445 !Do not modify the following lines by hand.
1446 #undef ABI_FUNC
1447 #define ABI_FUNC 'ifc_autocutoff'
1448 !End of the abilint section
1449 
1450  implicit none
1451 
1452 !Arguments ------------------------------------
1453 !scalars
1454  type(ifc_type),intent(inout) :: ifc
1455  type(crystal_t),intent(in) :: crystal
1456  integer,intent(in) :: comm
1457 
1458 !Local variables-------------------------------
1459 !scalars
1460  integer,parameter :: master=0
1461  integer :: iq_ibz,ierr,my_rank,nprocs,ii,nsphere,num_negw,jl,ju,jm,natom,nrpt
1462  real(dp),parameter :: rifcsph0=zero
1463  real(dp) :: adiff,qrad,min_negw,xval,rcut_min
1464  type(lebedev_t) :: lgrid
1465 !arrays
1466  real(dp) :: displ_cart(2*3*ifc%natom*3*ifc%natom)
1467  real(dp) :: qred(3),qred_vers(3),phfrqs(3*ifc%natom) !,dwdq(3,3*ifc%natom)
1468  real(dp),allocatable :: ref_phfrq(:,:),cut_phfrq(:,:)
1469  real(dp),allocatable :: save_wghatm(:,:,:),save_atmfrc(:,:,:,:,:)
1470 
1471 ! *********************************************************************
1472 
1473  my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
1474  natom = ifc%natom; nrpt = ifc%nrpt
1475 
1476  ! Compute frequencies on the ab-initio q-mesh without cutoff.
1477  ABI_CALLOC(ref_phfrq, (3*natom, ifc%nqibz))
1478  do iq_ibz=1,ifc%nqibz
1479    if (mod(iq_ibz, nprocs) /= my_rank) cycle ! mpi-parallelism
1480    call ifc_fourq(ifc, crystal, ifc%qibz(:,iq_ibz), ref_phfrq(:,iq_ibz), displ_cart)
1481  end do
1482  call xmpi_sum(ref_phfrq, comm, ierr)
1483 
1484  ABI_MALLOC(save_wghatm, (natom,natom,nrpt))
1485  ABI_MALLOC(save_atmfrc, (3,natom,3,natom,ifc%nrpt))
1486  save_wghatm = ifc%wghatm; save_atmfrc = ifc%atmfrc
1487 
1488  ABI_MALLOC(cut_phfrq, (3*natom, ifc%nqibz))
1489  qrad = 0.01; lgrid = lebedev_new(16)
1490 
1491  if (my_rank == master) then
1492    write(ab_out, "(a)")" Apply cutoff on IFCs. Using bisection algorithm to find initial guess for nsphere."
1493    write(ab_out, "(a,i0)")" Maximum nuber of atom-centered spheres: ",natom * nrpt
1494    write(ab_out, "(a,i0,a,f5.3)")" Using Lebedev-Laikov grid with npts: ",lgrid%npts, ", qrad: ",qrad
1495    write(ab_out, "(/,a)")" <adiff>: Average difference between ab-initio frequencies and frequencies with cutoff."
1496    write(ab_out, "(a)")" num_negw: Number of negative freqs detected in small sphere around Gamma."
1497    write(ab_out, "(a)")" min_negw: Min negative frequency on the small sphere."
1498    write(ab_out, "(a,/,/)")" rifcsph: Effective cutoff radius corresponding to nsphere."
1499    write(ab_out, "(a)")" nsphere   <adiff>[meV]   num_negw   min_negw[meV]   rifcsph"
1500  end if
1501 
1502  jl = 0; ju = natom * nrpt + 1 ! Initialize lower and upper limits.
1503  do
1504    if (ju - jl <= 1) then
1505      exit
1506    end if
1507    jm = (ju + jl) / 2  ! Compute a midpoint
1508    nsphere = jm
1509 
1510    ifc%wghatm = save_wghatm; ifc%atmfrc = save_atmfrc
1511    call corsifc9(ifc%acell,ifc%gprim, natom, nrpt,nsphere,rifcsph0,ifc%rcan,ifc%rprim,ifc%rpt,rcut_min,ifc%wghatm)
1512    if (ifc%asr > 0) call asrif9(ifc%asr,ifc%atmfrc,ifc%natom,ifc%nrpt,ifc%rpt,ifc%wghatm)
1513 
1514    cut_phfrq = zero
1515    do iq_ibz=1,ifc%nqibz
1516      if (mod(iq_ibz, nprocs) /= my_rank) cycle ! mpi-parallelism
1517      call ifc_fourq(ifc, crystal, ifc%qibz(:,iq_ibz), cut_phfrq(:,iq_ibz), displ_cart)
1518      !write(std_out,*)cut_phfrq(1,iq_ibz),ref_phfrq(1,iq_ibz)
1519    end do
1520    call xmpi_sum(cut_phfrq, comm, ierr)
1521 
1522    ! Test wether there are negative frequencies around gamma, including reciprocal lattice vectors.
1523    num_negw = 0; min_negw = zero
1524    do ii=1,lgrid%npts+3
1525      if (mod(ii, nprocs) /= my_rank) cycle ! mpi-parallelism
1526      if (ii <= 3) then
1527        qred = zero; qred(ii) = one
1528      else
1529        qred = matmul(crystal%rprimd, lgrid%versors(:, ii-3))
1530      end if
1531      qred_vers = (qred / normv(qred, crystal%gmet, "G"))
1532      qred = qrad * qred_vers
1533      call ifc_fourq(ifc, crystal, qred, phfrqs, displ_cart) !, dwdq=dwdq)
1534      if (any(phfrqs < +tol8)) then
1535        num_negw = num_negw + 1; min_negw = min(min_negw, minval(phfrqs))
1536      end if
1537      !do jj=1,3
1538      !  xval = dot_product(dwdq(:,jj), matmul(crystal%gprimd, qred_vers))
1539      !  if (xval < zero) num_negw = num_negw + 1
1540      !end do
1541    end do
1542    call xmpi_sum(num_negw, comm, ierr)
1543    xval = min_negw; call xmpi_min(xval, min_negw, comm, ierr)
1544 
1545    adiff = sum(abs(cut_phfrq - ref_phfrq)) / (ifc%nqibz * 3 * natom)
1546    if (my_rank == master) then
1547      write(ab_out,"(a,i7,1x,es13.4,4x,i8,1x,es13.4,2x,es13.4)") &
1548        "-",nsphere, adiff * Ha_meV, num_negw, min_negw * Ha_meV, rcut_min
1549    end if
1550 
1551    if (num_negw == 0) then
1552      jl = jm ! Replace lower limit
1553    else
1554      ju = jm ! Replace upper limit
1555    end if
1556  end do
1557 
1558  ABI_FREE(ref_phfrq)
1559  ABI_FREE(cut_phfrq)
1560  ABI_FREE(save_wghatm)
1561  ABI_FREE(save_atmfrc)
1562  call lebedev_free(lgrid)
1563 
1564 end subroutine ifc_autocutoff

m_ifc/ifc_build_phbspl [ Functions ]

[ Top ] [ m_ifc ] [ Functions ]

NAME

 ifc_build_phbspl

FUNCTION

 build the phbspl_t object to interpolate the phonon band structure.

INPUTS

  ifc<type(ifc_type)>=Object containing the dynamical matrix and the IFCs.
  crystal<type(crystal_t)> = Information on the crystalline structure.
  ngqpt(3)=Divisions of the q-mesh used to produce the B-spline.
  nshiftq=Number of shifts in Q-mesh
  shiftq(3,nshiftq)=Shifts of the q-mesh.
  ords(3)=order of the spline for the three directions. ord(1) must be in [0, nqx] where
    nqx is the number of points along the x-axis.
  comm=MPI communicator

OUTPUT

PARENTS

      m_ifc

CHILDREN

SOURCE

2849 type(phbspl_t) function ifc_build_phbspl(ifc, cryst, ngqpt, nshiftq, shiftq, ords, comm) result(new)
2850 
2851 
2852 !This section has been created automatically by the script Abilint (TD).
2853 !Do not modify the following lines by hand.
2854 #undef ABI_FUNC
2855 #define ABI_FUNC 'ifc_build_phbspl'
2856 !End of the abilint section
2857 
2858  implicit none
2859 
2860 !arguments ------------------------------------
2861 !scalars
2862  integer,intent(in) :: comm,nshiftq
2863  type(ifc_type),intent(in) :: ifc
2864  type(crystal_t),intent(in) :: cryst
2865 !arrays
2866  integer,intent(in) :: ngqpt(3),ords(3)
2867  real(dp),intent(in) :: shiftq(3,nshiftq)
2868 
2869 !local variables-------------------------------
2870 !scalars
2871  integer,parameter :: sppoldbl1=1,timrev1=1,qptopt1=1
2872  integer :: qxord,qyord,qzord,nxknot,nyknot,nzknot,nqibz,nqbz,nqfull,iqf,ierr
2873  integer :: nu,iq_ibz,ix,iy,iz,nqx,nqy,nqz,my_rank,nprocs,cnt
2874  real(dp) :: dksqmax
2875  character(len=500) :: msg
2876 !arrays
2877  integer :: qptrlatt(3,3)
2878  integer,allocatable :: bz2ibz(:,:)
2879  real(dp) :: displ_cart(2,3,cryst%natom,3*cryst%natom)
2880  real(dp),allocatable :: xvec(:),yvec(:),zvec(:),xyzdata(:,:,:)
2881  real(dp),allocatable :: ibz_freqs(:,:),ibzdata_qnu(:,:)
2882  real(dp),allocatable :: wtq(:),qbz(:,:),qfull(:,:),qibz(:,:)
2883 
2884 ! *********************************************************************
2885 
2886  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2887 
2888  ! Check input parameters
2889  ierr = 0
2890  if (nshiftq /= 1) then
2891    MSG_WARNING('Multiple shifts not allowed')
2892    ierr = ierr + 1
2893  end if
2894  if (ierr /= 0) then
2895    MSG_WARNING("bspline interpolation cannot be performed. See warnings above. Returning")
2896    return
2897  end if
2898 
2899  ! Setup IBZ, weights and BZ. Always use q --> -q symmetry for phonons even in systems wo inversion
2900  qptrlatt = 0; qptrlatt(1,1) = ngqpt(1); qptrlatt(2,2) = ngqpt(2); qptrlatt(3,3) = ngqpt(3)
2901  call kpts_ibz_from_kptrlatt(cryst, qptrlatt, qptopt1, nshiftq, shiftq, nqibz, qibz, wtq, nqbz, qbz)
2902 
2903  ABI_FREE(wtq)
2904  ABI_FREE(qbz)
2905 
2906  ! Build BZ mesh Note that:
2907  ! 1) q-point coordinates are in [0, 1]
2908  ! 2) The mesh is closed e.g. (0,0,0) and (1,1,1) are included
2909  nqx = ngqpt(1)+1; nqy = ngqpt(2)+1; nqz = ngqpt(3)+1
2910 
2911  ABI_MALLOC(xvec, (nqx))
2912  ABI_MALLOC(yvec, (nqy))
2913  ABI_MALLOC(zvec, (nqz))
2914 
2915  ! Multiple shifts are not supported here.
2916  do ix=1,nqx
2917    xvec(ix) = (ix-one+shiftq(1,1)) / ngqpt(1)
2918  end do
2919  do iy=1,nqy
2920    yvec(iy) = (iy-one+shiftq(2,1)) / ngqpt(2)
2921  end do
2922  do iz=1,nqz
2923    zvec(iz) = (iz-one+shiftq(3,1)) / ngqpt(3)
2924  end do
2925 
2926  ! Build list of q-points in full BZ (ordered as required by B-spline routines)
2927  nqfull = nqx*nqy*nqz
2928  ABI_MALLOC(qfull, (3,nqfull))
2929  iqf = 0
2930  do iz=1,nqz
2931    do iy=1,nqy
2932      do ix=1,nqx
2933        iqf = iqf + 1
2934        qfull(:,iqf) = [xvec(ix), yvec(iy), zvec(iz)]
2935      end do
2936    end do
2937  end do
2938 
2939  ! Build mapping qfull --> IBZ (q --> -q symmetry is always used)
2940  ABI_MALLOC(bz2ibz, (nqfull*sppoldbl1,6))
2941 
2942  call listkk(dksqmax,cryst%gmet,bz2ibz,qibz,qfull,nqibz,nqfull,cryst%nsym,&
2943    sppoldbl1,cryst%symafm,cryst%symrec,timrev1,use_symrec=.True.)
2944  ABI_FREE(qfull)
2945 
2946  if (dksqmax > tol12) then
2947    write(msg, '(3a,es16.6,4a)' )&
2948    'At least one of the q-points could not be generated from a symmetrical one.',ch10,&
2949    'dksqmax=',dksqmax,ch10,&
2950    'Action: check q-point input variables',ch10,&
2951    '        (e.g. kptopt or shiftk might be wrong in the present dataset or the preparatory one.'
2952    MSG_ERROR(msg)
2953  end if
2954 
2955  ! Generate knots (ords is input)
2956  qxord = ords(1); qyord = ords(2); qzord = ords(3)
2957  nxknot = nqx + qxord
2958  nyknot = nqy + qyord
2959  nzknot = nqz + qzord
2960 
2961  new%nqx = nqx; new%qxord = qxord
2962  new%nqy = nqy; new%qyord = qyord
2963  new%nqz = nqz; new%qzord = qzord
2964 
2965  ABI_MALLOC(new%xknot,(nxknot))
2966  ABI_MALLOC(new%yknot,(nyknot))
2967  ABI_MALLOC(new%zknot,(nzknot))
2968 
2969  call dbsnak(nqx, xvec, qxord, new%xknot)
2970  call dbsnak(nqy, yvec, qyord, new%yknot)
2971  call dbsnak(nqz, zvec, qzord, new%zknot)
2972 
2973  new%natom3 = 3 * cryst%natom
2974 
2975  ! Get phonon frequencies in IBZ
2976  ABI_CALLOC(ibz_freqs, (new%natom3, nqibz))
2977  cnt = 0
2978  do iq_ibz=1,nqibz
2979    cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! Mpi parallelism.
2980    call ifc_fourq(ifc, cryst, qibz(:,iq_ibz), ibz_freqs(:,iq_ibz), displ_cart)
2981  end do
2982  call xmpi_sum(ibz_freqs, comm, ierr)
2983 
2984  ABI_MALLOC(ibzdata_qnu, (nqibz, new%natom3))
2985  ibzdata_qnu = transpose(ibz_freqs)
2986  ABI_FREE(ibz_freqs)
2987 
2988  ABI_MALLOC(xyzdata,(nqx,nqy,nqz))
2989  ABI_DT_MALLOC(new%coeff, (new%natom3))
2990 
2991  do nu=1,new%natom3
2992 
2993    ABI_MALLOC(new%coeff(nu)%vals, (nqx,nqy,nqz))
2994 
2995    ! Build array in full bz to prepare call to dbs3in.
2996    iqf = 0
2997    do iz=1,nqz
2998      do iy=1,nqy
2999        do ix=1,nqx
3000          iqf = iqf + 1
3001          iq_ibz = bz2ibz(iqf,1)
3002          xyzdata(ix,iy,iz) = ibzdata_qnu(iq_ibz, nu)
3003        end do
3004      end do
3005    end do
3006 
3007    ! Construct 3D tensor for B-spline. Results in coeff(nu)%vals
3008    call dbs3in(nqx,xvec,nqy,yvec,nqz,zvec,xyzdata,nqx,nqy,qxord,qyord,qzord,new%xknot,new%yknot,new%zknot,&
3009       new%coeff(nu)%vals)
3010  end do
3011 
3012  ABI_FREE(xvec)
3013  ABI_FREE(yvec)
3014  ABI_FREE(zvec)
3015  ABI_FREE(bz2ibz)
3016  ABI_FREE(xyzdata)
3017  ABI_FREE(ibzdata_qnu)
3018  ABI_FREE(qibz)
3019 
3020 end function ifc_build_phbspl

m_ifc/ifc_build_skw [ Functions ]

[ Top ] [ m_ifc ] [ Functions ]

NAME

 ifc_build_skw

FUNCTION

INPUTS

  ifc<type(ifc_type)>=Object containing the dynamical matrix and the IFCs.
  crystal<type(crystal_t)> = Information on the crystalline structure.
  ngqpt(3)=Divisions of the q-mesh used to produce the B-spline.
  nshiftq=Number of shifts in Q-mesh
  shiftq(3,nshiftq)=Shifts of the q-mesh.
  ords(3)=order of the spline for the three directions. ord(1) must be in [0, nqx] where
    nqx is the number of points along the x-axis.
  comm=MPI communicator

OUTPUT

PARENTS

CHILDREN

SOURCE

3196 type(skw_t) function ifc_build_skw(ifc, cryst, ngqpt, nshiftq, shiftq, comm) result(new)
3197 
3198 
3199 !This section has been created automatically by the script Abilint (TD).
3200 !Do not modify the following lines by hand.
3201 #undef ABI_FUNC
3202 #define ABI_FUNC 'ifc_build_skw'
3203 !End of the abilint section
3204 
3205  implicit none
3206 
3207 !arguments ------------------------------------
3208 !scalars
3209  integer,intent(in) :: comm,nshiftq
3210  type(ifc_type),intent(in) :: ifc
3211  type(crystal_t),intent(in) :: cryst
3212 !arrays
3213  integer,intent(in) :: ngqpt(3)
3214  real(dp),intent(in) :: shiftq(3,nshiftq)
3215 
3216 !local variables-------------------------------
3217 !scalars
3218  integer,parameter :: sppoldbl1=1,timrev1=1,master=0,qptopt1=1
3219  integer :: nqibz,nqbz,iq_ibz,iq_bz,natom3,ierr,nu
3220  integer :: my_rank,nprocs,cnt
3221  real(dp) :: dksqmax
3222  character(len=500) :: msg
3223 !arrays
3224  integer :: qptrlatt(3,3)
3225  integer,allocatable :: bz2ibz(:,:)
3226  real(dp) :: params(4)
3227  real(dp) :: phfrq(3*cryst%natom),displ_cart(2,3,cryst%natom,3*cryst%natom)
3228  real(dp),allocatable :: ibz_freqs(:,:),wtq(:),qbz(:,:),qibz(:,:)
3229 
3230 ! *********************************************************************
3231 
3232  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
3233 
3234  natom3 = 3 * cryst%natom
3235 
3236  ! Setup IBZ, weights and BZ. Always use q --> -q symmetry for phonons even in systems wo inversion
3237  qptrlatt = 0; qptrlatt(1,1) = ngqpt(1); qptrlatt(2,2) = ngqpt(2); qptrlatt(3,3) = ngqpt(3)
3238  call kpts_ibz_from_kptrlatt(cryst, qptrlatt, qptopt1, nshiftq, shiftq, nqibz, qibz, wtq, nqbz, qbz)
3239 
3240  ! Get phonon frequencies in IBZ
3241  ABI_CALLOC(ibz_freqs, (natom3, nqibz))
3242  cnt = 0
3243  do iq_ibz=1,nqibz
3244    cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! Mpi parallelism.
3245    call ifc_fourq(ifc, cryst, qibz(:,iq_ibz), ibz_freqs(:,iq_ibz), displ_cart)
3246  end do
3247  call xmpi_sum(ibz_freqs, comm, ierr)
3248 
3249  params = zero; params(1) = 120
3250  new = skw_new(cryst, params, 1, natom3, nqibz, 1, qibz, ibz_freqs, [0,0], comm)
3251 
3252  if (.False. .and. my_rank == master) then
3253    ! Test whether SKW preserves symmetries.
3254    ! Build mapping qbz --> IBZ (q --> -q symmetry is always used)
3255    ABI_MALLOC(bz2ibz, (nqbz*sppoldbl1,6))
3256 
3257    call listkk(dksqmax,cryst%gmet,bz2ibz,qibz,qbz,nqibz,nqbz,cryst%nsym,&
3258      sppoldbl1,cryst%symafm,cryst%symrec,timrev1,use_symrec=.True.)
3259 
3260    if (dksqmax > tol12) then
3261      write(msg, '(3a,es16.6,4a)' )&
3262      'At least one of the q-points could not be generated from a symmetrical one.',ch10,&
3263      'dksqmax=',dksqmax,ch10,&
3264      'Action: check q-point input variables',ch10,&
3265      '        (e.g. kptopt or shiftk might be wrong in the present dataset or the preparatory one.'
3266      MSG_ERROR(msg)
3267    end if
3268 
3269    do iq_bz=1,nqbz
3270      iq_ibz = bz2ibz(iq_bz,1)
3271      do nu=1,natom3
3272        call skw_eval_bks(new, nu, qbz(:,iq_bz), 1, phfrq(nu))
3273      end do
3274      write(std_out,*)"BZ-IBZ:", maxval(abs(phfrq - ibz_freqs(:, iq_ibz)))
3275    end do
3276    ABI_FREE(bz2ibz)
3277  end if
3278 
3279  ABI_FREE(qibz)
3280  ABI_FREE(wtq)
3281  ABI_FREE(qbz)
3282  ABI_FREE(ibz_freqs)
3283 
3284 end function ifc_build_skw

m_ifc/ifc_calcnwrite_nana_terms [ Functions ]

[ Top ] [ m_ifc ] [ Functions ]

NAME

  ifc_calcnwrite_nana_terms

FUNCTION

  Compute frequencies and phonon displacement for q-->0 in the presence of non-analytical behaviour.

INPUTS

  nph2l=Number of qpoints.
  qph2l(3,nph2l)=List of phonon wavevector directions along which the non-analytical correction
    to the Gamma-point phonon frequencies will be calculated
    The direction is in CARTESIAN COORDINATES
  qnrml2(nph2l)=Normalization factor.

OUTPUT

  (Optional)
  phfrq2l(3*crystal%natom,nph2l)=List of phonon frequencies
  polarity2l(3,3*crystal%natom,nph2l)=List of mode-polarities
     (see Eq.(41) of Veithen et al, PRB71, 125107 (2005) [[cite:Veithen2005]])

 NOTES:
  This routine should be called by master node and when ifcflag == 1.

PARENTS

      m_gruneisen,m_phonons

CHILDREN

      dfpt_phfrq,gtdyn9,nctk_defwrite_nonana_terms

SOURCE

3529 subroutine ifc_calcnwrite_nana_terms(ifc, crystal, nph2l, qph2l, &
3530 &  qnrml2, ncid, phfrq2l, polarity2l) ! optional arguments
3531 
3532 
3533 !This section has been created automatically by the script Abilint (TD).
3534 !Do not modify the following lines by hand.
3535 #undef ABI_FUNC
3536 #define ABI_FUNC 'ifc_calcnwrite_nana_terms'
3537 !End of the abilint section
3538 
3539  implicit none
3540 
3541 !Arguments ------------------------------------
3542  integer,intent(in) :: nph2l
3543  integer,optional,intent(in) :: ncid
3544  type(ifc_type),intent(in) :: ifc
3545  type(crystal_t),intent(in) :: crystal
3546 !arrays
3547  real(dp),intent(in) :: qph2l(3, nph2l)
3548  real(dp),optional,intent(in) :: qnrml2(nph2l)
3549  real(dp),optional,intent(out) :: phfrq2l(3*crystal%natom,nph2l), polarity2l(3,3*crystal%natom,nph2l)
3550 
3551 !Local variables-------------------------------
3552 !scalars
3553  integer :: iatom,idir,imode,iphl2
3554  !character(len=500) :: msg
3555 !arrays
3556  real(dp) :: qphnrm(3),qphon(3,3)
3557  real(dp),allocatable :: displ_cart(:,:,:),phfrq(:),d2cart(:,:,:),eigvec(:,:,:),eigval(:)
3558 
3559 ! ************************************************************************
3560 
3561  if (nph2l == 0) return
3562 
3563  !Now treat the second list of vectors (only at the Gamma point, but can include non-analyticities)
3564  ABI_MALLOC(phfrq, (3*crystal%natom))
3565  ABI_MALLOC(displ_cart, (2, 3*crystal%natom, 3*crystal%natom))
3566  ABI_MALLOC(d2cart, (2, 3*ifc%mpert, 3*ifc%mpert))
3567  ABI_MALLOC(eigvec, (2, 3*crystal%natom, 3*crystal%natom))
3568  ABI_MALLOC(eigval, (3*crystal%natom))
3569 
3570  ! Before examining every direction or the dielectric tensor, generates the dynamical matrix at gamma
3571  qphon(:,1)=zero; qphnrm(1)=zero
3572 
3573  ! Generation of the dynamical matrix in cartesian coordinates
3574  ! Get d2cart using the interatomic forces and the
3575  ! long-range coulomb interaction through Ewald summation
3576  call gtdyn9(ifc%acell,ifc%atmfrc,ifc%dielt,ifc%dipdip, &
3577    ifc%dyewq0,d2cart,crystal%gmet,ifc%gprim,ifc%mpert,crystal%natom, &
3578    ifc%nrpt,qphnrm(1),qphon,crystal%rmet,ifc%rprim,ifc%rpt, &
3579    ifc%trans,crystal%ucvol,ifc%wghatm,crystal%xred,ifc%zeff)
3580 
3581 #ifdef HAVE_NETCDF
3582  if(present(ncid))then
3583    iphl2 = 0
3584    call nctk_defwrite_nonana_terms(ncid, iphl2, nph2l, qph2l, crystal%natom, phfrq, displ_cart, "define")
3585  endif
3586 #endif
3587 
3588  ! Examine every wavevector of this list
3589  do iphl2=1,nph2l
3590    ! Initialisation of the phonon wavevector
3591    qphon(:,1) = qph2l(:,iphl2)
3592    qphnrm(1)=zero
3593    if(present(qnrml2))then
3594      qphnrm(1) = qnrml2(iphl2)
3595    endif
3596 
3597 
3598    ! Calculation of the eigenvectors and eigenvalues of the dynamical matrix
3599    call dfpt_phfrq(ifc%amu,displ_cart,d2cart,eigval,eigvec,crystal%indsym, &
3600       ifc%mpert,crystal%nsym,crystal%natom,crystal%nsym,crystal%ntypat,phfrq,qphnrm(1),qphon, &
3601       crystal%rprimd,ifc%symdynmat,crystal%symrel,crystal%symafm,crystal%typat,crystal%ucvol)
3602 
3603    ! Write the phonon frequencies
3604    !call dfpt_prtph(displ_cart,inp%eivec,inp%enunit,ab_out,natom,phfrq,qphnrm(1),qphon)
3605 
3606    if(present(phfrq2l))then
3607      phfrq2l(:,iphl2)=phfrq(:)
3608    endif
3609 
3610    if(present(polarity2l))then
3611      polarity2l(:,:,iphl2)=zero
3612      do imode=1,3*crystal%natom
3613        do iatom=1,crystal%natom
3614          do idir=1,3
3615            polarity2l(:,imode,iphl2)=polarity2l(:,imode,iphl2)+ifc%zeff(:,idir,iatom)*displ_cart(1,idir+(iatom-1)*3,imode)
3616          enddo
3617        enddo
3618      enddo
3619    endif
3620 
3621 #ifdef HAVE_NETCDF
3622    if(present(ncid))then
3623      ! Loop is not MPI-parallelized --> no need for MPI-IO API.
3624      call nctk_defwrite_nonana_terms(ncid, iphl2, nph2l, qph2l, crystal%natom, phfrq, displ_cart, "write")
3625    endif
3626 #endif
3627  end do ! iphl2
3628 
3629  ABI_FREE(phfrq)
3630  ABI_FREE(displ_cart)
3631  ABI_FREE(d2cart)
3632  ABI_FREE(eigvec)
3633  ABI_FREE(eigval)
3634 
3635 end subroutine ifc_calcnwrite_nana_terms

m_ifc/ifc_fourq [ Functions ]

[ Top ] [ m_ifc ] [ Functions ]

NAME

  ifc_fourq

FUNCTION

  Compute the phonon frequencies and the group velocities at the specified q-point by performing
  a Fourier transform on the IFCs matrix in real space.

INPUTS

  Ifc<type(ifc_type)>=Object containing the dynamical matrix and the IFCs.
  Crystal<type(crystal_t)> = Information on the crystalline structure.
  qpt(3)=q-point in reduced coordinates (unless nanaqdir is specified)
  [nanaqdir]=If present, the qpt will be treated as a vector specifying the
    direction in q-space along which the non-analytic behaviour of the dynamical
    matrix will be treated. Possible values:
       "cart" if qpt defines a direction in Cartesian coordinates
       "reduced" if qpt defines a direction in reduced coordinates

OUTPUT

  phfrq(3*natom) = Phonon frequencies in Hartree
  displ_cart(2,3,natom,3*natom) = Phonon displacement in Cartesian coordinates
  [out_d2cart(2,3,3*natom,3,3*natom)] = The (interpolated) dynamical matrix for this q-point
  [out_eigvec(2*3*natom*3*natom) = The (interpolated) eigenvectors of the dynamical matrix in Cartesian coords..
  [out_displ_red(2*3*natom*3*natom) = The (interpolated) displacement in reduced coordinates.
  [dwdq(3,3*natom)] = Group velocities i.e. d(omega(q))/dq in Cartesian coordinates.

PARENTS

      get_nv_fs_en,get_tau_k,harmonic_thermo,interpolate_gkk,m_gruneisen
      m_ifc,m_phgamma,m_phonons,m_phpi,m_sigmaph,m_tdep_phdos,mka2f,mka2f_tr
      mka2f_tr_lova,mkph_linwid,read_gkk

CHILDREN

      dfpt_phfrq,gtdyn9,nctk_defwrite_nonana_terms

SOURCE

 966 subroutine ifc_fourq(ifc, crystal, qpt, phfrq, displ_cart, &
 967                      nanaqdir, out_d2cart, out_eigvec, out_displ_red, dwdq)   ! Optional [out]
 968 
 969 
 970 !This section has been created automatically by the script Abilint (TD).
 971 !Do not modify the following lines by hand.
 972 #undef ABI_FUNC
 973 #define ABI_FUNC 'ifc_fourq'
 974 !End of the abilint section
 975 
 976  implicit none
 977 
 978 !Arguments ------------------------------------
 979 !scalars
 980  character(len=*),optional,intent(in) :: nanaqdir
 981  type(ifc_type),intent(in) :: Ifc
 982  type(crystal_t),intent(in) :: Crystal
 983 !arrays
 984  real(dp),intent(in) :: qpt(3)
 985  real(dp),intent(out) :: displ_cart(2,3,Crystal%natom,3*Crystal%natom)
 986  real(dp),intent(out) :: phfrq(3*Crystal%natom)
 987  real(dp),optional,intent(out) :: out_d2cart(2,3,Crystal%natom,3,Crystal%natom)
 988  real(dp),optional,intent(out) :: out_eigvec(2,3,Crystal%natom,3*Crystal%natom)
 989  real(dp),optional,intent(out) :: out_displ_red(2,3,Crystal%natom,3*Crystal%natom)
 990  real(dp),optional,intent(out) :: dwdq(3,3*crystal%natom)
 991 
 992 !Local variables-------------------------------
 993 !scalars
 994  integer :: natom
 995  real(dp) :: qphnrm
 996 !arrays
 997  real(dp) :: my_qpt(3),eigvec(2,3,Crystal%natom,3*Crystal%natom),eigval(3*Crystal%natom)
 998  real(dp) :: d2cart(2,3,Ifc%mpert,3,Ifc%mpert)
 999 
1000 ! ************************************************************************
1001 
1002  natom = Crystal%natom
1003 
1004  ! Use my_qpt because dfpt_phfrq can change the q-point (very bad design)
1005  qphnrm = one; my_qpt = qpt
1006 
1007  if (present(nanaqdir)) then
1008    ! This will break backward compatibility because qpt is **always** in reduced coordinates.
1009    ! while dfpt_phfrq assume cartesian coordinates !!!!!!!!!!!
1010    ! It does not make sense to change API just to treat this particular case
1011    ! We should **alwayse use q-points in reduced coordinates.
1012    qphnrm = zero
1013    select case (nanaqdir)
1014    case ("reduced")
1015      ! Convert to Cartesian.
1016      my_qpt = matmul(Crystal%gprimd, qpt)
1017    case ("cart")
1018      continue
1019    case default
1020      MSG_ERROR(sjoin("Wrong value for nanaqdir:", nanaqdir))
1021    end select
1022  end if
1023 
1024  ! The dynamical matrix d2cart is calculated here:
1025  call gtdyn9(Ifc%acell,Ifc%atmfrc,Ifc%dielt,Ifc%dipdip,Ifc%dyewq0,d2cart,Crystal%gmet,Ifc%gprim,Ifc%mpert,natom,&
1026 &  Ifc%nrpt,qphnrm,my_qpt,Crystal%rmet,Ifc%rprim,Ifc%rpt,Ifc%trans,Crystal%ucvol,Ifc%wghatm,Crystal%xred,Ifc%zeff)
1027 
1028  ! Calculate the eigenvectors and eigenvalues of the dynamical matrix
1029  call dfpt_phfrq(Ifc%amu,displ_cart,d2cart,eigval,eigvec,Crystal%indsym,&
1030 &  Ifc%mpert,Crystal%nsym,natom,Crystal%nsym,Crystal%ntypat,phfrq,qphnrm,my_qpt,&
1031 &  Crystal%rprimd,Ifc%symdynmat,Crystal%symrel,Crystal%symafm,Crystal%typat,Crystal%ucvol)
1032 
1033  ! OmegaSRLR: Perform decomposition of dynamical matrix
1034  !if (srlr==1) call omega_decomp(amu,natom,ntypat,typat,dynmatfull,dynmatsr,dynmatlr,iqpt,nqpt,eigvec)
1035 
1036  ! Return the interpolated dynamical matrix and the eigenvector for this q-point
1037  if (present(out_d2cart)) out_d2cart = d2cart(:,:,:natom,:,:natom)
1038  if (present(out_eigvec)) out_eigvec = eigvec
1039 
1040  ! Return phonon displacement in reduced coordinates.
1041  if (present(out_displ_red)) call phdispl_cart2red(natom, crystal%gprimd, displ_cart, out_displ_red)
1042 
1043  ! Option to get vectors in reduced coordinates?
1044  !call phdispl_cart2red(natom, crystal%gprimd, out_eigvec, out_eigvec_red)
1045 
1046  ! Compute group velocities.
1047  if (present(dwdq)) call ifc_get_dwdq(ifc, crystal, my_qpt, phfrq, eigvec, dwdq)
1048 
1049 end subroutine ifc_fourq

m_ifc/ifc_free [ Functions ]

[ Top ] [ m_ifc ] [ Functions ]

NAME

 ifc_free

FUNCTION

  Deallocate memory for the ifc_type structure

PARENTS

      anaddb,compute_anharmonics,eph,m_anharmonics_terms
      m_effective_potential,m_effective_potential_file,m_gruneisen
      m_harmonics_terms,m_ifc

CHILDREN

      dfpt_phfrq,gtdyn9,nctk_defwrite_nonana_terms

SOURCE

272 subroutine ifc_free(ifc)
273 
274 
275 !This section has been created automatically by the script Abilint (TD).
276 !Do not modify the following lines by hand.
277 #undef ABI_FUNC
278 #define ABI_FUNC 'ifc_free'
279 !End of the abilint section
280 
281  implicit none
282 
283 !Arguments ------------------------------------
284 !array
285  type(ifc_type),intent(inout) :: ifc
286 
287 ! ************************************************************************
288  ! real
289  if (allocated(ifc%amu)) then
290    ABI_FREE(ifc%amu)
291  end if
292  if (allocated(ifc%atmfrc)) then
293    ABI_FREE(ifc%atmfrc)
294  end if
295  if (allocated(ifc%cell)) then
296    ABI_FREE(ifc%cell)
297  end if
298  if (allocated(ifc%ewald_atmfrc)) then
299    ABI_FREE(ifc%ewald_atmfrc)
300  end if
301  if (allocated(ifc%short_atmfrc)) then
302    ABI_FREE(ifc%short_atmfrc)
303  end if
304  if (allocated(ifc%qshft)) then
305    ABI_FREE(ifc%qshft)
306  end if
307  if (allocated(ifc%rpt)) then
308    ABI_FREE(ifc%rpt)
309  end if
310  if (allocated(ifc%wghatm)) then
311    ABI_FREE(ifc%wghatm)
312  end if
313  if (allocated(ifc%rcan)) then
314    ABI_FREE(ifc%rcan)
315  end if
316  if (allocated(ifc%trans)) then
317    ABI_FREE(ifc%trans)
318  end if
319  if (allocated(ifc%dyewq0)) then
320    ABI_FREE(ifc%dyewq0)
321  end if
322  if (allocated(ifc%qibz)) then
323    ABI_FREE(ifc%qibz)
324  end if
325  if (allocated(ifc%wtq)) then
326    ABI_FREE(ifc%wtq)
327  end if
328  if (allocated(ifc%qbz)) then
329    ABI_FREE(ifc%qbz)
330  end if
331  if (allocated(ifc%zeff))then
332    ABI_FREE(ifc%zeff)
333  end if
334  if (allocated(ifc%dynmat)) then
335    ABI_FREE(ifc%dynmat)
336  end if
337 
338  !if (allocated(ifc%dynmat_lr)) then
339  !  ABI_FREE(ifc%dynmat_lr)
340  !end if
341 
342 end subroutine ifc_free

m_ifc/ifc_get_dwdq [ Functions ]

[ Top ] [ m_ifc ] [ Functions ]

NAME

  ifc_get_dwdq

FUNCTION

  Compute phonon group velocities at an arbitrary q-point.

INPUTS

  ifc<ifc_type>=Object containing the dynamical matrix and the IFCs.
  crystal<crystal_t> = Information on the crystalline structure.
  qpt(3)=q-point in reduced coordinates.
  eigvec(2*3*natom*3*natom) = The eigenvectors of the dynamical matrix.

OUTPUT

  dwdq(3,3*natom) = Group velocities e.g. d(omega(q))/dq in Cartesian coordinates.

NOTES

  Using:

    D(q) u(q,nu) = w(q, nu)**2 and <u(q,nu) | u(q,nu')> = \delta_{nu, nu'}

  one can show, using the Hellman-Feynman theorem, that:

    \nabla_q w(q, nu) = 1/(2 w(q, nu))  <u(q, nu)| \nabla_q D(q) | u(q, nu)>

PARENTS

      m_ifc

CHILDREN

      dfpt_phfrq,gtdyn9,nctk_defwrite_nonana_terms

SOURCE

1085 subroutine ifc_get_dwdq(ifc, cryst, qpt, phfrq, eigvec, dwdq)
1086 
1087 
1088 !This section has been created automatically by the script Abilint (TD).
1089 !Do not modify the following lines by hand.
1090 #undef ABI_FUNC
1091 #define ABI_FUNC 'ifc_get_dwdq'
1092 !End of the abilint section
1093 
1094  implicit none
1095 
1096 !Arguments ------------------------------------
1097 !scalars
1098  type(ifc_type),intent(in) :: ifc
1099  type(crystal_t),intent(in) :: cryst
1100 !arrays
1101  real(dp),intent(in) :: qpt(3)
1102  real(dp),intent(in) :: phfrq(3*cryst%natom)
1103  real(dp),intent(in) :: eigvec(2,3*cryst%natom,3*cryst%natom)
1104  real(dp),intent(out) :: dwdq(3,3*cryst%natom)
1105 
1106 !Local variables-------------------------------
1107 !scalars
1108  !integer,save :: enough=0
1109  integer,parameter :: nqpt1=1,option2=2,sumg0=0
1110  integer :: ii,nu,natom3,jj
1111  real(dp) :: hh
1112 !arrays
1113  real(dp) :: dddq(2,3*cryst%natom,3*cryst%natom,3),dot(2),qfd(3)
1114  real(dp) :: omat(2,3*cryst%natom,3*cryst%natom)
1115  real(dp) :: dyew(2,3*cryst%natom,3*cryst%natom)
1116 
1117 ! ************************************************************************
1118 
1119  natom3 = cryst%natom * 3
1120 
1121  ! Generate the analytical part from the interatomic forces
1122  call dynmat_dq(qpt, cryst%natom, ifc%gprim, ifc%nrpt, ifc%rpt, ifc%atmfrc, ifc%wghatm, dddq)
1123 
1124  ! The analytical dynamical matrix dq has been generated
1125  ! in the normalized canonical coordinate system. Now, the
1126  ! phase is modified, in order to recover the usual (xred) coordinate of atoms.
1127  do ii=1,3
1128    call dymfz9(dddq(:,:,:,ii), cryst%natom, nqpt1, ifc%gprim, option2, qpt, ifc%trans)
1129    dddq(:,:,:,ii) = dddq(:,:,:,ii) * ifc%acell(ii)
1130  end do
1131 
1132  if (ifc%dipdip == 1) then
1133    ! Add the gradient of the non-analytical part.
1134    ! Note that we dddq is in cartesian cordinates.
1135    ! For the time being, the gradient is computed with finite difference and step hh.
1136    ! TODO: should generalize ewald9 to compute dq.
1137    !enough = enough + 1
1138    !if (enough <= 5)  MSG_WARNING("phonon velocities with dipdip==1 not yet tested.")
1139    hh = 0.01_dp
1140    do ii=1,3
1141      do jj=-1,1,2
1142        ! qcart --> qred
1143        qfd = zero; qfd(ii) = jj
1144        qfd = matmul(cryst%rprimd, qfd); qfd = qfd / normv(qfd, cryst%gmet, "G")
1145        !write(std_out,*)"normv:",normv(qfd, cryst%gmet, "G")
1146        qfd = qpt + hh * qfd
1147 
1148        call ewald9(ifc%acell,ifc%dielt,dyew,cryst%gmet,ifc%gprim,cryst%natom,qfd,&
1149           cryst%rmet,ifc%rprim,sumg0,cryst%ucvol,cryst%xred,ifc%zeff)
1150        call q0dy3_apply(cryst%natom,ifc%dyewq0,dyew)
1151        dddq(:,:,:,ii) = dddq(:,:,:,ii) + (jj * half / hh) * dyew
1152      end do
1153    end do
1154  end if
1155 
1156  do ii=1,3
1157    call massmult_and_breaksym(cryst%natom, cryst%ntypat, cryst%typat, ifc%amu, dddq(:,:,:,ii))
1158  end do
1159 
1160  ! Compute 1/(2w(q)) <u(q)|dD(q)/dq|u(q)>
1161  do ii=1,3
1162    call zgemm('N','N',natom3,natom3,natom3,cone,dddq(:,:,:,ii),natom3,eigvec,natom3,czero,omat,natom3)
1163    do nu=1,natom3
1164      if (abs(phfrq(nu)) > tol12) then
1165        dot = cg_zdotc(natom3, eigvec(1,1,nu), omat(1,1,nu))
1166        ! abs(w) is needed to get the correct derivative if we have a purely imaginary solution.
1167        dwdq(ii, nu) = dot(1) / (two * abs(phfrq(nu)))
1168      else
1169        dwdq(ii, nu) = zero
1170      end if
1171    end do
1172  end do
1173 
1174 end subroutine ifc_get_dwdq

m_ifc/ifc_getiaf [ Functions ]

[ Top ] [ m_ifc ] [ Functions ]

NAME

 ifc_getiaf

FUNCTION

 Extracts the IFCs needed for the output for one atom in the
 unit cell. Accumulates the results for writing in the NetCDF file.
 Prints to the output file

INPUTS

 Ifc<type(ifc_type)>=Object containing the dynamical matrix and the IFCs.
 ifcana= 0 => no analysis of ifc ; 1 => full analysis
 ifcout= Number of interatomic force constants written in the output file
 iout=unit number for nice output
 zeff(3,3,natom)=effective charge on each atom, versus electric field and atomic displacement
 ia=index of the atom in the unit cell for which the IFCs are being analyzed
 ra(3)=position of atom ia in cartesian coordinates
 list(ifcout)=index of permutation for distances from atom ia in ascending order
 dist(natom,natom,nrpt)=distance from atom ia to atom ib in unit cell irpt.
 invdlt(3,3)=inverse (transpose) of the dielectric tensor
 detdlt=determinant of the dielectric tensor

OUTPUT

 rsiaf(3,3,ifcout)=list of real space IFCs
 sriaf(3,3,ifcout)=list of the short range part of the real space IFCs
 vect(3,3,ifcout)=base vectors for local coordinates (longitudinal/transverse
 indngb(ifcout)=indices in the unit cell of the neighbouring atoms
 posngb(3,ifcout)=position of the neighbouring atoms in cartesian coordinates
 output file

NOTES

 This routine should be executed by one processor only

PARENTS

      m_ifc

CHILDREN

      dfpt_phfrq,gtdyn9,nctk_defwrite_nonana_terms

SOURCE

2104 subroutine ifc_getiaf(Ifc,ifcana,ifcout,iout,zeff,ia,ra,list,&
2105 & dist,invdlt,detdlt,rsiaf,sriaf,vect,indngb,posngb)
2106 
2107 
2108 !This section has been created automatically by the script Abilint (TD).
2109 !Do not modify the following lines by hand.
2110 #undef ABI_FUNC
2111 #define ABI_FUNC 'ifc_getiaf'
2112 !End of the abilint section
2113 
2114 implicit none
2115 
2116 !Arguments -------------------------------
2117 !scalars
2118  integer,intent(in) :: ia,ifcana,ifcout,iout
2119  real(dp), intent(in) :: detdlt
2120  type(ifc_type),intent(inout) :: Ifc
2121 !arrays
2122  real(dp),intent(in) :: invdlt(3,3),ra(3)
2123  real(dp),intent(in) :: dist(Ifc%natom,Ifc%natom,Ifc%nrpt)
2124  real(dp),intent(in) :: zeff(3,3,Ifc%natom)
2125  integer,intent(in) :: list(Ifc%natom*Ifc%nrpt)
2126  integer,intent(out) :: indngb(ifcout)
2127  real(dp),intent(out) :: rsiaf(3,3,ifcout),sriaf(3,3,ifcout),vect(3,3,ifcout),posngb(3,ifcout)
2128 
2129 !Local variables -------------------------
2130 !scalars
2131  integer :: flag,ib,ii,index,jj,kk,mu,nu,irpt
2132  real(dp) :: ew1,rsq,scprod,trace1,trace2,trace3
2133  real(dp) :: yy,dist1
2134  character(len=500) :: message
2135 !arrays
2136  real(dp) :: ewiaf0(3,3),ewiaf1(3,3),ewloc(3,3),ifcloc(3,3)
2137  real(dp) :: rcart(3),rdiff(3),rsloc(3,3)
2138  real(dp) :: srloc(3,3),vect1(3),vect2(3),vect3(3),work(3),xx(3)
2139 
2140 ! *********************************************************************
2141 
2142  if(ifcana==1)then
2143    ! Generate the local coordinate system for the atom ia
2144    index=list(2)
2145    write(std_out,*)index
2146    call canct9(Ifc%acell,Ifc%gprim,ib,index,irpt,Ifc%natom,Ifc%nrpt,Ifc%rcan,rcart,Ifc%rprim,Ifc%rpt)
2147    vect2(1)=rcart(1)-ra(1)
2148    vect2(2)=rcart(2)-ra(2)
2149    vect2(3)=rcart(3)-ra(3)
2150    flag=0
2151    do ii=3,Ifc%natom*Ifc%nrpt
2152      index=list(ii)
2153      call canct9(Ifc%acell,Ifc%gprim,ib,index,irpt,Ifc%natom,Ifc%nrpt,Ifc%rcan,rcart,Ifc%rprim,Ifc%rpt)
2154      vect1(1)=(rcart(1)-ra(1))-vect2(1)
2155      vect1(2)=(rcart(2)-ra(2))-vect2(2)
2156      vect1(3)=(rcart(3)-ra(3))-vect2(3)
2157      scprod=0.0_dp
2158      do jj=1,3
2159        scprod=scprod+vect1(jj)**2
2160      end do
2161      do jj=1,3
2162        vect1(jj)=vect1(jj)/scprod**0.5
2163      end do
2164      scprod=0.0_dp
2165      do jj=1,3
2166        scprod=scprod+vect2(jj)*vect1(jj)
2167      end do
2168      do jj=1,3
2169        work(jj)=vect2(jj)-vect1(jj)*scprod
2170      end do
2171      scprod=0.0_dp
2172      do jj=1,3
2173        scprod=scprod+work(jj)**2
2174      end do
2175      if(scprod>1.0d-10)then
2176        flag=1
2177      end if
2178      if(flag==1)exit
2179    end do
2180    if(flag==0)then
2181      write(message, '(3a)' )&
2182 &     'Unable to find a third atom not aligned',ch10,&
2183 &     'with the two selected ones.'
2184      MSG_BUG(message)
2185    end if
2186    vect2(1)=work(1)/scprod**0.5
2187    vect2(2)=work(2)/scprod**0.5
2188    vect2(3)=work(3)/scprod**0.5
2189    vect3(1)=vect1(2)*vect2(3)-vect1(3)*vect2(2)
2190    vect3(2)=vect1(3)*vect2(1)-vect1(1)*vect2(3)
2191    vect3(3)=vect1(1)*vect2(2)-vect1(2)*vect2(1)
2192    if (iout > 0) then
2193      write(iout, '(a)' )' Third atom defining local coordinates : '
2194      write(iout, '(a,i4,a,i4)' )'     ib = ',ib,'   irpt = ',irpt
2195    end if
2196  end if
2197 
2198  ! Analysis and output of force constants, ordered with respect to the distance from atom ia
2199  do ii=1,ifcout
2200    index=list(ii)
2201    call canct9(Ifc%acell,Ifc%gprim,ib,index,irpt,Ifc%natom,Ifc%nrpt,Ifc%rcan,posngb(:,ii),Ifc%rprim,Ifc%rpt)
2202    indngb(ii)=ib
2203    dist1=dist(ia,ib,irpt)
2204    if (iout > 0) then
2205      write(iout, '(a)' )
2206      write(iout, '(i4,a,i6,a,i8)' )ii,' interaction with atom',ib,' cell',irpt
2207      write(iout, '(a,3es16.6)' )' with coordinates ',posngb(1:3,ii)*(one+tol8)
2208      write(iout, '(a,es16.6)' )' and distance ',dist1
2209    end if
2210 
2211    if(ifcana==1.and.ii/=1)then
2212      vect1(1)=(posngb(1,ii)-ra(1))/dist1
2213      vect1(2)=(posngb(2,ii)-ra(2))/dist1
2214      vect1(3)=(posngb(3,ii)-ra(3))/dist1
2215    end if
2216 
2217    if(Ifc%dipdip==0)then
2218      ! Get the "total" force constants (=real space FC)
2219      ! without taking into account the dipole-dipole interaction
2220      do mu=1,3
2221        do nu=1,3
2222          rsiaf(mu,nu,ii)=Ifc%atmfrc(mu,ia,nu,ib,irpt) * Ifc%wghatm(ia,ib,irpt)
2223        end do
2224      end do
2225      ! Output of the ifcs in cartesian coordinates
2226      if (iout > 0) then
2227        do nu=1,3
2228          write(iout, '(1x,3f9.5)' )(rsiaf(mu,nu,ii)+tol10,mu=1,3)
2229 !       transfer short range and long range
2230          do mu=1,3
2231            Ifc%short_atmfrc(mu,ia,nu,ib,irpt) = rsiaf(mu,nu,ii) + tol10
2232          end do
2233 
2234        end do
2235      end if
2236 
2237      if(ifcana==1)then
2238        ! Further analysis
2239        trace1=rsiaf(1,1,ii)+rsiaf(2,2,ii)+rsiaf(3,3,ii)
2240        if (iout > 0) then
2241          write(iout, '(a,f9.5)' ) '  Trace         ',trace1+tol10
2242        end if
2243        if(ii/=1)then
2244          call axial9(rsiaf(:,:,ii),vect1,vect2,vect3)
2245        end if
2246        if (iout > 0) then
2247          write(iout, '(a)' )' Transformation to local coordinates '
2248          write(iout, '(a,3f16.6)' ) ' First  local vector :',vect1
2249          write(iout, '(a,3f16.6)' ) ' Second local vector :',vect2
2250          write(iout, '(a,3f16.6)' ) ' Third  local vector :',vect3
2251        end if
2252        call ifclo9(rsiaf(:,:,ii),ifcloc,vect1,vect2,vect3)
2253        if (iout > 0) then
2254          do nu=1,3
2255            write(iout, '(1x,3f9.5)' )(ifcloc(mu,nu)+tol10,mu=1,3)
2256          end do
2257        end if
2258 
2259        vect(:,1,ii) = vect1
2260        vect(:,2,ii) = vect2
2261        vect(:,3,ii) = vect3
2262 
2263      end if ! Further analysis finished
2264 
2265    else if(Ifc%dipdip==1)then
2266 
2267 !DEBUG
2268 !    write(iout,'(a)')
2269 !    write(iout,'(a)')' Enter dipdip section, for debugging'
2270 !    write(iout,'(a)')
2271 !ENDDEBUG
2272      ! Get the Coulomb part
2273      do jj=1,3
2274        rdiff(jj)=ra(jj)-posngb(jj,ii)
2275      end do
2276      rsq=0.0_dp
2277      xx(1:3)=0.0_dp
2278      do jj=1,3
2279        do kk=1,3
2280          ewiaf0(jj,kk)=0.0_dp
2281          rsq=rsq+rdiff(jj)*invdlt(kk,jj)*rdiff(kk)
2282          xx(kk)=xx(kk)+invdlt(kk,jj)*rdiff(jj)
2283        end do
2284      end do
2285      yy=sqrt(rsq)
2286      !  Avoid zero denominators in term:
2287      if (sqrt(rsq)>=tol12) then
2288        do mu=1,3
2289          do nu=1,3
2290            ewiaf0(mu,nu)=(-3*xx(nu)*xx(mu)+invdlt(nu,mu)*yy**2)/yy**5/sqrt(detdlt)
2291          end do
2292        end do
2293      else
2294        if (ia/=ib)then
2295          write(message, '(a,a,a,a,a,i5,a,i5,a)' )&
2296 &         'The distance between two atoms vanishes.',ch10,&
2297 &         'This is not allowed.',ch10,&
2298 &         'Action: check the input for the atoms number',ia,' and',ib,'.'
2299          MSG_ERROR(message)
2300        end if
2301      end if
2302 
2303      ! Take into account the effective charge tensor
2304      do mu=1,3
2305        do nu=1,3
2306          ew1=zero
2307          if(ii==1)then
2308            ew1=-Ifc%dyewq0(mu,nu,ia)
2309          end if
2310          do jj=1,3
2311            do kk=1,3
2312              ew1=ew1+zeff(jj,mu,ia)*(zeff(kk,nu,ib)* ewiaf0(jj,kk))
2313            end do
2314          end do
2315          ewiaf1(mu,nu)=ew1
2316        end do
2317      end do
2318      ! Get the short-range force constants and the
2319      ! "total" force constants (=real space FC)
2320      do mu=1,3
2321        do nu=1,3
2322          sriaf(mu,nu,ii)=Ifc%atmfrc(mu,ia,nu,ib,irpt)* Ifc%wghatm(ia,ib,irpt)
2323          rsiaf(mu,nu,ii)=ewiaf1(mu,nu)+sriaf(mu,nu,ii)
2324        end do
2325      end do
2326 
2327      ! Output of the results
2328      if (iout > 0) then
2329        do nu=1,3
2330          write(iout, '(1x,3(3f9.5,1x))' )&
2331 &         (rsiaf(mu,nu,ii) +tol10,mu=1,3),&
2332 &         (ewiaf1(mu,nu)+tol10,mu=1,3),&
2333 &         (sriaf(mu,nu,ii) +tol10,mu=1,3)
2334 
2335 !       transfer short range and long range
2336          do mu=1,3
2337            Ifc%short_atmfrc(mu,ia,nu,ib,irpt) = sriaf(mu,nu,ii) + tol10
2338            Ifc%ewald_atmfrc(mu,ia,nu,ib,irpt) = ewiaf1(mu,nu) + tol10
2339          end do
2340        end do
2341      end if
2342 
2343      if(ifcana==1)then
2344        ! Further analysis
2345        if (iout > 0) then
2346          write(iout, '(a)' )' Traces (and ratios) :'
2347        end if
2348        trace1=rsiaf(1,1,ii)+rsiaf(2,2,ii)+rsiaf(3,3,ii)
2349        trace2=ewiaf1(1,1)+ewiaf1(2,2)+ewiaf1(3,3)
2350        trace3=sriaf(1,1,ii)+sriaf(2,2,ii)+sriaf(3,3,ii)
2351        if (iout > 0) then
2352          write(iout,'(3(f9.5,17x))')trace1+tol10,trace2+tol10,trace3+tol10
2353          write(iout,'(3(f9.5,17x))')1.0,trace2/trace1+tol10,trace3/trace1+tol10 !
2354        end if
2355 
2356        if(ii/=1)then
2357          call axial9(rsiaf(:,:,ii),vect1,vect2,vect3)
2358        end if
2359        if (iout > 0) then
2360          write(iout, '(a)' )' Transformation to local coordinates '
2361          write(iout, '(a,3f16.6)' )' First  local vector :',vect1
2362          write(iout, '(a,3f16.6)' )' Second local vector :',vect2
2363          write(iout, '(a,3f16.6)' )' Third  local vector :',vect3
2364        end if
2365        call ifclo9(rsiaf(:,:,ii),rsloc,vect1,vect2,vect3)
2366        call ifclo9(ewiaf1,ewloc,vect1,vect2,vect3)
2367        call ifclo9(sriaf(:,:,ii),srloc,vect1,vect2,vect3)
2368        if (iout > 0) then
2369          do nu=1,3
2370            write(iout, '(1x,3(3f9.5,1x))' )&
2371 &           (rsloc(mu,nu)+tol10,mu=1,3),&
2372 &           (ewloc(mu,nu)+tol10,mu=1,3),&
2373 &           (srloc(mu,nu)+tol10,mu=1,3)
2374          end do
2375          if(ii/=1)then
2376            write(iout, '(a)' )' Ratio with respect to the longitudinal ifc'
2377          else
2378            write(iout, '(a)' )' Ratio with respect to the (1,1) element'
2379          end if
2380          do nu=1,3
2381            write(iout, '(1x,3(3f9.5,1x))' )&
2382 &           (rsloc(mu,nu)/rsloc(1,1)+tol10,mu=1,3),&
2383 &           (ewloc(mu,nu)/rsloc(1,1)+tol10,mu=1,3),&
2384 &           (srloc(mu,nu)/rsloc(1,1)+tol10,mu=1,3)
2385          end do
2386        end if
2387 
2388        vect(:,1,ii) = vect1
2389        vect(:,2,ii) = vect2
2390        vect(:,3,ii) = vect3
2391 
2392      end if ! Further analysis finished
2393    end if ! End the condition on dipdip
2394  end do ! End loop over all atoms in BigBox:
2395 
2396 end subroutine ifc_getiaf

m_ifc/ifc_init [ Functions ]

[ Top ] [ m_ifc ] [ Functions ]

NAME

  ifc_init

FUNCTION

  Initialize the dynamical matrix as well as the IFCs.
  taking into account the dipole-dipole interaction.

INPUTS

 crystal<type(crystal_t)> = Information on the crystalline structure.
 ddb<type(ddb_type)> = Database with derivatives.
 brav=bravais lattice (1 or -1=simple lattice,2=face centered lattice, 3=centered lattice,4=hexagonal lattice)
 asr= Option for the imposition of the ASR
   0 => no ASR,
   1 => modify "asymmetrically" the diagonal element
   2 => modify "symmetrically" the diagonal element
 symdynmat=if 1, (re)symmetrize the dynamical matrix, except if Gamma wavevector with electric field added.
 dipdip=
   if 0, no dipole-dipole interaction was subtracted in atmfrc
   if 1, atmfrc has been build without dipole-dipole part
 rfmeth =
   1 if non-stationary block
   2 if stationary block
   3 if third order derivatives
 dielt(3,3)=dielectric tensor.
 zeff(3,3,natom)=effective charge on each atom, versus electric field and atomic displacement
 prtsrlr: TODO: TO BE REMOVED
 enunit: TODO: TO BE REMOVED
 dielt(3,3)=dielectric tensor
 ngqpt_in = input values of ngqpt
 nqshft=Number of shifths in q-grid.
 q1shft(3,nqshft)=Shifts for q-grid
 nsphere=number of atoms to be included in the cut-off sphere for interatomic force constant.
    0: maximum extent allowed by the grid.
  > 0: Apply cutoff
   -1: Analyze the effect of different nsphere values on the phonon spectrum, in particular the
       frequencies around gamma.
 rifcsph=radius for cutoff of IFC.
 comm=MPI communicator.
 [Ifc_coarse]=Optional.

OUTPUT

 Ifc<ifc_type>=Object containing the dynamical matrix and the IFCs.

PARENTS

      anaddb,eph,m_effective_potential_file,m_gruneisen,m_tdep_abitypes

CHILDREN

      dfpt_phfrq,gtdyn9,nctk_defwrite_nonana_terms

SOURCE

399 subroutine ifc_init(ifc,crystal,ddb,brav,asr,symdynmat,dipdip,&
400   rfmeth,ngqpt_in,nqshft,q1shft,dielt,zeff,nsphere,rifcsph,&
401   prtsrlr,enunit, & ! TODO: TO BE REMOVED
402   comm, &
403   Ifc_coarse) ! Optional
404 
405 
406 !This section has been created automatically by the script Abilint (TD).
407 !Do not modify the following lines by hand.
408 #undef ABI_FUNC
409 #define ABI_FUNC 'ifc_init'
410 !End of the abilint section
411 
412  implicit none
413 
414 !Arguments ------------------------------------
415  integer,intent(in) :: asr,brav,dipdip,symdynmat,nqshft,rfmeth,nsphere,comm
416  real(dp),intent(in) :: rifcsph
417  type(ifc_type),intent(inout) :: Ifc
418  type(crystal_t),intent(in) :: Crystal
419  type(ddb_type),intent(in) :: ddb
420  type(ifc_type),optional,intent(in) :: Ifc_coarse
421 !arrays
422  integer,intent(in) :: ngqpt_in(3)
423  real(dp),intent(in) :: q1shft(3,nqshft)
424  real(dp),intent(in) :: dielt(3,3),zeff(3,3,Crystal%natom)
425 !anaddb variables (TO BE REMOVED)
426  integer,intent(in) :: prtsrlr,enunit
427 !end anaddb variables
428 
429 !Local variables -------------------------
430 !scalars
431  integer,parameter :: timrev1=1,iout0=0,chksymbreak0=0
432  integer :: mpert,iout,iqpt,mqpt,nsym,ntypat,iq_ibz,iq_bz,ii,natom
433  integer :: nqbz,option,plus,sumg0,irpt,irpt_new,new_wght
434  integer :: nprocs,my_rank,ierr
435  real(dp),parameter :: qphnrm=one
436  real(dp) :: xval,cpu,wall,gflops,rcut_min
437  real(dp) :: r_inscribed_sphere
438  character(len=500) :: message
439  type(ifc_type) :: ifc_tmp
440 !arrays
441  integer :: ngqpt(9),qptrlatt(3,3)
442  integer,allocatable :: qmissing(:),ibz2bz(:)
443  real(dp) :: gprim(3,3),rprim(3,3),qpt(3),rprimd(3,3)
444  real(dp):: rcan(3,Crystal%natom),trans(3,Crystal%natom),dyewq0(3,3,Crystal%natom)
445  real(dp) :: displ_cart(2*3*Crystal%natom*3*Crystal%natom)
446  real(dp) :: phfrq(3*Crystal%natom)
447  real(dp) :: eigvec(2,3,Crystal%natom,3,Crystal%natom)
448  real(dp),allocatable :: dyew(:,:,:,:,:),out_d2cart(:,:,:,:,:)
449  real(dp),allocatable :: dynmatfull(:,:,:,:,:,:),dynmat_sr(:,:,:,:,:,:),dynmat_lr(:,:,:,:,:,:) ! for OmegaSRLR
450  real(dp),allocatable :: wtq(:),wtq_folded(:),qbz(:,:)
451 
452 !******************************************************************
453 
454  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
455  call cwtime(cpu, wall, gflops, "start")
456 
457  ! This dimension should be encapsulated somewhere. We don't want to
458  ! change the entire code if someone adds a new kind of perturbation.
459  mpert = Crystal%natom + 6; iout = ab_out
460 
461  rprim = ddb%rprim; gprim = ddb%gprim
462 
463  nsym = Crystal%nsym
464  natom = Crystal%natom
465  ntypat = Crystal%ntypat
466  rprimd = Crystal%rprimd
467 
468  ngqpt=0; ngqpt(1:3)=ngqpt_in(1:3)
469 
470 ! Copy important parameters in Ifc
471  Ifc%natom = natom
472  Ifc%mpert = mpert
473  Ifc%asr = asr
474  Ifc%brav = brav
475  Ifc%dipdip = dipdip
476  Ifc%symdynmat = symdynmat
477  Ifc%nqshft = nqshft
478  call alloc_copy(q1shft(:,1:Ifc%nqshft),Ifc%qshft)
479  Ifc%ngqpt = ngqpt_in(1:3)
480  Ifc%rprim = ddb%rprim
481  Ifc%gprim = ddb%gprim
482  Ifc%acell = ddb%acell
483 
484  ! Check if the rprim are coherent with the choice used in the interatomic forces generation
485  call chkrp9(Ifc%brav,rprim)
486 
487  dyewq0 = zero
488  if (Ifc%dipdip==1 .and. (Ifc%asr==1.or.Ifc%asr==2)) then
489    ! Calculation of the non-analytical part for q=0
490    sumg0=0
491    qpt(:)=zero
492    ABI_MALLOC(dyew,(2,3,natom,3,natom))
493    call ewald9(ddb%acell,dielt,dyew,Crystal%gmet,gprim,natom,qpt,Crystal%rmet,rprim,sumg0,Crystal%ucvol,Crystal%xred,zeff)
494    call q0dy3_calc(natom,dyewq0,dyew,Ifc%asr)
495    ABI_FREE(dyew)
496  end if
497 
498  ! Sample the Brillouin zone
499  option=1
500  qptrlatt = 0; qptrlatt(1,1)=ngqpt(1); qptrlatt(2,2)=ngqpt(2); qptrlatt(3,3)=ngqpt(3)
501  mqpt=ngqpt(1)*ngqpt(2)*ngqpt(3)*nqshft
502  if (Ifc%brav==2) mqpt=mqpt/2
503  if (Ifc%brav==3) mqpt=mqpt/4
504 
505  ABI_MALLOC(qbz,(3,mqpt))
506  call smpbz(Ifc%brav,ab_out,qptrlatt,mqpt,nqbz,nqshft,option,q1shft,qbz)
507 
508  ! Find the irreducible zone (qibz)
509  ABI_MALLOC(ibz2bz, (nqbz))
510  ABI_MALLOC(wtq_folded, (nqbz))
511  ABI_MALLOC(wtq, (nqbz))
512  wtq = one / nqbz ! Weights sum up to one
513 
514  call symkpt(chksymbreak0,crystal%gmet,ibz2bz,iout0,qbz,nqbz,ifc%nqibz,crystal%nsym,&
515    crystal%symrec,timrev1,wtq,wtq_folded)
516 
517  ABI_MALLOC(ifc%qibz, (3,ifc%nqibz))
518  ABI_MALLOC(ifc%wtq, (ifc%nqibz))
519  do iq_ibz=1,ifc%nqibz
520    ifc%qibz(:,iq_ibz) = qbz(:, ibz2bz(iq_ibz))
521    ifc%wtq(iq_ibz) = wtq_folded(ibz2bz(iq_ibz))
522  end do
523  ABI_FREE(ibz2bz)
524  ABI_FREE(wtq_folded)
525  ABI_FREE(wtq)
526 
527  ABI_MALLOC(Ifc%dynmat,(2,3,natom,3,natom,nqbz))
528 
529 ! Find symmetrical dynamical matrices
530  if (.not.present(Ifc_coarse)) then
531 
532    ! Each q-point in the BZ mush be the symmetrical of one of the qpts in the ddb file.
533    call symdm9(ddb%flg,ddb%nrm,ddb%qpt,ddb%typ,ddb%val,&
534 &    Ifc%dynmat,gprim,Crystal%indsym,mpert,natom,ddb%nblok,nqbz,nsym,rfmeth,rprim,qbz,&
535 &    Crystal%symrec,Crystal%symrel)
536 
537  else
538    ! Symmetrize the qpts in the BZ using the q-points in the ddb.
539    ! Then use Ifc_coarse to fill the missing entries with Fourier interpolated matrices.
540    !
541    ! TODO: The previous version of refineblk was hacking the DDB database to add the q-points in the **IBZ**
542    ! Then D(q) was symmetrized in symdm9. This version avoids the symmetrization: the q-points
543    ! in the BZ that are not in the coarse q-mesh are obtained by an explicit FT.
544    ! This means that the final D(q) may break some symmetry in q-space if the FT does not preserve it.
545    ! The most elegant approach would be to get D(q_ibz) via FT if q_ibz is not in the coarse mesh and then
546    ! call symdm9 to get D(q) for each q point in the star of q_ibz.
547    call wrtout(std_out,"Will fill missing qpoints in the full BZ using the coarse q-mesh","COLL")
548 
549    call symdm9(ddb%flg,ddb%nrm,ddb%qpt,ddb%typ,ddb%val,&
550 &    Ifc%dynmat,gprim,Crystal%indsym,mpert,natom,ddb%nblok,nqbz,nsym,rfmeth,rprim,qbz,&
551 &    Crystal%symrec,Crystal%symrel,qmissing=qmissing)
552 
553    ! Compute dynamical matrix with Fourier interpolation on the coarse q-mesh.
554    write(message,"(a,i0,a)")"Will use Fourier interpolation to construct D(q) for ",size(qmissing)," q-points"
555    call wrtout(std_out,message,"COLL")
556 
557    ABI_MALLOC(out_d2cart, (2,3,natom,3,natom))
558    do ii=1,size(qmissing)
559      iq_bz = qmissing(ii)
560      qpt = qbz(:,iq_bz)
561      ! TODO: check dipdip option and phase, but I think this is correct!
562      call ifc_fourq(Ifc_coarse,Crystal,qpt,phfrq,displ_cart,out_d2cart=out_d2cart)
563      Ifc%dynmat(:,:,:,:,:,iq_bz) = out_d2cart
564 
565    end do
566 
567    ABI_FREE(qmissing)
568    ABI_FREE(out_d2cart)
569  end if
570 
571 !OmegaSRLR: Store full dynamical matrix for decomposition into short- and long-range parts
572  ABI_MALLOC(dynmatfull,(2,3,natom,3,natom,nqbz))
573  dynmatfull=Ifc%dynmat
574 
575  if (Ifc%dipdip==1) then
576    ! Take off the dipole-dipole part of the dynamical matrix
577    call wrtout(std_out, " Will extract the dipole-dipole part for every wavevector")
578    ABI_MALLOC(dyew,(2,3,natom,3,natom))
579 
580    do iqpt=1,nqbz
581      if (mod(iqpt, nprocs) /= my_rank) then ! mpi-parallelism
582        ifc%dynmat(:,:,:,:,:,iqpt) = zero; cycle
583      end if
584      qpt(:)=qbz(:,iqpt)
585      sumg0=0
586      call ewald9(ddb%acell,dielt,dyew,Crystal%gmet,gprim,natom,qpt,Crystal%rmet,rprim,sumg0,Crystal%ucvol,Crystal%xred,zeff)
587      call q0dy3_apply(natom,dyewq0,dyew)
588      plus=0
589      call nanal9(dyew,Ifc%dynmat,iqpt,natom,nqbz,plus)
590    end do
591 
592    call xmpi_sum(ifc%dynmat, comm, ierr)
593    ABI_FREE(dyew)
594  end if
595 
596 ! OmegaSRLR: Store the short-range dynmat and compute long-range as difference
597  ABI_MALLOC(dynmat_sr,(2,3,natom,3,natom,nqbz))
598  ABI_MALLOC(dynmat_lr,(2,3,natom,3,natom,nqbz))
599  dynmat_sr=Ifc%dynmat
600  dynmat_lr=dynmatfull-dynmat_sr
601 
602 ! Now, take care of the remaining part of the dynamical matrix
603 ! Move to canonical normalized coordinates
604  call canat9(Ifc%brav,natom,rcan,rprim,trans,Crystal%xred)
605 
606 ! Multiply the dynamical matrix by a phase shift
607  option=1
608  call dymfz9(Ifc%dynmat,natom,nqbz,gprim,option,qbz,trans)
609 
610 ! Create the Big Box of R vectors in real space and compute the number of points (cells) in real space
611  call make_bigbox(Ifc%brav,ifc_tmp%cell,ngqpt,nqshft,rprim,ifc_tmp%nrpt,ifc_tmp%rpt)
612 
613 ! Weights associated to these R points and to atomic pairs
614 ! MG FIXME: Why ngqpt is intent(inout)?
615  ABI_MALLOC(ifc_tmp%wghatm,(natom,natom,ifc_tmp%nrpt))
616  call wght9(Ifc%brav,gprim,natom,ngqpt,nqbz,nqshft,ifc_tmp%nrpt,q1shft,rcan,&
617 &     ifc_tmp%rpt,rprimd,r_inscribed_sphere,ifc_tmp%wghatm)
618 
619 ! Fourier transformation of the dynamical matrices (q-->R)
620  ABI_MALLOC(ifc_tmp%atmfrc,(3,natom,3,natom,ifc_tmp%nrpt))
621  call ftifc_q2r(ifc_tmp%atmfrc,Ifc%dynmat,gprim,natom,nqbz,ifc_tmp%nrpt,ifc_tmp%rpt,qbz)
622 
623 ! Eventually impose Acoustic Sum Rule to the interatomic forces
624  if (Ifc%asr>0) then
625    call asrif9(Ifc%asr,ifc_tmp%atmfrc,natom,ifc_tmp%nrpt,ifc_tmp%rpt,ifc_tmp%wghatm)
626  end if
627 
628 !*** The interatomic forces have been calculated ! ***
629  write(message, '(2a)')ch10,' The interatomic forces have been obtained '
630  call wrtout(ab_out,message,'COLL')
631  call wrtout(std_out,message,'COLL')
632 
633  ! Apply cutoff on ifc if needed
634  if (nsphere > 0 .or. abs(rifcsph) > tol10) then
635    call wrtout(std_out, ' Apply cutoff on IFCs.')
636    call wrtout(std_out, sjoin(" nsphere:", itoa(nsphere), ", rifcsph:", ftoa(rifcsph)))
637    call corsifc9(ddb%acell,gprim,natom,ifc_tmp%nrpt,nsphere,rifcsph,rcan,rprim,ifc_tmp%rpt,rcut_min,ifc_tmp%wghatm)
638    if (Ifc%asr > 0) then
639      call wrtout(std_out, ' Enforcing ASR on cutoffed IFCs.')
640      call asrif9(Ifc%asr,ifc_tmp%atmfrc,natom,ifc_tmp%nrpt,ifc_tmp%rpt,ifc_tmp%wghatm)
641    end if
642  end if
643 
644  ! Only conserve the necessary points in rpt: in the FT algorithm the order of the points is unimportant
645  ! In the case of effective potential, we need to keep all the points
646  Ifc%nrpt = 0
647  do irpt=1,ifc_tmp%nrpt
648    if (sum(ifc_tmp%wghatm(:,:,irpt)) /= 0) Ifc%nrpt = Ifc%nrpt+1
649  end do
650 
651  ABI_CALLOC(Ifc%atmfrc,(3,natom,3,natom,Ifc%nrpt))
652  ABI_CALLOC(Ifc%rpt,(3,Ifc%nrpt))
653  ABI_CALLOC(Ifc%cell,(3,Ifc%nrpt))
654  ABI_CALLOC(Ifc%wghatm,(natom,natom,Ifc%nrpt))
655  ABI_CALLOC(Ifc%short_atmfrc,(3,natom,3,natom,Ifc%nrpt))
656  ABI_CALLOC(Ifc%ewald_atmfrc,(3,natom,3,natom,Ifc%nrpt))
657 
658  irpt_new = 1
659  do irpt = 1, ifc_tmp%nrpt
660    if (sum(ifc_tmp%wghatm(:,:,irpt)) /= 0) then
661      Ifc%atmfrc(:,:,:,:,irpt_new) = ifc_tmp%atmfrc(:,:,:,:,irpt)
662      Ifc%rpt(:,irpt_new) = ifc_tmp%rpt(:,irpt)
663      Ifc%wghatm(:,:,irpt_new) = ifc_tmp%wghatm(:,:,irpt)
664      Ifc%cell(:,irpt_new) = ifc_tmp%cell(:,irpt)
665      Ifc%r_inscribed_sphere = r_inscribed_sphere
666      irpt_new = irpt_new + 1
667    end if
668  end do
669 
670  ! Copy other useful arrays.
671  Ifc%dielt = dielt
672  Ifc%nqbz = nqbz
673 
674  call alloc_copy(rcan, Ifc%rcan)
675  call alloc_copy(trans, Ifc%trans)
676  call alloc_copy(dyewq0, Ifc%dyewq0)
677  call alloc_copy(qbz(:,1:nqbz), Ifc%qbz)
678  call alloc_copy(zeff, Ifc%zeff)
679  call alloc_copy(ddb%amu, Ifc%amu)
680 
681  call ifc_free(ifc_tmp)
682 
683  ! Compute min/max ph frequency with ab-initio q-mesh.
684  ifc%omega_minmax(1) = huge(one); ifc%omega_minmax(2) = -huge(one)
685  do iq_ibz=1,ifc%nqibz
686    if (mod(iq_ibz, nprocs) /= my_rank) cycle ! mpi-parallelism
687    call ifc_fourq(ifc, crystal, ifc%qibz(:,iq_ibz), phfrq, displ_cart)
688    ifc%omega_minmax(1) = min(ifc%omega_minmax(1), minval(phfrq))
689    ifc%omega_minmax(2) = max(ifc%omega_minmax(2), maxval(phfrq))
690  end do
691  xval = ifc%omega_minmax(1); call xmpi_min(xval, ifc%omega_minmax(1), comm, ierr)
692  xval = ifc%omega_minmax(2); call xmpi_max(xval, ifc%omega_minmax(2), comm, ierr)
693  ! Enlarge boundaries by 30 cm-1
694  ifc%omega_minmax(1) = ifc%omega_minmax(1) - 30.0_dp/Ha_cmm1
695  ifc%omega_minmax(2) = ifc%omega_minmax(2) + 30.0_dp/Ha_cmm1
696 
697  ! TODO (This is to be suppressed in a future version)
698  if (prtsrlr == 1) then
699    ! Check that the starting values are well reproduced.
700    write(message, '(2a)' )' mkifc9 : now check that the starting values ',&
701      ' are reproduced after the use of interatomic forces '
702    call wrtout(std_out,message,"COLL")
703    do iqpt=1,nqbz
704      qpt(:)=Ifc%qbz(:,iqpt)
705      call ifc_fourq(Ifc,Crystal,qpt,phfrq,displ_cart,out_eigvec=eigvec)
706 
707      ! OmegaSRLR: Perform decomposition of dynamical matrix
708      ! MG: FIXME I don't think the implementation is correct when q !=0
709      if (prtsrlr==1) then
710        call omega_decomp(ddb%amu,natom,ntypat,Crystal%typat,dynmatfull,dynmat_sr,dynmat_lr,iqpt,nqbz,eigvec)
711      end if
712 
713      ! Write the phonon frequencies (this is for checking purposes).
714      ! Note: these phonon frequencies are not written on unit iout, only on unit std_out.
715      call dfpt_prtph(displ_cart,0,enunit,-1,natom,phfrq,qphnrm,qpt)
716    end do
717  end if
718 
719  ! OmegaSRLR: deallocate memory used by dynmat decomposition
720  ABI_FREE(dynmatfull)
721  ABI_FREE(dynmat_sr)
722  ABI_FREE(dynmat_lr)
723  ABI_FREE(qbz)
724 
725  if (nsphere == -1) call ifc_autocutoff(ifc, crystal, comm)
726 
727  call cwtime(cpu, wall, gflops, "stop")
728  write(message,"(2(a,f6.2))")" ifc_init: cpu: ",cpu,", wall: ",wall
729  call wrtout(std_out,message,'COLL')
730 
731 end subroutine ifc_init

m_ifc/ifc_init_fromFile [ Functions ]

[ Top ] [ m_ifc ] [ Functions ]

NAME

  ifc_init_fromFile

FUNCTION

  Need to be updated

INPUTS

OUTPUT

PARENTS

      anaddb,eph,m_effective_potential_file,m_gruneisen,m_tdep_abitypes

CHILDREN

      dfpt_phfrq,gtdyn9,nctk_defwrite_nonana_terms

SOURCE

755 subroutine ifc_init_fromFile(dielt,filename,Ifc,natom,ngqpt,nqshift,qshift,ucell_ddb,zeff,comm)
756 
757 
758 !This section has been created automatically by the script Abilint (TD).
759 !Do not modify the following lines by hand.
760 #undef ABI_FUNC
761 #define ABI_FUNC 'ifc_init_fromFile'
762 !End of the abilint section
763 
764  implicit none
765 
766 !Arguments ------------------------------------
767  !scalars
768  integer,intent(in) :: nqshift,comm
769  integer,intent(inout) :: natom
770  !arrays
771  integer,intent(in) :: ngqpt(3)
772  real(dp),intent(in) :: qshift(3,nqshift)
773  character(len=*),intent(in) :: filename
774  type(ifc_type),intent(out) :: Ifc
775  real(dp),intent(inout) :: dielt(3,3)
776  real(dp),allocatable,intent(inout) :: zeff(:,:,:)
777  type(crystal_t),intent(out) :: ucell_ddb
778 !Local variables -------------------------
779  !scalars
780  integer :: dipdip,i,iblok,iblok_tmp
781  logical :: file_exists
782  !arrays
783  integer,allocatable :: atifc(:)
784  type(ddb_type) :: ddb
785  type(ddb_hdr_type) :: ddb_hdr
786  character(len=500) :: msg
787 
788 !******************************************************************
789 
790  !check if ddb file exists
791  inquire(file=filename, exist=file_exists)
792 
793  if (file_exists .eqv. .true.)then
794    !Reading the ddb
795    call ddb_hdr_open_read(ddb_hdr,filename,2,DDB_VERSION,comm,dimonly=1)
796 
797    natom = ddb_hdr%natom
798    ABI_ALLOCATE(atifc,(ddb_hdr%natom))
799    do i=1,ddb_hdr%natom
800      atifc(i)=i
801    end do
802 
803    call ddb_from_file(ddb,filename,1,ddb_hdr%natom,ddb_hdr%natom,atifc,ucell_ddb,comm)
804 
805    else
806      write (msg, "(a,a,a)") "File ", filename , " is not present in the directory"
807      MSG_ERROR(msg)
808    end if
809 
810    ! Get Dielectric Tensor and Effective Charges
811    ABI_ALLOCATE(zeff,(3,3,natom))
812    iblok = ddb_get_dielt_zeff(ddb,ucell_ddb,1,1,0,dielt,zeff)
813 
814    ! Try to get dielt, in case just the DDE are present
815    if (iblok == 0) then
816      iblok_tmp = ddb_get_dielt(ddb,1,dielt)
817    end if
818 
819    ! ifc to be calculated for interpolation
820    write(msg, '(a,a,(80a),a,a,a,a)' ) ch10,('=',i=1,80),ch10,ch10,&
821      &   ' Calculation of the interatomic forces ',ch10
822    call wrtout(std_out,msg,'COLL')
823    call wrtout(ab_out,msg,'COLL')
824    if ((maxval(abs(zeff)) .lt. tol10) .OR. (maxval(dielt) .gt. 100000.0)) then
825      dipdip=0
826    else
827      dipdip=1
828    end if
829    call ifc_init(Ifc,ucell_ddb,ddb,1,1,1,dipdip,1,ngqpt,nqshift,&
830 &                qshift,dielt,zeff,0,0.0_dp,0,1,comm)
831 
832 
833 !  Free them all
834    ABI_DEALLOCATE(atifc)
835    call ddb_free(ddb)
836    call ddb_hdr_free(ddb_hdr)
837 
838  end subroutine ifc_init_fromFile

m_ifc/ifc_outphbtrap [ Functions ]

[ Top ] [ m_ifc ] [ Functions ]

NAME

 ifc_outphbtrap

FUNCTION

  Print out phonon frequencies on regular grid for BoltzTrap
  Flag in input file is outboltztrap=1

INPUTS

  ifc<ifc_type>=Stores data related to interatomic force constants.
  Crystal<crystal_t>=Info on the crystal structure
  basename = file name for output to disk
  ngqpt(3)=Divisions of the q-mesh
  nqshft=Number of shifts
  qshft(3,nqshft)=Shifts of the q-mesh.

OUTPUT

  only write to file. This routine should be called by a single processor.

PARENTS

      anaddb,eph

CHILDREN

      dfpt_phfrq,gtdyn9,nctk_defwrite_nonana_terms

SOURCE

2630 subroutine ifc_outphbtrap(ifc, cryst, ngqpt, nqshft, qshft, basename)
2631 
2632 
2633 !This section has been created automatically by the script Abilint (TD).
2634 !Do not modify the following lines by hand.
2635 #undef ABI_FUNC
2636 #define ABI_FUNC 'ifc_outphbtrap'
2637 !End of the abilint section
2638 
2639  implicit none
2640 
2641 !Arguments -------------------------------
2642 !scalars
2643  integer,intent(in) :: nqshft
2644  character(len=*),intent(in) :: basename
2645  type(ifc_type),intent(in) :: ifc
2646  type(crystal_t),intent(in) :: cryst
2647 !arrays
2648  integer,intent(in) :: ngqpt(3)
2649  real(dp),intent(in) :: qshft(3,nqshft)
2650 
2651 !Local variables -------------------------
2652 !scalars
2653  integer,parameter :: qptopt1=1
2654  integer :: natom,imode,iq_ibz,nqbz,nqibz, nreals,unit_btrap,iatom,idir
2655  character(len=500) :: msg,format_nreals,format_line_btrap
2656  character(len=fnlen) :: outfile
2657 !arrays
2658  integer :: qptrlatt(3,3)
2659  real(dp) :: d2cart(2,3,cryst%natom,3,cryst%natom),displ(2*3*cryst%natom*3*cryst%natom)
2660  real(dp) :: phfrq(3*cryst%natom),qphon(3)
2661  real(dp),allocatable :: qbz(:,:),qibz(:,:),wtq(:)
2662 
2663 ! *********************************************************************
2664 
2665  DBG_ENTER("COLL")
2666 
2667  natom = cryst%natom
2668 
2669  ! Setup IBZ, weights and BZ. Always use q --> -q symmetry for phonons even in systems wo inversion
2670  qptrlatt = 0; qptrlatt(1,1) = ngqpt(1); qptrlatt(2,2) = ngqpt(2); qptrlatt(3,3) = ngqpt(3)
2671  call kpts_ibz_from_kptrlatt(cryst, qptrlatt, qptopt1, nqshft, qshft, nqibz, qibz, wtq, nqbz, qbz)
2672 
2673  outfile = trim(basename) // '_BTRAP'
2674  write(msg, '(3a)')ch10,' Will write phonon FREQS in BoltzTrap format to file ',trim(outfile)
2675  call wrtout(ab_out,msg,'COLL')
2676  call wrtout(std_out,msg,'COLL')
2677 
2678  if (open_file(outfile,msg,newunit=unit_btrap,status="replace") /= 0) then
2679    MSG_ERROR(msg)
2680  end if
2681 
2682  write (unit_btrap,'(a)') '#'
2683  write (unit_btrap,'(a)') '# ABINIT package : Boltztrap phonon file. With old BT versions remove this header before feeding to BT'
2684  write (unit_btrap,'(a)') '#    for compatibility with PHON output the freq are in Ry (before the square)'
2685  write (unit_btrap,'(a)') '#'
2686  write (unit_btrap,'(a)') '#    nq, nband  '
2687  write (unit_btrap,'(a)') '#  qx, qy, qz   '
2688  write (unit_btrap,'(a)') '#  qpt weight   '
2689  write (unit_btrap,'(a)') '#  freq_1^2, dynmat column for mode 1 '
2690  write (unit_btrap,'(a)') '#  etc for mode 2,3,4... qpt 2,3,4... '
2691  write (unit_btrap,'(2I6)') nqibz, 3*natom
2692 
2693 ! Loop over irreducible q-points
2694  do iq_ibz=1,nqibz
2695    qphon(:)=qibz(:,iq_ibz)
2696 
2697    call ifc_fourq(ifc, cryst, qphon, phfrq, displ, out_d2cart=d2cart)
2698 
2699    write (unit_btrap,'(3E20.10)') qphon
2700    write (unit_btrap,'(E20.10)') wtq(iq_ibz)
2701    nreals=1+2*3*natom
2702    call appdig(nreals,'(',format_nreals)
2703    format_line_btrap=trim(format_nreals)//'E20.10)'
2704    do iatom = 1, natom
2705      do idir = 1, 3
2706        imode = idir + 3*(iatom-1)
2707 !      factor two for Ry output - this may change in definitive BT and abinit formats
2708        write (unit_btrap,trim(format_line_btrap))phfrq(imode)*two,d2cart(1:2,1:3,1:natom,idir,iatom)
2709      end do
2710    end do
2711 
2712  end do !irred q-points
2713  close (unit_btrap)
2714 
2715  ABI_DEALLOCATE(qibz)
2716  ABI_DEALLOCATE(qbz)
2717  ABI_DEALLOCATE(wtq)
2718 
2719  DBG_EXIT("COLL")
2720 
2721 end subroutine ifc_outphbtrap

m_ifc/ifc_print [ Functions ]

[ Top ] [ m_ifc ] [ Functions ]

NAME

  ifc_print

FUNCTION

  Print info on the object

INPUTS

  [unit]=Unit number for output. Defaults to std_out
  [prtvol]=Verbosity level.
  [header]=String to be printed as header for additional info.

OUTPUT

  Only printing

PARENTS

      anaddb,eph,m_tdep_abitypes

CHILDREN

      dfpt_phfrq,gtdyn9,nctk_defwrite_nonana_terms

SOURCE

867 subroutine ifc_print(ifc,header,unit,prtvol)
868 
869 
870 !This section has been created automatically by the script Abilint (TD).
871 !Do not modify the following lines by hand.
872 #undef ABI_FUNC
873 #define ABI_FUNC 'ifc_print'
874 !End of the abilint section
875 
876  implicit none
877 
878 !Arguments ------------------------------------
879 !scalars
880  integer,optional,intent(in) :: unit,prtvol
881  character(len=*),optional,intent(in) :: header
882  type(ifc_type),intent(in) :: ifc
883 
884 !Local variables-------------------------------
885  integer :: unt,my_prtvol,iatom,ii
886  character(len=500) :: msg
887 ! *********************************************************************
888 
889  unt = std_out; if (present(unit)) unt = unit
890  my_prtvol = 0; if (present(prtvol)) my_prtvol = prtvol
891 
892  msg = ' ==== Info on the ifc% object ==== '
893  if (present(header)) msg = ' ==== '//trim(adjustl(header))//' ==== '
894  call wrtout(unt, msg)
895 
896  call wrtout(unt,' Real(R)+Recip(G) space primitive vectors, cartesian coordinates (Bohr,Bohr^-1):')
897  do ii=1,3
898    write(msg,'(1x,a,i1,a,3f11.7,2x,a,i1,a,3f11.7)')&
899     'R(',ii,')=',ifc%rprim(:,ii),'G(',ii,')=',ifc%gprim(:,ii)
900    call wrtout(unt,msg)
901  end do
902  call wrtout(unt, sjoin(" acell:", ltoa(ifc%acell)))
903 
904  call wrtout(unt, sjoin(" Acoustic Sum Rule option (asr):", itoa(ifc%asr)))
905  call wrtout(unt, sjoin(" Option for the sampling of the BZ (brav):", itoa(ifc%brav)))
906  call wrtout(unt, sjoin(" Symmetrization flag (symdynmat):", itoa(ifc%symdynmat)))
907  call wrtout(unt, sjoin(" Dipole-dipole interaction flag (dipdip):", itoa(ifc%dipdip)))
908  call wrtout(unt, sjoin(" Dielectric tensor: ", ltoa(reshape(ifc%dielt, [9]), fmt="f10.2")))
909  call wrtout(unt, " Effective charges:")
910  do iatom=1,ifc%natom
911    call wrtout(unt, ltoa(reshape(ifc%zeff(:,:,iatom), [3*3]), fmt="f10.2"))
912  end do
913 
914  call wrtout(unt, sjoin(" Mass of the atoms (atomic mass unit): ", ltoa(ifc%amu)))
915  call wrtout(unt, sjoin(" Number of real-space points for IFC(R): ", itoa(ifc%nrpt)))
916  call wrtout(unt, " ")
917 
918  call wrtout(unt, " Q-mesh:")
919  call wrtout(unt, sjoin(" ngqpt:", ltoa(ifc%ngqpt),", nqshft:", itoa(ifc%nqshft)))
920  do ii=1,ifc%nqshft
921    call wrtout(unt, sjoin("  ", ktoa(ifc%qshft(:,ii))))
922  end do
923 
924 end subroutine ifc_print

m_ifc/ifc_printbxsf [ Functions ]

[ Top ] [ m_ifc ] [ Functions ]

NAME

 ifc_printbxsf

FUNCTION

  Output phonon isosurface in Xcrysden format.

INPUTS

  ifc<ifc_type>=Stores data related to interatomic force constants.
  crystal<crystal_t>=Info on the crystal structure
  ngqpt(3)=Divisions of the q-mesh
  nqshft=Number of shifts
  qshft(3,nqshft)=Shifts of the q-mesh.
  path=File name for output to disk
  comm=MPI communicator.

OUTPUT

  Only write to file

PARENTS

      eph

CHILDREN

      dfpt_phfrq,gtdyn9,nctk_defwrite_nonana_terms

SOURCE

2753 subroutine ifc_printbxsf(ifc, cryst, ngqpt, nqshft, qshft, path, comm)
2754 
2755 
2756 !This section has been created automatically by the script Abilint (TD).
2757 !Do not modify the following lines by hand.
2758 #undef ABI_FUNC
2759 #define ABI_FUNC 'ifc_printbxsf'
2760 !End of the abilint section
2761 
2762  implicit none
2763 
2764 !Arguments -------------------------------
2765 !scalars
2766  integer,intent(in) :: nqshft,comm
2767  character(len=*),intent(in) :: path
2768  type(ifc_type),intent(in) :: ifc
2769  type(crystal_t),intent(in) :: cryst
2770 !arrays
2771  integer,intent(in) :: ngqpt(3)
2772  real(dp),intent(in) :: qshft(3,nqshft)
2773 
2774 !Local variables -------------------------
2775 !scalars
2776  integer,parameter :: nsppol1=1,master=0,qptopt1=1
2777  integer :: my_rank,nprocs,iq_ibz,nqibz,nqbz,ierr
2778  character(len=500) :: msg
2779 !arrays
2780  integer :: qptrlatt(3,3),dummy_symafm(cryst%nsym)
2781  real(dp) :: displ_cart(2,3*cryst%natom,3*cryst%natom)
2782  real(dp),allocatable :: qibz(:,:),wtq(:),qbz(:,:),freqs_qibz(:,:)
2783 
2784 ! *********************************************************************
2785 
2786  my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
2787 
2788  ! Setup IBZ, weights and BZ. Always use q --> -q symmetry for phonons even in systems wo inversion
2789  qptrlatt = 0; qptrlatt(1,1) = ngqpt(1); qptrlatt(2,2) = ngqpt(2); qptrlatt(3,3) = ngqpt(3)
2790  call kpts_ibz_from_kptrlatt(cryst, qptrlatt, qptopt1, nqshft, qshft, nqibz, qibz, wtq, nqbz, qbz)
2791  ABI_FREE(qbz)
2792  ABI_FREE(wtq)
2793 
2794  ! Compute phonon frequencies in the irreducible wedge.
2795  ABI_CALLOC(freqs_qibz, (3*cryst%natom, nqibz))
2796 
2797  do iq_ibz=1,nqibz
2798    if (mod(iq_ibz, nprocs) /= my_rank) cycle ! mpi parallelism.
2799    call ifc_fourq(ifc, cryst, qibz(:,iq_ibz), freqs_qibz(:,iq_ibz), displ_cart)
2800  end do
2801  call xmpi_sum(freqs_qibz, comm, ierr)
2802 
2803  ! Output phonon isosurface.
2804  if (my_rank == master) then
2805    dummy_symafm = 1
2806    call printbxsf(freqs_qibz, zero, zero, cryst%gprimd, qptrlatt, 3*cryst%natom,&
2807      nqibz, qibz, cryst%nsym, .False., cryst%symrec, dummy_symafm, .True., nsppol1, qshft, nqshft, path, ierr)
2808    if (ierr /=0) then
2809      msg = "Cannot produce BXSF file with phonon isosurface, see log file for more info"
2810      MSG_WARNING(msg)
2811      call wrtout(ab_out, msg, 'COLL')
2812    end if
2813  end if
2814 
2815  ABI_FREE(freqs_qibz)
2816  ABI_FREE(qibz)
2817 
2818 end subroutine ifc_printbxsf

m_ifc/ifc_speedofsound [ Functions ]

[ Top ] [ m_ifc ] [ Functions ]

NAME

 ifc_speedofsound

FUNCTION

  Calculate the speed of sound by averaging the phonon group velocities of the
  three acoustic modes on a small sphere of radius qrad centered around Gamma.
  Perform spherical integration with Lebedev-Laikov grids

INPUTS

 ifc<ifc_type>=Object containing the dynamical matrix and the IFCs.
 crystal<crystal_t> = Information on the crystalline structure.
 qrad_tolkms(2):
   qrad=Radius of the sphere in reciprocal space
   atols_kms=Absolute tolerance in kilometer/second. The code generates spherical meshes
     until the results are converged twice within atols_kms.
 ncid=the id of the open NetCDF file. Use nctk_noid to disable netcdf output.
 comm=MPI communicator.

OUTPUT

PARENTS

      anaddb,m_gruneisen

CHILDREN

      dfpt_phfrq,gtdyn9,nctk_defwrite_nonana_terms

SOURCE

1209 subroutine ifc_speedofsound(ifc, crystal, qrad_tolkms, ncid, comm)
1210 
1211 
1212 !This section has been created automatically by the script Abilint (TD).
1213 !Do not modify the following lines by hand.
1214 #undef ABI_FUNC
1215 #define ABI_FUNC 'ifc_speedofsound'
1216 !End of the abilint section
1217 
1218  implicit none
1219 
1220 !Arguments -------------------------------
1221 !scalars
1222  integer,intent(in) :: comm,ncid
1223  type(ifc_type),intent(in) :: ifc
1224  type(crystal_t),intent(in) :: crystal
1225 !arrays
1226  real(dp),intent(in) :: qrad_tolkms(2)
1227 
1228 !Local variables -------------------------
1229 !scalars
1230  integer,parameter :: master=0
1231  integer :: ii,nu,igrid,my_rank,nprocs,ierr,converged,npts,num_negw,vs_ierr,ncerr
1232  integer :: iatom,iatref,num_acoustic,isacoustic
1233  real(dp) :: min_negw,cpu,wall,gflops
1234  real(dp) :: qrad,tolkms,diff
1235  character(len=500) :: msg
1236  type(lebedev_t) :: lgrid
1237 !arrays
1238  integer :: asnu(3)
1239  real(dp) :: qred(3),qvers_cart(3),qvers_red(3),quad(3),prev_quad(3),vs(7,3)
1240  real(dp) :: phfrqs(3*crystal%natom),dwdq(3,3*crystal%natom)
1241  real(dp) :: displ_cart(2,3*crystal%natom,3*crystal%natom),eigvec(2,3*crystal%natom,3*crystal%natom)
1242 
1243 ! *********************************************************************
1244 
1245  my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
1246 
1247  if (ifc%asr == 0) MSG_WARNING("Computing speed of sound with asr == 0! Use asr > 0")
1248  qrad = qrad_tolkms(1); tolkms = qrad_tolkms(2)
1249  ABI_CHECK(qrad > zero, "vs_qrad <= 0")
1250  ABI_CHECK(tolkms > zero, "vs_tolkms <= 0")
1251 
1252  call cwtime(cpu, wall, gflops, "start")
1253 
1254  ! Find the index of the first acoustic modes (needed to handle systems with unstable branches at Gamma
1255  ! In this case, indeed, we end up computing derivatives for the wrong branches, the integration becomes
1256  ! unstable and difficult to converge.
1257  qred = zero
1258  call ifc_fourq(ifc, crystal, qred, phfrqs, displ_cart,out_eigvec=eigvec)
1259  !write(std_out,*)"omega(q==Gamma): ",phfrqs
1260 
1261  num_acoustic = 0
1262  do nu = 1, 3*crystal%natom
1263 !  Check if this mode is acoustic like: scalar product of all displacement vectors are collinear
1264    isacoustic = 1
1265 !  Find reference atom with non-zero displacement
1266    do iatom=1,crystal%natom
1267      if(sum(displ_cart(:,(iatom-1)*3+1:(iatom-1)*3+3,nu)**2) >tol16)iatref=iatom
1268    enddo
1269 !  Now compute scalar product, and check they are all positive
1270    do iatom = 1, crystal%natom
1271      if (sum(eigvec(:,(iatom-1)*3+1:(iatom-1)*3+3, nu)*eigvec(:,(iatref-1)*3+1:(iatref-1)*3+3, nu)) < tol16 ) isacoustic = 0
1272    end do
1273    if (isacoustic == 1) then
1274      num_acoustic=num_acoustic+1
1275      asnu(num_acoustic)=nu
1276      if (num_acoustic==3) exit
1277    end if
1278  end do
1279 
1280  ABI_CHECK(num_acoustic == 3, sjoin("Wrong number of acoustic modes:", itoa(num_acoustic)))
1281 
1282  write(std_out,"(a,3i2,a)") "The bands with indices ",asnu(:)," will be used to calculate the sound velocities"
1283 
1284  ! Speed of sound along reduced directions.
1285  do ii=1,6
1286    qred = zero; qred(MOD(ii-1,3)+1) = one
1287    if (ii >= 4 .and. ii <= 6) qred = matmul(crystal%rprimd, qred) ! Cartesian directions.
1288    !if (ii >= 7 .and. ii <= 9) qred = matmul(crystal%rprimd, qred) ! Cartesian directions.
1289    qvers_red = (qred / normv(qred, crystal%gmet, "G"))
1290    qred = qrad * qvers_red
1291    !write(std_out,*)"dir",normv(qred, crystal%gmet, "G"), qrad
1292    call ifc_fourq(ifc, crystal, qred, phfrqs, displ_cart, dwdq=dwdq)
1293 
1294    do nu=1,3
1295      vs(ii, nu) = sqrt(sum(dwdq(1:3,asnu(nu)) ** 2)) * Bohr_meter * 0.001_dp / Time_Sec
1296    end do
1297    write(std_out,"(a,3es12.4,a)")" ||vs(nu)||:",vs(ii,:), " [km/s]"
1298 
1299    qvers_cart = matmul(crystal%gprimd, qvers_red) * two_pi
1300    do nu=1,3
1301      vs(ii, nu) = dot_product(dwdq(1:3,asnu(nu)), qvers_cart) * Bohr_meter * 0.001_dp / Time_Sec
1302    end do
1303    write(std_out,"(a,3es12.4,a)")" <q|vs(nu)>:",vs(ii,:), " [km/s]"
1304 
1305    !do nu=1,3
1306    !  write(std_out,"(a,3es12.4,a)")" vs(nu)_vect_red:",&
1307    !     matmul(crystal%gprimd, dwdq(1:3,asnu(nu))) * Bohr_meter * 0.001_dp / Time_Sec, " [km/s]"
1308    !end do
1309  end do
1310 
1311  ! Spherical average with Lebedev-Laikov grids.
1312  converged = 0
1313  do igrid=1,lebedev_ngrids
1314    lgrid = lebedev_new(igrid)
1315    npts = lgrid%npts; quad = zero; num_negw = 0; min_negw = zero
1316    do ii=1,npts
1317      if (mod(ii, nprocs) /= my_rank) cycle ! mpi-parallelism
1318 
1319      ! Build q-point on sphere of radius qrad. qcart --> qred
1320      qred = matmul(crystal%rprimd, lgrid%versors(:, ii))
1321      qred = qrad * (qred / normv(qred, crystal%gmet, "G"))
1322      !write(std_out,*)"lebe",normv(qred, crystal%gmet, "G"), qrad
1323      call ifc_fourq(ifc, crystal, qred, phfrqs, displ_cart, dwdq=dwdq)
1324      if (any(phfrqs(asnu) < -tol8)) then
1325        num_negw = num_negw + 1; min_negw = min(min_negw, minval(phfrqs(asnu)))
1326      end if
1327 
1328      do nu=1,3
1329        quad(nu) = quad(nu) + lgrid%weights(ii) * sqrt(sum(dwdq(1:3,asnu(nu)) ** 2))
1330        !quad(nu) = quad(nu) + lgrid%weights(ii) * abs(dot_product(lgrid%versors(:,ii), dwdq(:,asnu(nu))))
1331      end do
1332    end do
1333 
1334 !  Will use km/sec unit for echo purposes
1335    quad = quad * Bohr_meter * 0.001_dp / Time_Sec
1336    call xmpi_sum(quad, comm, ierr)
1337    call xmpi_sum(num_negw, comm, ierr)
1338    call lebedev_free(lgrid)
1339 
1340    write(std_out,'(2(a,i6),a,3es12.4,a,es12.4)') &
1341      " Lebedev-Laikov grid: ",igrid,", npts: ", npts, " vs_sphavg(ac_modes): ",quad, " <vs>: ",sum(quad)/3
1342 
1343    if (igrid > 1) then
1344      diff = zero
1345      do nu=1,3
1346        diff = diff + abs(quad(nu) - prev_quad(nu)) / 3
1347      end do
1348      !if (abs(sum(quad - prev_quad)/3) < tolkms) then
1349      if (diff < tolkms) then
1350         converged = converged + 1
1351      else
1352         converged = 0
1353      end if
1354    end if
1355    prev_quad = quad
1356    vs(7, :) = quad
1357    if (converged == 2) exit
1358  end do ! igrid
1359 
1360  if (my_rank == master) then
1361    ! vs_err: 1 if not converged, < 0 if negative freqs, == 0 if success.
1362    vs_ierr = 0
1363    do ii=1,3
1364      write(ab_out,"(a,3es12.4,a,i1)")" Speed of sound:",vs(ii,:)," [km/s] along reduced direction: ",ii
1365    end do
1366    write(ab_out,'(2(a,es12.4),a,i0)') &
1367      " Lebedev-Laikov integration with qradius: ", qrad, " tolkms: ",tolkms, " [km/s], npts: ", npts
1368    write(ab_out,"(a,3es12.4,a,es12.4)")" Spherical average:",vs(7,:)," [km/s], ",sum(vs(7,:))/3
1369    if (converged /= 2) then
1370      vs_ierr = 1
1371      write(msg,'(a,es12.4,a)')" WARNING: Results are not converged within: ",tolkms, " [km/s]"
1372      call wrtout(ab_out, msg)
1373      MSG_WARNING(msg)
1374    end if
1375    if (num_negw > 0) then
1376      vs_ierr = -num_negw
1377      write(msg,'(a,i0,a,es12.4,3a)') &
1378        " WARNING: Detected ",num_negw, " negative frequencies. Minimum was: ",min_negw * Ha_meV, "[meV]",ch10,&
1379        " Speed of sound could be wrong"
1380      call wrtout(ab_out, msg)
1381      MSG_WARNING(msg)
1382    end if
1383 
1384    ! Dump results to netcdf file.
1385    if (ncid /= nctk_noid) then
1386 #ifdef HAVE_NETCDF
1387      ncerr = nctk_def_arrays(ncid, [nctkarr_t("vsound", "dp", "seven, three")], defmode=.True.)
1388      NCF_CHECK(ncerr)
1389      ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "vsound_ierr"])
1390      NCF_CHECK(ncerr)
1391      ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "vsound_qrad", "vsound_tolkms"])
1392      NCF_CHECK(ncerr)
1393      ncerr = nctk_def_arrays(ncid, [nctkarr_t("asnu", "i", "three")], defmode=.True.)
1394      NCF_CHECK(ncerr)
1395      ! Write data.
1396      NCF_CHECK(nctk_set_datamode(ncid))
1397      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vsound_ierr"), vs_ierr))
1398      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vsound_qrad"), qrad))
1399      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vsound_tolkms"), tolkms))
1400      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vsound"), vs))
1401      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "asnu"), asnu))
1402 #endif
1403    end if
1404  end if
1405 
1406  call cwtime(cpu, wall, gflops, "stop")
1407  write(std_out,"(2(a,f6.2))")"ifc_speedofsound: cpu: ",cpu,", wall: ",wall
1408 
1409 end subroutine ifc_speedofsound

m_ifc/ifc_test_phinterp [ Functions ]

[ Top ] [ m_ifc ] [ Functions ]

NAME

 ifc_test_phinterp

INPUTS

  crystal<type(crystal_t)> = Information on the crystalline structure.
  ngqpt(3)=Divisions of the q-mesh used to produce the B-spline.
  nshiftq=Number of shifts in Q-mesh
  shiftq(3,nshiftq)=Shifts of the q-mesh.
  ords(3)=order of the spline for the three directions. ord(1) must be in [0, nqx] where
    nqx is the number of points along the x-axis.
  comm=MPI communicator

OUTPUT

  Only writing

PARENTS

      eph

CHILDREN

      dfpt_phfrq,gtdyn9,nctk_defwrite_nonana_terms

SOURCE

3313 subroutine ifc_test_phinterp(ifc, cryst, ngqpt, nshiftq, shiftq, ords, comm, test_dwdq)
3314 
3315  use m_bz_mesh,  only : kpath_t, kpath_new, kpath_free
3316 
3317 !This section has been created automatically by the script Abilint (TD).
3318 !Do not modify the following lines by hand.
3319 #undef ABI_FUNC
3320 #define ABI_FUNC 'ifc_test_phinterp'
3321 !End of the abilint section
3322 
3323  implicit none
3324 
3325 !arguments ------------------------------------
3326 !scalars
3327  integer,intent(in) :: comm,nshiftq
3328  type(ifc_type),intent(in) :: ifc
3329  type(crystal_t),intent(in) :: cryst
3330  logical,optional,intent(in) :: test_dwdq
3331 !arrays
3332  integer,intent(in) :: ngqpt(3),ords(3)
3333  real(dp),intent(in) :: shiftq(3,nshiftq)
3334 
3335 !local variables-------------------------------
3336 !scalars
3337  integer,parameter :: master=0
3338  integer :: iq,nq,natom3,my_rank,nprocs,ierr,nu
3339  real(dp) :: mare_bspl,mae_bspl,mare_skw,mae_skw,dq
3340  real(dp) :: cpu,wall,gflops,cpu_fourq,wall_fourq,gflops_fourq
3341  real(dp) :: cpu_bspl,wall_bspl,gflops_bspl,cpu_skw,wall_skw,gflops_skw
3342  logical :: do_dwdq
3343  type(phbspl_t) :: phbspl
3344  type(skw_t) :: skw
3345  type(kpath_t) :: qpath
3346 !arrays
3347  real(dp) :: phfrq(3*cryst%natom),ofreqs(3*cryst%natom),qpt(3)
3348  real(dp) :: adiff_mev(3*cryst%natom),rel_err(3*cryst%natom)
3349  real(dp) :: displ_cart(2,3,cryst%natom,3*cryst%natom),qvers_cart(3),qvers_red(3)
3350  real(dp) :: vals4(4),dwdq(3,3*cryst%natom)
3351  real(dp) :: q1(3),q2(3)
3352  real(dp) :: bounds(3,5)
3353  real(dp),allocatable :: winterp(:,:),wdata(:,:)
3354 
3355 ! *********************************************************************
3356 
3357  bounds = reshape([zero, zero, zero, half, zero, zero, zero, half, zero, zero, zero, zero, zero, zero, half], [3,5])
3358  qpath = kpath_new(bounds, cryst%gprimd, 10)
3359 
3360  my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
3361  natom3 = 3 * cryst%natom
3362 
3363  ! Test computation of group velocities
3364  do_dwdq = .False.; if (present(test_dwdq)) do_dwdq = test_dwdq
3365  if (do_dwdq) then
3366    call wrtout(std_out, "Testing computation of group velocities.")
3367    q1 = zero; q2 = [half, zero, zero]
3368    !q1 = -q2
3369    !q1 = [-0.1_dp, zero, zero]; q2 = [half, zero, zero]
3370    !q1 = [0.01_dp, zero, zero]; q2 = [half, half, half]
3371    !q2 = [0.01_dp, zero, zero]; q1 = [half, zero, zero]
3372    !q2 = [-0.1_dp, zero, zero]; q1 = [half, zero, zero]
3373    nq = 50
3374    dq = normv(q2 - q1, cryst%gmet, "G") / (nq - 1)
3375    qvers_red = (q2 -q1) / normv(q2 - q1, cryst%gmet, "G")
3376    qvers_cart = matmul(cryst%gprimd, qvers_red) * two_pi
3377    ABI_MALLOC(wdata, (natom3, nq))
3378    ABI_MALLOC(winterp, (natom3, nq))
3379 
3380    do iq=1,nq
3381       qpt = q1 + (iq - 1) * dq * qvers_red
3382       call ifc_fourq(ifc, cryst, qpt, phfrq, displ_cart, dwdq=dwdq)
3383       wdata(:,iq) = phfrq
3384       if (iq == 1) winterp(:,iq) = phfrq
3385       if (iq > 1) then
3386         do nu=1,natom3
3387           winterp(nu, iq) = winterp(nu, iq-1) + dq * dot_product(dwdq(:, nu), qvers_cart)
3388         end do
3389         !imax = maxloc(abs(phfrq - wnext)); ii = imax(1)
3390         !write(*,*)ii, phfrq(ii), wnext(ii), (phfrq(ii) - wnext(ii)) * Ha_meV
3391         !write(std_out,*)"wvel:", qpt, minval(abs(phfrq - wnext)) * Ha_meV, maxval(abs(phfrq - wnext)) * Ha_meV, imax
3392         !write(std_out,*)"dwdq:",dwdq
3393       end if
3394       write(100,*)wdata(:,iq) * Ha_meV
3395       write(101,*)winterp(:,iq) * Ha_meV
3396       !do nu=1,natom3
3397       !  wnext(nu) = phfrq(nu) + dq * dot_product(dwdq(:, nu), qvers_cart)
3398       !end do
3399    end do
3400    ABI_FREE(wdata)
3401    ABI_FREE(winterp)
3402    MSG_COMMENT("Test phonon group velocities done")
3403    return
3404  end if
3405 
3406  ! Build interpolator objects.
3407  phbspl = ifc_build_phbspl(ifc, cryst, ngqpt, nshiftq, shiftq, ords, comm)
3408  skw = ifc_build_skw(ifc, cryst, ngqpt, nshiftq, shiftq, comm)
3409 
3410  cpu_fourq = zero; wall_fourq = zero; gflops_fourq = zero
3411  cpu_bspl = zero; wall_bspl = zero; gflops_bspl = zero
3412  cpu_skw = zero; wall_skw = zero; gflops_skw = zero
3413  mae_bspl = zero; mare_bspl = zero
3414  mae_skw = zero; mare_skw = zero
3415 
3416  !nq = 100
3417  nq = qpath%npts
3418  do iq=1,nq
3419    !if (mod(iq, nprocs) /= my_rank) cycle ! mpi parallelism
3420    !call random_number(qred)
3421    !qred = qred + (0.1_dp * my_rank / nprocs)
3422    !call wrap2_pmhalf(qred, qpt, shift)
3423 
3424    if (my_rank /= master) cycle
3425    qpt = qpath%points(:,iq)
3426 
3427    ! Fourier interpolation
3428    call cwtime(cpu, wall, gflops, "start")
3429    call ifc_fourq(ifc, cryst, qpt, phfrq, displ_cart) !, dwdq=dwdq)
3430    call cwtime(cpu, wall, gflops, "stop")
3431    cpu_fourq = cpu_fourq + cpu; wall_fourq = wall_fourq + wall
3432 
3433    ! B-spline interpolation
3434    call cwtime(cpu, wall, gflops, "start")
3435    call phbspl_evalq(phbspl, qpt, ofreqs)
3436    call cwtime(cpu, wall, gflops, "stop")
3437    cpu_bspl = cpu_bspl + cpu; wall_bspl = wall_bspl + wall
3438 
3439    adiff_meV = abs(phfrq - ofreqs); rel_err = zero
3440    where (abs(phfrq) > tol16)
3441      rel_err = adiff_meV / abs(phfrq)
3442    end where
3443    rel_err = 100 * rel_err; adiff_meV = adiff_meV * Ha_meV
3444    mae_bspl = mae_bspl + sum(adiff_meV) / natom3
3445    mare_bspl = mare_bspl + sum(rel_err) / natom3
3446    if (my_rank == master) then
3447      !write(std_out,*)maxval(adiff_meV), maxval(rel_err), trim(ktoa(qpt))
3448      write(456, *)phfrq
3449      write(457, *)ofreqs
3450    end if
3451 
3452    ! SKW interpolation
3453    call cwtime(cpu, wall, gflops, "start")
3454    do nu=1,natom3
3455      call skw_eval_bks(skw, nu, qpt, 1, ofreqs(nu))
3456    end do
3457    call cwtime(cpu, wall, gflops, "stop")
3458    cpu_skw = cpu_skw + cpu; wall_skw = wall_skw + wall
3459 
3460    adiff_meV = abs(phfrq - ofreqs); rel_err = zero
3461    where (abs(phfrq) > tol16)
3462      rel_err = adiff_meV / abs(phfrq)
3463    end where
3464    rel_err = 100 * rel_err; adiff_meV = adiff_meV * Ha_meV
3465    mae_skw = mae_skw + sum(adiff_meV) / natom3
3466    mare_skw = mare_skw + sum(rel_err) / natom3
3467    if (my_rank == master) then
3468      !write(std_out,*)maxval(adiff_meV), maxval(rel_err), trim(ktoa(qpt))
3469      write(458, *)ofreqs
3470    endif
3471  end do
3472 
3473  vals4 = [mae_bspl, mare_bspl, mae_skw, mare_skw]
3474  call xmpi_sum(vals4, comm, ierr)
3475  mae_bspl = vals4(1); mare_bspl = vals4(2); mae_skw = vals4(3); mare_skw = vals4(4)
3476 
3477  mae_bspl = mae_bspl / nq; mare_bspl = mare_bspl / nq
3478  mae_skw = mae_skw / nq; mare_skw = mare_skw / nq
3479 
3480  if (my_rank == master) then
3481    write(std_out,"(2(a,f6.2),a,i0)")"B-spline MAE: ",mae_bspl," [meV], MARE: ",mare_bspl, "% with nqpt: ",nq
3482    write(std_out,"(2(a,f6.2),a,i0)")"SKW      MAE: ",mae_skw, " [meV], MARE: ",mare_skw, "% with nqpt: ",nq
3483    write(std_out,"(2(a,f6.2))")"fourq: cpu: ",cpu_fourq,", wall: ",wall_fourq
3484    write(std_out,"(2(a,f6.2))")"bspl:  cpu: ",cpu_bspl,", wall: ",wall_bspl
3485    write(std_out,"(2(a,f6.2))")"skw:   cpu: ",cpu_skw,", wall: ",wall_skw
3486  end if
3487 
3488  call phbspl_free(phbspl)
3489  call skw_free(skw)
3490  call kpath_free(qpath)
3491  MSG_COMMENT("ifc_test_phinterp done")
3492 
3493 end subroutine ifc_test_phinterp

m_ifc/ifc_type [ Types ]

[ Top ] [ m_ifc ] [ Types ]

NAME

 ifc_type

FUNCTION

  Contains the necessary data to interpolate the
  phonon bandstructure and eigenvectors in reciprocal space (ie.
  interatomic force constants and corresponding real space grid info).

SOURCE

 79  type,public :: ifc_type
 80 
 81    integer :: natom
 82      ! Number of atoms in the unit cell.
 83 
 84    integer :: mpert
 85      ! Maximum number of ipert.
 86 
 87    integer :: asr
 88      ! Option for the treatment of the Acoustic Sum Rule.
 89 
 90    integer :: brav
 91      ! Option for the sampling of the BZ (anaddb input variable)
 92 
 93    integer :: dipdip
 94      ! dipole dipole interaction flag.
 95 
 96    integer :: symdynmat
 97      ! If equal to 1, the dynamical matrix is symmetrized in dfpt_phfrq before the diagonalization.
 98 
 99    integer :: nqshft
100      ! Number of shifts in the q-mesh (usually 1 since the mesh is gamma-centered!)
101 
102    integer :: nqibz
103      ! Number of points in the IBZ
104 
105    integer :: nqbz
106      ! Number of points in the full BZ
107 
108    integer :: nrpt
109      ! Number of real space points used to integrate IFC (for interpolation of dynamical matrices)
110 
111    integer :: ngqpt(3)
112     ! Number of division in the Q mesh.
113 
114    real(dp) :: rprim(3,3),gprim(3,3),acell(3)
115      ! These values are used to call anaddb routines that don't use rprimd, gprimd
116 
117    real(dp) :: dielt(3,3)
118      ! Dielectric tensor (Cartesian coordinates)
119 
120    real(dp) :: omega_minmax(2)
121      ! Min and max frequency obtained on the initial ab-initio q-mesh (-+ 30 cmm1)
122      ! Used to generate frequency meshes for DOSes.
123 
124    real(dp) :: r_inscribed_sphere
125      ! radius of biggest sphere inscribed in the WS supercell
126 
127    real(dp),allocatable :: amu(:)
128      ! amu(ntypat)
129      ! mass of the atoms (atomic mass unit)
130 
131    real(dp),allocatable :: atmfrc(:,:,:,:,:)
132      ! atmfrc(3,natom,3,natom,nrpt)
133      ! Inter atomic forces in real space
134 
135    integer,allocatable :: cell(:,:)
136      ! cell(nrpt,3)
137      ! Give the index of the the cell and irpt
138 
139    real(dp),allocatable :: ewald_atmfrc(:,:,:,:,:)
140      ! Ewald_atmfrc(3,natom,3,natom,nrpt)
141      ! Ewald Inter atomic forces in real space
142 
143    real(dp),allocatable :: short_atmfrc(:,:,:,:,:)
144      ! short_atmfrc(3,natom,3,natom,nrpt)
145      ! Short range part of Inter atomic forces in real space
146 
147    real(dp),allocatable :: qshft(:,:)
148     ! qshft(3,nqshft)
149     ! The shifts of the q-mesh
150 
151    real(dp), allocatable :: rpt(:,:)
152      ! rpt(3,nrpt)
153      ! Real space points in canonical type coordinates.
154 
155    real(dp),allocatable :: wghatm(:,:,:)
156      ! wghatm(natom,natom,nrpt)
157      ! Weights for each point and atom in the Wigner Seitz supercell in real space.
158 
159    real(dp),allocatable :: rcan(:,:)
160      ! rcan(3,natom)
161      ! Atomic position in canonical coordinates.
162 
163    real(dp),allocatable :: trans(:,:)
164      ! trans(3,natom)
165      ! Atomic translations: xred = rcan + trans
166 
167    real(dp),allocatable :: dyewq0(:,:,:)
168      ! dyewq0(3,3,natom)
169      ! Atomic self-interaction correction to the dynamical matrix (only when dipdip=1).
170 
171    real(dp),allocatable :: zeff(:,:,:)
172      ! zeff(3,3,natom)
173      ! Effective charge on each atom, versus electric field and atomic displacement.
174      ! Cartesian coordinates
175 
176    real(dp),allocatable :: qibz(:,:)
177      ! qibz(3,nqibz))
178      ! List of q-points in the IBZ
179 
180    real(dp),allocatable :: wtq(:)
181      ! wtq(nqibz))
182      ! q-point Weights.
183 
184    real(dp),allocatable :: qbz(:,:)
185      ! qbz(3,nqbz))
186      ! List of q-points in the full BZ
187 
188    real(dp),allocatable :: dynmat(:,:,:,:,:,:)
189      ! dynmat(2,3,natom,3,natom,nqbz))
190      ! dynamical matrices relative to the q points of the B.Z. sampling
191      ! Note that the long-range dip-dip part has been removed if dipdip=1
192      ! Moreover the array is multiplied by a phase shift in mkifc9.
193 
194    !real(dp),allocatable :: dynmat_lr(:,:,:,:,:,:)
195     ! dynmat_lr(2,3,natom,3,natom,nqbz))
196     ! Long-range part of dynmat in q-space
197 
198  end type ifc_type
199 
200  public :: ifc_init          ! Constructor from DDB datatype
201  public :: ifc_init_fromFile ! Constructor from filename
202  public :: ifc_free          ! Release memory
203  public :: ifc_print         ! Print info on the object.
204  public :: ifc_fourq         ! Use Fourier interpolation to compute interpolated frequencies w(q) and eigenvectors e(q)
205  public :: ifc_speedofsound  ! Compute the speed of sound by averaging phonon group velocities.
206  public :: ifc_write         ! Print the ifc (output, netcdf and text file)
207  public :: ifc_outphbtrap    ! Print out phonon frequencies on regular grid for BoltzTrap code.
208  public :: ifc_printbxsf     ! Output phonon isosurface in Xcrysden format.
209  public :: ifc_calcnwrite_nana_terms   ! Compute phonons for q--> 0 with LO-TO

m_ifc/ifc_write [ Functions ]

[ Top ] [ m_ifc ] [ Functions ]

NAME

 ifc_write

FUNCTION

 Adds the real-space interatomic force constants to:
  the output file,
  a NetCDF file which is already open on ncid
  if prt_ifc==1, to the ifcinfo.out file
  to a TDEP file named outfile.forceconstants_ABINIT

INPUTS

 Ifc<type(ifc_type)>=Object containing the dynamical matrix and the IFCs.
 ifcana= 0 => no analysis of ifc ; 1 => full analysis
 atifc(natom) =  atifc(ia) equals 1 if the analysis of ifc has to be done for atom ia; otherwise 0.
 ifcout= Number of interatomic force constants written in the output file
 prt_ifc = flag to print out ifc information for dynamical matrix (AI2PS)
 ncid=the id of the open NetCDF file. Set to nctk_noid if netcdf output is not wanted.

OUTPUT

   written in the output file and in the NetCDF file

NOTES

 This routine should be executed by one processor only

 TODO:
  1) ifc_write should not have side effects

  2) the code is unreadable and horrible - 3/4 different file formats for the
  same stuff. We should make different subroutines, even if it duplicates some code

PARENTS

      anaddb

CHILDREN

      dfpt_phfrq,gtdyn9,nctk_defwrite_nonana_terms

SOURCE

1745 subroutine ifc_write(Ifc,ifcana,atifc,ifcout,prt_ifc,ncid)
1746 
1747 
1748 !This section has been created automatically by the script Abilint (TD).
1749 !Do not modify the following lines by hand.
1750 #undef ABI_FUNC
1751 #define ABI_FUNC 'ifc_write'
1752 !End of the abilint section
1753 
1754  implicit none
1755 
1756 !Arguments -------------------------------
1757 !scalars
1758  integer,intent(in) :: ifcout,ifcana,prt_ifc,ncid
1759  type(ifc_type),intent(inout) :: Ifc
1760 !arrays
1761  integer,intent(in) :: atifc(Ifc%natom)
1762 
1763 !Local variables -------------------------
1764 !scalars
1765  integer :: ia,ib,ii,ncerr,iatifc,ifcout1,mu,nu,iout, irpt
1766 ! unit number to print out ifc information for dynamical matrix (AI2PS)
1767  integer :: unit_ifc, unit_tdep
1768  real(dp) :: detdlt
1769  real(dp) :: maxdist_tdep
1770  character(len=500) :: message
1771  character(len=4) :: str1, str2
1772 !arrays
1773  integer,allocatable :: list(:),indngb(:)
1774  real(dp) :: invdlt(3,3),ra(3),xred(3),dielt(3,3)
1775  real(dp),allocatable :: dist(:,:,:),wkdist(:),rsiaf(:,:,:),sriaf(:,:,:),vect(:,:,:)
1776  real(dp),allocatable :: posngb(:,:)
1777  real(dp) :: gprimd(3,3),rprimd(3,3)
1778 
1779 ! *********************************************************************
1780 
1781  iout = ab_out
1782  dielt = ifc%dielt
1783 
1784  ! Compute the distances between atoms
1785  ABI_MALLOC(dist,(Ifc%natom,Ifc%natom,Ifc%nrpt))
1786  call dist9(Ifc%acell,dist,Ifc%gprim,Ifc%natom,Ifc%nrpt,Ifc%rcan,Ifc%rprim,Ifc%rpt)
1787  ! Now dist(ia,ib,irpt) contains the distance from atom ia to atom ib in unit cell irpt.
1788 
1789  ABI_MALLOC(list,(Ifc%natom*Ifc%nrpt))
1790  ABI_MALLOC(wkdist,(Ifc%natom*Ifc%nrpt))
1791 
1792  ! Calculating the inverse (transpose) of the dielectric tensor
1793  call matr3inv(dielt,invdlt)
1794 
1795  ! Calculating the determinant of the dielectric tensor
1796  detdlt=dielt(1,1)*dielt(2,2)*dielt(3,3)+dielt(1,3)*dielt(2,1)*&
1797 & dielt(3,2)+dielt(1,2)*dielt(2,3)*dielt(3,1)-dielt(1,3)*&
1798 & dielt(2,2)*dielt(3,1)-dielt(1,1)*dielt(2,3)*dielt(3,2)-&
1799 & dielt(1,2)*dielt(2,1)*dielt(3,3)
1800 
1801 ! echo to log file
1802  write(std_out,'(a)' )' ifc_write: analysis of interatomic force constants '
1803  call mkrdim(Ifc%acell,Ifc%rprim,rprimd)
1804  call matr3inv(rprimd,gprimd)
1805 
1806  if (iout > 0) then
1807    write(iout, '(/,a,/)' )' Analysis of interatomic force constants '
1808    if(Ifc%dipdip==1)then
1809      write(iout, '(a)' )' Are given : column(1-3), the total force constant'
1810      write(iout, '(a)' )'       then  column(4-6), the Ewald part'
1811      write(iout, '(a)' )'       then  column(7-9), the short-range part'
1812      write(iout, '(a)' )' Column 1, 4 and 7 are related to the displacement'
1813      write(iout, '(a)' )'       of the generic atom along x,               '
1814      write(iout, '(a)' )' column 2, 5 and 8 are related to the displacement'
1815      write(iout, '(a)' )'       of the generic atom along y,               '
1816      write(iout, '(a)' )' column 3, 6 and 9 are related to the displacement'
1817      write(iout, '(a)')'       of the generic atom along z.               '
1818    else if(Ifc%dipdip==0)then
1819      write(iout, '(a)' )' column 1 is related to the displacement'
1820      write(iout, '(a)' )'        of the generic atom along x,    '
1821      write(iout, '(a)' )' column 2 is related to the displacement'
1822      write(iout, '(a)' )'        of the generic atom along y,    '
1823      write(iout, '(a)' )' column 3 is related to the displacement'
1824      write(iout, '(a)' )'        of the generic atom along z,    '
1825    end if
1826  end if
1827 
1828  if (ifcout>Ifc%natom*Ifc%nrpt .or. ifcout == -1) then
1829    ifcout1=Ifc%natom*Ifc%nrpt
1830    write(message, '(3a,i0,a)' )&
1831 &   'The value of ifcout exceeds the number of atoms in the big box.', ch10, &
1832 &   'Output limited to ',Ifc%natom*Ifc%nrpt,' atoms.'
1833    MSG_WARNING(message)
1834  else
1835    ifcout1=ifcout
1836  end if
1837 
1838  ! set up file for real space ifc output, if required
1839  if (prt_ifc == 1) then
1840    if (open_file('ifcinfo.out', message, newunit=unit_ifc, status="replace") /= 0) then
1841      MSG_ERROR(message)
1842    end if
1843    write(iout, '(a,a)' )ch10,&
1844 &   '  NOTE: Open file ifcinfo.out, for the output of interatomic force constants. This is because prt_ifc==1. '
1845 
1846    if (open_file('outfile.forceconstants_ABINIT', message, newunit=unit_tdep, status="replace") /= 0) then
1847      MSG_ERROR(message)
1848    end if
1849    write(iout, '(a,a,a)' )ch10,&
1850 &   '  NOTE: Open file outfile.forceconstants_ABINIT, for the output of interatomic force',&
1851 &   ' constants in TDEP format. This is because prt_ifc==1. '
1852    ! Print necessary stuff for TDEP
1853    write(unit_tdep,"(1X,I10,15X,'How many atoms per unit cell')") Ifc%natom
1854 
1855    ! look at all pairs, find furthest one with weight 1
1856 !   do ia
1857 !Ifc%wghatm(ia,ib,irpt)
1858    maxdist_tdep = Ifc%r_inscribed_sphere !maxval(dist)*0.8_dp
1859    write(unit_tdep,"(1X,F20.15,5X,'Realspace cutoff (A)')") maxdist_tdep*Bohr_Ang
1860  end if
1861 
1862 #ifdef HAVE_NETCDF
1863  if (ncid /= nctk_noid) then
1864    ! initialize netcdf variables
1865    ncerr = nctk_def_dims(ncid, [nctkdim_t("natifc", SUM(atifc)), nctkdim_t("number_of_r_points_big_box", Ifc%nrpt), &
1866      nctkdim_t("number_of_atoms_big_box", Ifc%natom*Ifc%nrpt), nctkdim_t("ifcout", ifcout1)], defmode=.True.)
1867    NCF_CHECK(ncerr)
1868 
1869    ncerr = nctk_def_arrays(ncid, [&
1870      nctkarr_t('ifc_atoms_indices', "i", "natifc"),&
1871      nctkarr_t('ifc_neighbours_indices', "i", "ifcout, natifc"),&
1872      nctkarr_t('ifc_distances', "dp", "ifcout, natifc "),&
1873      nctkarr_t('ifc_matrix_cart_coord', "dp", "number_of_cartesian_directions,number_of_cartesian_directions, ifcout, natifc")])
1874    NCF_CHECK(ncerr)
1875 
1876    if (Ifc%dipdip==1) then
1877      ncerr = nctk_def_arrays(ncid, [&
1878        nctkarr_t('ifc_matrix_cart_coord_short_range', "dp", &
1879        "number_of_cartesian_directions, number_of_cartesian_directions, ifcout, natifc")])
1880      NCF_CHECK(ncerr)
1881    end if
1882 
1883    if (ifcana==1) then
1884      ncerr = nctk_def_arrays(ncid, [&
1885        nctkarr_t('ifc_local_vectors', "dp", "number_of_cartesian_directions, number_of_cartesian_directions, ifcout, natifc")])
1886      NCF_CHECK(ncerr)
1887    end if
1888 
1889    NCF_CHECK(nctk_set_datamode(ncid))
1890  end if
1891 #endif
1892 
1893  ABI_MALLOC(rsiaf,(3,3,ifcout1))
1894  ABI_MALLOC(sriaf,(3,3,ifcout1))
1895  ABI_MALLOC(vect,(3,3,ifcout1))
1896  ABI_MALLOC(indngb,(ifcout1))
1897  ABI_MALLOC(posngb,(3,ifcout1))
1898 
1899  iatifc=0
1900 
1901  ! BIG loop on all generic atoms
1902  do ia=1,Ifc%natom
1903 
1904    if(atifc(ia)==1)then
1905 
1906      iatifc=iatifc+1
1907 
1908      ! First transform canonical coordinates to reduced coordinates
1909      do ii=1,3
1910        xred(ii)=Ifc%gprim(1,ii)*Ifc%rcan(1,ia)+Ifc%gprim(2,ii)*Ifc%rcan(2,ia)+Ifc%gprim(3,ii)*Ifc%rcan(3,ia)
1911      end do
1912 
1913      ! Then to cartesian coordinates
1914      ra(:)=xred(1)*Ifc%acell(1)*Ifc%rprim(:,1)+ xred(2)*Ifc%acell(2)*Ifc%rprim(:,2)+ xred(3)*Ifc%acell(3)*Ifc%rprim(:,3)
1915 
1916      ! This sorting algorithm is slow ...
1917      wkdist(:)=reshape(dist(ia,:,:),(/Ifc%natom*Ifc%nrpt/))
1918      do ii=1,Ifc%natom*Ifc%nrpt
1919        list(ii)=ii
1920      end do
1921      call sort_dp(Ifc%natom*Ifc%nrpt,wkdist,list,tol14)
1922 
1923      if (iout > 0) then
1924        write(iout, '(a)' )
1925        write(std_out,'(a,i4)' )' generic atom number',ia
1926        write(iout, '(a,i4)' )' generic atom number',ia
1927        write(std_out,'(a,3es16.8)' ) ' with cartesian coordinates',ra(1:3)
1928        write(iout,'(a,3es16.8)' ) ' with cartesian coordinates',ra(1:3)
1929        write(iout, '(a)' )
1930      end if
1931 
1932      call ifc_getiaf(Ifc,ifcana,ifcout1,iout,ifc%zeff,ia,ra,list,dist,invdlt,&
1933 &                    detdlt,rsiaf,sriaf,vect,indngb,posngb)
1934 
1935 
1936      if (prt_ifc == 1) then
1937        do ii=1,ifcout1
1938          if (wkdist(ii) > maxdist_tdep) exit
1939        end do
1940        ii = ii - 1
1941        write(unit_tdep,"(1X,I10,15X,'How many neighbours does atom ',I3,' have')") ii, ia
1942 
1943        do ii=1,ifcout1
1944          ib = indngb(ii)
1945          irpt = (list(ii)-1)/Ifc%natom+1
1946          ! limit printing to maximum distance for tdep
1947          if (wkdist(ii) > maxdist_tdep) cycle
1948 
1949          !TDEP
1950          call int2char4(ii, str1)
1951          call int2char4(ia, str2)
1952          write(unit_tdep,"(1X,I10,15X,a,a,a,a)") ib, &
1953 &            'In the unit cell, what is the index of neighbour ', &
1954 &            trim(str1), " of atom ", trim(str2)
1955          ! The lattice vector needs to be in reduced coordinates?
1956          ! TODO: check whether this is correct for TDEP: might need just lattice
1957          ! vector part and not full vector, and could be in integers instead of
1958          ! cartesian vector...
1959          write (unit_tdep,'(3es28.16)') matmul(Ifc%rpt(1:3,irpt),Ifc%gprim)
1960 
1961          !AI2PS
1962          write(unit_ifc,'(i6,i6)') ia,ii
1963          write(unit_ifc,'(3es28.16)') posngb(1:3,ii)
1964          do nu=1,3
1965            !TDEp
1966            ! And the actual short ranged forceconstant: TODO: check if
1967            ! a transpose is needed or a swap between the nu and the mu
1968            !write(unit_tdep,'(3f28.16)') (sriaf(nu,mu,ii)*Ha_eV/amu_emass, mu=1, 3)
1969            write(unit_tdep,'(3f28.16)') (Ifc%short_atmfrc(mu,ia,nu,ib,irpt)*Ha_eV/Bohr_Ang**2, mu=1, 3)
1970 
1971            !AI2PS
1972            write(unit_ifc,'(3f28.16)')(rsiaf(nu,mu,ii),mu=1,3)
1973          end do
1974        end do
1975 
1976 #ifdef HAVE_NETCDF
1977        if (ncid /= nctk_noid) then
1978          NCF_CHECK(nf90_put_var(ncid, vid("ifc_atoms_indices"), ia, start=[iatifc]))
1979          NCF_CHECK(nf90_put_var(ncid, vid("ifc_neighbours_indices"), indngb, start=[1,iatifc], count=[ifcout1,1]))
1980          NCF_CHECK(nf90_put_var(ncid, vid("ifc_distances"), wkdist(:ifcout1), start=[1,iatifc],count=[ifcout1,1]))
1981          ncerr = nf90_put_var(ncid, vid("ifc_matrix_cart_coord"), rsiaf, start=[1,1,1,iatifc], count=[3,3,ifcout1,1])
1982          NCF_CHECK(ncerr)
1983          if (Ifc%dipdip==1) then
1984            ncerr = nf90_put_var(ncid, vid("ifc_matrix_cart_coord_short_range"), sriaf, &
1985              start=[1,1,1,iatifc], count=[3,3,ifcout1,1])
1986            NCF_CHECK(ncerr)
1987          end if
1988          if (ifcana==1) then
1989            ncerr = nf90_put_var(ncid, vid("ifc_local_vectors"), vect, start=[1,1,1,iatifc], count=[3,3,ifcout1,1])
1990            NCF_CHECK(ncerr)
1991          end if
1992        end if
1993 #endif
1994      end if
1995    end if ! End the condition on atifc
1996  end do ! End Big loop on atoms in the unit cell, and corresponding test
1997 
1998 
1999 ! NB for future use: in TDEP the following can also be provided.
2000 !        ! Print some auxiliary information, if it is there. Such as norm of
2001 !        ! forceconstant per shell, which shells there are and so on.
2002 !        if ( fc%npairshells .gt. 0 .and. allocated(fc%pairshell) ) then
2003 !            write(u,"(1X,I10,15X,'Number of irreducible coordination shells')") fc%npairshells
2004 !            do i=1,fc%npairshells
2005 !                write(u,"(1X,I10,1X,F16.10,1X,F16.10,15X,'number atoms in shell, radius, norm of forceconstant',I0)") fc%pairshell(i)%n,fc%pairshell(i)%rad,fc%pairshell(i)%norm,i
2006 !            enddo
2007 !            do i=1,fc%npairshells
2008 !                do j=1,fc%pairshell(i)%n
2009 !                    write(u,"(1X,3(1X,F18.12),2(1X,I0))") lo_chop(matmul(p%inv_latticevectors,fc%pairshell(i)%vec(:,j)),lo_sqtol),fc%pairshell(i)%atind(j),fc%pairshell(i)%pairind(j)
2010 !                enddo
2011 !            enddo
2012 !        endif
2013 
2014  if (prt_ifc == 1) then
2015    close(unit_ifc)
2016    close(unit_tdep)
2017 
2018    if (open_file('infile.lotosplitting_ABINIT', message, newunit=unit_tdep, status="replace") /= 0) then
2019      MSG_ERROR(message)
2020    end if
2021    write(unit_tdep,'(3es28.16)') dielt(:,1)
2022    write(unit_tdep,'(3es28.16)') dielt(:,2)
2023    write(unit_tdep,'(3es28.16)') dielt(:,3)
2024    do ia = 1, Ifc%natom
2025      do ii = 1, 3
2026        write(unit_tdep,'(3es28.16)') ifc%zeff(:,ii,ia)
2027      end do
2028    end do
2029    close(unit_tdep)
2030  end if
2031 
2032  ABI_FREE(rsiaf)
2033  ABI_FREE(sriaf)
2034  ABI_FREE(vect)
2035  ABI_FREE(indngb)
2036  ABI_FREE(posngb)
2037  ABI_FREE(dist)
2038  ABI_FREE(list)
2039  ABI_FREE(wkdist)
2040 
2041 #ifdef HAVE_NETCDF
2042 contains
2043  integer function vid(vname)
2044 
2045 
2046 !This section has been created automatically by the script Abilint (TD).
2047 !Do not modify the following lines by hand.
2048 #undef ABI_FUNC
2049 #define ABI_FUNC 'vid'
2050 !End of the abilint section
2051 
2052    character(len=*),intent(in) :: vname
2053    vid = nctk_idname(ncid, vname)
2054  end function vid
2055 #endif
2056 
2057 end subroutine ifc_write

m_ifc/omega_decomp [ Functions ]

[ Top ] [ m_ifc ] [ Functions ]

NAME

  omega_decomp

FUNCTION

 Compute and return the eigenvalues (frequencies) of the short-range and
 long-range part of the dynamical matrix  See Europhys. Lett. 33 p.713 (1996) for details.
 (included by U. Aschauer and EB)

INPUTS

  argin(sizein)=description

OUTPUT

  argout(sizeout)=description

SIDE EFFECTS

PARENTS

      m_ifc

CHILDREN

      dfpt_phfrq,gtdyn9,nctk_defwrite_nonana_terms

SOURCE

2427 subroutine omega_decomp(amu,natom,ntypat,typat,dynmatfl,dynmatsr,dynmatlr,iqpt,nqpt,eigenvec)
2428 
2429 
2430 !This section has been created automatically by the script Abilint (TD).
2431 !Do not modify the following lines by hand.
2432 #undef ABI_FUNC
2433 #define ABI_FUNC 'omega_decomp'
2434 !End of the abilint section
2435 
2436  implicit none
2437 
2438 !Arguments -------------------------------
2439 !scalars
2440  integer,intent(in) :: natom,ntypat
2441  integer,intent(in) :: iqpt,nqpt
2442 !arrays
2443  integer,intent(in) :: typat(natom)
2444  real(dp),intent(in) :: amu(ntypat)
2445  real(dp),intent(inout) :: dynmatfl(2,3,natom,3,natom,nqpt)
2446  real(dp),intent(inout) :: dynmatsr(2,3,natom,3,natom,nqpt)
2447  real(dp),intent(inout) :: dynmatlr(2,3,natom,3,natom,nqpt)
2448  real(dp),intent(in)    :: eigenvec(2*3*natom*3*natom)
2449 
2450 !Local variables -------------------------
2451 !scalars
2452  integer :: i1,i2,idir1,idir2,imode,ipert1,ipert2,index1,index2
2453  real(dp),parameter :: break_symm=1.0d-12
2454  real(dp) :: fac
2455 !arrays
2456  real(dp) :: nearidentity(3,3)
2457  real(dp) :: omegafl, omegasr, omegalr
2458  real(dp) :: sumfl,sumlr,sumsr,asr
2459 ! *********************************************************************
2460 
2461 !write(ab_out,*)''
2462 !write(std_out,*) 'SR/LR decomposition: enter for wavevector number :',iqpt
2463 
2464 !apply asr (note the lr part is already asred by construction in mkifc9)
2465  do ipert1=1,natom
2466    do idir1=1,3
2467      do idir2=1,3
2468        asr=0.0d0
2469        do ipert2=1,natom
2470          asr=asr+dynmatfl(1,idir1,ipert1,idir2,ipert2,iqpt)
2471        end do
2472        dynmatfl(1,idir1,ipert1,idir2,ipert1,iqpt)=dynmatfl(1,idir1,ipert1,idir2,ipert1,iqpt)-asr
2473        dynmatsr(1,idir1,ipert1,idir2,ipert1,iqpt)=dynmatsr(1,idir1,ipert1,idir2,ipert1,iqpt)-asr
2474      end do
2475    end do
2476  end do
2477 
2478 
2479 !This slight breaking of the symmetry allows the
2480 !results to be more portable between machines
2481  nearidentity(:,:)=1.0
2482  nearidentity(1,1)=1.0+break_symm
2483  nearidentity(3,3)=1.0-break_symm
2484 
2485 
2486 !Include Mass
2487  do ipert1=1,natom
2488    do ipert2=1,natom
2489 
2490      fac=1.0d0/sqrt(amu(typat(ipert1))*amu(typat(ipert2)))/amu_emass
2491 
2492      do idir1=1,3
2493        do idir2=1,3
2494 
2495          dynmatfl(1,idir1,ipert1,idir2,ipert2,iqpt)=&
2496 &         dynmatfl(1,idir1,ipert1,idir2,ipert2,iqpt)*&
2497 &         fac*nearidentity(idir1,idir2)
2498 
2499          dynmatsr(1,idir1,ipert1,idir2,ipert2,iqpt)=&
2500 &         dynmatsr(1,idir1,ipert1,idir2,ipert2,iqpt)*&
2501 &         fac*nearidentity(idir1,idir2)
2502 
2503          dynmatlr(1,idir1,ipert1,idir2,ipert2,iqpt)=&
2504 &         dynmatlr(1,idir1,ipert1,idir2,ipert2,iqpt)*&
2505 &         fac*nearidentity(idir1,idir2)
2506 
2507          ! This is to break slightly the translation invariance, and make
2508          ! the automatic tests more portable
2509          if(ipert1==ipert2 .and. idir1==idir2)then
2510            dynmatfl(1,idir1,ipert1,idir2,ipert2,iqpt)=&
2511 &           dynmatfl(1,idir1,ipert1,idir2,ipert2,iqpt)+&
2512 &           break_symm*natom/amu_emass/idir1*0.01d0
2513 
2514            dynmatsr(1,idir1,ipert1,idir2,ipert2,iqpt)=&
2515 &           dynmatsr(1,idir1,ipert1,idir2,ipert2,iqpt)+&
2516 &           break_symm*natom/amu_emass/idir1*0.01d0
2517 
2518            dynmatlr(1,idir1,ipert1,idir2,ipert2,iqpt)=&
2519 &           dynmatlr(1,idir1,ipert1,idir2,ipert2,iqpt)+&
2520 &           break_symm*natom/amu_emass/idir1*0.01d0
2521          end if
2522 
2523        end do
2524      end do
2525    end do
2526  end do
2527 
2528 !Calculation of <eigvec|Dyn_tot,Dyn_SR,Dyn_LR|eigenvec>=omega**2
2529 
2530 !write(ab_out,*)''
2531 !write(ab_out,*)'==============================================================================='
2532  write(ab_out,*)''
2533  write(ab_out,*) 'Long-Range/Short-Range decomposed phonon freq. (cm-1)**2'
2534  write(ab_out,*) 'at wavevector number:',iqpt
2535  write(ab_out,*)''
2536  write(ab_out,'(a13,1x,a16,2x,a16,2x,a16)') ' Mode number.','tot**2','SR**2','LR**2'
2537  write(std_out,'(a13,1x,a16,2x,a16,2x,a16)') ' Mode number.','tot**2','SR**2','LR**2'
2538 !write(ab_out,'(a12,2x,a10,2x,a10,2x,a10,2x,a16,2x,a16,2x,a16)') 'Mode number.','tot','SR','LR','tot**2','SR**2','LR**2'
2539 !write(std_out,'(a12,2x,a10,2x,a10,2x,a10,2x,a16,2x,a16,2x,a16)') 'Mode number.','tot','SR','LR','tot**2','SR**2','LR**2'
2540 
2541  do imode=1,3*natom
2542    sumfl=zero; sumlr=zero; sumsr=zero
2543 
2544    do ipert1=1,natom
2545      do ipert2=1,natom
2546        do i1=1,3
2547          do i2=1,3
2548 
2549            index1=i1+(ipert1-1)*3+3*natom*(imode-1)
2550            index2=i2+(ipert2-1)*3+3*natom*(imode-1)
2551            ! MG FIXME: I don't think these expressions are correct when q != 0
2552            ! We should also include the imaginary part
2553 
2554            sumfl = sumfl + eigenvec(2*index1-1) * dynmatfl(1,i1,ipert1,i2,ipert2,iqpt) * eigenvec(2*index2-1)
2555            sumlr = sumlr + eigenvec(2*index1-1) * dynmatlr(1,i1,ipert1,i2,ipert2,iqpt) * eigenvec(2*index2-1)
2556            sumsr = sumsr + eigenvec(2*index1-1) * dynmatsr(1,i1,ipert1,i2,ipert2,iqpt) * eigenvec(2*index2-1)
2557          end do
2558        end do
2559      end do
2560    end do
2561 
2562    sumfl = sumfl * Ha_cmm1 * Ha_cmm1
2563    sumsr = sumsr * Ha_cmm1 * Ha_cmm1
2564    sumlr = sumlr * Ha_cmm1 * Ha_cmm1
2565 
2566 !  Compute omega=sqrt(omega**2)
2567    if(sumfl>=1.0d-16)then
2568      omegafl=sqrt(sumfl)
2569    else if(sumfl>=-1.0d-16)then
2570      omegafl=zero
2571    else
2572      omegafl=-sqrt(-sumfl)
2573    end if
2574 
2575    if(sumsr>=1.0d-16)then
2576      omegasr=sqrt(sumsr)
2577    else if(sumsr>=-1.0d-16)then
2578      omegasr=zero
2579    else
2580      omegasr=-sqrt(-sumsr)
2581    end if
2582 
2583    if(sumlr>=1.0d-16)then
2584      omegalr=sqrt(sumlr)
2585    else if(sumlr>=-1.0d-16)then
2586      omegalr=zero
2587    else
2588      omegalr=-sqrt(-sumlr)
2589    end if
2590 
2591 !  Output
2592    write(ab_out,'(i4,10x,s,f16.4,2x,f16.4,2x,f16.4)') imode,sumfl,sumsr,sumlr  !vz_d
2593    write(std_out,'(i4,10x,s,f16.4,2x,f16.4,2x,f16.4)') imode,sumfl,sumsr,sumlr  !vz_d
2594 !  write(ab_out,'(i4,8x,f10.4,2x,f10.4,2x,f10.4,2x,s,f16.6,2x,f16.6,2x,f16.6)') imode,omegafl,omegasr,omegalr,sumfl,sumsr,sumlr
2595 !  write(std_out,'(i4,8x,f10.4,2x,f10.4,2x,f10.4,2x,s,f16.6,2x,f16.6,2x,f16.6)') imode,omegafl,omegasr,omegalr,sumfl,sumsr,sumlr
2596  end do
2597 
2598 end subroutine omega_decomp

m_ifc/phbspl_evalq [ Functions ]

[ Top ] [ m_ifc ] [ Functions ]

NAME

 phbspl_evalq

FUNCTION

  Interpolate phonon frequencies at an arbitrary q-point.

INPUTS

  qpt(3)=Q-point in reduced coordinate (will be wrapped in the interval [0,1[

OUTPUT

  ofreqs(%natom3)=Interpolated phonon frequencies.
    Note that ofreqs is not necessarily sorted in ascending order.
    The routine does not reorder the interpolated frequencies
    to be consistent with the interpolation of the derivatives.
  [oder1(3,%natom3)]=First order derivatives.

PARENTS

      m_ifc

CHILDREN

      dfpt_phfrq,gtdyn9,nctk_defwrite_nonana_terms

SOURCE

3050 subroutine phbspl_evalq(phbspl, qpt, ofreqs, oder1)
3051 
3052 
3053 !This section has been created automatically by the script Abilint (TD).
3054 !Do not modify the following lines by hand.
3055 #undef ABI_FUNC
3056 #define ABI_FUNC 'phbspl_evalq'
3057 !End of the abilint section
3058 
3059  implicit none
3060 
3061 !Arguments ------------------------------------
3062 !scalars
3063  type(phbspl_t),intent(in) :: phbspl
3064 !arrays
3065  real(dp),intent(in) :: qpt(3)
3066  real(dp),intent(out) :: ofreqs(phbspl%natom3)
3067  real(dp),optional,intent(out) :: oder1(3,phbspl%natom3)
3068 
3069 !Local variables-------------------------------
3070 !scalars
3071  integer :: nu,ii
3072 !arrays
3073  integer :: iders(3)!, iperm(phbspl%natom3)
3074  real(dp) :: qred(3),shift(3)
3075 
3076 ! *********************************************************************
3077 
3078  ! Wrap k-point in the interval [0,1[ where 1 is not included (tol12)
3079  ! This is required because the spline has been constructed in this region.
3080  call wrap2_zero_one(qpt, qred, shift)
3081 
3082  do nu=1,phbspl%natom3
3083    ! B-spline interpolation.
3084    ofreqs(nu) = dbs3vl(qred(1), qred(2), qred(3), phbspl%qxord, phbspl%qyord, phbspl%qzord,&
3085                        phbspl%xknot, phbspl%yknot, phbspl%zknot, phbspl%nqx, phbspl%nqy, phbspl%nqz,&
3086                        phbspl%coeff(nu)%vals)
3087  end do
3088 
3089  ! Sort frequencies.
3090  !iperm = [(nu, nu=1, phbspl%natom3)]; call sort_dp(phbspl%natom3, ofreqs, iperm, tol14)
3091 
3092  if (present(oder1)) then
3093    ! Compute first-order derivatives.
3094    do nu=1,phbspl%natom3
3095      do ii=1,3
3096        iders = 0; iders(ii) = 1
3097        oder1(ii,nu) = dbs3dr(iders(1), iders(2), iders(3), &
3098                              qred(1), qred(2), qred(3), phbspl%qxord, phbspl%qyord, phbspl%qzord,&
3099                              phbspl%xknot, phbspl%yknot, phbspl%zknot, phbspl%nqx, phbspl%nqy, phbspl%nqz,&
3100                              phbspl%coeff(nu)%vals)
3101      end do
3102    end do
3103 
3104  end if
3105 
3106 end subroutine phbspl_evalq

m_ifc/phbspl_free [ Functions ]

[ Top ] [ m_ifc ] [ Functions ]

NAME

 phbspl_free

FUNCTION

  Free dynamic memory.

PARENTS

      m_ifc

CHILDREN

      dfpt_phfrq,gtdyn9,nctk_defwrite_nonana_terms

SOURCE

3126 subroutine phbspl_free(phbspl)
3127 
3128 
3129 !This section has been created automatically by the script Abilint (TD).
3130 !Do not modify the following lines by hand.
3131 #undef ABI_FUNC
3132 #define ABI_FUNC 'phbspl_free'
3133 !End of the abilint section
3134 
3135  implicit none
3136 
3137 !Arguments ------------------------------------
3138 !scalars
3139  type(phbspl_t),intent(inout) :: phbspl
3140 
3141 !Local variables-------------------------------
3142 !scalars
3143  integer :: ii
3144 
3145 ! *********************************************************************
3146 
3147  if (allocated(phbspl%xknot)) then
3148    ABI_FREE(phbspl%xknot)
3149  end if
3150  if (allocated(phbspl%yknot)) then
3151    ABI_FREE(phbspl%yknot)
3152  end if
3153  if (allocated(phbspl%zknot)) then
3154    ABI_FREE(phbspl%zknot)
3155  end if
3156 
3157  ! Free B-spline coefficients.
3158  if (allocated(phbspl%coeff)) then
3159    do ii=1,size(phbspl%coeff, dim=1)
3160      if (allocated(phbspl%coeff(ii)%vals)) then
3161        ABI_FREE(phbspl%coeff(ii)%vals)
3162      end if
3163    end do
3164    ABI_DT_FREE(phbspl%coeff)
3165  end if
3166 
3167 end subroutine phbspl_free

m_ifc/phbspl_t [ Types ]

[ Top ] [ m_ifc ] [ Types ]

NAME

 phbspl_t

FUNCTION

SOURCE

222  type :: bcoeff_t
223    real(dp),allocatable :: vals(:,:,:)
224  end type bcoeff_t
225 
226  type,public :: phbspl_t
227 
228    integer :: natom3
229 
230    integer :: nqx,nqy,nqz
231    ! Number of input data points
232 
233    integer :: qxord,qyord,qzord
234    ! Order of the spline.
235 
236    !real(dp),allocatable :: xvec(:),yvec(:),zvec(:)
237    real(dp),allocatable :: xknot(:),yknot(:),zknot(:)
238    ! Array of length ndata+korder containing the knot
239 
240    type(bcoeff_t),allocatable :: coeff(:)
241    ! coeff(natom3)
242 
243  end type phbspl_t
244 
245  public :: ifc_build_phbspl     ! Build B-spline object.
246  public :: phbspl_evalq         ! Interpolate frequencies at arbitrary q-point.
247  public :: phbspl_free          ! Free memory.
248  public :: ifc_test_phinterp
249 
250 CONTAINS  !===========================================================