TABLE OF CONTENTS


ABINIT/m_dvdb [ Modules ]

[ Top ] [ Modules ]

NAME

  m_dvdb

FUNCTION

  Objects and methods to extract data from the DVDB file.
  The DVDB file is Fortran binary file with a collection of DFPT potentials
  associated to the different phonon perturbations (idir, ipert, qpt).
  DVDB files are produced with the `mrgdv` utility and are used in the EPH code
  to compute the matrix elements <k+q| dvscf_{idir, ipert, qpt} |k>.

COPYRIGHT

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

PARENTS

SOURCE

24 #if defined HAVE_CONFIG_H
25 #include "config.h"
26 #endif
27 
28 #include "abi_common.h"
29 
30 MODULE m_dvdb
31 
32  use defs_basis
33  use m_abicore
34  use m_errors
35  use m_xmpi
36  use m_distribfft
37  use m_nctk
38  use m_sort
39 #ifdef HAVE_NETCDF
40  use netcdf
41 #endif
42  use m_hdr
43 
44  use defs_abitypes,   only : hdr_type, mpi_type, datafiles_type
45  use m_fstrings,      only : strcat, sjoin, itoa, ktoa, ltoa, ftoa, yesno, endswith
46  use m_io_tools,      only : open_file, file_exists
47  use m_numeric_tools, only : wrap2_pmhalf, vdiff_eval, vdiff_print
48  use m_symtk,         only : mati3inv, littlegroup_q
49  use m_geometry,      only : littlegroup_pert, irreducible_set_pert
50  use m_copy,          only : alloc_copy
51  use m_mpinfo,        only : destroy_mpi_enreg, initmpi_seq
52  use m_fftcore,       only : ngfft_seq
53  use m_fft_mesh,      only : rotate_fft_mesh, times_eigr, times_eikr, ig2gfft, get_gftt, calc_ceikr, calc_eigr
54  use m_fft,           only : fourdp
55  use m_crystal,       only : crystal_t, crystal_free, crystal_print
56  use m_crystal_io,    only : crystal_from_hdr
57  use m_kpts,          only : kpts_ibz_from_kptrlatt, listkk
58  use m_spacepar,      only : symrhg, setsym
59  use m_fourier_interpol,only : fourier_interpol
60 
61  implicit none
62 
63  private

