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