TABLE OF CONTENTS
ABINIT/ddb_interpolate [ 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.
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (GA) This file is distributed under the terms of the GNU General Public Licence, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
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
33 #if defined HAVE_CONFIG_H 34 #include "config.h" 35 #endif 36 37 #include "abi_common.h" 38 39 40 subroutine ddb_interpolate(ifc, crystal, inp, ddb, ddb_hdr, & 41 & asrq0, prefix, comm) 42 43 44 use defs_basis 45 use m_errors 46 use m_xmpi 47 use m_profiling_abi 48 use m_ddb 49 use m_ddb_hdr 50 use m_ifc 51 use m_nctk 52 #ifdef HAVE_NETCDF 53 use netcdf 54 #endif 55 56 use m_anaddb_dataset, only : anaddb_dataset_type 57 use m_crystal, only : crystal_t 58 use m_io_tools, only : get_unit 59 use m_fstrings, only : strcat 60 use m_dynmat, only : gtdyn9, d2cart_to_red 61 62 !This section has been created automatically by the script Abilint (TD). 63 !Do not modify the following lines by hand. 64 #undef ABI_FUNC 65 #define ABI_FUNC 'ddb_interpolate' 66 use interfaces_14_hidewrite 67 use interfaces_72_response 68 !End of the abilint section 69 70 implicit none 71 72 !Arguments ------------------------------- 73 !scalars 74 type(ifc_type),intent(in) :: ifc 75 type(anaddb_dataset_type),target,intent(inout) :: inp 76 type(crystal_t),intent(in) :: crystal 77 type(ddb_type),intent(inout) :: ddb 78 type(ddb_hdr_type),intent(inout) :: ddb_hdr 79 type(asrq0_t),intent(inout) :: asrq0 80 integer,intent(in) :: comm 81 character(len=*),intent(in) :: prefix 82 !arrays 83 84 !Local variables ------------------------- 85 !scalars 86 integer,parameter :: master=0 87 integer :: unddb 88 integer :: nsym,natom,ntypat,mband,nqpt_fine 89 integer :: msize,nsize,mpert,nblok,mtyp 90 integer :: rftyp,choice 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), qred(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 = natom + 6 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_malloc(ddb_new,msize,nblok,natom,ntypat) 140 ddb_new%flg = zero 141 ddb_new%amu = ddb%amu 142 ddb_new%typ = one 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 = one 151 152 rfphon(1:2)=1; rfelfd(1:2)=0; rfstrs(1:2)=0 153 qpt_padded = zero 154 155 ddb_out_filename = strcat(prefix, "_DDB") 156 157 ddb_hdr%dscrpt = 'Interpolated DDB using interatomic force constants' 158 ddb_hdr%nblok = nblok 159 160 ! ================================================ 161 ! Interpolate the dynamical matrix at each q-point 162 ! ================================================ 163 164 qptnrm = one 165 166 do iqpt=1,nqpt_fine 167 168 ! Initialisation of the phonon wavevector 169 qpt(:)=qpt_fine(:,iqpt) 170 171 if (inp%nph1l /= 0) qptnrm(1) = inp%qnrml1(iqpt) 172 173 ! Look for the information in the DDB 174 qpt_padded(:,1) = qpt 175 call gtblk9(ddb,iblok,qpt_padded,qptnrm,rfphon,rfelfd,rfstrs,rftyp) 176 177 if (iblok /= 0) then 178 ! DEBUG 179 write(*,*) 'DDB found in file for qpt=', qpt 180 ! END DEBUG 181 182 ! q-point is present in the ddb. No interpolation needed. 183 184 !d2cart(:,1:msize) = ddb%val(:,:,iblok) 185 d2cart(1,:,:,:,:) = reshape(ddb%val(1,:,iblok), & 186 & shape = (/3,mpert,3,mpert/)) 187 d2cart(2,:,:,:,:) = reshape(ddb%val(2,:,iblok), & 188 & shape = (/3,mpert,3,mpert/)) 189 else 190 191 ! Get d2cart using the interatomic forces and the 192 ! long-range coulomb interaction through Ewald summation 193 call gtdyn9(ddb%acell,Ifc%atmfrc,Ifc%dielt,Ifc%dipdip,Ifc%dyewq0,d2cart, & 194 & crystal%gmet,ddb%gprim,ddb%mpert,natom,Ifc%nrpt,qptnrm(1), & 195 & qpt, crystal%rmet,ddb%rprim,Ifc%rpt,Ifc%trans,crystal%ucvol, & 196 & Ifc%wghatm,crystal%xred,ifc%zeff) 197 198 end if 199 200 ! Eventually impose the acoustic sum rule based on previously calculated d2asr 201 call asrq0_apply(asrq0, natom, ddb%mpert, ddb%msize, crystal%xcart, d2cart) 202 203 ! Transform d2cart into reduced coordinates. 204 call d2cart_to_red(d2cart,d2red,crystal%gprimd,crystal%rprimd,mpert, & 205 & natom,ntypat,crystal%typat,crystal%ucvol,crystal%zion) 206 207 208 ! Store the dynamical matrix into a block of the new ddb 209 jblok = iqpt 210 ddb_new%val(1,1:nsize,jblok) = reshape(d2red(1,:,:,:,:), shape = (/nsize/)) 211 ddb_new%val(2,1:nsize,jblok) = reshape(d2red(2,:,:,:,:), shape = (/nsize/)) 212 213 ! Store the q-point 214 ddb_new%qpt(1:3,jblok) = qpt 215 ddb_new%nrm(1,jblok) = qptnrm(1) 216 217 ! Set the flags 218 ii=0 219 do ipert2=1,mpert 220 do idir2=1,3 221 do ipert1=1,mpert 222 do idir1=1,3 223 ii=ii+1 224 if (ipert1<=natom.and.ipert2<=natom) then 225 ddb_new%flg(ii,jblok) = one 226 end if 227 end do 228 end do 229 end do 230 end do 231 232 end do ! iqpt 233 234 ! Copy the flags for Gamma 235 qpt_padded(:,1) = zero 236 qptnrm = one 237 call gtblk9(ddb,iblok,qpt_padded,qptnrm,rfphon,rfelfd,rfstrs,rftyp) 238 call gtblk9(ddb_new,jblok,qpt_padded,qptnrm,rfphon,rfelfd,rfstrs,rftyp) 239 240 if (iblok /= 0 .and. jblok /= 0) then 241 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 ddb_new%flg(ii,jblok) = ddb%flg(ii,iblok) 249 end do 250 end do 251 end do 252 end do 253 254 end if 255 256 257 ! ======================= 258 ! Write the new DDB files 259 ! ======================= 260 261 if (my_rank == master) then 262 263 unddb = get_unit() 264 265 ! Write the DDB header 266 call ddb_hdr_open_write(ddb_hdr, ddb_out_filename, unddb) 267 268 ! Write the whole database 269 call wrtout(std_out,' write the DDB ','COLL') 270 choice=2 271 do iblok=1,nblok 272 call ddb_write_blok(ddb_new,iblok,choice,ddb_hdr%mband,mpert,msize,& 273 & ddb_hdr%nkpt,unddb) 274 end do 275 276 ! Also write summary of bloks at the end 277 write(unddb, '(/,a)' )' List of bloks and their characteristics ' 278 choice=3 279 do iblok=1,nblok 280 call ddb_write_blok(ddb_new,iblok,choice,ddb_hdr%mband,mpert,msize,& 281 & ddb_hdr%nkpt,unddb) 282 end do 283 284 close (unddb) 285 286 #ifdef HAVE_NETCDF 287 ! Write the DDB.nc files (one for each q-point) 288 do jblok=1,nblok 289 290 write(ddb_out_nc_filename,'(2a,i5.5,a)') trim(prefix),'_qpt_',jblok,'_DDB.nc' 291 292 d2red(1,:,:,:,:) = reshape(ddb_new%val(1,1:nsize,jblok), & 293 & shape = (/3,mpert,3,mpert/)) 294 d2red(2,:,:,:,:) = reshape(ddb_new%val(2,1:nsize,jblok), & 295 & shape = (/3,mpert,3,mpert/)) 296 blkflg(:,:,:,:) = reshape(ddb_new%flg(1:nsize,jblok), & 297 & shape = (/3,mpert,3,mpert/)) 298 299 do ii=1,3 300 qred(ii) = ddb_new%qpt(ii,jblok) / ddb_new%nrm(1,jblok) 301 end do 302 303 call outddbnc(ddb_out_nc_filename,mpert,d2red,blkflg,qred,crystal) 304 305 end do 306 #endif 307 308 end if 309 310 311 ! =========== 312 ! Free memory 313 ! =========== 314 315 call ddb_free(ddb_new) 316 ABI_FREE(d2cart) 317 ABI_FREE(d2red) 318 ABI_FREE(blkflg) 319 320 end subroutine ddb_interpolate