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

PARENTS

      anaddb

CHILDREN

      d2cart_to_red,ddb_free,ddb_hdr_open_write,ddb_malloc,ddb_write_blok
      gtblk9,gtdyn9,outddbnc,wrtout

SOURCE

 85 subroutine ddb_interpolate(ifc, crystal, inp, ddb, ddb_hdr, asrq0, prefix, comm)
 86 
 87 
 88 !This section has been created automatically by the script Abilint (TD).
 89 !Do not modify the following lines by hand.
 90 #undef ABI_FUNC
 91 #define ABI_FUNC 'ddb_interpolate'
 92 !End of the abilint section
 93 
 94  implicit none
 95 
 96 !Arguments -------------------------------
 97 !scalars
 98  type(ifc_type),intent(in) :: ifc
 99  type(anaddb_dataset_type),target,intent(inout) :: inp
100  type(crystal_t),intent(in) :: crystal
101  type(ddb_type),intent(inout) :: ddb
102  type(ddb_hdr_type),intent(inout) :: ddb_hdr
103  type(asrq0_t),intent(inout) :: asrq0
104  integer,intent(in) :: comm
105  character(len=*),intent(in) :: prefix
106 !arrays
107 
108 !Local variables -------------------------
109 !scalars
110  integer,parameter :: master=0
111  integer :: unddb
112  integer :: nsym,natom,ntypat,mband,nqpt_fine
113  integer :: msize,nsize,mpert,nblok,mtyp
114  integer :: rftyp,choice
115  integer :: ii,iblok,jblok,iqpt,ipert1,ipert2,idir1,idir2
116  integer :: nprocs,my_rank
117  character(len=500) :: msg
118  character(len=fnlen) :: ddb_out_filename, ddb_out_nc_filename
119  type(ddb_type) :: ddb_new
120 !arrays
121  integer :: rfphon(4),rfelfd(4),rfstrs(4)
122  integer,allocatable :: blkflg(:,:,:,:)
123  real(dp) :: qpt(3), qptnrm(3), qpt_padded(3,3), qred(3)
124  real(dp),allocatable :: d2cart(:,:,:,:,:),d2red(:,:,:,:,:)
125  real(dp),pointer :: qpt_fine(:,:)
126 
127 ! *********************************************************************
128 
129 
130  ! Only master works for the time being
131  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
132  if (my_rank /= master) return
133 
134  ! ===================================
135  ! Copy dimensions and allocate arrays
136  ! ===================================
137 
138  write(msg, '(a,(80a),a,a,a,a)' ) ch10,('=',ii=1,80),ch10,ch10,' Treat the first list of vectors ',ch10
139  call wrtout(std_out,msg,'COLL')
140  call wrtout(ab_out,msg,'COLL')
141 
142  nullify(qpt_fine)
143  nqpt_fine = inp%nph1l
144  qpt_fine => inp%qph1l
145 
146  rftyp=inp%rfmeth
147 
148  nsym = Crystal%nsym
149  natom = Crystal%natom
150  ntypat = Crystal%ntypat
151 
152  mband = ddb_hdr%mband
153 
154  mtyp = max(ddb_hdr%mblktyp, 2)  ! Limited to 2nd derivatives of total energy
155  ddb_hdr%mblktyp = mtyp
156 
157  mpert = natom + 6
158  msize = 3 * mpert * 3 * mpert  !; if (mtyp==3) msize=msize*3*mpert
159  nsize = 3 * mpert * 3 * mpert
160  nblok = nqpt_fine
161 
162  ddb_new%nblok = nblok
163  call ddb_malloc(ddb_new,msize,nblok,natom,ntypat)
164  ddb_new%flg = 0
165  ddb_new%amu = ddb%amu
166  ddb_new%typ = 1
167  ddb_new%qpt = zero
168  ddb_new%nrm = one
169 
170  ABI_MALLOC(d2cart,(2,3,mpert,3,mpert))
171  ABI_MALLOC(d2red,(2,3,mpert,3,mpert))
172  ABI_MALLOC(blkflg,(3,mpert,3,mpert))
173 
174  blkflg = 1
175 
176  rfphon(1:2)=1; rfelfd(1:2)=0; rfstrs(1:2)=0
177  qpt_padded = zero
178 
179  ddb_out_filename = strcat(prefix, "_DDB")
180 
181  ddb_hdr%dscrpt = 'Interpolated DDB using interatomic force constants'
182  ddb_hdr%nblok = nblok
183 
184  ! ================================================
185  ! Interpolate the dynamical matrix at each q-point
186  ! ================================================
187 
188  qptnrm = one
189 
190  do iqpt=1,nqpt_fine
191 
192    ! Initialisation of the phonon wavevector
193    qpt(:)=qpt_fine(:,iqpt)
194 
195    if (inp%nph1l /= 0) qptnrm(1) = inp%qnrml1(iqpt)
196 
197    ! Look for the information in the DDB
198    qpt_padded(:,1) = qpt
199    call gtblk9(ddb,iblok,qpt_padded,qptnrm,rfphon,rfelfd,rfstrs,rftyp)
200 
201    if (iblok /= 0) then
202      ! DEBUG
203      write(*,*) 'DDB found in file for qpt=', qpt
204      ! END DEBUG
205 
206      ! q-point is present in the ddb. No interpolation needed.
207 
208      !d2cart(:,1:msize) = ddb%val(:,:,iblok)
209      d2cart(1,:,:,:,:) = reshape(ddb%val(1,:,iblok), &
210 &     shape = (/3,mpert,3,mpert/))
211      d2cart(2,:,:,:,:) = reshape(ddb%val(2,:,iblok), &
212 &     shape = (/3,mpert,3,mpert/))
213    else
214 
215      ! Get d2cart using the interatomic forces and the
216      ! long-range coulomb interaction through Ewald summation
217      call gtdyn9(ddb%acell,Ifc%atmfrc,Ifc%dielt,Ifc%dipdip,Ifc%dyewq0,d2cart, &
218 &     crystal%gmet,ddb%gprim,ddb%mpert,natom,Ifc%nrpt,qptnrm(1), &
219 &     qpt, crystal%rmet,ddb%rprim,Ifc%rpt,Ifc%trans,crystal%ucvol, &
220 &     Ifc%wghatm,crystal%xred,ifc%zeff)
221 
222    end if
223 
224    ! Eventually impose the acoustic sum rule based on previously calculated d2asr
225    call asrq0_apply(asrq0, natom, ddb%mpert, ddb%msize, crystal%xcart, d2cart)
226 
227    ! Transform d2cart into reduced coordinates.
228    call d2cart_to_red(d2cart,d2red,crystal%gprimd,crystal%rprimd,mpert, &
229 &   natom,ntypat,crystal%typat,crystal%ucvol,crystal%zion)
230 
231 
232    ! Store the dynamical matrix into a block of the new ddb
233    jblok = iqpt
234    ddb_new%val(1,1:nsize,jblok) = reshape(d2red(1,:,:,:,:), shape = (/nsize/))
235    ddb_new%val(2,1:nsize,jblok) = reshape(d2red(2,:,:,:,:), shape = (/nsize/))
236 
237    ! Store the q-point
238    ddb_new%qpt(1:3,jblok) = qpt
239    ddb_new%nrm(1,jblok) = qptnrm(1)
240 
241    ! Set the flags
242    ii=0
243    do ipert2=1,mpert
244      do idir2=1,3
245        do ipert1=1,mpert
246          do idir1=1,3
247            ii=ii+1
248            if (ipert1<=natom.and.ipert2<=natom) then
249              ddb_new%flg(ii,jblok) = 1
250            end if
251          end do
252        end do
253      end do
254    end do
255 
256  end do ! iqpt
257 
258  ! Copy the flags for Gamma
259  qpt_padded(:,1) = zero
260  qptnrm = one
261  call gtblk9(ddb,iblok,qpt_padded,qptnrm,rfphon,rfelfd,rfstrs,rftyp)
262  call gtblk9(ddb_new,jblok,qpt_padded,qptnrm,rfphon,rfelfd,rfstrs,rftyp)
263 
264  if (iblok /= 0 .and. jblok /= 0) then
265 
266    ii=0
267    do ipert2=1,mpert
268      do idir2=1,3
269        do ipert1=1,mpert
270          do idir1=1,3
271            ii=ii+1
272            ddb_new%flg(ii,jblok) = ddb%flg(ii,iblok)
273          end do
274        end do
275      end do
276    end do
277 
278  end if
279 
280 
281  ! =======================
282  ! Write the new DDB files
283  ! =======================
284 
285  if (my_rank == master) then
286 
287    unddb = get_unit()
288 
289    ! Write the DDB header
290    call ddb_hdr_open_write(ddb_hdr, ddb_out_filename, unddb)
291 
292    ! Write the whole database
293    call wrtout(std_out,' write the DDB ','COLL')
294    choice=2
295    do iblok=1,nblok
296      call ddb_write_blok(ddb_new,iblok,choice,ddb_hdr%mband,mpert,msize,&
297 &     ddb_hdr%nkpt,unddb)
298    end do
299 
300    ! Also write summary of bloks at the end
301    write(unddb, '(/,a)' )' List of bloks and their characteristics '
302    choice=3
303    do iblok=1,nblok
304      call ddb_write_blok(ddb_new,iblok,choice,ddb_hdr%mband,mpert,msize,&
305 &     ddb_hdr%nkpt,unddb)
306    end do
307 
308    close (unddb)
309 
310 #ifdef HAVE_NETCDF
311    ! Write the DDB.nc files (one for each q-point)
312    do jblok=1,nblok
313 
314      write(ddb_out_nc_filename,'(2a,i5.5,a)') trim(prefix),'_qpt_',jblok,'_DDB.nc'
315 
316      d2red(1,:,:,:,:) = reshape(ddb_new%val(1,1:nsize,jblok), &
317 &     shape = (/3,mpert,3,mpert/))
318      d2red(2,:,:,:,:) = reshape(ddb_new%val(2,1:nsize,jblok), &
319 &     shape = (/3,mpert,3,mpert/))
320      blkflg(:,:,:,:)  = reshape(ddb_new%flg(1:nsize,jblok), &
321 &     shape = (/3,mpert,3,mpert/))
322 
323      do ii=1,3
324        qred(ii) = ddb_new%qpt(ii,jblok) / ddb_new%nrm(1,jblok)
325      end do
326 
327      call outddbnc(ddb_out_nc_filename,mpert,d2red,blkflg,qred,crystal)
328 
329    end do
330 #endif
331  end if
332 
333  ! ===========
334  ! Free memory
335  ! ===========
336 
337  call ddb_free(ddb_new)
338  ABI_FREE(d2cart)
339  ABI_FREE(d2red)
340  ABI_FREE(blkflg)
341 
342 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-2018 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 .

