TABLE OF CONTENTS


ABINIT/ddb_interpolate [ Functions ]

[ Top ] [ Functions ]

NAME

 ddb_interpolate

FUNCTION

 Interpolate the ddb onto a fine set of q-points using
 the interatomic force constants and write the result in a DDB file.

INPUTS

OUTPUT

NOTES

SOURCE

 71 subroutine ddb_interpolate(ifc, crystal, inp, ddb, ddb_hdr, asrq0, prefix, comm)
 72 
 73 !Arguments -------------------------------
 74 !scalars
 75  type(ifc_type),intent(in) :: ifc
 76  type(anaddb_dataset_type),target,intent(inout) :: inp
 77  type(crystal_t),intent(in) :: crystal
 78  type(ddb_type),intent(inout) :: ddb
 79  type(ddb_hdr_type),intent(inout) :: ddb_hdr
 80  type(asrq0_t),intent(inout) :: asrq0
 81  integer,intent(in) :: comm
 82  character(len=*),intent(in) :: prefix
 83 !arrays
 84 
 85 !Local variables -------------------------
 86 !scalars
 87  integer,parameter :: master=0
 88  integer :: nsym,natom,ntypat,mband,nqpt_fine
 89  integer :: msize,nsize,mpert,nblok,mtyp
 90  integer :: rftyp
 91  integer :: ii,iblok,jblok,iqpt,ipert1,ipert2,idir1,idir2
 92  integer :: nprocs,my_rank
 93  character(len=500) :: msg
 94  character(len=fnlen) :: ddb_out_filename, ddb_out_nc_filename
 95  type(ddb_type) :: ddb_new
 96 !arrays
 97  integer :: rfphon(4),rfelfd(4),rfstrs(4)
 98  integer,allocatable :: blkflg(:,:,:,:)
 99  real(dp) :: qpt(3), qptnrm(3), qpt_padded(3,3)