m_dvdb/calc_eiqr [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  calc_eiqr

FUNCTION

INPUTS

OUTPUT

PARENTS

      m_dvdb

CHILDREN

      cwtime,dvdb_get_v1scf_qpt,dvdb_get_v1scf_rpt,dvdb_open_read
      hdr_fort_read,hdr_fort_write,hdr_free,irreducible_set_pert
      kpts_ibz_from_kptrlatt,littlegroup_q,wrtout,xmpi_barrier

SOURCE

3559 subroutine calc_eiqr(qpt,nrpt,rpt,eiqr)
3560 
3561 
3562 !This section has been created automatically by the script Abilint (TD).
3563 !Do not modify the following lines by hand.
3564 #undef ABI_FUNC
3565 #define ABI_FUNC 'calc_eiqr'
3566 !End of the abilint section
3567 
3568  implicit none
3569 
3570 !Arguments -------------------------------
3571 !scalars
3572  integer,intent(in) :: nrpt
3573 !arrays
3574  real(dp),intent(in) :: qpt(3),rpt(3,nrpt)
3575  real(dp),intent(out) :: eiqr(2,nrpt)
3576 
3577 !Local variables -------------------------
3578 !scalars
3579  integer :: ir
3580  real(dp) :: qr
3581 
3582 ! *********************************************************************
3583 
3584  do ir=1,nrpt
3585    qr = two_pi * dot_product(qpt, rpt(:,ir))
3586    eiqr(1,ir) = cos(qr); eiqr(2,ir) = sin(qr)
3587  end do
3588 
3589 end subroutine calc_eiqr

m_dvdb/dvdb_check_fform [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  dvdb_check_fform

FUNCTION

  Check the value of fform. Return exit status and error message.

INPUTS

   fform=Value read from the header
   mode="merge_dvdb" to check the value of fform when we are merging POT1 files
        "read_dvdb" when we are reading POT1 files from a DVDB file.

OUTPUT

   errmsg=String with error message if ierr /= 0

PARENTS

CHILDREN

SOURCE

3615 integer function dvdb_check_fform(fform, mode, errmsg) result(ierr)
3616 
3617 
3618 !This section has been created automatically by the script Abilint (TD).
3619 !Do not modify the following lines by hand.
3620 #undef ABI_FUNC
3621 #define ABI_FUNC 'dvdb_check_fform'
3622 !End of the abilint section
3623 
3624  implicit none
3625 
3626 !Arguments ------------------------------------
3627  integer,intent(in) :: fform
3628  character(len=*),intent(in) :: mode
3629  character(len=*),intent(out) :: errmsg
3630 
3631 ! *************************************************************************
3632  ierr = 0
3633 
3634 ! Here I made a mistake because 102 corresponds to GS potentials
3635 ! as a consequence DVDB files generated with version <= 8.1.6
3636 ! contain list of potentials with fform = 102.
3637  !integer :: fform_pot=102
3638  !integer :: fform_pot=109
3639  !integer :: fform_pot=111
3640 
3641  if (fform == 0) then
3642    errmsg = "fform == 0! Either corrupted/nonexistent file or IO error"
3643    ierr = 42; return
3644  end if
3645 
3646  select case (mode)
3647  case ("merge_dvdb")
3648     if (all(fform /= [109, 111])) then
3649       errmsg = sjoin("fform:", itoa(fform), "is not supported in `merge_dvdb` mode")
3650       ierr = 1; return
3651     end if
3652 
3653  case ("read_dvdb")
3654     if (all(fform /= [102, 109, 111])) then
3655       errmsg = sjoin("fform:", itoa(fform), "is not supported in `read_dvdb` mode")
3656       ierr = 1; return
3657     end if
3658 
3659  case default
3660    errmsg = sjoin("Invalid mode:", mode)
3661    ierr = -1; return
3662  end select
3663 
3664 end function dvdb_check_fform

m_dvdb/dvdb_close [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  dvdb_close

FUNCTION

 Close the file

PARENTS

CHILDREN

SOURCE

581 subroutine dvdb_close(db)
582 
583 
584 !This section has been created automatically by the script Abilint (TD).
585 !Do not modify the following lines by hand.
586 #undef ABI_FUNC
587 #define ABI_FUNC 'dvdb_close'
588 !End of the abilint section
589 
590  implicit none
591 
592 !Arguments ------------------------------------
593 !scalars
594  type(dvdb_t),intent(inout) :: db
595 
596 !************************************************************************
597 
598  select case (db%iomode)
599  case (IO_MODE_FORTRAN)
600    close(db%fh)
601  case default
602    MSG_ERROR(sjoin("Unsupported iomode:", itoa(db%iomode)))
603  end select
604 
605  db%rw_mode = DVDB_NOMODE
606 
607 end subroutine dvdb_close

m_dvdb/dvdb_findq [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  dvdb_findq

FUNCTION

  Find the index of the q-point in db%qpts. Non zero umklapp vectors are not allowed.
  Returns -1 if not found.

INPUTS

  qpt(3)=q-point in reduced coordinates.
  [qtol]=Optional tolerance for q-point comparison.
         For each reduced direction the absolute difference between the coordinates must be less that qtol

PARENTS

CHILDREN

SOURCE

2927 integer pure function dvdb_findq(db, qpt, qtol) result(iqpt)
2928 
2929 
2930 !This section has been created automatically by the script Abilint (TD).
2931 !Do not modify the following lines by hand.
2932 #undef ABI_FUNC
2933 #define ABI_FUNC 'dvdb_findq'
2934 !End of the abilint section
2935 
2936  implicit none
2937 
2938 !Arguments ------------------------------------
2939 !scalars
2940  real(dp),optional,intent(in) :: qtol
2941  type(dvdb_t),intent(in) :: db
2942 !arrays
2943  real(dp),intent(in) :: qpt(3)
2944 
2945 !Local variables-------------------------------
2946 !scalars
2947  integer :: iq
2948  real(dp) :: my_qtol
2949 
2950 ! *************************************************************************
2951 
2952  my_qtol = 0.0001_dp; if (present(qtol)) my_qtol = qtol
2953 
2954  iqpt = -1
2955  do iq=1,db%nqpt
2956    if (all(abs(db%qpts(:, iq) - qpt) < my_qtol)) then
2957       iqpt = iq; exit
2958    end if
2959  end do
2960 
2961 end function dvdb_findq

m_dvdb/dvdb_free [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  dvdb_free

FUNCTION

 Close the file and release the memory allocated.

PARENTS

      eph,m_dvdb,mrgdv

CHILDREN

      cwtime,dvdb_get_v1scf_qpt,dvdb_get_v1scf_rpt,dvdb_open_read
      hdr_fort_read,hdr_fort_write,hdr_free,irreducible_set_pert
      kpts_ibz_from_kptrlatt,littlegroup_q,wrtout,xmpi_barrier

SOURCE

629 subroutine dvdb_free(db)
630 
631 
632 !This section has been created automatically by the script Abilint (TD).
633 !Do not modify the following lines by hand.
634 #undef ABI_FUNC
635 #define ABI_FUNC 'dvdb_free'
636 !End of the abilint section
637 
638  implicit none
639 
640 !Arguments ------------------------------------
641 !scalars
642  type(dvdb_t),intent(inout) :: db
643 
644 !************************************************************************
645 
646  ! integer arrays
647  if (allocated(db%pos_dpq)) then
648    ABI_FREE(db%pos_dpq)
649  end if
650  if (allocated(db%cplex_v1)) then
651    ABI_FREE(db%cplex_v1)
652  end if
653  if (allocated(db%symq_table)) then
654    ABI_FREE(db%symq_table)
655  end if
656  if (allocated(db%iv_pinfoq)) then
657    ABI_FREE(db%iv_pinfoq)
658  end if
659  if (allocated(db%ngfft3_v1)) then
660    ABI_FREE(db%ngfft3_v1)
661  end if
662 
663  ! real arrays
664  if (allocated(db%qpts)) then
665    ABI_FREE(db%qpts)
666  end if
667  if (allocated(db%rpt)) then
668    ABI_FREE(db%rpt)
669  end if
670  if (allocated(db%v1scf_rpt)) then
671    ABI_FREE(db%v1scf_rpt)
672  end if
673  if (allocated(db%rhog1_g0)) then
674    ABI_FREE(db%rhog1_g0)
675  end if
676  if (allocated(db%zeff)) then
677    ABI_FREE(db%zeff)
678  end if
679 
680  ! types
681  call crystal_free(db%cryst)
682  call destroy_mpi_enreg(db%mpi_enreg)
683 
684  ! Close the file but only if we have performed IO.
685  if (db%rw_mode == DVDB_NOMODE) return
686  call dvdb_close(db)
687 
688 end subroutine dvdb_free

m_dvdb/dvdb_ftinterp_qpt [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  dvdb_ftinterp_qpt

FUNCTION

  Fourier interpolation of potentials for a given q-point
  Internal tables must be prepared in advance by calling `dvdb_ftinterp_setup`.

INPUTS

  qpt(3)=q-point in reduced coordinates.
  nfft=Number of fft-points treated by this processors
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  comm=MPI communicator

OUTPUT

  ov1r(2*nfft, nspden, 3*natom)=Interpolated DFPT potentials at the given q-point.

PARENTS

      m_dvdb,m_phgamma,m_phpi,m_sigmaph

CHILDREN

      cwtime,dvdb_get_v1scf_qpt,dvdb_get_v1scf_rpt,dvdb_open_read
      hdr_fort_read,hdr_fort_write,hdr_free,irreducible_set_pert
      kpts_ibz_from_kptrlatt,littlegroup_q,wrtout,xmpi_barrier

SOURCE

2224 subroutine dvdb_ftinterp_qpt(db, qpt, nfft, ngfft, ov1r, comm)
2225 
2226 
2227 !This section has been created automatically by the script Abilint (TD).
2228 !Do not modify the following lines by hand.
2229 #undef ABI_FUNC
2230 #define ABI_FUNC 'dvdb_ftinterp_qpt'
2231 !End of the abilint section
2232 
2233  implicit none
2234 
2235 !Arguments ------------------------------------
2236 !scalars
2237  integer,intent(in) :: nfft,comm
2238  type(dvdb_t),intent(in) :: db
2239 !arrays
2240  integer,intent(in) :: ngfft(18)
2241  real(dp),intent(in) :: qpt(3)
2242  real(dp),intent(out) :: ov1r(2,nfft,db%nspden,db%natom3)
2243 
2244 !Local variables-------------------------------
2245 !scalars
2246  integer,parameter :: cplex2=2
2247  integer :: ir,ispden,ifft,mu,idir,ipert,timerev_q,nproc,my_rank,cnt,ierr
2248  real(dp) :: wr,wi
2249 !arrays
2250  integer :: symq(4,2,db%cryst%nsym)
2251  real(dp),allocatable :: eiqr(:,:), v1r_lr(:,:,:)
2252 
2253 ! *************************************************************************
2254 
2255  ABI_CHECK(allocated(db%v1scf_rpt), "v1scf_rpt is not allocated (call dvdb_ftinterp_setup)")
2256 
2257  my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
2258 
2259  ABI_MALLOC(v1r_lr, (2,nfft,db%natom3))
2260 
2261  ! Examine the symmetries of the q wavevector
2262  call littlegroup_q(db%cryst%nsym,qpt,symq,db%cryst%symrec,db%cryst%symafm,timerev_q,prtvol=db%prtvol)
2263 
2264  ! Compute FT phases for this q-point.
2265  ABI_MALLOC(eiqr, (2,db%nrpt))
2266  call calc_eiqr(qpt,db%nrpt,db%rpt,eiqr)
2267 
2268  ! Compute long-range part of the coupling potential
2269  v1r_lr = zero; cnt = 0
2270  if (db%has_dielt_zeff) then
2271    do mu=1,db%natom3
2272      cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI-parallelism
2273      idir = mod(mu-1, 3) + 1; ipert = (mu - idir) / 3 + 1
2274      call dvdb_v1r_long_range(db,qpt,ipert,idir,nfft,ngfft,v1r_lr(:,:,mu))
2275    end do
2276    call xmpi_sum(v1r_lr, comm, ierr)
2277  end if
2278 
2279  ! TODO: If high-symmetry q-points, one could save flops by FFT interpolating the independent
2280  ! perturbations and then rotate ...
2281  ov1r = zero; cnt = 0
2282  do mu=1,db%natom3
2283    idir = mod(mu-1, 3) + 1; ipert = (mu - idir) / 3 + 1
2284 
2285    do ispden=1,db%nspden
2286      do ifft=1,nfft
2287        cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI-parallelism
2288 
2289        do ir=1,db%nrpt
2290          wr = db%v1scf_rpt(1, ir,ifft,ispden,mu)
2291          wi = db%v1scf_rpt(2, ir,ifft,ispden,mu)
2292          ov1r(1,ifft,ispden,mu) = ov1r(1,ifft,ispden,mu) + wr*eiqr(1,ir) - wi * eiqr(2,ir)
2293          ov1r(2,ifft,ispden,mu) = ov1r(2,ifft,ispden,mu) + wr*eiqr(2,ir) + wi * eiqr(1,ir)
2294        end do
2295 
2296        ! Add the long-range part of the potential
2297        ov1r(1,ifft,ispden,mu) = ov1r(1,ifft,ispden,mu) + v1r_lr(1,ifft,mu)
2298        ov1r(2,ifft,ispden,mu) = ov1r(2,ifft,ispden,mu) + v1r_lr(2,ifft,mu)
2299 
2300      end do
2301 
2302      ! Remove the phase.
2303      call times_eikr(-qpt, ngfft, nfft, 1, ov1r(:,:,ispden,mu))
2304    end do
2305 
2306    ! Be careful with gamma and cplex!
2307    if (db%symv1) then
2308      call v1phq_symmetrize(db%cryst,idir,ipert,symq,ngfft,cplex2,nfft,db%nspden,db%nsppol,db%mpi_enreg,ov1r(:,:,:,mu))
2309    end if
2310  end do ! mu
2311 
2312  call xmpi_sum(ov1r, comm, ierr)
2313 
2314  ABI_FREE(eiqr)
2315  ABI_FREE(v1r_lr)
2316 
2317 
2318 end subroutine dvdb_ftinterp_qpt

m_dvdb/dvdb_ftinterp_setup [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  dvdb_ftinterp_setup

FUNCTION

  Prepare the internal tables used for the Fourier interpolation of the DFPT potentials.

INPUTS

  ngqpt(3)=Divisions of the ab-initio q-mesh.
  nqshift=Number of shifts used to generated the ab-initio q-mesh.
  qshift(3,nqshift)=The shifts of the ab-initio q-mesh.
  nfft=Number of fft-points treated by this processors
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  comm=MPI communicator

PARENTS

      m_dvdb,m_phgamma,m_sigmaph

CHILDREN

      cwtime,dvdb_get_v1scf_qpt,dvdb_get_v1scf_rpt,dvdb_open_read
      hdr_fort_read,hdr_fort_write,hdr_free,irreducible_set_pert
      kpts_ibz_from_kptrlatt,littlegroup_q,wrtout,xmpi_barrier

SOURCE

1852 subroutine dvdb_ftinterp_setup(db,ngqpt,nqshift,qshift,nfft,ngfft,comm,cryst_op)
1853 
1854 
1855 !This section has been created automatically by the script Abilint (TD).
1856 !Do not modify the following lines by hand.
1857 #undef ABI_FUNC
1858 #define ABI_FUNC 'dvdb_ftinterp_setup'
1859 !End of the abilint section
1860 
1861  implicit none
1862 
1863 !Arguments ------------------------------------
1864 !scalars
1865  integer,intent(in) :: nqshift,nfft,comm
1866  type(dvdb_t),target,intent(inout) :: db
1867 !arrays
1868  integer,intent(in) :: ngqpt(3),ngfft(18)
1869  real(dp),intent(in) :: qshift(3,nqshift)
1870  type(crystal_t),intent(in),target,optional :: cryst_op
1871 
1872 !Local variables-------------------------------
1873 !scalars
1874  integer,parameter :: sppoldbl1=1,timrev1=1,tim_fourdp0=0
1875  integer :: my_qptopt,iq_ibz,nqibz,iq_bz,nqbz
1876  integer :: ii,iq_dvdb,cplex_qibz,ispden,mu,irpt,idir,ipert
1877  integer :: iqst,nqst,itimrev,tsign,isym,ix,iy,iz,nq1,nq2,nq3,r1,r2,r3
1878  integer :: nproc,my_rank,ifft,cnt,ierr
1879  real(dp) :: dksqmax,phre,phim
1880  logical :: isirr_q, found
1881  type(crystal_t),pointer :: cryst
1882  type(mpi_type) :: mpi_enreg_seq
1883 !arrays
1884  integer :: qptrlatt(3,3),g0q(3),ngfft_qspace(18)
1885  integer,allocatable :: indqq(:,:),iperm(:),bz2ibz_sort(:),nqsts(:),iqs_dvdb(:)
1886  real(dp) :: qpt_bz(3),shift(3) !,qpt_ibz(3)
1887  real(dp),allocatable :: qibz(:,:),qbz(:,:),wtq(:),emiqr(:,:)
1888  real(dp),allocatable :: v1r_qibz(:,:,:,:),v1r_qbz(:,:,:,:), v1r_lr(:,:,:) !,all_v1qr(:,:,:,:)
1889 
1890 ! *************************************************************************
1891 
1892  if (allocated(db%v1scf_rpt)) then
1893    if (db%debug) then
1894      call wrtout(std_out, "v1scf_rpt is already computed. Returning")
1895    end if
1896    return
1897  end if
1898 
1899  nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
1900 
1901  if (present(cryst_op)) then
1902    cryst => cryst_op
1903  else
1904    cryst => db%cryst
1905  end if
1906  nq1 = ngqpt(1); nq2 = ngqpt(2); nq3 = ngqpt(3); my_qptopt = 1 !; if (present(qptopt)) my_qptopt = qptopt
1907 
1908  ! Generate the q-mesh by finding the IBZ and the corresponding weights.
1909  qptrlatt = 0
1910  do ii=1,3
1911    qptrlatt(ii,ii) = ngqpt(ii)
1912  end do
1913 
1914  ! Get IBZ and BZ.
1915  call kpts_ibz_from_kptrlatt(cryst, qptrlatt, my_qptopt, nqshift, qshift, &
1916    nqibz, qibz, wtq, nqbz, qbz) ! new_kptrlatt, new_shiftk)
1917 
1918  !write(std_out,*)"Irreducible q-points:"
1919  !do iq_ibz=1,nqibz; write(std_out,*)trim(ktoa(qibz(:,iq_ibz))),wtq(iq_ibz)*nqbz; end do
1920  ABI_CHECK(nqbz == product(ngqpt) * nqshift, "nqbz /= product(ngqpt) * nqshift")
1921 
1922  ! We want a gamma centered q-mesh for FFT.
1923  ABI_CHECK(nqshift == 1, "nshift > 1 not supported")
1924  ABI_CHECK(all(qshift(:,1) == zero), "qshift != 0 not supported")
1925 
1926  ABI_FREE(qbz)
1927  ABI_MALLOC(qbz, (3, nqbz))
1928  ii = 0
1929  do iz=0,nq3-1
1930    do iy=0,nq2-1
1931      do ix=0,nq1-1
1932        ii = ii + 1
1933        qbz(:, ii) = [ix/dble(nq1), iy/dble(nq2), iz/dble(nq3)]
1934        call wrap2_pmhalf([ix/dble(nq1), iy/dble(nq2), iz/dble(nq3)], qbz(:,ii), shift)
1935      end do
1936    end do
1937  end do
1938 
1939  ! Compute real-space points.
1940  ! Use the following indexing (N means ngfft of the adequate direction)
1941  ! 0 1 2 3 ... N/2    -(N-1)/2 ... -1    <= gc
1942  ! 1 2 3 4 ....N/2+1  N/2+2    ...  N    <= index ig
1943  db%nrpt = nqbz
1944  ABI_MALLOC(db%rpt, (3, db%nrpt))
1945  ii = 0
1946  do iz=1,nq3
1947    r3 = ig2gfft(iz,nq3) !* nq3
1948    do iy=1,nq2
1949      r2 = ig2gfft(iy,nq2) !* nq2
1950      do ix=1,nq1
1951        r1 = ig2gfft(ix,nq1) !* nq1
1952        ii = ii + 1
1953        db%rpt(:,ii) = [r1, r2, r3]
1954      end do
1955    end do
1956  end do
1957 
1958  ! Find correspondence BZ --> IBZ. Note:
1959  !   * q --> -q symmetry is always used for phonons.
1960  !   * we use symrec instead of symrel
1961  ABI_MALLOC(indqq, (nqbz*sppoldbl1,6))
1962  call listkk(dksqmax,cryst%gmet,indqq,qibz,qbz,nqibz,nqbz,cryst%nsym,&
1963    sppoldbl1,cryst%symafm,cryst%symrec,timrev1,use_symrec=.True.)
1964 
1965  if (dksqmax > tol12) then
1966    MSG_BUG("Something wrong in the generation of the q-points in the BZ! Cannot map BZ --> IBZ")
1967  end if
1968 
1969  ! Construct sorted mapping BZ --> IBZ to speedup qbz search below.
1970  ABI_MALLOC(iperm, (nqbz))
1971  ABI_MALLOC(bz2ibz_sort, (nqbz))
1972  iperm = [(ii, ii=1,nqbz)]
1973  bz2ibz_sort = indqq(:,1)
1974  call sort_int(nqbz,bz2ibz_sort,iperm)
1975 
1976  ! Reconstruct the IBZ according to what is present in the DVDB .
1977  ABI_MALLOC(nqsts, (nqibz))
1978  ABI_MALLOC(iqs_dvdb, (nqibz))
1979 
1980  ABI_MALLOC(v1r_lr, (2,nfft,db%natom3))
1981 
1982  iqst = 0
1983  do iq_ibz=1,nqibz
1984 
1985    ! In each q-point star, count the number of q-points
1986    ! and find the one present in DVDB.
1987    nqst = 0
1988    found = .false.
1989    do ii=iqst+1,nqbz
1990      if (bz2ibz_sort(ii) /= iq_ibz) exit
1991      nqst = nqst + 1
1992 
1993      iq_bz = iperm(ii)
1994      if (.not. found) then
1995        iq_dvdb = dvdb_findq(db, qbz(:,iq_bz))
1996        if (iq_dvdb /= -1) then
1997          qibz(:,iq_ibz) = qbz(:,iq_bz)
1998          iqs_dvdb(iq_ibz) = iq_dvdb
1999          found = .true.
2000        end if
2001      end if
2002    end do
2003 
2004    ! Check that nqst has been counted properly.
2005    ABI_CHECK(nqst > 0 .and. bz2ibz_sort(iqst+1) == iq_ibz, "Wrong iqst")
2006    if (abs(nqst - wtq(iq_ibz) * nqbz) > tol12) then
2007      write(std_out,*)nqst, wtq(iq_ibz) * nqbz
2008      MSG_ERROR("Error in counting q-point star or in the weights.")
2009    end if
2010 
2011    ! Check that the q-point has been found in DVDB.
2012    if (.not. found) then
2013      MSG_ERROR(sjoin("Cannot find symmetric q-point of:", ktoa(qibz(:,iq_ibz)), "in DVDB file"))
2014    end if
2015    !write(std_out,*)sjoin("qpt irred:",ktoa(qibz(:,iq_ibz)))
2016 
2017    iqst = iqst + nqst
2018    nqsts(iq_ibz) = nqst
2019 
2020  end do
2021 
2022  ! Redo the mapping with the new IBZ
2023  call listkk(dksqmax,cryst%gmet,indqq,qibz,qbz,nqibz,nqbz,cryst%nsym,&
2024    sppoldbl1,cryst%symafm,cryst%symrec,timrev1,use_symrec=.True.)
2025 
2026  ABI_MALLOC(emiqr, (2,db%nrpt))
2027  ABI_STAT_MALLOC(db%v1scf_rpt,(2,db%nrpt, nfft, db%nspden, db%natom3), ierr)
2028  ABI_CHECK(ierr==0, 'out of memory in db%v1scf_rpt')
2029  db%v1scf_rpt = zero
2030 
2031  ABI_MALLOC(v1r_qbz, (2, nfft, db%nspden, db%natom3))
2032  !v1r_qbz = huge(one)
2033 
2034  iqst = 0
2035  do iq_ibz=1,nqibz
2036 
2037    ! Get potentials for this IBZ q-point on the real-space FFT mesh.
2038    ! This call allocates v1r_qibz(cplex_qibz, nfft, nspden, 3*natom)
2039    call dvdb_readsym_allv1(db, iqs_dvdb(iq_ibz), cplex_qibz, nfft, ngfft, v1r_qibz, comm)
2040 
2041    ! Reconstruct by symmetry the potentials for the star of this q-point, perform FT and accumulate
2042    ! Be careful with the gamma point.
2043    do ii=1,nqsts(iq_ibz)
2044      iqst = iqst + 1
2045      iq_bz = iperm(iqst)
2046      ABI_CHECK(iq_ibz == indqq(iq_bz,1), "iq_ibz !/ ind qq(1)")
2047      isym = indqq(iq_bz,2); itimrev = indqq(iq_bz,6) + 1; g0q = indqq(iq_bz,3:5) ! IS(q_ibz) + g0q = q_bz
2048      tsign = 3-2*itimrev
2049 
2050      qpt_bz = qbz(:, iq_bz)
2051      !write(std_out,*)"  treating:",trim(ktoa(qpt_bz))
2052      isirr_q = (isym == 1 .and. itimrev == 1 .and. all(g0q == 0))
2053      !ABI_CHECK(all(g0q == 0), "g0q /= 0")
2054 
2055      ! Compute long-range part of the coupling potential
2056      v1r_lr = zero; cnt = 0
2057      if (db%has_dielt_zeff) then
2058        do mu=1,db%natom3
2059          cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI parallelism.
2060          idir = mod(mu-1, 3) + 1; ipert = (mu - idir) / 3 + 1
2061          call dvdb_v1r_long_range(db,qpt_bz,ipert,idir,nfft,ngfft,v1r_lr(:,:,mu))
2062        end do
2063        call xmpi_sum(v1r_lr, comm, ierr)
2064      end if
2065 
2066      if (cplex_qibz == 1) then
2067        ! Gamma point.
2068        ABI_CHECK(nqsts(iq_ibz) == 1, "cplex_qibz == 1 and nq nqst /= 1 (should be gamma)")
2069        ABI_CHECK(all(g0q == 0), "gamma point with g0q /= 0")
2070 
2071        ! Substract the long-range part of the potential
2072        do mu=1,db%natom3
2073          do ispden=1,db%nspden
2074            v1r_qibz(1, :, ispden, mu) = v1r_qibz(1, :, ispden, mu) - v1r_lr(1, :, mu)
2075          end do
2076        end do
2077 
2078        ! Slow FT.
2079        cnt = 0
2080        do mu=1,db%natom3
2081          do ispden=1,db%nspden
2082            do irpt=1,db%nrpt
2083              ! MPI-parallelism
2084              cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle
2085              do ifft=1,nfft
2086                db%v1scf_rpt(1, irpt, ifft, ispden, mu) = db%v1scf_rpt(1, irpt, ifft, ispden, mu) + &
2087                                                          v1r_qibz(1, ifft, ispden, mu)
2088              end do
2089            end do
2090          end do
2091        end do
2092 
2093      else
2094        ! Get the periodic part of the potential in BZ (v1r_qbz)
2095        if (isirr_q) then
2096          !write(std_out,*)sjoin("qpt irred:",ktoa(qpt_bz))
2097          v1r_qbz = v1r_qibz
2098        else
2099          call v1phq_rotate(cryst, qibz(:,iq_ibz), isym, itimrev, g0q,&
2100            ngfft, cplex_qibz, nfft, db%nspden, db%nsppol, db%mpi_enreg, v1r_qibz, v1r_qbz)
2101          !v1r_qbz = zero; v1r_qbz = v1r_qibz
2102 
2103          !call times_eigr(-tsign * g0q, ngfft, nfft, db%nspden*db%natom3, v1r_qbz)
2104          !call times_eigr(+tsign * g0q, ngfft, nfft, db%nspden*db%natom3, v1r_qbz)
2105          !if (itimrev == 2) v1r_qbz(2,:,:,:) = -v1r_qbz(2,:,:,:)
2106          !call times_eigr(-tsign * g0q, ngfft, nfft, db%nspden*db%natom3, v1r_qbz)
2107          !call times_eigr(+tsign * g0q, ngfft, nfft, db%nspden*db%natom3, v1r_qbz)
2108        end if
2109 
2110        ! Multiply by e^{iqpt_bz.r}
2111        call times_eikr(qpt_bz, ngfft, nfft, db%nspden*db%natom3, v1r_qbz)
2112 
2113        ! Substract the long-range part of the potential
2114        do mu=1,db%natom3
2115          do ispden=1,db%nspden
2116            v1r_qbz(1, :, ispden, mu) = v1r_qbz(1, :, ispden, mu) - v1r_lr(1, :, mu)
2117            v1r_qbz(2, :, ispden, mu) = v1r_qbz(2, :, ispden, mu) - v1r_lr(2, :, mu)
2118          end do
2119        end do
2120 
2121        ! Compute FT phases for this qpt_bz.
2122        call calc_eiqr(-qpt_bz,db%nrpt,db%rpt,emiqr)
2123 
2124        ! Slow FT.
2125        cnt = 0
2126        do mu=1,db%natom3
2127 
2128          do ispden=1,db%nspden
2129            do irpt=1,db%nrpt
2130              cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI parallelism.
2131              phre = emiqr(1,irpt); phim = emiqr(2,irpt)
2132 
2133              do ifft=1,nfft
2134                db%v1scf_rpt(1,irpt, ifft, ispden, mu) = db%v1scf_rpt(1, irpt, ifft, ispden, mu) + &
2135                   phre * v1r_qbz(1, ifft, ispden, mu) - phim * v1r_qbz(2, ifft, ispden, mu)
2136 
2137                db%v1scf_rpt(2, irpt, ifft, ispden, mu) = db%v1scf_rpt(2, irpt, ifft, ispden, mu) + &
2138                   phre * v1r_qbz(2, ifft, ispden, mu) + phim * v1r_qbz(1, ifft, ispden, mu)
2139              end do
2140 
2141            end do
2142          end do
2143        end do
2144      end if
2145 
2146    end do ! iqst
2147 
2148    ABI_FREE(v1r_qibz)
2149  end do ! iq_ibz
2150  ABI_CHECK(iqst == nqbz, "iqst /= nqbz")
2151 
2152  db%v1scf_rpt = db%v1scf_rpt / nqbz
2153  call xmpi_sum(db%v1scf_rpt, comm, ierr)
2154 
2155  ! Build mpi_type for executing fourdp in sequential.
2156  call ngfft_seq(ngfft_qspace, [nq1, nq2, nq3])
2157  call initmpi_seq(mpi_enreg_seq)
2158  call init_distribfft_seq(mpi_enreg_seq%distribfft,'c',ngfft_qspace(2),ngfft_qspace(3),'all')
2159  call init_distribfft_seq(mpi_enreg_seq%distribfft,'f',ngfft_qspace(2),ngfft_qspace(3),'all')
2160 
2161  !ABI_STAT_MALLOC(all_v1qr, (nqbz, 2, nfft, db%nspden, db%natom3), ierr)
2162  !ABI_CHECK(ierr==0, "oom in all_v1qr")
2163  !all_v1qr = zero
2164 
2165  cnt = 0
2166  do mu=1,db%natom3
2167    do ispden=1,db%nspden
2168      do ifft=1,nfft
2169        !cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle
2170        !call fourdp(1,all_v1qr(:,:,ifft,ispden,mu),db%v1scf_rpt(:,ifft,ispden,mu),+1,&
2171        ! mpi_enreg_seq,nqbz,ngfft_qspace,paral_kgb0,tim_fourdp0)
2172      end do
2173    end do
2174  end do
2175 
2176  !ABI_FREE(all_v1qr)
2177  call destroy_mpi_enreg(mpi_enreg_seq)
2178 
2179  ABI_FREE(emiqr)
2180  ABI_FREE(qibz)
2181  ABI_FREE(wtq)
2182  ABI_FREE(qbz)
2183  ABI_FREE(indqq)
2184  ABI_FREE(iperm)
2185  ABI_FREE(bz2ibz_sort)
2186  ABI_FREE(iqs_dvdb)
2187  ABI_FREE(nqsts)
2188  ABI_FREE(v1r_qbz)
2189  ABI_FREE(v1r_lr)
2190 
2191 
2192 end subroutine dvdb_ftinterp_setup

m_dvdb/dvdb_get_pinfo [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  dvdb_get_pinfo

FUNCTION

  Return information on the perturbations available for a given q-point index.

INPUTS

  iqpt=Index of the q-point

OUTPUT

  nperts=Number of perturbations found.
  cplex=2 if potentials are complex, 1 for real
  pinfo(3,3*db%mpert)=Array with info on the perturbations present on file
     pinfo(1, ip) gives the `idir` index of the ip-th perturbation.
     pinfo(2, ip) gives the `ipert` index of the ip-th perturbation.
     pinfo(3, ip) gives `pertcase`=idir + (ipert-1)*3

PARENTS

CHILDREN

SOURCE

824 integer function dvdb_get_pinfo(db, iqpt, cplex, pinfo) result(nperts)
825 
826 
827 !This section has been created automatically by the script Abilint (TD).
828 !Do not modify the following lines by hand.
829 #undef ABI_FUNC
830 #define ABI_FUNC 'dvdb_get_pinfo'
831 !End of the abilint section
832 
833  implicit none
834 
835 !Arguments ------------------------------------
836 !scalars
837  type(dvdb_t),intent(in) :: db
838  integer,intent(in) :: iqpt
839  integer,intent(out) :: cplex
840 !arrays
841  integer,intent(out) :: pinfo(3,3*db%mpert)
842 
843 !Local variables-------------------------------
844 !scalars
845  integer :: idir,ipert,iv1
846 
847 ! *************************************************************************
848 
849  ! Get the number of perturbations computed for this iqpt
850  pinfo = 0; cplex = 0; nperts = 0
851  do ipert=1,db%natom ! selects atomic perturbations only.
852     do idir=1,3
853       iv1 = db%pos_dpq(idir,ipert,iqpt)
854       if (iv1 /= 0) then
855         nperts = nperts + 1
856         pinfo(:, nperts) = [idir, ipert, idir + (ipert-1)*3]
857         if (cplex == 0) cplex = db%cplex_v1(iv1)
858         ABI_CHECK(cplex == db%cplex_v1(iv1), "cplex should be constant for given q!")
859       end if
860     end do
861  end do
862 
863 end function dvdb_get_pinfo

m_dvdb/dvdb_get_v1scf_qpt [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  dvdb_get_v1scf_qpt

FUNCTION

  Fourier interpolation of potentials for a given q-point
  This routine is meant to replace dvdb_ftinterp_qpt
  by performing the interpolation one perturbation at a time.

INPUTS

  qpt(3)=q-point in reduced coordinates.
  nfft=Number of fft-points treated by this processors
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  nrpt=Number of R-points = number of q-points in the full BZ
  nspden=Number of spin densities.
  ipert=index of the perturbation to be treated [1,natom3]
  v1scf_rpt(2,nrpt,nfft,nspden)=phonon perturbation potential
                                in real space lattice representation.
  comm=MPI communicator

OUTPUT

  v1scf_qpt(2*nfft, nspden)=Interpolated DFPT potentials at the given q-point.

PARENTS

      m_dvdb

CHILDREN

      cwtime,dvdb_get_v1scf_qpt,dvdb_get_v1scf_rpt,dvdb_open_read
      hdr_fort_read,hdr_fort_write,hdr_free,irreducible_set_pert
      kpts_ibz_from_kptrlatt,littlegroup_q,wrtout,xmpi_barrier

SOURCE

2715 subroutine dvdb_get_v1scf_qpt(db, cryst, qpt, nfft, ngfft, nrpt, nspden, &
2716 &                             ipert, v1scf_rpt, v1scf_qpt, comm)
2717 
2718 
2719 !This section has been created automatically by the script Abilint (TD).
2720 !Do not modify the following lines by hand.
2721 #undef ABI_FUNC
2722 #define ABI_FUNC 'dvdb_get_v1scf_qpt'
2723 !End of the abilint section
2724 
2725  implicit none
2726 
2727 !Arguments ------------------------------------
2728 !scalars
2729  integer,intent(in) :: nfft,nrpt,nspden,ipert,comm
2730  type(dvdb_t),intent(in) :: db
2731  type(crystal_t),intent(in) :: cryst
2732 !arrays
2733  integer,intent(in) :: ngfft(18)
2734  real(dp),intent(in) :: qpt(3)
2735  real(dp),intent(in) :: v1scf_rpt(2,nrpt,nfft,db%nspden)
2736  real(dp),intent(out) :: v1scf_qpt(2,nfft,db%nspden)
2737 
2738 !Local variables-------------------------------
2739 !scalars
2740  integer,parameter :: cplex2=2
2741  integer :: ir,ispden,ifft,idir,iat,timerev_q,nproc,my_rank,cnt,ierr
2742  real(dp) :: wr,wi
2743 !arrays
2744  integer :: symq(4,2,db%cryst%nsym)
2745  real(dp),allocatable :: eiqr(:,:), v1r_lr(:,:)
2746 
2747 ! *************************************************************************
2748  ABI_UNUSED(cryst%natom)
2749  ABI_UNUSED(nspden)
2750 
2751  my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
2752 
2753  ABI_MALLOC(v1r_lr, (2,nfft))
2754 
2755  ! Examine the symmetries of the q wavevector
2756  call littlegroup_q(db%cryst%nsym,qpt,symq,db%cryst%symrec,db%cryst%symafm,timerev_q,prtvol=db%prtvol)
2757 
2758  ! Compute FT phases for this q-point.
2759  ABI_MALLOC(eiqr, (2,db%nrpt))
2760  call calc_eiqr(qpt,db%nrpt,db%rpt,eiqr)
2761 
2762  idir = mod(ipert-1, 3) + 1; iat = (ipert - idir) / 3 + 1
2763 
2764  ! Compute long-range part of the coupling potential
2765  v1r_lr = zero; cnt = 0
2766  if (db%has_dielt_zeff) then
2767    call dvdb_v1r_long_range(db,qpt,iat,idir,nfft,ngfft,v1r_lr)
2768  end if
2769 
2770  ! TODO: If high-symmetry q-points, one could save flops by FFT interpolating the independent
2771  ! perturbations and then rotate ...
2772  v1scf_qpt = zero; cnt = 0
2773  do ispden=1,db%nspden
2774    do ifft=1,nfft
2775      cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI-parallelism
2776 
2777      do ir=1,db%nrpt
2778        wr = v1scf_rpt(1, ir,ifft,ispden)
2779        wi = v1scf_rpt(2, ir,ifft,ispden)
2780        v1scf_qpt(1,ifft,ispden) = v1scf_qpt(1,ifft,ispden) + wr*eiqr(1,ir) - wi * eiqr(2,ir)
2781        v1scf_qpt(2,ifft,ispden) = v1scf_qpt(2,ifft,ispden) + wr*eiqr(2,ir) + wi * eiqr(1,ir)
2782      end do
2783 
2784      ! Add the long-range part of the potential
2785      v1scf_qpt(1,ifft,ispden) = v1scf_qpt(1,ifft,ispden) + v1r_lr(1,ifft)
2786      v1scf_qpt(2,ifft,ispden) = v1scf_qpt(2,ifft,ispden) + v1r_lr(2,ifft)
2787 
2788    end do
2789 
2790    ! Remove the phase.
2791    call times_eikr(-qpt, ngfft, nfft, 1, v1scf_qpt(:,:,ispden))
2792  end do
2793 
2794  ! Be careful with gamma and cplex!
2795  if (db%symv1) then
2796    call v1phq_symmetrize(db%cryst,idir,iat,symq,ngfft,cplex2,nfft,db%nspden,db%nsppol,db%mpi_enreg,v1scf_qpt(:,:,:))
2797  end if
2798 
2799  call xmpi_sum(v1scf_qpt, comm, ierr)
2800 
2801  ABI_FREE(eiqr)
2802  ABI_FREE(v1r_lr)
2803 
2804 
2805 end subroutine dvdb_get_v1scf_qpt

m_dvdb/dvdb_get_v1scf_rpt [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  dvdb_get_v1scf_rpt

FUNCTION

  Compute the phonon perturbation potential in real space lattice
  representation.
  This routine is meant to replace dvdb_ftinterp_setup
  and performs the potential interpolation one perturbation at a time.

INPUTS

  ngqpt(3)=Divisions of the ab-initio q-mesh.
  nqshift=Number of shifts used to generated the ab-initio q-mesh.
  qshift(3,nqshift)=The shifts of the ab-initio q-mesh.
  nfft=Number of fft-points treated by this processors
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  nrpt=Number of R-points = number of q-points in the full BZ
  nspden=Number of spin densities.
  ipert=index of the perturbation to be treated [1,natom3]
  comm=MPI communicator

OUTPUT

  v1scf_rpt(2,nrpt,nfft,nspden)

PARENTS

      m_dvdb

CHILDREN

      cwtime,dvdb_get_v1scf_qpt,dvdb_get_v1scf_rpt,dvdb_open_read
      hdr_fort_read,hdr_fort_write,hdr_free,irreducible_set_pert
      kpts_ibz_from_kptrlatt,littlegroup_q,wrtout,xmpi_barrier

SOURCE

2357 subroutine dvdb_get_v1scf_rpt(db, cryst, ngqpt, nqshift, qshift, nfft, ngfft, &
2358 &                             nrpt, nspden, ipert, v1scf_rpt, comm)
2359 
2360 
2361 !This section has been created automatically by the script Abilint (TD).
2362 !Do not modify the following lines by hand.
2363 #undef ABI_FUNC
2364 #define ABI_FUNC 'dvdb_get_v1scf_rpt'
2365 !End of the abilint section
2366 
2367  implicit none
2368 
2369 !Arguments ------------------------------------
2370 !scalars
2371  integer,intent(in) :: nqshift,nfft,nrpt,nspden,ipert,comm
2372  type(dvdb_t),target,intent(inout) :: db
2373 !arrays
2374  integer,intent(in) :: ngqpt(3),ngfft(18)
2375  real(dp),intent(in) :: qshift(3,nqshift)
2376  real(dp),intent(out) :: v1scf_rpt(2,nrpt,nfft,nspden)
2377  type(crystal_t),intent(in) :: cryst
2378 
2379 !Local variables-------------------------------
2380 !scalars
2381  integer,parameter :: sppoldbl1=1,timrev1=1,tim_fourdp0=0
2382  integer :: my_qptopt,iq_ibz,nqibz,iq_bz,nqbz
2383  integer :: ii,iq_dvdb,cplex_qibz,ispden,irpt,idir,iat
2384  integer :: iqst,nqst,itimrev,tsign,isym,ix,iy,iz,nq1,nq2,nq3,r1,r2,r3
2385  integer :: nproc,my_rank,ifft,cnt,ierr
2386  real(dp) :: dksqmax,phre,phim
2387  logical :: isirr_q, found
2388 !arrays
2389  integer :: qptrlatt(3,3),g0q(3)
2390  integer,allocatable :: indqq(:,:),iperm(:),bz2ibz_sort(:),nqsts(:),iqs_dvdb(:)
2391  real(dp) :: qpt_bz(3),shift(3) !,qpt_ibz(3)
2392  real(dp),allocatable :: qibz(:,:),qbz(:,:),wtq(:),emiqr(:,:)
2393  real(dp),allocatable :: v1r_qibz(:,:,:,:),v1r_qbz(:,:,:,:), v1r_lr(:,:)
2394 
2395 ! *************************************************************************
2396 
2397  nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2398 
2399  nq1 = ngqpt(1); nq2 = ngqpt(2); nq3 = ngqpt(3); my_qptopt = 1 !; if (present(qptopt)) my_qptopt = qptopt
2400 
2401  ! Generate the q-mesh by finding the IBZ and the corresponding weights.
2402  qptrlatt = 0
2403  do ii=1,3
2404    qptrlatt(ii,ii) = ngqpt(ii)
2405  end do
2406 
2407  ! Get IBZ and BZ.
2408  call kpts_ibz_from_kptrlatt(cryst, qptrlatt, my_qptopt, nqshift, qshift, &
2409    nqibz, qibz, wtq, nqbz, qbz) ! new_kptrlatt, new_shiftk)
2410 
2411  !write(std_out,*)"Irreducible q-points:"
2412  !do iq_ibz=1,nqibz; write(std_out,*)trim(ktoa(qibz(:,iq_ibz))),wtq(iq_ibz)*nqbz; end do
2413  ABI_CHECK(nqbz == product(ngqpt) * nqshift, "nqbz /= product(ngqpt) * nqshift")
2414  ABI_CHECK(nqbz == nrpt, "nqbz /= nrpt")
2415 
2416  db%nrpt = nqbz
2417 
2418  ABI_CHECK(nspden == db%nspden, "nspden /= db%nspden")
2419 
2420  ! We want a gamma centered q-mesh for FFT.
2421  ABI_CHECK(nqshift == 1, "nshift > 1 not supported")
2422  ABI_CHECK(all(qshift(:,1) == zero), "qshift != 0 not supported")
2423 
2424  ABI_FREE(qbz)
2425  ABI_MALLOC(qbz, (3, nqbz))
2426  ii = 0
2427  do iz=0,nq3-1
2428    do iy=0,nq2-1
2429      do ix=0,nq1-1
2430        ii = ii + 1
2431        qbz(:, ii) = [ix/dble(nq1), iy/dble(nq2), iz/dble(nq3)]
2432        call wrap2_pmhalf([ix/dble(nq1), iy/dble(nq2), iz/dble(nq3)], qbz(:,ii), shift)
2433      end do
2434    end do
2435  end do
2436 
2437 
2438  ! Compute real-space points.
2439  ! Use the following indexing (N means ngfft of the adequate direction)
2440  ! 0 1 2 3 ... N/2    -(N-1)/2 ... -1    <= gc
2441  ! 1 2 3 4 ....N/2+1  N/2+2    ...  N    <= index ig
2442  ABI_MALLOC(db%rpt, (3, db%nrpt))
2443  ii = 0
2444  do iz=1,nq3
2445    r3 = ig2gfft(iz,nq3) !* nq3
2446    do iy=1,nq2
2447      r2 = ig2gfft(iy,nq2) !* nq2
2448      do ix=1,nq1
2449        r1 = ig2gfft(ix,nq1) !* nq1
2450        ii = ii + 1
2451        db%rpt(:,ii) = [r1, r2, r3]
2452      end do
2453    end do
2454  end do
2455 
2456  ! Find correspondence BZ --> IBZ. Note:
2457  !   * q --> -q symmetry is always used for phonons.
2458  !   * we use symrec instead of symrel
2459  ABI_MALLOC(indqq, (nqbz*sppoldbl1,6))
2460  call listkk(dksqmax,cryst%gmet,indqq,qibz,qbz,nqibz,nqbz,cryst%nsym,&
2461    sppoldbl1,cryst%symafm,cryst%symrec,timrev1,use_symrec=.True.)
2462 
2463  if (dksqmax > tol12) then
2464    MSG_BUG("Something wrong in the generation of the q-points in the BZ! Cannot map BZ --> IBZ")
2465  end if
2466 
2467  ! Construct sorted mapping BZ --> IBZ to speedup qbz search below.
2468  ABI_MALLOC(iperm, (nqbz))
2469  ABI_MALLOC(bz2ibz_sort, (nqbz))
2470  iperm = [(ii, ii=1,nqbz)]
2471  bz2ibz_sort = indqq(:,1)
2472  call sort_int(nqbz,bz2ibz_sort,iperm)
2473 
2474  ! Reconstruct the IBZ according to what is present in the DVDB .
2475  ABI_MALLOC(nqsts, (nqibz))
2476  ABI_MALLOC(iqs_dvdb, (nqibz))
2477 
2478  ABI_MALLOC(v1r_lr, (2,nfft))
2479 
2480  iqst = 0
2481  do iq_ibz=1,nqibz
2482 
2483    ! In each q-point star, count the number of q-points
2484    ! and find the one present in DVDB.
2485    nqst = 0
2486    found = .false.
2487    do ii=iqst+1,nqbz
2488      if (bz2ibz_sort(ii) /= iq_ibz) exit
2489      nqst = nqst + 1
2490 
2491      iq_bz = iperm(ii)
2492      if (.not. found) then
2493        iq_dvdb = dvdb_findq(db, qbz(:,iq_bz))
2494        if (iq_dvdb /= -1) then
2495          qibz(:,iq_ibz) = qbz(:,iq_bz)
2496          iqs_dvdb(iq_ibz) = iq_dvdb
2497          found = .true.
2498        end if
2499      end if
2500    end do
2501 
2502    ! Check that nqst has been counted properly.
2503    ABI_CHECK(nqst > 0 .and. bz2ibz_sort(iqst+1) == iq_ibz, "Wrong iqst")
2504    if (abs(nqst - wtq(iq_ibz) * nqbz) > tol12) then
2505      write(std_out,*)nqst, wtq(iq_ibz) * nqbz
2506      MSG_ERROR("Error in counting q-point star or in the weights.")
2507    end if
2508 
2509    ! Check that the q-point has been found in DVDB.
2510    if (.not. found) then
2511      MSG_ERROR(sjoin("Cannot find symmetric q-point of:", ktoa(qibz(:,iq_ibz)), "in DVDB file"))
2512    end if
2513    !write(std_out,*)sjoin("qpt irred:",ktoa(qibz(:,iq_ibz)))
2514 
2515    iqst = iqst + nqst
2516    nqsts(iq_ibz) = nqst
2517 
2518  end do
2519 
2520  ! Redo the mapping with the new IBZ
2521  call listkk(dksqmax,cryst%gmet,indqq,qibz,qbz,nqibz,nqbz,cryst%nsym,&
2522    sppoldbl1,cryst%symafm,cryst%symrec,timrev1,use_symrec=.True.)
2523 
2524  ABI_MALLOC(emiqr, (2,db%nrpt))
2525  v1scf_rpt = zero
2526 
2527  ABI_MALLOC(v1r_qbz, (2, nfft, db%nspden, db%natom3))
2528  !v1r_qbz = huge(one)
2529 
2530  iqst = 0
2531  do iq_ibz=1,nqibz
2532 
2533    ! Get potentials for this IBZ q-point on the real-space FFT mesh.
2534    ! This call allocates v1r_qibz(cplex_qibz, nfft, nspden, 3*natom)
2535    call dvdb_readsym_allv1(db, iqs_dvdb(iq_ibz), cplex_qibz, nfft, ngfft, v1r_qibz, comm)
2536 
2537    ! Reconstruct by symmetry the potentials for the star of this q-point, perform FT and accumulate
2538    ! Be careful with the gamma point.
2539    do ii=1,nqsts(iq_ibz)
2540      iqst = iqst + 1
2541      iq_bz = iperm(iqst)
2542      ABI_CHECK(iq_ibz == indqq(iq_bz,1), "iq_ibz !/ ind qq(1)")
2543      isym = indqq(iq_bz,2); itimrev = indqq(iq_bz,6) + 1; g0q = indqq(iq_bz,3:5) ! IS(q_ibz) + g0q = q_bz
2544      tsign = 3-2*itimrev
2545 
2546      qpt_bz = qbz(:, iq_bz)
2547      !write(std_out,*)"  treating:",trim(ktoa(qpt_bz))
2548      isirr_q = (isym == 1 .and. itimrev == 1 .and. all(g0q == 0))
2549      !ABI_CHECK(all(g0q == 0), "g0q /= 0")
2550 
2551      ! Compute long-range part of the coupling potential
2552      v1r_lr = zero; cnt = 0
2553      if (db%has_dielt_zeff) then
2554        idir = mod(ipert-1, 3) + 1; iat = (ipert - idir) / 3 + 1
2555        call dvdb_v1r_long_range(db,qpt_bz,iat,idir,nfft,ngfft,v1r_lr)
2556      end if
2557 
2558      if (cplex_qibz == 1) then
2559        ! Gamma point.
2560        ABI_CHECK(nqsts(iq_ibz) == 1, "cplex_qibz == 1 and nq nqst /= 1 (should be gamma)")
2561        ABI_CHECK(all(g0q == 0), "gamma point with g0q /= 0")
2562 
2563        ! Substract the long-range part of the potential
2564        do ispden=1,db%nspden
2565          v1r_qibz(1,:,ispden,ipert) = v1r_qibz(1,:,ispden,ipert) - v1r_lr(1,:)
2566        end do
2567 
2568        ! Slow FT.
2569        cnt = 0
2570        do ispden=1,db%nspden
2571          do irpt=1,db%nrpt
2572            ! MPI-parallelism
2573            cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle
2574            do ifft=1,nfft
2575              v1scf_rpt(1,irpt,ifft,ispden) = v1scf_rpt(1,irpt,ifft,ispden) + &
2576                                              v1r_qibz(1, ifft, ispden, ipert)
2577            end do
2578          end do
2579        end do
2580 
2581      else
2582        ! Get the periodic part of the potential in BZ (v1r_qbz)
2583        if (isirr_q) then
2584          !write(std_out,*)sjoin("qpt irred:",ktoa(qpt_bz))
2585          v1r_qbz = v1r_qibz
2586        else
2587          call v1phq_rotate(cryst, qibz(:,iq_ibz), isym, itimrev, g0q,&
2588            ngfft, cplex_qibz, nfft, db%nspden, db%nsppol, db%mpi_enreg, v1r_qibz, v1r_qbz)
2589          !v1r_qbz = zero; v1r_qbz = v1r_qibz
2590 
2591          !call times_eigr(-tsign * g0q, ngfft, nfft, db%nspden*db%natom3, v1r_qbz)
2592          !call times_eigr(+tsign * g0q, ngfft, nfft, db%nspden*db%natom3, v1r_qbz)
2593          !if (itimrev == 2) v1r_qbz(2,:,:,:) = -v1r_qbz(2,:,:,:)
2594          !call times_eigr(-tsign * g0q, ngfft, nfft, db%nspden*db%natom3, v1r_qbz)
2595          !call times_eigr(+tsign * g0q, ngfft, nfft, db%nspden*db%natom3, v1r_qbz)
2596        end if
2597 
2598        ! Multiply by e^{iqpt_bz.r}
2599        call times_eikr(qpt_bz, ngfft, nfft, db%nspden*db%natom3, v1r_qbz)
2600 
2601        ! Substract the long-range part of the potential
2602        do ispden=1,db%nspden
2603          v1r_qbz(1,:,ispden,ipert) = v1r_qbz(1,:,ispden,ipert) - v1r_lr(1,:)
2604          v1r_qbz(2,:,ispden,ipert) = v1r_qbz(2,:,ispden,ipert) - v1r_lr(2,:)
2605        end do
2606 
2607        ! Compute FT phases for this qpt_bz.
2608        call calc_eiqr(-qpt_bz,db%nrpt,db%rpt,emiqr)
2609 
2610        ! Slow FT.
2611        cnt = 0
2612        do ispden=1,db%nspden
2613          do irpt=1,db%nrpt
2614            cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI parallelism.
2615            phre = emiqr(1,irpt); phim = emiqr(2,irpt)
2616 
2617            do ifft=1,nfft
2618              v1scf_rpt(1,irpt,ifft,ispden) = v1scf_rpt(1,irpt,ifft,ispden) &
2619              &                      + phre * v1r_qbz(1, ifft, ispden, ipert) &
2620              &                      - phim * v1r_qbz(2, ifft, ispden, ipert)
2621 
2622              v1scf_rpt(2,irpt,ifft,ispden) = v1scf_rpt(2,irpt,ifft,ispden) &
2623              &                      + phre * v1r_qbz(2, ifft, ispden, ipert) &
2624              &                      + phim * v1r_qbz(1, ifft, ispden, ipert)
2625            end do
2626 
2627          end do
2628        end do
2629      end if
2630 
2631    end do ! iqst
2632 
2633    ABI_FREE(v1r_qibz)
2634  end do ! iq_ibz
2635  ABI_CHECK(iqst == nqbz, "iqst /= nqbz")
2636 
2637  v1scf_rpt = v1scf_rpt / nqbz
2638  call xmpi_sum(v1scf_rpt, comm, ierr)
2639 
2640  !! Build mpi_type for executing fourdp in sequential.
2641  !call ngfft_seq(ngfft_qspace, [nq1, nq2, nq3])
2642  !call initmpi_seq(mpi_enreg_seq)
2643  !call init_distribfft_seq(mpi_enreg_seq%distribfft,'c',ngfft_qspace(2),ngfft_qspace(3),'all')
2644  !call init_distribfft_seq(mpi_enreg_seq%distribfft,'f',ngfft_qspace(2),ngfft_qspace(3),'all')
2645 
2646  !!ABI_STAT_MALLOC(all_v1qr, (nqbz, 2, nfft, db%nspden, db%natom3), ierr)
2647  !!ABI_CHECK(ierr==0, "oom in all_v1qr")
2648  !!all_v1qr = zero
2649 
2650  !cnt = 0
2651  !do mu=1,db%natom3
2652  !  do ispden=1,db%nspden
2653  !    do ifft=1,nfft
2654  !      !cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle
2655  !      !call fourdp(1,all_v1qr(:,:,ifft,ispden,mu),db%v1scf_rpt(:,ifft,ispden,mu),+1,&
2656  !      ! mpi_enreg_seq,nqbz,ngfft_qspace,paral_kgb0,tim_fourdp0)
2657  !    end do
2658  !  end do
2659  !end do
2660 
2661  !!ABI_FREE(all_v1qr)
2662  !call destroy_mpi_enreg(mpi_enreg_seq)
2663 
2664  ABI_FREE(emiqr)
2665  ABI_FREE(qibz)
2666  ABI_FREE(wtq)
2667  ABI_FREE(qbz)
2668  ABI_FREE(indqq)
2669  ABI_FREE(iperm)
2670  ABI_FREE(bz2ibz_sort)
2671  ABI_FREE(iqs_dvdb)
2672  ABI_FREE(nqsts)
2673  ABI_FREE(v1r_qbz)
2674  ABI_FREE(v1r_lr)
2675 
2676 
2677 end subroutine dvdb_get_v1scf_rpt

m_dvdb/dvdb_init [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  dvdb_init

FUNCTION

  Initialize the object from file. This is a COLLECTIVE procedure that must be called
  by each process in the MPI communicator comm.

INPUTS

   path=DVDB Filename.
   comm=MPI communicator.

PARENTS

      eph,m_dvdb,mrgdv

CHILDREN

      cwtime,dvdb_get_v1scf_qpt,dvdb_get_v1scf_rpt,dvdb_open_read
      hdr_fort_read,hdr_fort_write,hdr_free,irreducible_set_pert
      kpts_ibz_from_kptrlatt,littlegroup_q,wrtout,xmpi_barrier

SOURCE

289 subroutine dvdb_init(db, path, comm)
290 
291 
292 !This section has been created automatically by the script Abilint (TD).
293 !Do not modify the following lines by hand.
294 #undef ABI_FUNC
295 #define ABI_FUNC 'dvdb_init'
296 !End of the abilint section
297 
298  implicit none
299 
300 !Arguments ------------------------------------
301 !scalars
302  character(len=*),intent(in) :: path
303  integer,intent(in) :: comm
304  type(dvdb_t),intent(inout) :: db
305 
306 !Local variables-------------------------------
307 !scalars
308  integer,parameter :: master=0,timrev2=2
309  integer :: iv1,ii,ierr,unt,fform,nqpt,iq,iq_found,cplex,trev_q
310  integer :: idir,ipert,my_rank
311  character(len=500) :: msg
312  type(hdr_type) :: hdr1,hdr_ref
313 !arrays
314  integer,allocatable :: tmp_pos(:,:,:)
315  real(dp),allocatable :: tmp_qpts(:,:)
316 
317 !************************************************************************
318 
319  my_rank = xmpi_comm_rank(comm)
320  db%path = path; db%comm = comm; db%iomode = IO_MODE_FORTRAN
321 
322  ! Master reads the header and builds useful tables
323  if (my_rank == master) then
324 
325    if (open_file(path, msg, newunit=unt, form="unformatted", status="old", action="read") /= 0) then
326      MSG_ERROR(msg)
327    end if
328    read(unt, err=10, iomsg=msg) db%version
329    read(unt, err=10, iomsg=msg) db%numv1
330 
331    ! Get important dimensions from the first header and rewind the file.
332    call hdr_fort_read(hdr_ref, unt, fform)
333    if (dvdb_check_fform(fform, "read_dvdb", msg) /= 0) then
334      MSG_ERROR(sjoin("While reading:", path, ch10, msg))
335    end if
336    if (db%debug) call hdr_echo(hdr_ref,fform,4,unit=std_out)
337 
338    rewind(unt)
339    read(unt, err=10, iomsg=msg)
340    read(unt, err=10, iomsg=msg)
341 
342    ! The code below must be executed by the other procs if MPI.
343    db%natom = hdr_ref%natom
344    db%natom3 = 3 * hdr_ref%natom
345    db%nspden = hdr_ref%nspden
346    db%nsppol = hdr_ref%nsppol
347    db%nspinor = hdr_ref%nspinor
348    db%usepaw = hdr_ref%usepaw
349    ABI_CHECK(db%usepaw == 0, "PAW not yet supported")
350 
351    ! TODO: Write function that returns mpert from natom!
352    db%mpert = db%natom + 6
353 
354    ABI_MALLOC(tmp_qpts, (3, db%numv1))
355    ABI_MALLOC(tmp_pos, (3, db%mpert, db%numv1))
356    tmp_pos = 0
357 
358    ABI_MALLOC(db%cplex_v1, (db%numv1))
359    ABI_MALLOC(db%ngfft3_v1, (3, db%numv1))
360    ABI_MALLOC(db%iv_pinfoq, (4, db%numv1))
361    ABI_MALLOC(db%rhog1_g0, (2, db%numv1))
362 
363    nqpt = 0
364    do iv1=1,db%numv1
365      call hdr_fort_read(hdr1, unt, fform)
366      if (dvdb_check_fform(fform, "read_dvdb", msg) /= 0) then
367        MSG_ERROR(sjoin("While reading hdr of v1 potential of index:", itoa(iv1), ch10, msg))
368      end if
369 
370      ! Save cplex and FFT mesh associated to this perturbation.
371      cplex = 2; if (hdr1%qptn(1)**2+hdr1%qptn(2)**2+hdr1%qptn(3)**2<1.d-14) cplex = 1
372      db%cplex_v1(iv1) = cplex
373      db%ngfft3_v1(:, iv1) = hdr1%ngfft(:3)
374 
375      ! Skip the records with v1.
376      do ii=1,hdr1%nspden
377        read(unt, err=10, iomsg=msg)
378      end do
379      ! Read rhog1_g0 (if available)
380      db%rhog1_g0(:, iv1) = zero
381      if (db%version > 1) read(unt, err=10, iomsg=msg) db%rhog1_g0(:, iv1)
382 
383      ! Check whether this q-point is already in the list.
384      iq_found = 0
385      do iq=1,nqpt
386        if (all(abs(hdr1%qptn - tmp_qpts(:,iq)) < tol14)) then
387          iq_found = iq; exit
388        end if
389      end do
390 
391      ! pertcase = idir + (ipert-1)*3 where ipert=iatom in the interesting cases
392      idir = mod(hdr1%pertcase-1, 3) + 1
393      ipert = (hdr1%pertcase - idir) / 3 + 1
394 
395      ! Increment nqpt is new q-points and update tmp_pos
396      if (iq_found == 0) then
397        nqpt = nqpt + 1
398        tmp_qpts(:, nqpt) = hdr1%qptn
399        iq_found = nqpt
400      end if
401      tmp_pos(idir, ipert, iq_found) = iv1
402      db%iv_pinfoq(:,iv1) = [idir, ipert, hdr1%pertcase, iq_found]
403 
404      call hdr_free(hdr1)
405    end do
406 
407    ! Allocate arrays with correct nqpt dimension
408    db%nqpt = nqpt
409    ABI_MALLOC(db%qpts, (3, nqpt))
410    db%qpts = tmp_qpts(:,1:nqpt)
411    ABI_FREE(tmp_qpts)
412 
413    ABI_MALLOC(db%pos_dpq, (3, db%mpert, nqpt))
414    db%pos_dpq = tmp_pos(:, :, 1:nqpt)
415    ABI_FREE(tmp_pos)
416 
417    close(unt)
418  end if
419 
420  ! Master broadcasts data.
421  if (xmpi_comm_size(comm) > 1) then
422    call xmpi_bcast(db%version, master, comm, ierr)
423    call xmpi_bcast(db%numv1, master, comm, ierr)
424    call xmpi_bcast(db%nqpt, master, comm, ierr)
425    call hdr_bcast(hdr_ref, master, my_rank, comm)
426 
427    db%natom = hdr_ref%natom
428    db%natom3 = 3 * hdr_ref%natom
429    db%nspden = hdr_ref%nspden
430    db%nsppol = hdr_ref%nsppol
431    db%nspinor = hdr_ref%nspinor
432    db%usepaw = hdr_ref%usepaw
433    db%mpert = db%natom + 6
434 
435    if (my_rank /= master) then
436      ABI_MALLOC(db%cplex_v1, (db%numv1))
437      ABI_MALLOC(db%ngfft3_v1, (3, db%numv1))
438      ABI_MALLOC(db%iv_pinfoq, (4, db%numv1))
439      ABI_MALLOC(db%qpts, (3, db%nqpt))
440      ABI_MALLOC(db%pos_dpq, (3, db%mpert, db%nqpt))
441      ABI_MALLOC(db%rhog1_g0, (2, db%numv1))
442    end if
443 
444    call xmpi_bcast(db%cplex_v1, master, comm, ierr)
445    call xmpi_bcast(db%ngfft3_v1, master, comm, ierr)
446    call xmpi_bcast(db%iv_pinfoq, master, comm, ierr)
447    call xmpi_bcast(db%qpts, master, comm, ierr)
448    call xmpi_bcast(db%pos_dpq, master, comm, ierr)
449    call xmpi_bcast(db%rhog1_g0, master, comm, ierr)
450  end if
451 
452  ! Init crystal_t from the hdr read from file.
453  call crystal_from_hdr(db%cryst,hdr_ref,timrev2)
454  call hdr_free(hdr_ref)
455 
456  ! Init Born effective charges
457  ABI_MALLOC(db%zeff, (3, 3, db%natom))
458 
459  ! Internal MPI_type needed for calling fourdp!
460  call initmpi_seq(db%mpi_enreg)
461 
462  ! Precompute symq_table for all q-points in the DVDB.
463  ABI_MALLOC(db%symq_table, (4, 2, db%cryst%nsym, db%nqpt))
464  do iq=1,db%nqpt
465    call littlegroup_q(db%cryst%nsym, db%qpts(:,iq), db%symq_table(:,:,:,iq), &
466      db%cryst%symrec, db%cryst%symafm, trev_q, prtvol=0)
467  end do
468 
469  return
470 
471  ! Handle Fortran IO error
472 10 continue
473  MSG_ERROR(sjoin("Error while reading:", path, ch10, msg))
474 
475 end subroutine dvdb_init

m_dvdb/dvdb_interpolate_and_write [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  dvdb_interpolate_and_write

FUNCTION

  Interpolate the phonon potential onto a fine q-point grid
  and write the data in a new DVDB file.

INPUTS

OUTPUT

PARENTS

      eph

CHILDREN

      cwtime,dvdb_get_v1scf_qpt,dvdb_get_v1scf_rpt,dvdb_open_read
      hdr_fort_read,hdr_fort_write,hdr_free,irreducible_set_pert
      kpts_ibz_from_kptrlatt,littlegroup_q,wrtout,xmpi_barrier

SOURCE

4275 subroutine dvdb_interpolate_and_write(dtfil, ngfft, ngfftf, cryst, dvdb, &
4276 &          ngqpt_coarse, nqshift_coarse, qshift_coarse, &
4277 &          ngqpt, qptopt, mpi_enreg, comm)
4278 
4279  use m_time,           only : cwtime
4280 
4281 !This section has been created automatically by the script Abilint (TD).
4282 !Do not modify the following lines by hand.
4283 #undef ABI_FUNC
4284 #define ABI_FUNC 'dvdb_interpolate_and_write'
4285 !End of the abilint section
4286 
4287  implicit none
4288 
4289 !Arguments ------------------------------------
4290 !scalars
4291  integer,intent(in) :: comm
4292  integer,intent(in) :: qptopt,nqshift_coarse
4293  type(datafiles_type),intent(in) :: dtfil
4294  type(crystal_t),intent(in) :: cryst
4295  type(dvdb_t),target,intent(inout) :: dvdb
4296  type(mpi_type),intent(inout) :: mpi_enreg
4297 !arrays
4298  integer,intent(in) :: ngfft(18), ngfftf(18)
4299  integer,intent(in) :: ngqpt(3), ngqpt_coarse(3)
4300  real(dp),intent(in) :: qshift_coarse(3,nqshift_coarse)
4301 
4302 !Local variables ------------------------------
4303 !scalars
4304  integer,parameter :: master=0
4305  integer :: fform_pot=111
4306  integer :: ierr
4307  integer :: my_rank,nproc,idir,ipert,iat,ipc,ispden
4308  integer :: cplex,db_iqpt,natom,natom3,npc,trev_q,nspden
4309  integer :: nqbz, nqibz, iq, ifft, nqbz_coarse
4310  integer :: nperts_read, nperts_interpolate, nperts
4311  integer :: nqpt_read, nqpt_interpolate
4312  integer :: nfft,nfftf
4313  integer :: ount, unt, fform
4314  real(dp) :: cpu,wall,gflops
4315  logical :: i_am_master
4316  character(len=500) :: msg
4317  character(len=fnlen) :: new_ddb_fname
4318 !arrays
4319  integer :: qptrlatt(3,3)
4320  integer :: symq(4,2,cryst%nsym)
4321  integer :: rfdir(3)
4322  integer,allocatable :: pinfo(:,:),rfpert(:),pertsy(:,:,:),iq_read(:)
4323  real(dp) :: qpt(3)
4324  real(dp) :: rhog1_g0(2)
4325  real(dp),allocatable :: v1scf(:,:,:), v1scf_rpt(:,:,:,:),v1(:)
4326  real(dp),allocatable :: wtq(:),qibz(:,:),qbz(:,:),q_interp(:,:),q_read(:,:)
4327  type(hdr_type) :: hdr_ref
4328 
4329 !************************************************************************
4330  ABI_UNUSED(mpi_enreg%nproc)
4331 
4332  call cwtime(cpu,wall,gflops,"start")
4333 
4334  write(msg, '(2a)') "Interpolation of the electron-phonon coupling potential", ch10
4335  call wrtout(ab_out, msg, "COLL", do_flush=.True.)
4336  call wrtout(std_out, msg, "COLL", do_flush=.True.)
4337 
4338  my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
4339  i_am_master = (my_rank == master)
4340 
4341  ! ======================= !
4342  ! Setup fine q-point grid
4343  ! ======================= !
4344 
4345  nfftf = product(ngfftf(1:3))
4346  nfft = product(ngfft(1:3))
4347 
4348  ! Generate the list of irreducible q-points in the grid
4349  qptrlatt = zero
4350  qptrlatt(1,1) = ngqpt(1); qptrlatt(2,2) = ngqpt(2); qptrlatt(3,3) = ngqpt(3)
4351  call kpts_ibz_from_kptrlatt(cryst, qptrlatt, qptopt, 1, &
4352  &                           [zero, zero, zero], nqibz, qibz, wtq, nqbz, qbz)
4353 
4354  nqbz_coarse = product(ngqpt_coarse) * nqshift_coarse
4355 
4356  ! ========================================== !
4357  ! Prepare the header to write the potentials
4358  ! ========================================== !
4359 
4360  ! Read the first header
4361  if (my_rank == master) then
4362    if (open_file(dvdb%path, msg, newunit=unt, form="unformatted", status="old", action="read") /= 0) then
4363      MSG_ERROR(msg)
4364    end if
4365    read(unt, err=10, iomsg=msg) dvdb%version
4366    read(unt, err=10, iomsg=msg) dvdb%numv1
4367 
4368    call hdr_fort_read(hdr_ref, unt, fform)
4369    if (dvdb_check_fform(fform, "read_dvdb", msg) /= 0) then
4370      MSG_ERROR(sjoin("While reading:", dvdb%path, ch10, msg))
4371    end if
4372    close(unt)
4373  end if
4374 
4375  ! Reset the symmetries of the header
4376  ! One might have disable the symmetries in the response function calculation
4377  ! that produced the initial set of potentials present in the DVDB.
4378  ! This is because the symmetry features are not used in all parts
4379  ! of the response function driver.
4380  if (allocated(hdr_ref%symrel)) ABI_FREE(hdr_ref%symrel)
4381  if (allocated(hdr_ref%tnons)) ABI_FREE(hdr_ref%tnons)
4382  if (allocated(hdr_ref%symafm)) ABI_FREE(hdr_ref%symafm)
4383  hdr_ref%nsym = cryst%nsym
4384  ABI_MALLOC(hdr_ref%symrel, (3,3,hdr_ref%nsym))
4385  ABI_MALLOC(hdr_ref%tnons, (3,hdr_ref%nsym))
4386  ABI_MALLOC(hdr_ref%symafm, (hdr_ref%nsym))
4387  hdr_ref%symrel(:,:,:) = cryst%symrel(:,:,:)
4388  hdr_ref%tnons(:,:) = cryst%tnons(:,:)
4389  hdr_ref%symafm(:) = cryst%symafm(:)
4390 
4391  hdr_ref%ngfft = ngfftf(1:3)
4392 
4393 
4394  ! ======================================= !
4395  ! Open DVDB and copy important dimensions
4396  ! ======================================= !
4397 
4398  call dvdb_open_read(dvdb, ngfftf, xmpi_comm_self)
4399 
4400  cplex = 2
4401  natom = cryst%natom
4402  natom3 = 3 * natom
4403  nspden = dvdb%nspden
4404 
4405  ! ================================================== !
4406  ! Sort the q-points to read and those to interpolate
4407  ! and find the irreducible perturbations
4408  ! ================================================== !
4409  ABI_MALLOC(iq_read, (nqibz))
4410  ABI_MALLOC(q_read, (3,nqibz))
4411  ABI_MALLOC(q_interp, (3,nqibz))
4412 
4413  ABI_MALLOC(pertsy, (nqibz,3,dvdb%mpert))
4414  ABI_MALLOC(rfpert, (dvdb%mpert))
4415  ABI_MALLOC(pinfo, (3,3*dvdb%mpert))
4416  rfpert = 0; rfpert(1:cryst%natom) = 1; rfdir = 1
4417 
4418  pertsy = zero
4419 
4420  nqpt_read = 0
4421  nperts_read = 0
4422  nqpt_interpolate = 0
4423  nperts_interpolate = 0
4424 
4425  do iq=1,nqibz
4426 
4427    qpt = qibz(:,iq)
4428 
4429    ! Find the index of the q-point in the DVDB.
4430    db_iqpt = dvdb_findq(dvdb, qpt)
4431 
4432    if (db_iqpt /= -1) then
4433 
4434      call wrtout(std_out, sjoin("Q-point: ",ktoa(qpt)," found in DVDB with index ",itoa(db_iqpt)))
4435      nqpt_read = nqpt_read + 1
4436      q_read(:,nqpt_read) = qpt(:)
4437      iq_read(nqpt_read) = db_iqpt
4438 
4439      ! Count the perturbations
4440      npc = dvdb_get_pinfo(dvdb, db_iqpt, cplex, pinfo)
4441      do ipc=1,npc
4442        idir = pinfo(1,ipc); iat = pinfo(2,ipc); ipert = pinfo(3, ipc)
4443        if (iat .le. natom) nperts_read = nperts_read + 1
4444      end do
4445 
4446    else
4447 
4448      call wrtout(std_out, sjoin("Q-point: ",ktoa(qpt), "not found in DVDB. Will interpolate."))
4449      nqpt_interpolate = nqpt_interpolate + 1
4450      q_interp(:,nqpt_interpolate) = qpt(:)
4451 
4452 
4453      ! Examine the symmetries of the q wavevector
4454      call littlegroup_q(cryst%nsym,qpt,symq,cryst%symrec,cryst%symafm,trev_q,prtvol=0)
4455 
4456      ! Find the list of irreducible perturbations for this q-point.
4457      call irreducible_set_pert(cryst%indsym,dvdb%mpert,cryst%natom,cryst%nsym,&
4458      &    pertsy(nqpt_interpolate,:,:),rfdir,rfpert,symq,cryst%symrec,cryst%symrel)
4459 
4460      do iat=1,natom
4461        do idir=1,3
4462          ipert = (iat-1) * 3 + idir
4463          if (pertsy(nqpt_interpolate,idir,iat) == 1) then
4464            nperts_interpolate = nperts_interpolate + 1
4465          end if
4466 
4467        end do
4468      end do
4469 
4470    end if
4471  end do
4472 
4473 
4474  ! ================================================= !
4475  ! Open the new DVDB file and write preliminary info
4476  ! ================================================= !
4477  nperts = nperts_read + nperts_interpolate
4478 
4479  if (my_rank == master) then
4480    new_ddb_fname = strcat(dtfil%filnam_ds(4), '_DVDB')
4481    if (open_file(new_ddb_fname, msg, newunit=ount, form="unformatted", action="write", status="unknown") /= 0) then
4482      MSG_ERROR(msg)
4483    end if
4484    write(ount, err=10, iomsg=msg) dvdb_last_version
4485    write(ount, err=10, iomsg=msg) nperts
4486  end if
4487 
4488 
4489  ! ================================================================= !
4490  ! Master reads all available perturbations and copy in the new DVDB
4491  ! ================================================================= !
4492 
4493  ABI_MALLOC(v1scf, (cplex,nfftf,nspden))
4494  ABI_MALLOC(v1, (cplex*nfftf))
4495 
4496  rhog1_g0 = zero
4497 
4498  if (my_rank == master) then
4499    do iq=1,nqpt_read
4500 
4501      qpt = q_read(:,iq)
4502      db_iqpt = iq_read(iq)
4503 
4504      ! Read each irreducible perturbation potentials
4505      npc = dvdb_get_pinfo(dvdb, db_iqpt, cplex, pinfo)
4506      ABI_CHECK(npc /= 0, "npc == 0!")
4507 
4508      do ipc=1,npc
4509        idir = pinfo(1,ipc); iat = pinfo(2,ipc); ipert = pinfo(3, ipc)
4510        if (dvdb_read_onev1(dvdb, idir, iat, db_iqpt, cplex, nfftf, ngfftf, v1scf, msg) /= 0) then
4511          MSG_ERROR(msg)
4512        end if
4513 
4514        hdr_ref%qptn = qpt
4515        hdr_ref%pertcase = ipert
4516 
4517        ! Write header
4518        call hdr_fort_write(hdr_ref, ount, fform_pot, ierr)
4519        ABI_CHECK(ierr == 0, "hdr_fort_write returned ierr = 0")
4520 
4521        do ispden=1,nspden
4522          v1 = reshape(v1scf(:,:,ispden), (/cplex*nfftf/))
4523          write(ount, err=10, iomsg=msg) (v1(ifft), ifft=1,cplex*nfftf)
4524        end do
4525 
4526        if (dvdb_last_version > 1) write(ount, err=10, iomsg=msg) rhog1_g0
4527 
4528      end do
4529 
4530    end do
4531  end if
4532 
4533  call xmpi_barrier(comm)
4534 
4535 
4536  ! ================================================================ !
4537  ! Interpolate the potential for q-points not in the original DVDB
4538  ! ================================================================ !
4539 
4540  dvdb%nrpt = nqbz_coarse
4541  ABI_MALLOC(v1scf_rpt, (2,dvdb%nrpt,nfftf,dvdb%nspden))
4542 
4543  do iat=1,natom
4544    do idir=1,3
4545      ipert = (iat-1) * 3 + idir
4546 
4547      if (sum(pertsy(:,idir,iat)) == -nqpt_interpolate) cycle
4548 
4549      call wrtout(std_out, sjoin("Interpolating perturbation iat,idir = ",itoa(iat), itoa(idir)))
4550 
4551      ! Compute phonon potential in real space lattice representation
4552      call dvdb_get_v1scf_rpt(dvdb, cryst, ngqpt_coarse, nqshift_coarse, &
4553      &                       qshift_coarse, nfftf, ngfftf, &
4554      &                       dvdb%nrpt, dvdb%nspden, ipert, v1scf_rpt, comm)
4555 
4556 
4557      do iq=1,nqpt_interpolate
4558 
4559        if (pertsy(iq,idir,iat) == -1) cycle
4560 
4561        qpt = q_interp(:,iq)
4562 
4563        ! Interpoalte the phonon potential
4564        call dvdb_get_v1scf_qpt(dvdb, cryst, qpt, nfftf, ngfftf, dvdb%nrpt, &
4565        &                       dvdb%nspden, ipert, v1scf_rpt, v1scf, comm)
4566 
4567        ! Master writes the file
4568        if (my_rank == master) then
4569 
4570          hdr_ref%qptn = qpt
4571          hdr_ref%pertcase = ipert
4572          ! Write header
4573          call hdr_fort_write(hdr_ref, ount, fform_pot, ierr)
4574          ABI_CHECK(ierr == 0, "hdr_fort_write returned ierr = 0")
4575 
4576          do ispden=1,nspden
4577            v1 = reshape(v1scf(:,:,ispden), (/cplex*nfftf/))
4578            write(ount, err=10, iomsg=msg) (v1(ifft), ifft=1,cplex*nfftf)
4579          end do
4580 
4581          if (dvdb_last_version > 1) write(ount, err=10, iomsg=msg) rhog1_g0
4582 
4583        end if
4584 
4585      end do
4586 
4587      ABI_FREE(dvdb%rpt)
4588 
4589    end do
4590  end do
4591 
4592  if (my_rank == master) close(ount)
4593 
4594  ! ========================================================================== !
4595  ! Free memory
4596  ! ========================================================================== !
4597 
4598  ABI_FREE(v1scf)
4599  ABI_FREE(v1)
4600  ABI_FREE(v1scf_rpt)
4601  ABI_FREE(qbz)
4602  ABI_FREE(qibz)
4603  ABI_FREE(q_interp)
4604  ABI_FREE(q_read)
4605  ABI_FREE(wtq)
4606  ABI_FREE(iq_read)
4607  ABI_FREE(pertsy)
4608  ABI_FREE(rfpert)
4609  ABI_FREE(pinfo)
4610 
4611  call hdr_free(hdr_ref)
4612 
4613  call cwtime(cpu,wall,gflops,"stop")
4614 
4615  write(msg, '(2a)') "Interpolation of the electron-phonon coupling potential completed", ch10
4616  call wrtout(ab_out, msg, "COLL", do_flush=.True.)
4617  call wrtout(std_out, msg, "COLL", do_flush=.True.)
4618 
4619 
4620  return
4621 
4622  ! Handle Fortran IO error
4623 10 continue
4624  MSG_ERROR(msg)
4625 
4626 end subroutine dvdb_interpolate_and_write

m_dvdb/dvdb_interpolate_v1scf [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  dvdb_interpolate_v1scf

FUNCTION

  Interpolate the phonon perturbation potential.
  This routine is meant to replace dvdb_ftinterp_setup and dvdb_ftinterp_qpt.
  It performs the interpolation one perturbation at a time.

INPUTS

  ngqpt(3)=Divisions of the ab-initio q-mesh.
  nqshift=Number of shifts used to generated the ab-initio q-mesh.
  qshift(3,nqshift)=The shifts of the ab-initio q-mesh.
  nfft=Number of fft-points treated by this processors
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  nfftf=Number of fft-points on the fine grid for interpolated potential
  ngfftf(18)=information on 3D FFT for interpolated potential
  comm=MPI communicator

OUTPUT

  v1scf(2, nfft, nspden, 3*natom)= v1scf potentials on the real-space FFT mesh for the 3*natom perturbations.

PARENTS

      m_gkk

CHILDREN

      cwtime,dvdb_get_v1scf_qpt,dvdb_get_v1scf_rpt,dvdb_open_read
      hdr_fort_read,hdr_fort_write,hdr_free,irreducible_set_pert
      kpts_ibz_from_kptrlatt,littlegroup_q,wrtout,xmpi_barrier

SOURCE

2842 subroutine dvdb_interpolate_v1scf(db, cryst, qpt, ngqpt, nqshift, qshift, &
2843 &                                 nfft, ngfft, nfftf, ngfftf, v1scf, comm)
2844 
2845 
2846 !This section has been created automatically by the script Abilint (TD).
2847 !Do not modify the following lines by hand.
2848 #undef ABI_FUNC
2849 #define ABI_FUNC 'dvdb_interpolate_v1scf'
2850 !End of the abilint section
2851 
2852  implicit none
2853 
2854 !Arguments ------------------------------------
2855 !scalars
2856  integer,intent(in) :: nqshift,nfft,nfftf,comm
2857  type(dvdb_t),target,intent(inout) :: db
2858 !arrays
2859  real(dp),intent(in) :: qpt(3)
2860  integer,intent(in) :: ngqpt(3),ngfft(18),ngfftf(18)
2861  real(dp),intent(in) :: qshift(3,nqshift)
2862  real(dp),allocatable,intent(out) :: v1scf(:,:,:,:)
2863  type(crystal_t),intent(in) :: cryst
2864 
2865 !Local variables-------------------------------
2866 !scalars
2867  integer :: ipert, nqbz
2868  integer :: nproc,my_rank
2869 !arrays
2870  real(dp),allocatable :: v1scf_rpt(:,:,:,:)
2871 
2872 ! *************************************************************************
2873 
2874  nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2875 
2876  nqbz = product(ngqpt) * nqshift
2877  db%nrpt = nqbz
2878 
2879  ABI_MALLOC(v1scf, (2,nfftf,db%nspden,db%natom3))
2880 
2881  ABI_MALLOC(v1scf_rpt, (2,db%nrpt,nfft,db%nspden))
2882 
2883  do ipert=1,db%natom3
2884 
2885    write(std_out, "(a,i4,a,i4,a)") "Interpolating potential for perturbation ", &
2886    &                           ipert, " / ", db%natom3, ch10
2887 
2888    ! FIXME I think this should be ngfftf and not ngfft
2889    !       Also, other calls to dvdb_ftinterp_setup should use ngfftf.
2890    call dvdb_get_v1scf_rpt(db, cryst, ngqpt, nqshift, qshift, nfft, ngfft, &
2891    &                       db%nrpt, db%nspden, ipert, v1scf_rpt, comm)
2892 
2893    call dvdb_get_v1scf_qpt(db, cryst, qpt, nfftf, ngfftf, db%nrpt, db%nspden, &
2894    &                       ipert, v1scf_rpt, v1scf(:,:,:,ipert), comm)
2895 
2896    ABI_FREE(db%rpt)
2897 
2898  end do
2899 
2900  ABI_FREE(v1scf_rpt)
2901 
2902 
2903 end subroutine dvdb_interpolate_v1scf

m_dvdb/dvdb_list_perts [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  dvdb_list_perts

FUNCTION

  Given a q-point mesh, this routine checks if all the (phonon) perturbations
  are available taking into account symmetries.

INPUTS

  ngqpt(3)=Q-mesh divisions. If all(ngqpt == -1), the list of q-points in the DVDB
    (i.e. db%qpts) is analyzed instead of the q-points generated from ngqpt.
  [unit]=Unit number for output. Default `std_out`.

OUTPUT

  Only writing.

PARENTS

      eph,m_dvdb,mrgdv

CHILDREN

      cwtime,dvdb_get_v1scf_qpt,dvdb_get_v1scf_rpt,dvdb_open_read
      hdr_fort_read,hdr_fort_write,hdr_free,irreducible_set_pert
      kpts_ibz_from_kptrlatt,littlegroup_q,wrtout,xmpi_barrier

SOURCE

3191 subroutine dvdb_list_perts(db, ngqpt, unit)
3192 
3193 
3194 !This section has been created automatically by the script Abilint (TD).
3195 !Do not modify the following lines by hand.
3196 #undef ABI_FUNC
3197 #define ABI_FUNC 'dvdb_list_perts'
3198 !End of the abilint section
3199 
3200  implicit none
3201 
3202 !Arguments ------------------------------------
3203  type(dvdb_t),target,intent(in) :: db
3204  integer,optional,intent(in) :: unit
3205 !arrays
3206  integer,intent(in) :: ngqpt(3)
3207 
3208 !Local variables-------------------------------
3209 !scalars
3210  integer :: tot_miss,tot_weird,miss_q,idir,ipert,iv1,psy,weird_q
3211  integer :: iq_ibz,nqibz,iq_file,qptopt,nshiftq,ii,timerev_q,unt,nqbz
3212  character(len=500) :: msg,ptype,found
3213  type(crystal_t),pointer :: cryst
3214 !arrays
3215  integer :: rfdir(3),qptrlatt(3,3)
3216  integer,allocatable :: pertsy(:,:),symq(:,:,:),rfpert(:)
3217  real(dp) :: qq(3),shiftq(3,1)
3218  real(dp),allocatable :: qibz(:,:),wtq(:),qbz(:,:)
3219 
3220 ! *************************************************************************
3221 
3222  unt = std_out; if (present(unit)) unt = unit
3223  cryst => db%cryst
3224 
3225  if (all(ngqpt == -1)) then
3226    ! Will test the q-points in db
3227    call alloc_copy(db%qpts, qibz)
3228    nqibz = db%nqpt
3229  else
3230    ! Will test the q-points in the IBZ associated to ngqpt hence build IBZ and BZ from ngqpt.
3231    qptopt = 1; shiftq = zero; nshiftq = 1; qptrlatt = 0
3232    do ii=1,3
3233      qptrlatt(ii, ii) = ngqpt(ii)
3234    end do
3235 
3236    call kpts_ibz_from_kptrlatt(cryst, qptrlatt, qptopt, nshiftq, shiftq, &
3237      nqibz, qibz, wtq, nqbz, qbz)
3238 
3239    ABI_FREE(qbz)
3240    ABI_FREE(wtq)
3241  end if
3242 
3243  ! Initialize the list of perturbations rfpert and rdfir
3244  ! WARNING: Only phonon perturbations are considered for the time being.
3245  ABI_MALLOC(rfpert,(db%mpert))
3246  rfpert = 0; rfpert(1:cryst%natom) = 1; rfdir = 1
3247 
3248  ABI_MALLOC(symq, (4,2,cryst%nsym))
3249  ABI_MALLOC(pertsy, (3,db%mpert))
3250 
3251  ! Loop over the q-points in the IBZ and test whether the q-point is present
3252  ! and if all the independent perturbations are available.
3253  !   `tot_miss` is the number of irreducible perturbations not found in the DVDB (critical)
3254  !   `tot_weird` is the number of redundant perturbations found in the DVDB (not critical)
3255  tot_miss = 0; tot_weird = 0
3256  do iq_ibz=1,nqibz
3257    qq = qibz(:,iq_ibz)
3258    iq_file = dvdb_findq(db, qq)
3259 
3260    ! Examine the symmetries of the q wavevector
3261    call littlegroup_q(cryst%nsym,qq,symq,cryst%symrec,cryst%symafm,timerev_q,prtvol=db%prtvol)
3262 
3263    ! Determine the symmetrical perturbations. Meaning of pertsy:
3264    !    0 for non-target perturbations
3265    !    1 for basis perturbations
3266    !   -1 for perturbations that can be found from basis perturbations
3267    call irreducible_set_pert(cryst%indsym,db%mpert,cryst%natom,cryst%nsym,&
3268      pertsy,rfdir,rfpert,symq,cryst%symrec,cryst%symrel)
3269 
3270    if (iq_file /= -1) then
3271      ! This q-point is in the DVDB. Test if all the independent perturbations are available.
3272      call wrtout(unt, sjoin("qpoint:", ktoa(qq), "is present in the DVDB file"))
3273      call wrtout(unt,' The list of irreducible perturbations for this q vector is:')
3274      ii = 0; weird_q = 0; miss_q = 0
3275      do ipert=1,db%mpert
3276        do idir=1,3
3277          psy = pertsy(idir,ipert)
3278          if (psy == 0) cycle
3279          iv1 = db%pos_dpq(idir,ipert,iq_file)
3280          ptype = "independent"; if (psy == -1) ptype = "symmetric"
3281          found = "Yes"; if (iv1 == 0) found = "No"
3282 
3283          if (psy == 1 .and. iv1 == 0) miss_q = miss_q + 1
3284          if (psy == -1 .and. iv1 /= 0) weird_q = weird_q + 1
3285 
3286          ii=ii+1
3287          write(msg,'(i5,a,i2,a,i4,4a)')ii,')  idir=',idir,', ipert=',ipert,", type=",trim(ptype),", found=",trim(found)
3288          call wrtout(unt, msg)
3289        end do
3290      end do
3291 
3292      if (weird_q /= 0) then
3293        write(msg,"(a,i0,a)")"DVDB is overcomplete. ",weird_q, " perturbation(s) can be reconstructed by symmetry."
3294        call wrtout(unt, msg)
3295      end if
3296 
3297      tot_weird = tot_weird + weird_q
3298      tot_miss = tot_miss + miss_q
3299      if (miss_q /=0) then
3300        call wrtout(unt, sjoin("WARNING:", itoa(miss_q), "independent perturbation(s) are missing!."))
3301      end if
3302 
3303    else
3304      ! This q-point is not present in dvdb. Print the list of independent perturbations.
3305      call wrtout(unt, sjoin("qpoint:", ktoa(qq), "is NOT present in the DVDB file"))
3306      call wrtout(unt,' The list of irreducible perturbations for this q vector is:')
3307      ii = 0
3308      do ipert=1,db%mpert
3309        do idir=1,3
3310          if (pertsy(idir,ipert) == 1) then
3311            ii=ii+1
3312            write(msg,'(i5,a,i2,a,i4,a)')ii,')  idir=',idir,', ipert=',ipert,", type=independent, found=No"
3313            call wrtout(unt, msg)
3314            tot_miss = tot_miss + 1
3315          end if
3316        end do
3317      end do
3318    end if
3319 
3320    call wrtout(unt," ")
3321  end do ! iq_ibz
3322 
3323  if (tot_miss /= 0) then
3324     call wrtout(unt, sjoin(ch10, "There are ",itoa(tot_miss), "independent perturbations missing!"))
3325  else
3326     call wrtout(unt, "All the independent perturbations are available")
3327     if (tot_weird /= 0) then
3328       call wrtout(unt, "Note however that the DVDB is overcomplete as symmetric perturbations are present.")
3329     end if
3330  end if
3331 
3332  ABI_FREE(qibz)
3333  ABI_FREE(rfpert)
3334  ABI_FREE(symq)
3335  ABI_FREE(pertsy)
3336 
3337 end subroutine dvdb_list_perts

m_dvdb/dvdb_merge_files [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  dvdb_merge_files

FUNCTION

  Merge a list of POT1 files.

 INPUT
  nfiles=Number of files to be merged.
  dvdb_path=Name of output DVDB file.
  prtvol=Verbosity level.

SIDE EFFECTS

   v1files=List of file names to merge. This list could be changed if POT1 files in netcdf format are found.

PARENTS

      mrgdv

CHILDREN

      cwtime,dvdb_get_v1scf_qpt,dvdb_get_v1scf_rpt,dvdb_open_read
      hdr_fort_read,hdr_fort_write,hdr_free,irreducible_set_pert
      kpts_ibz_from_kptrlatt,littlegroup_q,wrtout,xmpi_barrier

SOURCE

3367 subroutine dvdb_merge_files(nfiles, v1files, dvdb_path, prtvol)
3368 
3369 
3370 !This section has been created automatically by the script Abilint (TD).
3371 !Do not modify the following lines by hand.
3372 #undef ABI_FUNC
3373 #define ABI_FUNC 'dvdb_merge_files'
3374 !End of the abilint section
3375 
3376  implicit none
3377 
3378 !Arguments ------------------------------------
3379 !scalars
3380  integer,intent(in) :: nfiles,prtvol
3381  character(len=*),intent(in) :: dvdb_path
3382  character(len=*),intent(inout) :: v1files(nfiles)
3383 
3384 !Local variables-------------------------------
3385 !scalars
3386 ! Here I made a mistake because 102 corresponds to GS potentials
3387 ! as a consequence DVDB files generated with version <= 8.1.6
3388 ! contain list of potentials with fform = 102.
3389  !integer :: fform_pot=102
3390  integer :: fform_pot=111
3391  integer :: ii,jj,fform,ount,cplex,nfft,ifft,ispden,nperts
3392  integer :: n1,n2,n3,v1_varid,ierr
3393  logical :: qeq0
3394  character(len=500) :: msg
3395  type(hdr_type),pointer :: hdr1
3396  type(dvdb_t) :: dvdb
3397 !arrays
3398  integer :: units(nfiles)
3399  real(dp) :: rhog1_g0(2)
3400  real(dp),allocatable :: v1(:)
3401  logical :: has_rhog1_g0(nfiles)
3402  type(hdr_type),target,allocatable :: hdr1_list(:)
3403 
3404 !************************************************************************
3405 
3406  if (file_exists(dvdb_path)) then
3407    MSG_ERROR(sjoin("Cannot overwrite existing file:", dvdb_path))
3408  end if
3409 
3410  ! If a file is not found, try the netcdf version and change v1files accordingly.
3411  do ii=1,nfiles
3412    if (nctk_try_fort_or_ncfile(v1files(ii), msg) /= 0) then
3413      MSG_ERROR(msg)
3414    end if
3415  end do
3416 
3417  ! Read the headers
3418  ABI_DT_MALLOC(hdr1_list, (nfiles))
3419 
3420  nperts = size(hdr1_list)
3421 
3422  ! Write dvdb file (we only support fortran binary format)
3423  if (open_file(dvdb_path, msg, newunit=ount, form="unformatted", action="write", status="unknown") /= 0) then
3424    MSG_ERROR(msg)
3425  end if
3426  write(ount, err=10, iomsg=msg) dvdb_last_version
3427  write(ount, err=10, iomsg=msg) nperts
3428 
3429  ! Validate headers.
3430  ! TODO: Should perform consistency check on the headers, rearrange them in blocks of q-points.
3431  ! ignore POT1 files that do not correspond to atomic perturbations.
3432 
3433  do ii=1,nfiles
3434    write(std_out,"(a,i0,2a)")"- Reading header of file [",ii,"]: ",trim(v1files(ii))
3435 
3436    if (endswith(v1files(ii), ".nc")) then
3437 #ifdef HAVE_NETCDF
3438       NCF_CHECK(nctk_open_read(units(ii), v1files(ii), xmpi_comm_self))
3439       call hdr_ncread(hdr1_list(ii),units(ii),fform)
3440 #endif
3441    else
3442      if (open_file(v1files(ii), msg, newunit=units(ii), form="unformatted", action="read", status="old") /= 0) then
3443        MSG_ERROR(msg)
3444      end if
3445      call hdr_fort_read(hdr1_list(ii), units(ii), fform)
3446    end if
3447 
3448    if (dvdb_check_fform(fform, "merge_dvdb", msg) /= 0) then
3449      MSG_ERROR(sjoin("While reading:", v1files(ii), msg))
3450    end if
3451    if (prtvol > 0) call hdr_echo(hdr1_list(ii), fform, 3, unit=std_out)
3452    if (hdr1_list(ii)%pertcase == 0) then
3453      MSG_ERROR(sjoin("Found GS potential:", v1files(ii)))
3454    end if
3455    !write(std_out,*)"done", trim(v1files(ii))
3456 
3457    ! Supported fform:
3458    ! 109  POT1 files without vh1(G=0)
3459    ! 111  POT1 files with extra record with vh1(G=0) after FFT data.
3460    has_rhog1_g0(ii) = .True.
3461    if (fform == 109) has_rhog1_g0(ii) = .False.
3462 
3463    write(std_out,"(a,i0,2a)")"- Merging file [",ii,"]: ",trim(v1files(ii))
3464    jj = ii
3465    hdr1 => hdr1_list(jj)
3466    call hdr_fort_write(hdr1, ount, fform_pot, ierr)
3467    ABI_CHECK(ierr == 0, "hdr_fort_write returned ierr = 0")
3468 
3469    qeq0 = (hdr1%qptn(1)**2+hdr1%qptn(2)**2+hdr1%qptn(3)**2<1.d-14)
3470    cplex = 2; if (qeq0) cplex = 1
3471    nfft = product(hdr1%ngfft(1:3))
3472    n1 = hdr1%ngfft(1); n2 = hdr1%ngfft(2); n3 = hdr1%ngfft(3)
3473 
3474    ABI_MALLOC(v1, (cplex*nfft))
3475 
3476    if (.not. endswith(v1files(ii), ".nc")) then
3477       ! Fortran IO
3478       do ispden=1,hdr1%nspden
3479         read(units(jj), err=10, iomsg=msg) (v1(ifft), ifft=1,cplex*nfft)
3480         write(ount, err=10, iomsg=msg) (v1(ifft), ifft=1,cplex*nfft)
3481       end do
3482       ! Add rhog1(G=0)
3483       rhog1_g0 = zero
3484       if (has_rhog1_g0(jj)) read(units(jj), err=10, iomsg=msg) rhog1_g0
3485       if (dvdb_last_version > 1) write(ount, err=10, iomsg=msg) rhog1_g0
3486    else
3487 #ifdef HAVE_NETCDF
3488       ! Netcdf IO
3489       ! netcdf array has shape [cplex, n1, n2, n3, nspden]
3490       NCF_CHECK(nf90_inq_varid(units(ii), "first_order_potential", v1_varid))
3491       do ispden=1,hdr1%nspden
3492         NCF_CHECK(nf90_get_var(units(ii), v1_varid, v1, start=[1,1,1,ispden], count=[cplex, n1, n2, n3, 1]))
3493         write(ount, err=10, iomsg=msg) (v1(ifft), ifft=1,cplex*nfft)
3494       end do
3495       ! Add rhog1(G=0)
3496       rhog1_g0 = zero
3497       if (has_rhog1_g0(jj)) then
3498         NCF_CHECK(nf90_get_var(units(ii), nctk_idname(units(ii), "rhog1_g0"), rhog1_g0))
3499       end if
3500       if (dvdb_last_version > 1) write(ount, err=10, iomsg=msg) rhog1_g0
3501 #endif
3502    end if
3503 
3504    if (.not. endswith(v1files(ii), ".nc")) then
3505      close(units(ii))
3506    else
3507 #ifdef HAVE_NETCDF
3508      NCF_CHECK(nf90_close(units(ii)))
3509 #endif
3510    end if
3511 
3512    ABI_FREE(v1)
3513  end do ! nperts
3514 
3515  close(ount)
3516 
3517  do ii=1,size(hdr1_list)
3518    call hdr_free(hdr1_list(ii))
3519  end do
3520  ABI_DT_FREE(hdr1_list)
3521 
3522  write(std_out,"(a,i0,a)")"Merged successfully ",nfiles," files"
3523 
3524  ! List available perturbations.
3525  call dvdb_init(dvdb, dvdb_path, xmpi_comm_self)
3526  call dvdb_print(dvdb)
3527  call dvdb_list_perts(dvdb, [-1,-1,-1])
3528  call dvdb_free(dvdb)
3529 
3530  return
3531 
3532  ! Handle Fortran IO error
3533 10 continue
3534  MSG_ERROR(sjoin("Error while merging files", ch10, msg))
3535 
3536 end subroutine dvdb_merge_files

m_dvdb/dvdb_open_read [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  dvdb_open_read

FUNCTION

  Open the file in read-only mode.

INPUTS

   ngfft(18)=Info on the FFT mesh used for the DFPT potentials. Note that ngfft
     is the mesh used by the parent. In principle, it can differ from the one
     found in the file. In this case a Fourier interpolation is required.
   comm=MPI communicator

PARENTS

      m_dvdb,m_gkk,m_phgamma,m_phpi,m_sigmaph

CHILDREN

      cwtime,dvdb_get_v1scf_qpt,dvdb_get_v1scf_rpt,dvdb_open_read
      hdr_fort_read,hdr_fort_write,hdr_free,irreducible_set_pert
      kpts_ibz_from_kptrlatt,littlegroup_q,wrtout,xmpi_barrier

SOURCE

503 subroutine dvdb_open_read(db, ngfft, comm)
504 
505 
506 !This section has been created automatically by the script Abilint (TD).
507 !Do not modify the following lines by hand.
508 #undef ABI_FUNC
509 #define ABI_FUNC 'dvdb_open_read'
510 !End of the abilint section
511 
512  implicit none
513 
514 !Arguments ------------------------------------
515 !scalars
516  integer,intent(in) :: comm
517  type(dvdb_t),intent(inout) :: db
518 !arrays
519  integer,intent(in) :: ngfft(18)
520 
521 !Local variables-------------------------------
522 !scalars
523  integer :: nprocs
524  character(len=500) :: msg
525 
526 !************************************************************************
527 
528  if (db%rw_mode /= DVDB_NOMODE) then
529    MSG_ERROR("DVDB should be in DVDB_NOMODE when open_read is called.")
530  end if
531  db%rw_mode = DVDB_READMODE
532 
533  nprocs = xmpi_comm_size(comm)
534 
535  ! Initialize tables to call fourdp in sequential
536  db%ngfft = ngfft
537  call init_distribfft_seq(db%mpi_enreg%distribfft,'c',ngfft(2),ngfft(3),'all')
538  call init_distribfft_seq(db%mpi_enreg%distribfft,'f',ngfft(2),ngfft(3),'all')
539 
540  ! Open the file.
541  select case (db%iomode)
542  case (IO_MODE_FORTRAN)
543    if (open_file(db%path, msg, newunit=db%fh, form="unformatted", status="old", action="read") /= 0) then
544      MSG_ERROR(msg)
545    end if
546    read(db%fh, err=10, iomsg=msg)
547    read(db%fh, err=10, iomsg=msg)
548    db%current_fpos = 1
549 
550  case (IO_MODE_MPI)
551    MSG_ERROR("MPI not coded")
552 
553  case default
554    MSG_ERROR(sjoin("Unsupported iomode:", itoa(db%iomode)))
555  end select
556 
557  return
558 
559  ! Handle Fortran IO error
560 10 continue
561  MSG_ERROR(sjoin("Error while reading", db%path, ch10, msg))
562 
563 end subroutine dvdb_open_read

m_dvdb/dvdb_print [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  dvdb_print

FUNCTION

  Print info on the object.

INPUTS

 [unit]=the unit number for output
 [prtvol]=verbosity level
 [mode_paral]=either "COLL" or "PERS"

OUTPUT

  Only printing.

PARENTS

      eph,m_dvdb,m_sigmaph,mrgdv

CHILDREN

      cwtime,dvdb_get_v1scf_qpt,dvdb_get_v1scf_rpt,dvdb_open_read
      hdr_fort_read,hdr_fort_write,hdr_free,irreducible_set_pert
      kpts_ibz_from_kptrlatt,littlegroup_q,wrtout,xmpi_barrier

SOURCE

718 subroutine dvdb_print(db, header, unit, prtvol, mode_paral)
719 
720 
721 !This section has been created automatically by the script Abilint (TD).
722 !Do not modify the following lines by hand.
723 #undef ABI_FUNC
724 #define ABI_FUNC 'dvdb_print'
725 !End of the abilint section
726 
727  implicit none
728 
729 !Arguments ------------------------------------
730 !scalars
731  integer,optional,intent(in) :: prtvol,unit
732  character(len=4),optional,intent(in) :: mode_paral
733  character(len=*),optional,intent(in) :: header
734  type(dvdb_t),intent(in) :: db
735 
736 !Local variables-------------------------------
737 !scalars
738  integer :: my_unt,my_prtvol,iv1,iq,idir,ipert,iatom
739  character(len=4) :: my_mode
740  character(len=500) :: msg
741 
742 ! *************************************************************************
743 
744  my_unt = std_out; if (PRESENT(unit)) my_unt = unit
745  my_prtvol = 0   ; if (PRESENT(prtvol)) my_prtvol = prtvol
746  my_mode = 'COLL'; if (PRESENT(mode_paral)) my_mode = mode_paral
747 
748  msg=' ==== Info on the dvdb% object ==== '
749  if (PRESENT(header)) msg=' ==== '//TRIM(ADJUSTL(header))//' ==== '
750  call wrtout(my_unt,msg,my_mode)
751 
752  write(std_out,"(a)")sjoin("DVDB version:", itoa(db%version))
753  write(std_out,"(a)")sjoin("File path:", db%path)
754  write(std_out,"(a)")sjoin("Number of v1scf potentials:", itoa(db%numv1))
755  write(std_out,"(a)")sjoin("Number of q-points in DVDB: ", itoa(db%nqpt))
756  write(std_out,"(a)")sjoin("Activate symmetrization of v1scf(r):", yesno(db%symv1))
757  write(std_out,"(a)")"List of q-points: min(10, nqpt)"
758  do iq=1,min(db%nqpt, 10)
759    write(std_out,"(a)")sjoin("[", itoa(iq),"]", ktoa(db%qpts(:,iq)))
760  end do
761  if (db%nqpt > 10) write(std_out,"(a)")"..."
762 
763  write(std_out,"(a)")sjoin("Have dielectric tensor and Born effective charges:", yesno(db%has_dielt_zeff))
764  if (db%has_dielt_zeff) then
765    write(std_out, '(a,3(/,3es16.6))') ' Dielectric Tensor:', &
766      db%dielt(1,1), db%dielt(1,2), db%dielt(1,3), &
767      db%dielt(2,1), db%dielt(2,2), db%dielt(2,3), &
768      db%dielt(3,1), db%dielt(3,2), db%dielt(3,3)
769    write(std_out, '(a)') ' Effectives Charges: '
770    do iatom=1,db%natom
771      write(std_out,'(a,i0,3(/,3es16.6))')' iatom: ',iatom, &
772        db%zeff(1,1,iatom), db%zeff(1,2,iatom), db%zeff(1,3,iatom), &
773        db%zeff(2,1,iatom), db%zeff(2,2,iatom), db%zeff(2,3,iatom), &
774        db%zeff(3,1,iatom), db%zeff(3,2,iatom), db%zeff(3,3,iatom)
775    end do
776  end if
777 
778  if (my_prtvol > 0) then
779    call crystal_print(db%cryst, header="Crystal structure in DVDB file")
780    write(std_out,"(a)")"FFT mesh for potentials on file:"
781    write(std_out,"(a)")"q-point, idir, ipert, ngfft(:3)"
782    do iv1=1,db%numv1
783      idir = db%iv_pinfoq(1, iv1); ipert = db%iv_pinfoq(2, iv1); iq = db%iv_pinfoq(4, iv1)
784      write(std_out,"(a)")sjoin(ktoa(db%qpts(:,iq)), itoa(idir), itoa(ipert), ltoa(db%ngfft3_v1(:,iv1)))
785    end do
786 
787    write(std_out, "(a)")"q-point, idir, ipert, rhog1(q,G=0)"
788    do iv1=1,db%numv1
789      idir = db%iv_pinfoq(1, iv1); ipert = db%iv_pinfoq(2, iv1); iq = db%iv_pinfoq(4, iv1)
790      write(std_out,"(a)")sjoin(ktoa(db%qpts(:,iq)), itoa(idir), itoa(ipert), &
791        ftoa(db%rhog1_g0(1, iv1)), ftoa(db%rhog1_g0(2, iv1)))
792    end do
793  end if
794 
795 end subroutine dvdb_print

m_dvdb/dvdb_read_onev1 [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  dvdb_read_onev1

FUNCTION

  Read the DFPT potential for the specified (idir, ipert, iqpt).
  Note that iqpt is the index in dvdb%qpts. Use dvdb_findq to
  get the index from the q-point in reduced coordinates.

INPUTS

  idir=Direction of the perturbation
  ipert=Perturbation type.
  iqpt=Index of the q-point in dvdb%qpts
  cplex=1 if real, 2 if complex potentials.
  nfft=Number of fft-points treated by this processors
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft

OUTPUT

  ierr=Non-zero if error.
  v1scf(cplex*nfft, nspden)=DFT potential associated to (idir, ipert, iqpt).
  msg=String with error message if ierr /= 0.

PARENTS

CHILDREN

SOURCE

896 integer function dvdb_read_onev1(db, idir, ipert, iqpt, cplex, nfft, ngfft, v1scf, msg) result(ierr)
897 
898 
899 !This section has been created automatically by the script Abilint (TD).
900 !Do not modify the following lines by hand.
901 #undef ABI_FUNC
902 #define ABI_FUNC 'dvdb_read_onev1'
903 !End of the abilint section
904 
905  implicit none
906 
907 !Arguments ------------------------------------
908 !scalars
909  integer,intent(in) :: idir,ipert,iqpt,cplex,nfft
910  character(len=*),intent(out) :: msg
911  type(dvdb_t),intent(inout) :: db
912 !arrays
913  integer,intent(in) :: ngfft(18)
914  real(dp),intent(out) :: v1scf(cplex*nfft,db%nspden)
915 
916 !Local variables-------------------------------
917 !scalars
918  integer,parameter :: paral_kgb0=0
919  integer :: iv1,ispden,nfftot_file,nfftot_out,ifft
920  type(MPI_type) :: MPI_enreg_seq
921 !arrays
922  integer :: ngfft_in(18),ngfft_out(18)
923  real(dp),allocatable :: v1r_file(:,:),v1g_in(:,:),v1g_out(:,:)
924 
925 ! *************************************************************************
926 
927  ! Consistency checks
928  ierr = 1
929  iv1 = db%pos_dpq(idir,ipert,iqpt)
930 
931  if (iv1 == 0) then
932    write(msg,"(3(a,i0))")"Cannot find idir=",idir,", ipert=",ipert,", iqpt=",iqpt
933    return
934  end if
935 
936  if (cplex /= db%cplex_v1(iv1)) then
937    write(msg,"(2(a,i0))")"Wrong cplex. Expecting=",db%cplex_v1(iv1),", received= ",cplex
938    return
939  end if
940 
941  ! Find (idir, ipert, iqpt) and skip the header.
942  call dvdb_seek(db, idir, ipert, iqpt)
943  ierr = my_hdr_skip(db%fh, idir, ipert, db%qpts(:,iqpt), msg)
944  if (ierr /= 0) return
945 
946  ! Read v1 from file.
947  nfftot_out = product(ngfft(:3)); nfftot_file = product(db%ngfft3_v1(:3, iv1))
948 
949  if (all(ngfft(:3) == db%ngfft3_v1(:3, iv1))) then
950    do ispden=1,db%nspden
951      read(db%fh, err=10, iomsg=msg) (v1scf(ifft, ispden), ifft=1,cplex*nfftot_file)
952    end do
953  else
954    ! The FFT mesh used in the caller differ from the one found in the DBDB --> Fourier interpolation
955    ! TODO: Add linear interpolation as well.
956    !MSG_ERROR("FFT interpolation of DFPT potentials must be tested.")
957    ABI_MALLOC(v1r_file, (cplex*nfftot_file, db%nspden))
958    do ispden=1,db%nspden
959      read(db%fh, err=10, iomsg=msg) (v1r_file(ifft, ispden), ifft=1,cplex*nfftot_file)
960    end do
961 
962    ! Call fourier_interpol to get v1scf on ngfft mesh.
963    ngfft_in = ngfft; ngfft_out = ngfft
964    ngfft_in(1:3) = db%ngfft3_v1(1:3, iv1); ngfft_out(1:3) = ngfft(1:3)
965    ngfft_in(4:6) = ngfft_in(1:3); ngfft_out(4:6) = ngfft_out(1:3)
966    ngfft_in(9:18) = 0; ngfft_out(9:18) = 0
967    ngfft_in(10) = 1; ngfft_out(10) = 1
968 
969    call initmpi_seq(MPI_enreg_seq)
970    ! Which one is coarse? Note that this part is not very robust and can fail!
971    if (ngfft_in(2) * ngfft_in(3) < ngfft_out(2) * ngfft_out(3)) then
972      call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfft_in(2),ngfft_in(3),'all')
973      call init_distribfft_seq(MPI_enreg_seq%distribfft,'f',ngfft_out(2),ngfft_out(3),'all')
974    else
975      call init_distribfft_seq(MPI_enreg_seq%distribfft,'f',ngfft_in(2),ngfft_in(3),'all')
976      call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfft_out(2),ngfft_out(3),'all')
977    end if
978 
979    ABI_MALLOC(v1g_in,  (2, nfftot_file))
980    ABI_MALLOC(v1g_out, (2, nfftot_out))
981 
982    call fourier_interpol(cplex,db%nspden,0,0,nfftot_file,ngfft_in,nfftot_out,ngfft_out,&
983      paral_kgb0,MPI_enreg_seq,v1r_file,v1scf,v1g_in,v1g_out)
984 
985    ABI_FREE(v1g_in)
986    ABI_FREE(v1g_out)
987    ABI_FREE(v1r_file)
988    call destroy_mpi_enreg(MPI_enreg_seq)
989  end if
990 
991  return
992 
993  ! Handle Fortran IO error
994 10 continue
995  ierr = 1
996  msg = sjoin("Error while reading", db%path, ch10, msg)
997 
998 end function dvdb_read_onev1

m_dvdb/dvdb_readsym_allv1 [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  dvdb_readsym_allv1

FUNCTION

  Read all 3*natom DFPT potentials for the given iqpt (only atomic perturbations).

  The routine will:

     1) Reconstruct the potentials by symmetry if the DVDB contains less than 3*natom potentials.
     2) interpolate the data if the input FFT mesh defined by `ngfft` differs
        from the one used to store data in the file.

  Note that iqpt is the index in dvdb%qpts. Use dvdb_findq to
  get the index from the q-point in reduced coordinates.

INPUTS

  iqpt=Index of the q-point in dvdb%qpts
  nfft=Number of fft-points treated by this processors
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  comm=MPI communicator

OUTPUT

  cplex=1 if real, 2 if complex.
  v1scf(cplex, nfft, nspden, 3*natom)= v1scf potentials on the real-space FFT mesh for the 3*natom perturbations.

PARENTS

      m_dvdb,m_gkk,m_phgamma,m_phpi,m_sigmaph

CHILDREN

      cwtime,dvdb_get_v1scf_qpt,dvdb_get_v1scf_rpt,dvdb_open_read
      hdr_fort_read,hdr_fort_write,hdr_free,irreducible_set_pert
      kpts_ibz_from_kptrlatt,littlegroup_q,wrtout,xmpi_barrier

SOURCE

1039 subroutine dvdb_readsym_allv1(db, iqpt, cplex, nfft, ngfft, v1scf, comm)
1040 
1041 
1042 !This section has been created automatically by the script Abilint (TD).
1043 !Do not modify the following lines by hand.
1044 #undef ABI_FUNC
1045 #define ABI_FUNC 'dvdb_readsym_allv1'
1046 !End of the abilint section
1047 
1048  implicit none
1049 
1050 !Arguments ------------------------------------
1051 !scalars
1052  integer,intent(in) :: iqpt,nfft,comm
1053  integer,intent(out) :: cplex
1054  type(dvdb_t),target,intent(inout) :: db
1055 !arrays
1056  integer,intent(in) :: ngfft(18)
1057  real(dp),allocatable,intent(out) :: v1scf(:,:,:,:)
1058 
1059 !Local variables-------------------------------
1060 !scalars
1061  integer,parameter :: master=0
1062  integer :: ipc,npc,idir,ipert,pcase,my_rank,nproc,ierr,mu
1063  character(len=500) :: msg
1064 !arrays
1065  integer :: pinfo(3,3*db%mpert),pflag(3, db%natom)
1066 
1067 ! *************************************************************************
1068 
1069  my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
1070 
1071  ! Get number of perturbations computed for this iqpt as well as cplex.
1072  npc = dvdb_get_pinfo(db, iqpt, cplex, pinfo)
1073  ABI_CHECK(npc /= 0, "npc == 0!")
1074 
1075  ABI_STAT_MALLOC(v1scf, (cplex,nfft,db%nspden,3*db%natom), ierr)
1076  ABI_CHECK(ierr==0, "OOM in v1scf")
1077 
1078  ! Master read all available perturbations and broadcasts data.
1079  if (my_rank == master) then
1080    do ipc=1,npc
1081      idir = pinfo(1,ipc); ipert = pinfo(2,ipc); pcase = pinfo(3, ipc)
1082      if (dvdb_read_onev1(db, idir, ipert, iqpt, cplex, nfft, ngfft, v1scf(:,:,:,pcase), msg) /= 0) then
1083        MSG_ERROR(msg)
1084      end if
1085    end do
1086  end if
1087  call xmpi_bcast(v1scf, master, comm, ierr)
1088 
1089  ! Return if all perts are available.
1090  if (npc == 3*db%natom) then
1091    if (db%symv1) then
1092      if (db%debug) write(std_out,*)"Potentials are available but will call v1phq_symmetrize because of symv1"
1093      do mu=1,db%natom3
1094        idir = mod(mu-1, 3) + 1; ipert = (mu - idir) / 3 + 1
1095        call v1phq_symmetrize(db%cryst,idir,ipert,db%symq_table(:,:,:,iqpt),ngfft,cplex,nfft,&
1096          db%nspden,db%nsppol,db%mpi_enreg,v1scf(:,:,:,mu))
1097      end do
1098    end if
1099    if (db%debug) write(std_out,*)ABI_FUNC,": All perts available. Returning"
1100    return
1101  end if
1102 
1103  ! Perturbation are missing and we have to reconstruct them by symmetry.
1104  ! This is the common case when DFPT calculations are done for independent perturbations only.
1105  if (db%debug) then
1106    write(std_out,*)sjoin("Will use symmetries to recostruct:", itoa(3*db%natom - npc), "perturbations")
1107  end if
1108 
1109  ! 0 if pert is not available.
1110  ! 1 if pert is on file.
1111  ! 2 if pert has been reconstructed by symmetry.
1112  pflag = 0
1113  do ipc=1,npc
1114    pflag(pinfo(1,ipc), pinfo(2,ipc)) = 1
1115  end do
1116 
1117  call v1phq_complete(db%cryst,db%qpts(:,iqpt),ngfft,cplex,nfft,db%nspden,db%nsppol,db%mpi_enreg,db%symv1,pflag,v1scf)
1118 
1119 end subroutine dvdb_readsym_allv1

m_dvdb/dvdb_rewind [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  dvdb_rewind

FUNCTION

   Rewind the file and move to the first header. Needed only if dvdb%iomode==IO_MODE_FORTRAN
   Return exit code and error message in msg if ierr != 0

PARENTS

      m_dvdb

CHILDREN

SOURCE

3065 integer function dvdb_rewind(db, msg) result(ierr)
3066 
3067 
3068 !This section has been created automatically by the script Abilint (TD).
3069 !Do not modify the following lines by hand.
3070 #undef ABI_FUNC
3071 #define ABI_FUNC 'dvdb_rewind'
3072 !End of the abilint section
3073 
3074  implicit none
3075 
3076 !Arguments ------------------------------------
3077  type(dvdb_t),intent(inout) :: db
3078  character(len=*),intent(out) :: msg
3079 
3080 ! *************************************************************************
3081 
3082  ierr = 0
3083  if (db%iomode == IO_MODE_FORTRAN) then
3084    rewind(db%fh, err=10, iomsg=msg)
3085    read(db%fh, err=10, iomsg=msg)  ! version
3086    read(db%fh, err=10, iomsg=msg)  ! numv1
3087    db%current_fpos = 1
3088 
3089  else
3090    ierr = -1
3091    msg = "should not be called when iomode /= IO_MODE_FORTRAN"
3092  end if
3093 
3094  return
3095 
3096  ! Handle Fortran IO error
3097 10 continue
3098  ierr = 1
3099  msg = sjoin("Error while reading", db%path, ch10, msg)
3100 
3101 end function dvdb_rewind

m_dvdb/dvdb_seek [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  dvdb_seek

FUNCTION

   Move the internal file pointer so that it points to the
   block with (idir, ipert, iqpt). Needed only if dvdb%iomode==IO_MODE_FORTRAN

INPUTS

   idir,ipert,iqpt = (direction, perturbation, q-point) indices

SIDE EFFECTS

   db<type(dvdb_t)> : modifies db%f90_fptr and the internal F90 file pointer.

PARENTS

      m_dvdb

CHILDREN

      cwtime,dvdb_get_v1scf_qpt,dvdb_get_v1scf_rpt,dvdb_open_read
      hdr_fort_read,hdr_fort_write,hdr_free,irreducible_set_pert
      kpts_ibz_from_kptrlatt,littlegroup_q,wrtout,xmpi_barrier

SOURCE

2990 subroutine dvdb_seek(db, idir, ipert, iqpt)
2991 
2992 
2993 !This section has been created automatically by the script Abilint (TD).
2994 !Do not modify the following lines by hand.
2995 #undef ABI_FUNC
2996 #define ABI_FUNC 'dvdb_seek'
2997 !End of the abilint section
2998 
2999  implicit none
3000 
3001 !Arguments ------------------------------------
3002  integer,intent(in)  :: idir, ipert, iqpt
3003  type(dvdb_t),intent(inout) :: db
3004 
3005 !Local variables-------------------------------
3006  integer :: pos_wanted,ii,ispden
3007  real(dp),parameter :: fake_qpt(3)=zero
3008  character(len=500) :: msg
3009 
3010 ! *************************************************************************
3011 
3012  if (db%iomode == IO_MODE_FORTRAN) then
3013    ! Numb code: rewind the file and read it from the beginning
3014    ! TODO: Here I need a routine to backspace the hdr!
3015    if (dvdb_rewind(db, msg) /= 0) then
3016      MSG_ERROR(msg)
3017    end if
3018    pos_wanted = db%pos_dpq(idir,ipert,iqpt)
3019 
3020    do ii=1,pos_wanted-1
3021      !write(std_out,*)"in seek with ii: ",ii,"pos_wanted: ",pos_wanted
3022      if (my_hdr_skip(db%fh, -1, -1, fake_qpt, msg) /= 0) then
3023        MSG_ERROR(msg)
3024      end if
3025      ! Skip the records with v1.
3026      do ispden=1,db%nspden
3027        read(db%fh, err=10, iomsg=msg)
3028      end do
3029      ! Skip record with rhog1_g0 (if present)
3030      if (db%version > 1) read(db%fh, err=10, iomsg=msg)
3031    end do
3032    db%current_fpos = pos_wanted
3033 
3034  else
3035    MSG_ERROR("Should not be called when iomode /= IO_MODE_FORTRAN")
3036  end if
3037 
3038  return
3039 
3040  ! Handle Fortran IO error
3041 10 continue
3042  !ierr = 1
3043  msg = sjoin("Error while reading", db%path, ch10, msg)
3044 
3045 end subroutine dvdb_seek

m_dvdb/dvdb_t [ Types ]

[ Top ] [ m_dvdb ] [ Types ]

NAME

  dvdb_t

FUNCTION

  Database of DFPT results. The database contains `numv1` perturbations
  and the corresponding first order local potential in real space on the FFT mesh.
  Note that one can have different FFT meshes for the different perturbations.

NOTES

  natom, nspden, nspinor, and usepaw are global variables in the sense that it's not possible to add
  new entries to the database if these dimension differ from the global ones.

TODO

  1) Add option to store potentials in memory to reduce IO pressure
  2) handle q_bz = S q_ibz case by symmetrizing the potentials already available in the DVDB.
     without performing FT interpolation.
  3) Implement treatment of long-range behaviour in FT interpolation in polar semiconductors.

SOURCE

104  type,public :: dvdb_t
105 
106   integer :: fh
107    ! file handler
108    !  Fortran unit number if iomode==IO_MODE_FORTRAN
109    !  MPI file handler if iomode==IO_MODE_MPI
110 
111   integer :: comm
112   ! MPI communicator used for IO.
113 
114   integer :: version
115   ! File format version read from file.
116 
117   integer :: iomode=IO_MODE_FORTRAN
118   ! Method used to access the DVDB file:
119   !   IO_MODE_FORTRAN for usual Fortran IO routines
120   !   IO_MODE_MPI if MPI/IO routines.
121 
122   integer :: rw_mode = DVDB_NOMODE
123    ! (Read|Write) mode
124 
125   integer :: current_fpos
126   ! The current position of the file pointer used for sequential access with Fortran-IO
127   !  FPOS_EOF signals the end of file
128 
129   integer :: numv1
130   ! Number of v1 potentials present in file.
131 
132   integer :: nqpt
133   ! Number of q-points
134 
135   integer :: natom
136    ! Number of atoms
137 
138   integer :: natom3
139    ! 3 * natom
140 
141   integer :: nspden
142    ! Number of spin density components
143 
144   integer :: nsppol
145    ! Number of spin polarizations.
146 
147   integer :: nspinor
148    ! Number of spinor components.
149 
150   integer :: usepaw
151    ! 1 if PAW calculation, 0 otherwise
152 
153   integer :: mpert
154    ! Maximum number of perturbations
155 
156   integer :: nrpt=0
157   ! Number of real space points used for Fourier interpolation
158 
159   integer :: prtvol=0
160    ! Verbosity level
161 
162   logical :: debug=.False.
163    ! Debug flag
164 
165   logical :: has_dielt_zeff=.False.
166    ! Does the dvdb have the dielectric tensor and Born effective charges
167 
168   logical :: symv1=.False.
169    ! Activate symmetrization of v1 potentials.
170 
171   character(len=fnlen) :: path = ABI_NOFILE
172    ! File name
173 
174   real(dp) :: dielt(3,3)
175    ! Dielectric tensor
176 
177   integer,allocatable :: pos_dpq(:,:,:)
178    ! pos_dpq(3, mpert, nqpt)
179    ! The position of the (idir, ipert, iqpt) potential in the file (in units for POT1 blocks)
180    ! 0 if the corresponding entry is not available.
181 
182   integer,allocatable :: cplex_v1(:)
183   ! cplex_v1(numv1)
184   ! The value of cplex for each v1(cplex*nfft, nspden) potential
185   ! 2 if the potential is complex, 1 if real (q==Gamma)
186 
187   integer,allocatable :: symq_table(:,:,:,:)
188   ! symq(4,2,nsym,nqpt)
189   ! Table computed by littlegroup_q for all q-points found in the DVDB.
190   !   three first numbers define the G vector;
191   !   fourth number is zero if the q-vector is not preserved, is 1 otherwise
192   !   second index is one without time-reversal symmetry, two with time-reversal symmetry
193 
194   integer :: ngfft(18)=-1
195    ! Info on the FFT to be used for the potentials.
196 
197   integer,allocatable :: iv_pinfoq(:,:)
198    !iv_pinfoq(4, numv1)
199    !  iv_pinfoq(1, iv1) gives the `idir` index of the iv1 potential
200    !  iv_pinfoq(2, iv1) gives the `ipert` index of the iv1 potential
201    !  iv_pinfoq(3, iv1) gives `pertcase`=idir + (ipert-1)*3
202    !  iv_pinfoq(4, iv1) gives the `iqpt` index of the iv1 potential
203 
204   ! FIXME: (3) or (18)
205   integer,allocatable :: ngfft3_v1(:,:)
206    ! ngfft3_v1(3, numv1)
207    ! The FFT mesh used for each v1 potential (the one used to store data in the file).
208 
209   real(dp),allocatable :: qpts(:,:)
210    ! qpts(3,nqpt)
211    ! List of q-points in reduced coordinates.
212 
213   real(dp),allocatable :: rpt(:,:)
214   ! rpt(3,nrpt)
215   ! Real space points for Fourier interpolation.
216 
217   real(dp),allocatable :: v1scf_rpt(:,:,:,:,:)
218   ! Wannier representation (real values)
219   ! v1scf_rpt(nrpt, nfft , nspden, 3*natom)
220 
221   real(dp),allocatable :: rhog1_g0(:,:)
222   ! rhog1_g0(2, numv1)
223   ! G=0 component of rhog1. Used to treat the long range component
224   ! in (polar) semiconductors
225 
226   real(dp),allocatable :: zeff(:,:,:)
227   ! zeff(3,3,natom)
228   ! Effective charge on each atom, versus electric field and atomic displacement.
229 
230   type(crystal_t) :: cryst
231   ! Crystalline structure read from the the DVDB file.
232 
233   type(mpi_type) :: mpi_enreg
234   ! Internal object used to call fourdp
235 
236  end type dvdb_t
237 
238  public :: dvdb_init              ! Initialize the object.
239  public :: dvdb_open_read         ! Open the file in read-only mode.
240  public :: dvdb_close             ! Close file.
241  public :: dvdb_free              ! Release the memory allocated and close the file.
242  public :: dvdb_print             ! Release memory.
243  public :: dvdb_findq             ! Returns the index of the q-point.
244  public :: dvdb_read_onev1        ! Read and return the DFPT potential for given (idir, ipert, iqpt).
245  public :: dvdb_readsym_allv1     ! Read and return all the 3*natom DFPT potentials (either from file or symmetrized)
246  public :: dvdb_list_perts        ! Check if all the (phonon) perts are available taking into account symmetries.
247  public :: dvdb_ftinterp_setup    ! Prepare the internal tables for Fourier interpolation.
248  public :: dvdb_ftinterp_qpt      ! Fourier interpolation of potentials for given q-point
249  public :: dvdb_merge_files       ! Merge a list of POT1 files.
250  public :: dvdb_v1r_long_range    ! Long-range part of the phonon potential
251  public :: dvdb_interpolate_v1scf ! Fourier interpolation of the phonon potentials
252  public :: dvdb_get_v1scf_rpt     ! Fourier transform of the phonon potential from qpt to R
253  public :: dvdb_get_v1scf_qpt     ! Fourier transform of the phonon potential from R to qpt
254  public :: dvdb_interpolate_and_write ! Interpolate the phonon potentials and write a new DVDB file.
255 
256 ! Debugging tools.
257  public :: dvdb_test_v1rsym       ! Check symmetries of the DFPT potentials.
258  public :: dvdb_test_v1complete   ! Debugging tool used to test the symmetrization of the DFPT potentials.
259  public :: dvdb_test_ftinterp     ! Test Fourier interpolation of DFPT potentials.

m_dvdb/dvdb_test_ftinterp [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  dvdb_test_ftinterp

FUNCTION

  Debugging tool used to test the Fourier interpolation of the DFPT potentials.

INPUTS

  db_path=Filename
  ngqpt(3)=Divisions of the Q-mesh reported in the DVDB file
  comm=MPI communicator.

OUTPUT

  Only writing.

PARENTS

      mrgdv

CHILDREN

      cwtime,dvdb_get_v1scf_qpt,dvdb_get_v1scf_rpt,dvdb_open_read
      hdr_fort_read,hdr_fort_write,hdr_free,irreducible_set_pert
      kpts_ibz_from_kptrlatt,littlegroup_q,wrtout,xmpi_barrier

SOURCE

4015 subroutine dvdb_test_ftinterp(db_path, ngqpt, comm)
4016 
4017 
4018 !This section has been created automatically by the script Abilint (TD).
4019 !Do not modify the following lines by hand.
4020 #undef ABI_FUNC
4021 #define ABI_FUNC 'dvdb_test_ftinterp'
4022 !End of the abilint section
4023 
4024  implicit none
4025 
4026 !Arguments ------------------------------------
4027  character(len=*),intent(in) :: db_path
4028  integer,intent(in) :: comm
4029  integer,intent(in) :: ngqpt(3)
4030 
4031 !Local variables-------------------------------
4032 !scalars
4033  integer :: nfft,iq,cplex,mu,ispden
4034  type(dvdb_t) :: db
4035 !arrays
4036  integer :: ngfft(18)
4037  real(dp),allocatable :: file_v1r(:,:,:,:),intp_v1r(:,:,:,:),tmp_v1r(:,:,:,:)
4038 
4039 ! *************************************************************************
4040 
4041  call dvdb_init(db, db_path, comm)
4042  db%debug = .True.
4043  db%symv1 = .True.
4044  db%symv1 = .False.
4045  call dvdb_print(db)
4046 
4047  ! Define FFT mesh for real space representation.
4048  call ngfft_seq(ngfft, db%ngfft3_v1(:,1))
4049  nfft = product(ngfft(1:3))
4050  call dvdb_open_read(db, ngfft, xmpi_comm_self)
4051 
4052  ABI_MALLOC(intp_v1r, (2,nfft,db%nspden,db%natom3))
4053  ABI_MALLOC(file_v1r, (2,nfft,db%nspden,db%natom3))
4054 
4055  call dvdb_ftinterp_setup(db,ngqpt,1,[zero,zero,zero],nfft,ngfft,comm)
4056 
4057  do iq=1,db%nqpt
4058    ! Read data from file
4059    call dvdb_readsym_allv1(db, dvdb_findq(db, db%qpts(:,iq)), cplex, nfft, ngfft, tmp_v1r, comm)
4060    if (cplex == 1) then
4061      file_v1r(1,:,:,:) = tmp_v1r(1,:,:,:)
4062      file_v1r(2,:,:,:) = zero
4063    else
4064      file_v1r = tmp_v1r
4065    end if
4066    ABI_FREE(tmp_v1r)
4067 
4068    ! Interpolate data on the same q-point
4069    call dvdb_ftinterp_qpt(db, db%qpts(:,iq), nfft, ngfft, intp_v1r, comm)
4070 
4071    write(std_out,"(a)")sjoin("=== For q-point:", ktoa(db%qpts(:,iq)), "===")
4072    do mu=1,db%natom3
4073      do ispden=1,db%nspden
4074        write(std_out,"(2(a,i0))")"  iatom3: ", mu, ", ispden: ", ispden
4075        call vdiff_print(vdiff_eval(2,nfft,file_v1r(:,:,ispden,mu),intp_v1r(:,:,ispden,mu),db%cryst%ucvol))
4076 
4077        !do ifft=1,nfft
4078        !  write(std_out,*)file_v1r(1,ifft,ispden,mu),intp_v1r(1,ifft,ispden,mu),&
4079        !  file_v1r(2,ifft,ispden,mu),intp_v1r(2,ifft,ispden,mu)
4080        !end do
4081        write(std_out,*)" "
4082      end do
4083    end do
4084    write(std_out,*)" "
4085  end do ! iq
4086 
4087  ABI_FREE(intp_v1r)
4088  ABI_FREE(file_v1r)
4089 
4090  call dvdb_free(db)
4091 
4092 end subroutine dvdb_test_ftinterp

m_dvdb/dvdb_test_v1complete [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  dvdb_test_v1complete

FUNCTION

  Debugging tool used to test the symmetrization of the DFPT potentials.

INPUTS

  db_path=Filename of the DVDB file.
  dump_path=File used to dump potentials (empty string to disable output)
  comm=MPI communicator.

OUTPUT

  Only writing.

PARENTS

      mrgdv

CHILDREN

      cwtime,dvdb_get_v1scf_qpt,dvdb_get_v1scf_rpt,dvdb_open_read
      hdr_fort_read,hdr_fort_write,hdr_free,irreducible_set_pert
      kpts_ibz_from_kptrlatt,littlegroup_q,wrtout,xmpi_barrier

SOURCE

3840 subroutine dvdb_test_v1complete(db_path, dump_path, comm)
3841 
3842 
3843 !This section has been created automatically by the script Abilint (TD).
3844 !Do not modify the following lines by hand.
3845 #undef ABI_FUNC
3846 #define ABI_FUNC 'dvdb_test_v1complete'
3847 !End of the abilint section
3848 
3849  implicit none
3850 
3851 !Arguments ------------------------------------
3852  character(len=*),intent(in) :: db_path,dump_path
3853  integer,intent(in) :: comm
3854 
3855 !Local variables-------------------------------
3856 !scalars
3857  integer,parameter :: master=0
3858  integer :: iqpt,pcase,idir,ipert,cplex,nfft,ispden,timerev_q,ifft,unt,my_rank
3859  character(len=500) :: msg
3860  type(crystal_t),pointer :: cryst
3861  type(dvdb_t),target :: db
3862 !arrays
3863  integer :: ngfft(18),rfdir(3)
3864  integer,allocatable :: pflag(:,:)
3865  real(dp) :: qpt(3)
3866  integer,allocatable :: pertsy(:,:),rfpert(:),symq(:,:,:)
3867  real(dp),allocatable :: file_v1scf(:,:,:,:),symm_v1scf(:,:,:,:)
3868 
3869 ! *************************************************************************
3870 
3871  my_rank = xmpi_comm_rank(comm)
3872 
3873  call dvdb_init(db, db_path, comm)
3874  db%debug = .True.
3875  db%symv1 = .False. !; db%symv1 = .True.
3876  call dvdb_print(db)
3877  call dvdb_list_perts(db, [-1,-1,-1])
3878 
3879  call ngfft_seq(ngfft, db%ngfft3_v1(:,1))
3880  nfft = product(ngfft(1:3))
3881  call dvdb_open_read(db, ngfft, comm)
3882 
3883  cryst => db%cryst
3884  call ngfft_seq(ngfft, db%ngfft3_v1(:,1))
3885  nfft = product(ngfft(:3))
3886 
3887  ABI_MALLOC(pflag, (3, db%natom))
3888 
3889  ! Initialize the list of perturbations rfpert and rdfir
3890  ! WARNING: Only phonon perturbations are considered for the time being.
3891  ABI_MALLOC(rfpert,(db%mpert))
3892  rfpert = 0; rfpert(1:cryst%natom) = 1; rfdir = 1
3893 
3894  ABI_MALLOC(symq, (4,2,cryst%nsym))
3895  ABI_MALLOC(pertsy, (3,db%mpert))
3896 
3897  unt = -1
3898  if (len_trim(dump_path) /= 0 .and. my_rank == master) then
3899    if (open_file(dump_path, msg, newunit=unt, action="write", status="unknown", form="formatted") /= 0) then
3900      MSG_ERROR(msg)
3901    end if
3902    write(std_out,"(a)")sjoin("Will write potentials to:", dump_path)
3903  end if
3904 
3905  do iqpt=1,db%nqpt
3906    qpt = db%qpts(:,iqpt)
3907 
3908    ! Examine the symmetries of the q wavevector
3909    call littlegroup_q(cryst%nsym,qpt,symq,cryst%symrec,cryst%symafm,timerev_q,prtvol=db%prtvol)
3910 
3911    ! Determine the symmetrical perturbations. Meaning of pertsy:
3912    !    0 for non-target perturbations
3913    !    1 for basis perturbations
3914    !   -1 for perturbations that can be found from basis perturbations
3915    call irreducible_set_pert(cryst%indsym,db%mpert,cryst%natom,cryst%nsym,&
3916      pertsy,rfdir,rfpert,symq,cryst%symrec,cryst%symrel)
3917 
3918    ! Read all potentials (here I assume that all perturbations are available)
3919    call dvdb_readsym_allv1(db, iqpt, cplex, nfft, ngfft, file_v1scf, db%comm)
3920 
3921    ! Copy basis perturbations in symm_v1scf and set pflag
3922    ABI_MALLOC(symm_v1scf, (cplex,nfft,db%nspden,db%natom3))
3923    symm_v1scf = huge(one); pflag = 0
3924    do pcase=1,3*db%cryst%natom
3925      idir = mod(pcase-1, 3) + 1; ipert = (pcase - idir) / 3 + 1
3926      if (pertsy(idir, ipert) == 1) then
3927        symm_v1scf(:,:,:,pcase) = file_v1scf(:,:,:,pcase)
3928        pflag(idir,ipert) = 1
3929      end if
3930    end do
3931 
3932    ! Complete potentials
3933    call v1phq_complete(cryst,qpt,ngfft,cplex,nfft,db%nspden,db%nsppol,db%mpi_enreg,db%symv1,pflag,symm_v1scf)
3934 
3935    ! Compare values.
3936    do pcase=1,3*cryst%natom
3937      idir = mod(pcase-1, 3) + 1; ipert = (pcase - idir) / 3 + 1
3938      if (pflag(idir,ipert) == 1) cycle
3939 
3940      write(std_out,"(2(a,i0),2a)")"For idir: ",idir, ", ipert: ", ipert, ", qpt: ",trim(ktoa(qpt))
3941      do ispden=1,db%nspden
3942        !write(std_out,"(a,es10.3)")"  max(abs(f1-f2))", &
3943        ! maxval(abs(file_v1scf(:,:,ispden,pcase) - symm_v1scf(:,:,ispden,pcase)))
3944        call vdiff_print(vdiff_eval(cplex,nfft,file_v1scf(:,:,ispden,pcase),symm_v1scf(:,:,ispden,pcase),cryst%ucvol))
3945 
3946        ! Debug: Write potentials to file.
3947        if (unt /= -1) then
3948          write(unt,*)"# q-point:", trim(ktoa(qpt))
3949          write(unt,*)"# idir: ",idir,", ipert: ",ipert,", ispden:", ispden
3950          write(unt,*)"# file_v1scf, symmetrized_v1scf, diff"
3951          if (cplex == 1) then
3952            do ifft=1,nfft
3953              write(unt,"(3(es12.4,2x))") &
3954                file_v1scf(1,ifft,ispden,pcase), symm_v1scf(1,ifft,ispden,pcase), &
3955                file_v1scf(1,ifft,ispden,pcase) - symm_v1scf(1,ifft,ispden,pcase)
3956            end do
3957          else
3958            do ifft=1,nfft
3959              write(unt, "(6(es12.4,2x))") &
3960                file_v1scf(1,ifft,ispden,pcase), symm_v1scf(1,ifft,ispden,pcase),   &
3961                file_v1scf(1,ifft,ispden,pcase) - symm_v1scf(1,ifft,ispden,pcase), &
3962                file_v1scf(2,ifft,ispden,pcase), symm_v1scf(2,ifft,ispden,pcase),   &
3963                file_v1scf(2,ifft,ispden,pcase) - symm_v1scf(2,ifft,ispden,pcase)
3964            end do
3965          end if
3966        end if
3967 
3968      end do
3969      write(std_out,*)""
3970    end do
3971 
3972    ABI_FREE(symm_v1scf)
3973    ABI_FREE(file_v1scf)
3974  end do
3975 
3976  ABI_FREE(pflag)
3977  ABI_FREE(rfpert)
3978  ABI_FREE(symq)
3979  ABI_FREE(pertsy)
3980 
3981  call dvdb_free(db)
3982 
3983  if (unt /= -1) close(unt)
3984 
3985 end subroutine dvdb_test_v1complete

m_dvdb/dvdb_test_v1rsym [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  dvdb_test_v1rsym

FUNCTION

  Debugging tool used to check whether the DFPT potentials in real space fulfill
  the correct symmetries on the real space FFT mesh.

INPUTS

  db_path=Filename
  comm=MPI communicator.

OUTPUT

  Only writing.

PARENTS

      mrgdv

CHILDREN

      cwtime,dvdb_get_v1scf_qpt,dvdb_get_v1scf_rpt,dvdb_open_read
      hdr_fort_read,hdr_fort_write,hdr_free,irreducible_set_pert
      kpts_ibz_from_kptrlatt,littlegroup_q,wrtout,xmpi_barrier

SOURCE

3692 subroutine dvdb_test_v1rsym(db_path, comm)
3693 
3694 
3695 !This section has been created automatically by the script Abilint (TD).
3696 !Do not modify the following lines by hand.
3697 #undef ABI_FUNC
3698 #define ABI_FUNC 'dvdb_test_v1rsym'
3699 !End of the abilint section
3700 
3701  implicit none
3702 
3703 !Arguments ------------------------------------
3704  character(len=*),intent(in) :: db_path
3705  integer,intent(in) :: comm
3706 
3707 !Local variables-------------------------------
3708 !scalars
3709  integer,parameter :: rfmeth2=2,syuse0=0
3710  integer :: iqpt,idir,ipert,nsym1,cplex,v1pos
3711  integer :: isym,nfft,ifft,ifft_rot,ispden
3712  real(dp) :: max_err,re,im,vre,vim !,pre,pim
3713  character(len=500) :: msg
3714  logical :: isok
3715  type(dvdb_t),target :: db
3716  type(crystal_t),pointer :: cryst
3717 !arrays
3718  integer :: ngfft(18)
3719  integer,allocatable :: symafm1(:),symrel1(:,:,:),irottb(:,:)
3720  real(dp) :: qpt(3)
3721  real(dp),allocatable :: tnons1(:,:),v1scf(:,:)
3722 
3723 ! *************************************************************************
3724 
3725  call dvdb_init(db, db_path, comm)
3726  db%debug = .True.
3727  !db%symv1 = .True.
3728  db%symv1 = .False.
3729  call dvdb_print(db)
3730 
3731  call ngfft_seq(ngfft, db%ngfft3_v1(:,1))
3732  nfft = product(ngfft(1:3))
3733  call dvdb_open_read(db, ngfft, comm)
3734 
3735  ABI_CHECK(db%nspinor==1, "nspinor == 2 not coded")
3736 
3737  cryst => db%cryst
3738  ABI_MALLOC(symafm1, (cryst%nsym))
3739  ABI_MALLOC(symrel1, (3,3,cryst%nsym))
3740  ABI_MALLOC(tnons1, (3,cryst%nsym))
3741 
3742  do iqpt=1,db%nqpt
3743    qpt = db%qpts(:, iqpt)
3744 
3745    do ipert=1,db%natom
3746      do idir=1,3
3747        v1pos = db%pos_dpq(idir, ipert, iqpt); if (v1pos == 0) cycle
3748 
3749        ! Determines the set of symmetries that leaves the perturbation invariant.
3750        call littlegroup_pert(cryst%gprimd,idir,cryst%indsym,dev_null,ipert,cryst%natom,cryst%nsym,nsym1,rfmeth2,&
3751          cryst%symafm,symafm1,db%symq_table(:,:,:,iqpt),cryst%symrec,cryst%symrel,symrel1,syuse0,cryst%tnons,&
3752          tnons1,unit=dev_null)
3753 
3754        ! TODO: (3) or (18)
3755        cplex = db%cplex_v1(v1pos)
3756        ngfft(1:3) = db%ngfft3_v1(:, v1pos)
3757        nfft = product(ngfft(:3))
3758        ABI_MALLOC(v1scf, (cplex*nfft, db%nspden))
3759 
3760        if (dvdb_read_onev1(db, idir, ipert, iqpt, cplex, nfft, ngfft, v1scf, msg) /= 0) then
3761          MSG_ERROR(msg)
3762        end if
3763 
3764        ABI_MALLOC(irottb, (nfft,nsym1))
3765        call rotate_fft_mesh(nsym1,symrel1,tnons1,ngfft,irottb,isok)
3766        if (.not. isok) then
3767          MSG_WARNING("Real space FFT mesh is not compatible with symmetries!")
3768        end if
3769 
3770        max_err = zero
3771        do isym=1,nsym1
3772          do ispden=1,db%nspden
3773            do ifft=1,nfft
3774              ifft_rot = irottb(ifft, isym)
3775              !pre =  cos(two_pi * dot_product(qpt, tnons1(:,isym)))
3776              !pim = -sin(two_pi * dot_product(qpt, tnons1(:,isym)))
3777              if (cplex == 2) then
3778                re = v1scf(2*ifft_rot-1, ispden)
3779                im = v1scf(2*ifft_rot  , ispden)
3780                !vre = re * pre - im * pim
3781                !vim = re * pim + im * pre
3782                vre = re; vim = im
3783 
3784                re = v1scf(2*ifft-1, ispden) - vre
3785                im = v1scf(2*ifft  , ispden) - vim
3786              else
3787                re = v1scf(ifft, ispden) - v1scf(ifft_rot, ispden)
3788                im = zero
3789              end if
3790              !if (sqrt(re**2 + im**2) > tol6) write(std_out,*)"ifft,isym,err: ",ifft,isym,sqrt(re**2 + im**2)
3791              max_err = max(max_err, sqrt(re**2 + im**2))
3792            end do
3793          end do
3794        end do
3795        write(std_out,"(3(a,i0),a,es16.8)")"For iqpt= ",iqpt,", idir= ",idir,", ipert= ",ipert,", max_err= ",max_err
3796 
3797        ABI_FREE(irottb)
3798        ABI_FREE(v1scf)
3799      end do
3800    end do
3801 
3802  end do ! iqpt
3803 
3804  ABI_FREE(symafm1)
3805  ABI_FREE(symrel1)
3806  ABI_FREE(tnons1)
3807 
3808  call dvdb_free(db)
3809 
3810 end subroutine dvdb_test_v1rsym

m_dvdb/dvdb_v1r_long_range [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  dvdb_v1r_long_range

FUNCTION

  Compute the long-range part of the phonon potential
  due to the Born effective charges, PRL 115, 176401 (2015) [[cite:Verdi2015]].

    V^L_{iatom,idir}(r) = i (4pi/vol) sum_G (q+G) . Zeff_{iatom,idir}
                           e^{i (q + G) . (r - tau_{iatom})} / ((q + G) . dielt . (q + G))

  where Zeff and dielt are the Born effective charge tensor and the dielectric tensor,
  tau is the atom position, and vol is the volume of the unit cell.

INPUTS

  db = the DVDB object.
  qpt = the q-point in reduced coordinates.
  iatom = atom index.
  idir = direction index.
  nfft = number of fft points.
  ngfft(18) = FFT mesh.

OUTPUT

  v1r_lr = dipole potential

PARENTS

      m_dvdb

CHILDREN

      cwtime,dvdb_get_v1scf_qpt,dvdb_get_v1scf_rpt,dvdb_open_read
      hdr_fort_read,hdr_fort_write,hdr_free,irreducible_set_pert
      kpts_ibz_from_kptrlatt,littlegroup_q,wrtout,xmpi_barrier

SOURCE

4133 subroutine dvdb_v1r_long_range(db,qpt,iatom,idir,nfft,ngfft,v1r_lr)
4134 
4135 
4136 !This section has been created automatically by the script Abilint (TD).
4137 !Do not modify the following lines by hand.
4138 #undef ABI_FUNC
4139 #define ABI_FUNC 'dvdb_v1r_long_range'
4140 !End of the abilint section
4141 
4142  implicit none
4143 
4144 !Arguments ------------------------------------
4145 !scalars
4146  type(dvdb_t),intent(in) :: db
4147  integer,intent(in) :: iatom, idir
4148  integer,intent(in) :: nfft
4149 !arrays
4150  integer,intent(in) :: ngfft(18)
4151  real(dp),intent(in) :: qpt(3)
4152  real(dp),intent(out) :: v1r_lr(2,nfft)
4153 
4154 !Local variables-------------------------------
4155 !scalars
4156  integer :: n1, n2, n3, nfftot, ig
4157  real(dp) :: fac, qGZ, denom, denom_inv, qtau
4158  real(dp) :: re, im, phre, phim
4159  real(dp) :: dummy
4160  type(MPI_type) :: MPI_enreg_seq
4161 !arrays
4162  integer, allocatable :: gfft(:,:)
4163  real(dp) :: gprimd(3,3), rprimd(3,3)
4164  real(dp) :: dielt(3,3), zeff(3,3), dielt_red(3,3)
4165  real(dp) :: qG(3), Zstar(3), tau(3)
4166  real(dp), allocatable :: v1G_lr(:,:)
4167 
4168 
4169 ! *************************************************************************
4170 
4171  ! Make sure FFT parallelism is not used
4172  n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3); nfftot = product(ngfft(1:3))
4173  ABI_CHECK(nfftot == nfft, "FFT parallelism not supported")
4174 
4175  ! Fake MPI_type for sequential fft execution
4176  call initmpi_seq(MPI_enreg_seq)
4177  call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfft(2),ngfft(3),'all')
4178 
4179  ! Allocate memory
4180  ABI_MALLOC(gfft, (3, nfft))
4181  ABI_MALLOC(v1G_lr, (2,nfft))
4182 
4183  ! Reciprocal and real space primitive vectors
4184  gprimd = db%cryst%gprimd
4185  rprimd = db%cryst%rprimd
4186 
4187  !Prefactor
4188  fac = four_pi / db%cryst%ucvol
4189 
4190  ! Transform the Born effective charge tensor from Cartesian to reduced coordinates
4191  ! and select the relevant direction.
4192  zeff = db%zeff(:,:,iatom)
4193  Zstar = matmul(transpose(gprimd), matmul(zeff, rprimd(:,idir))) * two_pi
4194 
4195  ! Transform the dielectric tensor from Cartesian to reduced coordinates.
4196  dielt = db%dielt
4197  dielt_red = matmul(transpose(gprimd), matmul(dielt, gprimd)) * two_pi ** 2
4198 
4199  ! Atom position
4200  tau = db%cryst%xred(:,iatom)
4201 
4202  ! Get the set of G vectors
4203  call get_gftt(ngfft,qpt,db%cryst%gmet,dummy,gfft)
4204 
4205  ! Compute the long-range potential in G-space
4206  v1G_lr = zero
4207  do ig=1,nfft
4208 
4209    ! (q + G)
4210    qG(:) = qpt(:) + gfft(:,ig)
4211 
4212    ! (q + G) . Zeff(:,idir,iatom)
4213    qGZ = dot_product(qG, Zstar)
4214 
4215    ! (q + G) . dielt . (q + G)
4216    denom = dot_product(qG, matmul(dielt_red, qG))
4217 
4218    ! Avoid (q+G) = 0
4219    denom_inv = denom / (denom ** 2 + tol10 ** 2)
4220 
4221    ! Phase factor exp(-i (q+G) . tau)
4222    qtau = - two_pi * dot_product(qG, tau)
4223    phre = cos(qtau); phim = sin(qtau)
4224 
4225    re = zero
4226    im = fac * qGZ * denom_inv
4227 
4228    v1G_lr(1,ig) = phre * re - phim * im
4229    v1G_lr(2,ig) = phim * re + phre * im
4230 
4231  end do
4232 
4233  ! Free memory
4234  ABI_FREE(gfft)
4235 
4236  ! FFT to get the long-range potential in r space
4237  v1r_lr = zero
4238  call fourdp(2,v1G_lr,v1r_lr,1,MPI_enreg_seq,nfft,ngfft,1,0)
4239 
4240  ! Multiply by  exp(i q . r)
4241  call times_eikr(qpt,ngfft,nfft,1,v1r_lr)
4242 
4243  ! Free memory
4244  ABI_FREE(v1G_lr)
4245 
4246  call destroy_mpi_enreg(MPI_enreg_seq)
4247 
4248 end subroutine dvdb_v1r_long_range

m_dvdb/find_symeq [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

 find_symeq

FUNCTION

  Find symmetry which links to the perturbation specified by (idir, ipert)

INPUTS

  cryst<crystal_t>=crystal structure parameters
  idir=Direction of the perturbation
  ipert=Perturbation type.
  symq(4,2,nsym)=Table produced by littlegroup_q
  pflag(3,natom)= For each atomic perturbation:
     0 if pert is not available. 1 if pert is available. 2 if pert can been reconstructed by symmetry.

OUTPUT

  ipert_eq
  isym_eq
  itirev_eq
  g0_qpt(3)

PARENTS

      m_dvdb

CHILDREN

      cwtime,dvdb_get_v1scf_qpt,dvdb_get_v1scf_rpt,dvdb_open_read
      hdr_fort_read,hdr_fort_write,hdr_free,irreducible_set_pert
      kpts_ibz_from_kptrlatt,littlegroup_q,wrtout,xmpi_barrier

SOURCE

1404 subroutine find_symeq(cryst, idir, ipert, symq, pflag, ipert_eq, isym_eq, itirev_eq, g0_qpt)
1405 
1406 
1407 !This section has been created automatically by the script Abilint (TD).
1408 !Do not modify the following lines by hand.
1409 #undef ABI_FUNC
1410 #define ABI_FUNC 'find_symeq'
1411 !End of the abilint section
1412 
1413  implicit none
1414 
1415 !Arguments ------------------------------------
1416 !scalars
1417  integer,intent(in) :: idir,ipert
1418  integer,intent(out) :: ipert_eq,isym_eq,itirev_eq
1419  type(crystal_t),intent(in) :: cryst
1420 !arrays
1421  integer,intent(in) :: symq(4,2,cryst%nsym),pflag(3,cryst%natom)
1422  integer,intent(out) :: g0_qpt(3)
1423 
1424 !Local variables-------------------------------
1425 !scalars
1426  integer :: isym,idir_eq,ip,itirev
1427 
1428 ! *************************************************************************
1429 
1430  isym_eq = -1; ipert_eq = -1
1431 
1432 symloop: &
1433  do itirev=1,2
1434    itirev_eq = itirev
1435    do isym=1,cryst%nsym
1436 
1437      ! Check that isym preserves the q-point
1438      !if (symq(4,itirev,isym) /= 1 .or. any(symq(1:3,itirev,isym) /= 0)) cycle
1439      if (symq(4,itirev,isym) /= 1) cycle ! .or. any(symq(1:3,itirev,isym) /= 0)) cycle
1440      g0_qpt = symq(1:3,itirev,isym)
1441 
1442      do ip=1,cryst%natom
1443        !if (.not. cryst%indsym(4,isym,ip) == ipert) cycle
1444        if (.not. cryst%indsym(4,isym,ipert) == ip) cycle
1445        isym_eq = isym; ipert_eq = ip
1446        do idir_eq=1,3
1447          if (idir_eq == idir .and. ip == ipert .and. cryst%symrec(idir,idir_eq,isym) /= 0) isym_eq = -1
1448          if (cryst%symrec(idir,idir_eq,isym) /= 0 .and. pflag(idir_eq, ip) == 0) then
1449          !if (cryst%symrel(idir,idir_eq,isym) /= 0 .and. pflag(idir_eq, ip) == 0) then
1450            !if (idir_eq == idir .and. ip == ipert) cycle
1451            isym_eq = -1
1452          end if
1453        end do
1454        if (isym_eq /= -1) exit symloop
1455      end do
1456    end do
1457  end do symloop
1458 
1459  if (isym_eq == -1) then
1460    ipert_eq = -1; itirev_eq = -1
1461  end if
1462 
1463 end subroutine find_symeq

m_dvdb/my_hdr_skip [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  my_hdr_skip

FUNCTION

  Skip the header without rewinding the file. Return exit code.

NOTES

  Because hdr_skip rewinds the file and I'm not gonna change that ugly code.

PARENTS

CHILDREN

SOURCE

3122 integer function my_hdr_skip(unit, idir, ipert, qpt, msg) result(ierr)
3123 
3124 
3125 !This section has been created automatically by the script Abilint (TD).
3126 !Do not modify the following lines by hand.
3127 #undef ABI_FUNC
3128 #define ABI_FUNC 'my_hdr_skip'
3129 !End of the abilint section
3130 
3131  implicit none
3132 
3133 !Arguments ------------------------------------
3134 !scalars
3135  integer,intent(in) :: unit,idir,ipert
3136  real(dp),intent(in) :: qpt(3)
3137  character(len=500),intent(out) :: msg
3138 
3139 !Local variables-------------------------------
3140  integer :: fform
3141  type(hdr_type) :: tmp_hdr
3142 !************************************************************************
3143 
3144  ierr = 0; msg = ""
3145  call hdr_fort_read(tmp_hdr, unit, fform)
3146  ierr = dvdb_check_fform(fform, "read_dvdb", msg)
3147  if (ierr /= 0) return
3148 
3149  if (idir /= -1 .and. ipert /= -1) then
3150    if (idir /= mod(tmp_hdr%pertcase-1, 3) + 1 .or. &
3151        ipert /= (tmp_hdr%pertcase - idir) / 3 + 1 .or. &
3152        any(abs(qpt - tmp_hdr%qptn) > tol14)) then
3153          msg = "Perturbation index on file does not match the one expected by caller"
3154          ierr = -1
3155    end if
3156  end if
3157 
3158  call hdr_free(tmp_hdr)
3159 
3160 end function my_hdr_skip

m_dvdb/rotate_fqg [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

  rotate_fqg

FUNCTION

INPUTS

PARENTS

      m_dvdb

CHILDREN

      cwtime,dvdb_get_v1scf_qpt,dvdb_get_v1scf_rpt,dvdb_open_read
      hdr_fort_read,hdr_fort_write,hdr_free,irreducible_set_pert
      kpts_ibz_from_kptrlatt,littlegroup_q,wrtout,xmpi_barrier

SOURCE

1726 subroutine rotate_fqg(itirev, symm, qpt, tnon, ngfft, nfft, nspden, infg, outfg)
1727 
1728 
1729 !This section has been created automatically by the script Abilint (TD).
1730 !Do not modify the following lines by hand.
1731 #undef ABI_FUNC
1732 #define ABI_FUNC 'rotate_fqg'
1733 !End of the abilint section
1734 
1735  implicit none
1736 
1737 !Arguments ------------------------------------
1738 !scalars
1739  integer,intent(in) :: itirev,nfft,nspden
1740 !arrays
1741  integer,intent(in) :: symm(3,3),ngfft(18)
1742  real(dp),intent(in) :: qpt(3),tnon(3)
1743  real(dp),intent(in) :: infg(2,nfft,nspden)
1744  real(dp),intent(out) :: outfg(2,nfft,nspden)
1745 
1746 !Local variables-------------------------------
1747 !scalars
1748  integer :: i1,i2,i3,id1,id2,id3,n1,n2,n3,ind1,ind2,j1,j2,j3,l1,l2,l3,k1,k2,k3,nfftot,isp,tsign
1749  real(dp) :: arg
1750  logical :: has_phase
1751 !arrays
1752  integer :: tsg(3)
1753  real(dp) :: phnon1(2)
1754 
1755 ! *************************************************************************
1756 
1757  n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3); nfftot = product(ngfft(1:3))
1758  ABI_CHECK(nfftot == nfft, "FFT parallelism not supported")
1759  id1 = n1/2+2; id2 = n2/2+2; id3=n3/2+2
1760 
1761  ABI_CHECK(any(itirev == [1,2]), "Wrong itirev")
1762  tsign = 3-2*itirev; has_phase = any(abs(tnon) > tol12)
1763 
1764  do isp=1,nspden
1765    ind1 = 0
1766    do i3=1,n3
1767      do i2=1,n2
1768        do i1=1,n1
1769          ind1 = ind1 + 1
1770 
1771          ! Get location of G vector (grid point) centered at 0 0 0
1772          l1 = i1-(i1/id1)*n1-1
1773          l2 = i2-(i2/id2)*n2-1
1774          l3 = i3-(i3/id3)*n3-1
1775 
1776          ! Get rotated G vector. IS(G)
1777          j1 = tsign * (symm(1,1)*l1+symm(1,2)*l2+symm(1,3)*l3)
1778          j2 = tsign * (symm(2,1)*l1+symm(2,2)*l2+symm(2,3)*l3)
1779          j3 = tsign * (symm(3,1)*l1+symm(3,2)*l2+symm(3,3)*l3)
1780 
1781          ! FIXME :TO BE CLARIFIED:
1782          ! We are not working on the G-sphere thus SG may be outside
1783          ! of the box. This check is not done in irrzg!!!
1784          if ( (j1 > n1/2 .or. j1 < -(n1-1)/2) .or. &
1785               (j2 > n2/2 .or. j1 < -(n2-1)/2) .or. &
1786               (j3 > n3/2 .or. j3 < -(n3-1)/2) ) then
1787            !write(std_out,*)"got it"
1788            outfg(:,ind1,isp) = zero
1789            cycle
1790          end if
1791 
1792          tsg = [j1,j2,j3] ! +- S^{-1} G
1793 
1794          ! Map into [0,n-1] and then add 1 for array index in [1,n]
1795          k1=1+mod(n1+mod(j1,n1),n1)
1796          k2=1+mod(n2+mod(j2,n2),n2)
1797          k3=1+mod(n3+mod(j3,n3),n3)
1798 
1799          ! Get linear index of rotated point Gj
1800          ind2 = k1+n1*((k2-1)+n2*(k3-1))
1801 
1802          if (has_phase) then
1803            ! compute exp(-2*Pi*I*G dot tau) using original G
1804            arg = two_pi * dot_product(qpt + tsg, tnon)
1805            phnon1(1) = cos(arg); phnon1(2) =-sin(arg)
1806 
1807            ! rho(Strans*G)=exp(2*Pi*I*(G) dot tau_S) rho(G)
1808            outfg(1, ind1, isp) = phnon1(1) * infg(1, ind2, isp) - phnon1(2) * infg(2, ind2, isp)
1809            outfg(2, ind1, isp) = phnon1(1) * infg(2, ind2, isp) + phnon1(2) * infg(1, ind2, isp)
1810          else
1811            outfg(1, ind1, isp) = infg(1, ind2, isp)
1812            outfg(2, ind1, isp) = infg(2, ind2, isp)
1813          end if
1814 
1815          ! Take complex conjugate if time-reversal is used.
1816          if (tsign == -1) outfg(2, ind1, isp) = -outfg(2, ind1, isp)
1817        end do
1818      end do
1819    end do
1820  end do ! isp
1821 
1822 end subroutine rotate_fqg

m_dvdb/v1phq_complete [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

 v1phq_complete

FUNCTION

  Use the symmetries of the little group of the q-point to reconstruct
  the first order potentials starting from an initial irreducible set.

INPUTS

  cryst<crystal_t>=crystal structure parameters
  qpt(3)=q-point in reduced coordinates.
  ngfft(18)=Info of FFT grid.
  cplex=1 if real potentials (qpt==gamma), 2 if complex
  nfft=(effective) number of FFT grid points (for this proc).
  nspden=number of spin-density components
  nsppol=Number of independent spin polarizations
  mpi_enreg=information about MPI parallelization
  symv1=If True, the new potentials are symmetrized using the set of symmetries that leaves the
    perturbation invariant.

SIDE EFFECTS

  pflag(3,natom)= For each atomic perturbation:
     0 if pert is not available. 1 if pert is available. 2 if pert has been reconstructed by symmetry.
     Initialized by the caller. Changed in output.
  v1scf(cplex*nfft,nspden,3*cryst%natom)=Array with first order potentials.
    in input: filled with the irreducible potentials (corresponding pflag set to 1)
    output: Contains full set of perturbations.

PARENTS

      m_dvdb

CHILDREN

      cwtime,dvdb_get_v1scf_qpt,dvdb_get_v1scf_rpt,dvdb_open_read
      hdr_fort_read,hdr_fort_write,hdr_free,irreducible_set_pert
      kpts_ibz_from_kptrlatt,littlegroup_q,wrtout,xmpi_barrier

SOURCE

1162 subroutine v1phq_complete(cryst,qpt,ngfft,cplex,nfft,nspden,nsppol,mpi_enreg,symv1,pflag,v1scf)
1163 
1164 
1165 !This section has been created automatically by the script Abilint (TD).
1166 !Do not modify the following lines by hand.
1167 #undef ABI_FUNC
1168 #define ABI_FUNC 'v1phq_complete'
1169 !End of the abilint section
1170 
1171  implicit none
1172 
1173 !Arguments ------------------------------------
1174 !scalars
1175  integer,intent(in) :: cplex,nfft,nspden,nsppol
1176  logical,intent(in) :: symv1
1177  type(crystal_t),intent(in) :: cryst
1178  type(MPI_type),intent(in) :: mpi_enreg
1179 !arrays
1180  integer,intent(in) :: ngfft(18)
1181  integer,intent(inout) :: pflag(3, cryst%natom)
1182  real(dp),intent(in) :: qpt(3)
1183  real(dp),intent(inout) :: v1scf(cplex*nfft,nspden,3*cryst%natom)
1184 
1185 !Local variables-------------------------------
1186 !scalars
1187  integer,parameter :: syuse0=0,rfmeth2=2,tim_fourdp0=0
1188  integer :: idir,ipert,tsign,isym_eq,itirev_eq,ipert_eq !,itirev
1189  integer :: pcase,trev_q,idir_eq,pcase_eq,ispden,cnt
1190  integer :: i1,i2,i3,id1,id2,id3,n1,n2,n3,ind1,ind2,j1,j2,j3,l1,l2,l3,k1,k2,k3,nfftot
1191  real(dp) :: arg
1192  logical :: has_phase
1193  logical,parameter :: debug=.False.
1194  character(len=500) :: msg
1195  integer,save :: enough=0
1196 !arrays
1197  integer :: symrel_eq(3,3),symrec_eq(3,3),symm(3,3),g0_qpt(3),l0(3),tsm1g(3)
1198  integer :: symq(4,2,cryst%nsym)
1199  real(dp) :: phnon1(2),tnon(3)
1200  real(dp),allocatable :: workg(:,:), workg_eq(:,:),v1g(:,:,:)
1201 
1202 ! *************************************************************************
1203 
1204  n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3); nfftot = product(ngfft(1:3))
1205  ABI_CHECK(nfftot == nfft, "FFT parallelism not supported")
1206  id1 = n1/2+2; id2 = n2/2+2; id3 = n3/2+2
1207 
1208  ABI_MALLOC(v1g, (2,nfft,nspden))
1209  ABI_MALLOC(workg_eq, (2, nfft))
1210  ABI_MALLOC(workg, (2, nfft))
1211 
1212  ! Examine the symmetries of the q wavevector
1213  call littlegroup_q(cryst%nsym,qpt,symq,cryst%symrec,cryst%symafm,trev_q,prtvol=0)
1214 
1215 pcase_loop: &
1216  do pcase=1,3*cryst%natom
1217    idir = mod(pcase-1, 3) + 1; ipert = (pcase - idir) / 3 + 1
1218    if (pflag(idir, ipert) /= 0) cycle ! This pcase is available
1219 
1220    ! Find symmetry which links to the perturbation requested (pcase)
1221    call find_symeq(cryst, idir, ipert, symq, pflag, ipert_eq, isym_eq, itirev_eq, g0_qpt)
1222    if (isym_eq == -1) then
1223      if (debug) write(std_out,*)"Cannot find isym eq for idir, ipert:", idir,ipert
1224      cycle pcase_loop
1225    end if
1226 
1227    ! set flag since we will reconstruct pcase from isym_eq.
1228    pflag(idir, ipert) = 2
1229 
1230    symrel_eq = cryst%symrel(:,:,isym_eq)
1231    symrec_eq = cryst%symrec(:,:,isym_eq)
1232    tsign = 3-2*itirev_eq
1233 
1234    ! Phase due to L0 + R^{-1}tau
1235    l0 = cryst%indsym(1:3,isym_eq,ipert)
1236    call mati3inv(symrel_eq, symm); symm = transpose(symm)
1237 
1238    tnon = l0 + matmul(transpose(symrec_eq), cryst%tnons(:,isym_eq))
1239    has_phase = any(abs(tnon) > tol12)
1240    ! FIXME
1241    !ABI_CHECK(.not. has_phase, "has phase must be tested")
1242    if (has_phase) then
1243      enough = enough + 1
1244      if (enough == 1) MSG_WARNING("has phase must be tested")
1245    end if
1246 
1247    workg = zero
1248 
1249    ! Reconstruct DFPT potential. Final results stored in v1g.
1250    !write(std_out,*)"Reconstructing idir, ipert",idir,ipert
1251    v1g = zero; cnt = 0
1252    do idir_eq=1,3
1253      if (symrec_eq(idir, idir_eq) == 0) cycle
1254      cnt = cnt + 1
1255      pcase_eq = idir_eq + (ipert_eq-1)*3
1256      if (debug) write(std_out,*)"idir_eq, ipert_eq, tsign",idir_eq, ipert_eq, tsign
1257 
1258      do ispden=1,nspden
1259        ! Get symmetric perturbation in G-space in workg_eq array.
1260        call fourdp(cplex,workg_eq,v1scf(:,ispden,pcase_eq),-1,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb,tim_fourdp0)
1261 
1262        !call rotate_fqg(itirev_eq,symrec_eq,qpt,tnon,ngfft,nfft,nspden,workg_eq,workg)
1263        ind1=0
1264        do i3=1,n3
1265          do i2=1,n2
1266            do i1=1,n1
1267              ind1=ind1+1
1268 
1269              ! Get location of G vector (grid point) centered at 0 0 0
1270              l1 = i1-(i1/id1)*n1-1
1271              l2 = i2-(i2/id2)*n2-1
1272              l3 = i3-(i3/id3)*n3-1
1273 
1274              ! Get rotated G vector Gj for each symmetry element
1275              ! -- here we use the TRANSPOSE of symrel_eq; assuming symrel_eq expresses
1276              ! the rotation in real space, the transpose is then appropriate
1277              ! for G space symmetrization (p. 1172d,e of notes, 2 June 1995).
1278              j1 = tsign * (symrel_eq(1,1)*l1+symrel_eq(2,1)*l2+symrel_eq(3,1)*l3)
1279              j2 = tsign * (symrel_eq(1,2)*l1+symrel_eq(2,2)*l2+symrel_eq(3,2)*l3)
1280              j3 = tsign * (symrel_eq(1,3)*l1+symrel_eq(2,3)*l2+symrel_eq(3,3)*l3)
1281 
1282              ! FIXME :TO BE CLARIFIED:
1283              ! We are not working on the G-sphere thus SG may be outside
1284              ! of the box. This check is not done in irrzg!!!
1285              if ( (j1 > n1/2 .or. j1 < -(n1-1)/2) .or. &
1286                   (j2 > n2/2 .or. j1 < -(n2-1)/2) .or. &
1287                   (j3 > n3/2 .or. j3 < -(n3-1)/2) ) then
1288                !write(std_out,*)"got it"
1289                workg(:, ind1) = zero; cycle
1290              end if
1291 
1292              tsm1g = [j1,j2,j3] ! +- S^{-1} G
1293 
1294              ! Map into [0,n-1] and then add 1 for array index in [1,n]
1295              k1=1+mod(n1+mod(j1,n1),n1)
1296              k2=1+mod(n2+mod(j2,n2),n2)
1297              k3=1+mod(n3+mod(j3,n3),n3)
1298 
1299              ! Get linear index of rotated point Gj
1300              ind2 = k1+n1*((k2-1)+n2*(k3-1))
1301 
1302              if (has_phase) then
1303                ! compute exp(-2*Pi*I*G dot tau) using original G
1304                ! NB: this phase is same as that in irrzg and phnons1, and corresponds
1305                ! to complex conjugate of phase from G to Gj;
1306                ! we use it immediately below, to go _to_ workg_eq(ind1)
1307                arg = two_pi * dot_product(qpt + tsm1g, tnon)
1308                phnon1(1) = cos(arg); phnon1(2) = -sin(arg)
1309 
1310                ! rho(Strans*G)=exp(2*Pi*I*(G) dot tau_S) rho(G)
1311                workg(1, ind1) = phnon1(1) * workg_eq(1, ind2) - phnon1(2) * workg_eq(2, ind2)
1312                workg(2, ind1) = phnon1(1) * workg_eq(2, ind2) + phnon1(2) * workg_eq(1, ind2)
1313              else
1314                workg(1, ind1) = workg_eq(1, ind2)
1315                workg(2, ind1) = workg_eq(2, ind2)
1316              end if
1317 
1318              ! Take complex conjugate if time-reversal is used.
1319              if (tsign == -1) workg(2, ind1) = -workg(2, ind1)
1320            end do
1321          end do
1322        end do
1323 
1324        v1g(:,:,ispden) = v1g(:,:,ispden) + workg * symrec_eq(idir, idir_eq)
1325      end do ! ispden
1326    end do ! idir_eq
1327    if (debug) write(std_out,*)"Used ",cnt," equivalent perturbations"
1328 
1329    ! Get potential in real space (results in v1scf)
1330    do ispden=1,nspden
1331      call fourdp(cplex,v1g(:,:,ispden),v1scf(:,ispden,pcase),+1,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb,tim_fourdp0)
1332 
1333      ! IS(q) = q + G0
1334      ! we want q so we have to multiply by exp(iG0r) in real space.
1335      if (any(g0_qpt /= 0)) then
1336        ABI_CHECK(cplex==2, "cplex == 1")
1337        if (debug) write(std_out,*)"Found not zero g0_qpt for idir, ipert:",idir,ipert
1338        call times_eigr(g0_qpt, ngfft, nfft, 1, v1scf(:,ispden,pcase))
1339      end if
1340    end do
1341 
1342    if (symv1) then
1343      if (debug) write(std_out,*)"Calling v1phq_symmetrize"
1344      call v1phq_symmetrize(cryst,idir,ipert,symq,ngfft,cplex,nfft,nspden,nsppol,mpi_enreg,v1scf(:,:,pcase))
1345    end if
1346  end do pcase_loop
1347 
1348  ABI_FREE(v1g)
1349  ABI_FREE(workg)
1350  ABI_FREE(workg_eq)
1351 
1352  ! Handle possible error.
1353  if (any(pflag == 0)) then
1354    write(std_out,"(2a)")&
1355      "The following perturbations cannot be recostructed by symmetry for q-point: ",trim(ktoa(qpt))
1356    do ipert=1,cryst%natom
1357      do idir=1,3
1358         if (pflag(idir, ipert) == 0) write(std_out,"(2(a,i0))")"idir= ",idir,", ipert= ",ipert
1359      end do
1360    end do
1361    write(msg,"(5a)")&
1362      "Cannot recostruct all 3*natom atomic perturbations from file",ch10,&
1363      "This usually happens when the DVDB does not contain all the independent perturbations for this q-point",ch10,&
1364      "See message above for further information."
1365    MSG_ERROR(msg)
1366  end if
1367 
1368 end subroutine v1phq_complete

m_dvdb/v1phq_rotate [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

 v1phq_rotate

FUNCTION

  Reconstruct the DFPT potential for a q-point in the BZ starting from its symmetrical image in the IBZ.

INPUTS

  cryst<crystal_t>=crystal structure parameters
  qpt_ibz(3)=q-point in the IBZ in reduced coordinates.
  ngfft=array of dimensions for different FFT grids
  cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX
  nfft=(effective) number of FFT grid points (for this proc) for the "fine" grid (see NOTES in respfn.F90)
  nspden=number of spin-density components
  nsppol=Number of independent spin polarizations
  mpi_enreg=information about MPI parallelization
  v1r_qibz(cplex*nfft,nspden,3*cryst%natom)=Array with first order potentials in real space
   for the irreducible q-point `qpt_ibz`

OUTPUT

  v1r_qbz(cplex*nfft,nspden,3*cryst%natom)=Array with first order potentials in real space for the q-point in the BZ

PARENTS

      m_dvdb

CHILDREN

      cwtime,dvdb_get_v1scf_qpt,dvdb_get_v1scf_rpt,dvdb_open_read
      hdr_fort_read,hdr_fort_write,hdr_free,irreducible_set_pert
      kpts_ibz_from_kptrlatt,littlegroup_q,wrtout,xmpi_barrier

SOURCE

1500 subroutine v1phq_rotate(cryst,qpt_ibz,isym,itimrev,g0q,ngfft,cplex,nfft,nspden,nsppol,mpi_enreg,v1r_qibz,v1r_qbz)
1501 
1502 
1503 !This section has been created automatically by the script Abilint (TD).
1504 !Do not modify the following lines by hand.
1505 #undef ABI_FUNC
1506 #define ABI_FUNC 'v1phq_rotate'
1507 !End of the abilint section
1508 
1509  implicit none
1510 
1511 !Arguments ------------------------------------
1512 !scalars
1513  integer,intent(in) :: isym,itimrev,cplex,nfft,nspden,nsppol
1514  type(crystal_t),intent(in) :: cryst
1515  type(MPI_type),intent(in) :: mpi_enreg
1516 !arrays
1517  integer,intent(in) :: g0q(3),ngfft(18)
1518  real(dp),intent(in) :: qpt_ibz(3)
1519  real(dp),intent(inout) :: v1r_qibz(cplex*nfft,nspden,3*cryst%natom)
1520  real(dp),intent(out) :: v1r_qbz(cplex*nfft,nspden,3*cryst%natom)
1521 
1522 !Local variables-------------------------------
1523 !scalars
1524  integer,parameter :: tim_fourdp0=0
1525  integer,save :: enough=0
1526  integer :: natom3,mu,ispden,idir,ipert,idir_eq,ipert_eq,mu_eq,cnt,tsign
1527 !arrays
1528  integer :: symrec_eq(3,3),sm1(3,3),l0(3) !g0_qpt(3), symrel_eq(3,3),
1529  real(dp) :: tnon(3)
1530  real(dp),allocatable :: v1g_qibz(:,:,:),workg(:,:),v1g_mu(:,:)
1531 
1532 ! *************************************************************************
1533 
1534  ABI_UNUSED(nsppol)
1535  ABI_CHECK(cplex == 2, "cplex != 2")
1536 
1537  natom3 = 3 * cryst%natom; tsign = 3-2*itimrev
1538 
1539  ! Compute IBZ potentials in G-space (results in v1g_qibz)
1540  ABI_MALLOC(v1g_qibz, (2*nfft,nspden,natom3))
1541  do mu=1,natom3
1542    do ispden=1,nspden
1543      call fourdp(cplex,v1g_qibz(:,ispden,mu),v1r_qibz(:,ispden,mu),-1,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb,tim_fourdp0)
1544    end do
1545  end do
1546 
1547  ABI_MALLOC(workg, (2*nfft,nspden))
1548  ABI_MALLOC(v1g_mu, (2*nfft,nspden))
1549 
1550  ! For each perturbation:
1551  ! FIXME: This is wrong if natom > 1
1552  symrec_eq = cryst%symrec(:,:,isym)
1553  call mati3inv(symrec_eq, sm1); sm1 = transpose(sm1)
1554 
1555  do mu=1,natom3
1556    idir = mod(mu-1, 3) + 1; ipert = (mu - idir) / 3 + 1
1557 
1558    ! Phase due to L0 + R^{-1}tau
1559    l0 = cryst%indsym(1:3,isym,ipert)
1560    tnon = l0 + matmul(transpose(symrec_eq), cryst%tnons(:,isym))
1561    ! FIXME
1562    !ABI_CHECK(all(abs(tnon) < tol12), "tnon!")
1563    if (.not.all(abs(tnon) < tol12)) then
1564      enough = enough + 1
1565      if (enough == 1) MSG_WARNING("tnon must be tested!")
1566    end if
1567 
1568    ipert_eq = cryst%indsym(4,isym,ipert)
1569 
1570    v1g_mu = zero; cnt = 0
1571    do idir_eq=1,3
1572      if (symrec_eq(idir, idir_eq) == 0) cycle
1573      mu_eq = idir_eq + (ipert_eq-1)*3
1574      cnt = cnt + 1
1575 
1576      ! Rotate in G-space and accumulate in workg
1577      call rotate_fqg(itimrev,sm1,qpt_ibz,tnon,ngfft,nfft,nspden,v1g_qibz(:,:,mu_eq),workg)
1578      v1g_mu = v1g_mu + workg * symrec_eq(idir, idir_eq)
1579    end do ! idir_eq
1580    ABI_CHECK(cnt /= 0, "cnt should not be zero!")
1581 
1582    ! Transform to real space and take into account a possible shift.
1583    ! (results are stored in v1r_qbz)
1584    do ispden=1,nspden
1585      call fourdp(cplex,v1g_mu(:,ispden),v1r_qbz(:,ispden,mu),+1,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb,tim_fourdp0)
1586      call times_eigr(-g0q, ngfft, nfft, 1, v1r_qbz(:,ispden,mu))
1587      !call times_eigr(tsign * g0q, ngfft, nfft, 1, v1r_qbz(:,ispden,mu))
1588    end do
1589 
1590    !call v1phq_symmetrize(cryst,idir,ipert,symq,ngfft,cplex,nfft,nspden,nsppol,mpi_enreg,v1r)
1591  end do ! mu
1592 
1593  ABI_FREE(workg)
1594  ABI_FREE(v1g_mu)
1595  ABI_FREE(v1g_qibz)
1596 
1597 end subroutine v1phq_rotate

m_dvdb/v1phq_symmetrize [ Functions ]

[ Top ] [ m_dvdb ] [ Functions ]

NAME

 v1phq_symmetrize

FUNCTION

  Enforce spatial-symmetry on the DFPT potential.

INPUTS

 cryst<crystal_t>=crystal structure parameters
 idir=Direction of the perturbation
 ipert=Perturbation type.
 symq(4,2,nsym)= Table computed by littlegroup_q.
   three first numbers define the G vector;
   fourth number is zero if the q-vector is not preserved, is 1 otherwise
   second index is one without time-reversal symmetry, two with time-reversal symmetry
 ngfft=array of dimensions for different FFT grids
 cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX
 nfft=(effective) number of FFT grid points (for this proc) for the "fine" grid (see NOTES in respfn.F90)
 nspden=number of spin-density components
 mpi_enreg=information about MPI parallelization

SIDE EFFECTS

  v1r(cplex*nfft,nspden,3*cryst%natom)=Array with first order potentials in real space
   Symmetrized in output.

PARENTS

      m_dvdb

CHILDREN

      cwtime,dvdb_get_v1scf_qpt,dvdb_get_v1scf_rpt,dvdb_open_read
      hdr_fort_read,hdr_fort_write,hdr_free,irreducible_set_pert
      kpts_ibz_from_kptrlatt,littlegroup_q,wrtout,xmpi_barrier

SOURCE

1635 subroutine v1phq_symmetrize(cryst,idir,ipert,symq,ngfft,cplex,nfft,nspden,nsppol,mpi_enreg,v1r)
1636 
1637 
1638 !This section has been created automatically by the script Abilint (TD).
1639 !Do not modify the following lines by hand.
1640 #undef ABI_FUNC
1641 #define ABI_FUNC 'v1phq_symmetrize'
1642 !End of the abilint section
1643 
1644  implicit none
1645 
1646 !Arguments ------------------------------------
1647 !scalars
1648  integer,intent(in) :: idir,ipert,cplex,nfft,nspden,nsppol
1649  type(crystal_t),intent(in) :: cryst
1650  type(MPI_type),intent(in) :: mpi_enreg
1651 !arrays
1652  integer,intent(in) :: symq(4,2,cryst%nsym),ngfft(18)
1653  real(dp),intent(inout) :: v1r(cplex*nfft,nspden)
1654 
1655 !Local variables-------------------------------
1656  integer,parameter :: syuse0=0,rfmeth2=2,iscf1=1
1657  integer :: nsym1,nfftot
1658 !arrays
1659  integer :: symafm1(cryst%nsym),symrel1(3,3,cryst%nsym),symrc1(3,3,cryst%nsym)
1660  integer,allocatable :: irrzon1(:,:,:),indsy1(:,:,:)
1661  real(dp) :: tnons1(3,cryst%nsym)
1662  real(dp),allocatable :: phnons1(:,:,:),v1g(:,:)
1663 
1664 ! *************************************************************************
1665 
1666  if (cryst%nsym == 1) return
1667 
1668  nfftot = product(ngfft(1:3))
1669  ABI_CHECK(nfft == nfftot, "MPI-FFT not coded")
1670 
1671  ! Symmetrize (copied from dfpt_looppert)
1672  ! Determines the set of symmetries that leaves the perturbation invariant.
1673  call littlegroup_pert(cryst%gprimd,idir,cryst%indsym,dev_null,ipert,cryst%natom,cryst%nsym,nsym1,rfmeth2,&
1674    cryst%symafm,symafm1,symq,cryst%symrec,cryst%symrel,symrel1,syuse0,cryst%tnons,tnons1,unit=dev_null)
1675 
1676  ! Set up corresponding symmetry data
1677  ABI_MALLOC(irrzon1, (nfft**(1-1/nsym1),2,(nspden/nsppol)-3*(nspden/4)))
1678  ABI_MALLOC(phnons1, (2,nfft**(1-1/nsym1),(nspden/nsppol)-3*(nspden/4)))
1679  ABI_MALLOC(indsy1,(4,nsym1,cryst%natom))
1680 
1681 
1682  call setsym(indsy1,irrzon1,iscf1,cryst%natom,nfft,ngfft,nspden,nsppol,&
1683    nsym1,phnons1,symafm1,symrc1,symrel1,tnons1,cryst%typat,cryst%xred)
1684 
1685  !if (psps%usepaw==1) then
1686  !  ! Allocate/initialize only zarot in pawang1 datastructure
1687  !  call pawang_init(pawang1,0,pawang%l_max-1,0,nsym1,0,1,0,0,0)
1688  !  call setsym_ylm(gprimd,pawang1%l_max-1,pawang1%nsym,0,rprimd,symrc1,pawang1%zarot)
1689  !end if
1690 
1691  ! FIXME Be careful here because symrhg was written for densities!
1692  ABI_CHECK(nsppol == 1 .and. nspden == 1, "symrhg was written for densities, not for potentials")
1693 
1694  ABI_MALLOC(v1g, (2,nfft))
1695  call symrhg(cplex,cryst%gprimd,irrzon1,mpi_enreg,nfft,nfftot,ngfft,nspden,nsppol,nsym1,&
1696     mpi_enreg%paral_kgb,phnons1,v1g,v1r,cryst%rprimd,symafm1,symrel1)
1697 
1698  ABI_FREE(irrzon1)
1699  ABI_FREE(phnons1)
1700  ABI_FREE(indsy1)
1701  ABI_FREE(v1g)
1702 
1703 end subroutine v1phq_symmetrize