PARENTS

CHILDREN

SOURCE

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

ABINIT/outddbnc [ Functions ]

[ Top ] [ Functions ]

NAME

 outddbnc

FUNCTION

 Write Derivative DataBase in NetCDF format.
 See ~abinit/scripts/post_processing/merge_ddb_nc.py
 for a merging utility.

INPUTS

  filename=name of the DDB.nc file to be written.
  Crystal <type(crystal_t)>=info on the crystal
  natom = number of atoms in the unit cell.
  mpert = maximum number of perturbations.
  qpt(3)=curret q-point, in reduced coordinates.
  d2matr(2,3,mpert,3,mpert)=second-order derivative of the total energy
  blkflg(3,mpert,3,mpert)=mask telling whether each element is computed (1) or not (0).

OUTPUT

  Only writing

PARENTS

      ddb_interpolate,respfn

CHILDREN

SOURCE

375 subroutine outddbnc (filename, mpert, d2matr, blkflg, qpt, Crystal)
376 
377  !use defs_datatypes
378  !use defs_abitypes
379 
380 !This section has been created automatically by the script Abilint (TD).
381 !Do not modify the following lines by hand.
382 #undef ABI_FUNC
383 #define ABI_FUNC 'outddbnc'
384 !End of the abilint section
385 
386  implicit none
387 
388 !Arguments -------------------------------
389 !scalars
390  character(len=*),intent(in) :: filename
391  integer,intent(in) :: mpert
392  integer,intent(in) :: blkflg(3,mpert,3,mpert)
393  real(dp),intent(in) :: d2matr(2,3,mpert,3,mpert)
394  real(dp),intent(in) :: qpt(3)
395  type(crystal_t), intent(in) :: Crystal
396 
397 !Local variables -------------------------
398 !scalars
399  integer :: natom
400  integer :: ncid, ncerr
401  integer :: cplex, cart_dir, one_dim
402  integer :: ipert1, ipert2, idir1, idir2
403  integer,allocatable :: dynmat_mask(:,:,:,:)
404  integer,allocatable :: born_effective_charge_tensor_mask(:,:,:)
405  real(dp),allocatable :: dynmat(:,:,:,:,:)
406  real(dp),allocatable :: born_effective_charge_tensor(:,:,:)
407 
408 ! *********************************************************************
409 
410 #ifdef HAVE_NETCDF
411 
412  natom = Crystal%natom
413 
414  ABI_ALLOCATE(dynmat, (2,3,natom,3,natom))
415  ABI_ALLOCATE(dynmat_mask, (3,natom,3,natom))
416  ABI_ALLOCATE(born_effective_charge_tensor, (3,natom,3))
417  ABI_ALLOCATE(born_effective_charge_tensor_mask, (3,natom,3))
418 
419  ! Initialize NetCDF file.
420  NCF_CHECK(nctk_open_create(ncid, filename, xmpi_comm_self))
421 
422  ! ------------------------------
423  ! Construct local DDB
424  ! ------------------------------
425 
426  ! Construct the dynamical matrix
427  do ipert1=1,natom
428    do idir1=1,3
429      do ipert2=1,natom
430        do idir2=1,3
431          dynmat_mask(idir1,ipert1,idir2,ipert2) = blkflg(idir1,ipert1,idir2,ipert2)
432          if (blkflg(idir1,ipert1,idir2,ipert2)==1) then
433            dynmat(1,idir1,ipert1,idir2,ipert2) = d2matr(1,idir1,ipert1,idir2,ipert2)
434            dynmat(2,idir1,ipert1,idir2,ipert2) = d2matr(2,idir1,ipert1,idir2,ipert2)
435          else
436            dynmat(1,idir1,ipert1,idir2,ipert2) = zero
437            dynmat(2,idir1,ipert1,idir2,ipert2) = zero
438          end if
439        end do
440      end do
441    end do
442  end do
443 
444  ! Construct the Born effective charge tensor
445  ipert2 = natom + 2
446  do ipert1=1,natom
447    do idir1=1,3
448      do idir2=1,3
449        born_effective_charge_tensor_mask(idir1,ipert1,idir2) = blkflg(idir1,ipert1,idir2,ipert2)
450        if (blkflg(idir1,ipert1,idir2,ipert2)==1) then
451             ! This is a real quantity
452          born_effective_charge_tensor(idir1,ipert1,idir2) = d2matr(1,idir1,ipert1,idir2,ipert2)
453        else
454          born_effective_charge_tensor(idir1,ipert1,idir2) = zero
455        end if
456      end do
457    end do
458  end do
459 
460  ! TODO: also store the dielectric matrix
461 
462  ! ------------------------------
463  ! Write crystal info
464  ! ------------------------------
465  ncerr = crystal_ncwrite(Crystal, ncid)
466  NCF_CHECK(ncerr)
467 
468 
469  ! ------------------------------
470  ! Write DDB
471  ! ------------------------------
472 
473  ! Write the dimensions specified by ETSF
474  one_dim = 1
475  cplex = 2
476  cart_dir = 3
477 
478  ncerr = nctk_def_dims(ncid, [&
479 & nctkdim_t('current_one_dim', one_dim), &
480 & nctkdim_t('number_of_atoms', natom), &
481 & nctkdim_t('number_of_cartesian_directions', cart_dir), &
482 & nctkdim_t('number_of_perturbations', mpert), &
483 & nctkdim_t('cplex',cplex)], defmode=.True.)
484  NCF_CHECK(ncerr)
485 
486  ! Create the arrays
487  ncerr = nctk_def_arrays(ncid, [&
488  nctkarr_t('atomic_masses_amu', "dp", 'number_of_atom_species'),&
489  nctkarr_t('q_point_reduced_coord', "dp", 'number_of_cartesian_directions'),&
490  nctkarr_t('second_derivative_of_energy', "dp", 'cplex, &
491 & number_of_cartesian_directions, number_of_atoms, &
492 & number_of_cartesian_directions, number_of_atoms'), &
493  nctkarr_t('second_derivative_of_energy_mask', "i", '&
494 & number_of_cartesian_directions, number_of_atoms, &
495 & number_of_cartesian_directions, number_of_atoms'), &
496  nctkarr_t('born_effective_charge_tensor', "dp", '&
497 & number_of_cartesian_directions, number_of_atoms, &
498 & number_of_cartesian_directions'), &
499  nctkarr_t('born_effective_charge_tensor_mask', "i", ' &
500 & number_of_cartesian_directions, number_of_atoms, &
501 & number_of_cartesian_directions')])
502  NCF_CHECK(ncerr)
503 
504 ! Write data
505  NCF_CHECK(nctk_set_datamode(ncid))
506  NCF_CHECK(nf90_put_var(ncid, vid('atomic_masses_amu'), Crystal%amu))
507  NCF_CHECK(nf90_put_var(ncid, vid('q_point_reduced_coord'), qpt))
508  NCF_CHECK(nf90_put_var(ncid, vid('second_derivative_of_energy'), dynmat))
509  NCF_CHECK(nf90_put_var(ncid, vid('second_derivative_of_energy_mask'), dynmat_mask))
510  NCF_CHECK(nf90_put_var(ncid, vid('born_effective_charge_tensor'), born_effective_charge_tensor))
511  NCF_CHECK(nf90_put_var(ncid, vid('born_effective_charge_tensor_mask'), born_effective_charge_tensor_mask))
512 
513  ! Close file
514  NCF_CHECK(nf90_close(ncid))
515 
516  ! Deallocate stuff
517  ABI_FREE(dynmat)
518  ABI_FREE(dynmat_mask)
519  ABI_FREE(born_effective_charge_tensor)
520  ABI_FREE(born_effective_charge_tensor_mask)
521 
522 #else
523  MSG_ERROR("NETCDF support required to write DDB.nc file.")
524 #endif
525 
526  contains
527    integer function vid(vname)
528 
529 
530 !This section has been created automatically by the script Abilint (TD).
531 !Do not modify the following lines by hand.
532 #undef ABI_FUNC
533 #define ABI_FUNC 'vid'
534 !End of the abilint section
535 
536    character(len=*),intent(in) :: vname
537    vid = nctk_idname(ncid, vname)
538  end function vid
539 
540 end subroutine outddbnc