100  real(dp),allocatable :: d2cart(:,:,:,:,:),d2red(:,:,:,:,:)
101  real(dp),pointer :: qpt_fine(:,:)
102 
103 ! *********************************************************************
104 
105 
106  ! Only master works for the time being
107  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
108  if (my_rank /= master) return
109 
110  ! ===================================
111  ! Copy dimensions and allocate arrays
112  ! ===================================
113 
114  write(msg, '(a,(80a),a,a,a,a)' ) ch10,('=',ii=1,80),ch10,ch10,' Treat the first list of vectors ',ch10
115  call wrtout(std_out,msg,'COLL')
116  call wrtout(ab_out,msg,'COLL')
117 
118  nullify(qpt_fine)
119  nqpt_fine = inp%nph1l
120  qpt_fine => inp%qph1l
121 
122  rftyp=inp%rfmeth
123 
124  nsym = Crystal%nsym
125  natom = Crystal%natom
126  ntypat = Crystal%ntypat
127 
128  mband = ddb_hdr%mband
129 
130  mtyp = max(ddb_hdr%mblktyp, 2)  ! Limited to 2nd derivatives of total energy
131  ddb_hdr%mblktyp = mtyp
132 
133  mpert = ddb%mpert
134  msize = 3 * mpert * 3 * mpert  !; if (mtyp==3) msize=msize*3*mpert
135  nsize = 3 * mpert * 3 * mpert
136  nblok = nqpt_fine
137 
138  ddb_new%nblok = nblok
139  call ddb_new%malloc(msize,nblok,natom,ntypat,mpert)
140  ddb_new%flg = 0
141  ddb_new%amu = ddb%amu
142  ddb_new%typ = 1
143  ddb_new%qpt = zero
144  ddb_new%nrm = one
145 
146  ABI_MALLOC(d2cart,(2,3,mpert,3,mpert))
147  ABI_MALLOC(d2red,(2,3,mpert,3,mpert))
148  ABI_MALLOC(blkflg,(3,mpert,3,mpert))
149 
150  blkflg = 1
151 
152  rfphon(1:2)=1; rfelfd(1:2)=0; rfstrs(1:2)=0
153  qpt_padded = zero
154 
155 
156  ddb_hdr%dscrpt = 'Interpolated DDB using interatomic force constants'
157  ddb_hdr%nblok = nblok
158 
159  ! ================================================
160  ! Interpolate the dynamical matrix at each q-point
161  ! ================================================
162 
163  qptnrm = one
164 
165  do iqpt=1,nqpt_fine
166 
167    ! Initialisation of the phonon wavevector
168    qpt(:)=qpt_fine(:,iqpt)
169 
170    if (inp%nph1l /= 0) qptnrm(1) = inp%qnrml1(iqpt)
171 
172    ! Look for the information in the DDB
173    qpt_padded(:,1) = qpt
174    call ddb%get_block(iblok,qpt_padded,qptnrm,rfphon,rfelfd,rfstrs,rftyp)
175 
176    if (iblok /= 0) then
177 
178      ! q-point is present in the ddb. No interpolation needed.
179 
180      d2cart(1,:,:,:,:) = reshape(ddb%val(1,:,iblok), shape = (/3,mpert,3,mpert/))
181      d2cart(2,:,:,:,:) = reshape(ddb%val(2,:,iblok), shape = (/3,mpert,3,mpert/))
182    else
183 
184      ! Get d2cart using the interatomic forces and the
185      ! long-range coulomb interaction through Ewald summation
186      call gtdyn9(ddb%acell,Ifc%atmfrc,Ifc%dielt,Ifc%dipdip,Ifc%dyewq0,d2cart, &
187       crystal%gmet,ddb%gprim,ddb%mpert,natom,Ifc%nrpt,qptnrm(1), &
188       qpt, crystal%rmet,ddb%rprim,Ifc%rpt,Ifc%trans,crystal%ucvol, &
189       Ifc%wghatm,crystal%xred,ifc%zeff,ifc%qdrp_cart,ifc%ewald_option,xmpi_comm_self, &
190       dipquad=Ifc%dipquad,quadquad=Ifc%quadquad)
191 
192    end if
193 
194    ! Eventually impose the acoustic sum rule based on previously calculated d2asr
195    call asrq0%apply(natom, ddb%mpert, ddb%msize, crystal%xcart, d2cart)
196 
197    ! Transform d2cart into reduced coordinates.
198    call d2cart_to_red(d2cart,d2red,crystal%gprimd,crystal%rprimd,mpert, &
199 &   natom,ntypat,crystal%typat,crystal%ucvol,crystal%zion)
200 
201    ! Store the dynamical matrix into a block of the new ddb
202    jblok = iqpt
203    ddb_new%val(1,1:nsize,jblok) = reshape(d2red(1,:,:,:,:), shape = (/nsize/))
204    ddb_new%val(2,1:nsize,jblok) = reshape(d2red(2,:,:,:,:), shape = (/nsize/))
205 
206    ! Store the q-point
207    ddb_new%qpt(1:3,jblok) = qpt
208    ddb_new%nrm(1,jblok) = qptnrm(1)
209 
210    ! Set the flags
211    ii=0
212    do ipert2=1,mpert
213      do idir2=1,3
214        do ipert1=1,mpert
215          do idir1=1,3
216            ii=ii+1
217            if (ipert1<=natom.and.ipert2<=natom) then
218              ddb_new%flg(ii,jblok) = 1
219            end if
220          end do
221        end do
222      end do
223    end do
224 
225  end do ! iqpt
226 
227  ! Copy the flags for Gamma
228  qpt_padded(:,1) = zero
229  qptnrm = one
230  call ddb%get_block(iblok,qpt_padded,qptnrm,rfphon,rfelfd,rfstrs,rftyp)
231  call ddb_new%get_block(jblok,qpt_padded,qptnrm,rfphon,rfelfd,rfstrs,rftyp)
232 
233  if (iblok /= 0 .and. jblok /= 0) then
234 
235    ii=0
236    do ipert2=1,mpert
237      do idir2=1,3
238        do ipert1=1,mpert
239          do idir1=1,3
240            ii=ii+1
241            ddb_new%flg(ii,jblok) = ddb%flg(ii,iblok)
242          end do
243        end do
244      end do
245    end do
246 
247  end if
248 
249 
250  ! =======================
251  ! Write the new DDB files
252  ! =======================
253 
254  if (my_rank == master) then
255 
256    ddb_out_filename = strcat(prefix, "_DDB")
257 
258    call ddb_new%write_txt(ddb_hdr, ddb_out_filename)
259 
260    ! Write one separate nc file for each q-point
261    do jblok=1,nblok
262      write(ddb_out_nc_filename,'(2a,i5.5,a)') trim(prefix),'_qpt_',jblok,'_DDB.nc'
263      call ddb_new%write_nc(ddb_hdr, ddb_out_nc_filename, jblok)
264    end do
265 
266  end if
267 
268  ! ===========
269  ! Free memory
270  ! ===========
271 
272  call ddb_new%free()
273  ABI_FREE(d2cart)
274  ABI_FREE(d2red)
275  ABI_FREE(blkflg)
276 
277 end subroutine ddb_interpolate

ABINIT/m_ddb_interpolate [ Modules ]

[ Top ] [ Modules ]

NAME

  m_ddb_interpolate

FUNCTION

 Interpolate the ddb onto a fine set of q-points using
 the interatomic force constants and write the result in a DDB file.

COPYRIGHT

  Copyright (C) 2008-2022 ABINIT group (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 .

SOURCE

17 #if defined HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20 
21 #include "abi_common.h"
22 
23 module m_ddb_interpolate
24 
25  use defs_basis
26  use m_errors
27  use m_xmpi
28  use m_abicore
29  use m_ddb
30  use m_ddb_hdr
31  use m_ifc
32  use m_nctk
33 #ifdef HAVE_NETCDF
34  use netcdf
35 #endif
36 
37  use m_anaddb_dataset, only : anaddb_dataset_type
38  use m_crystal,        only : crystal_t
39  use m_io_tools,       only : get_unit
40  use m_fstrings,       only : strcat
41  use m_dynmat,         only : gtdyn9, d2cart_to_red
42 
43  implicit none
44 
45  private