TABLE OF CONTENTS
ABINIT/m_gkk [ Modules ]
NAME
FUNCTION
Tools for the computation of electron-phonon coupling matrix elements (gkk)
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (GKA, MG) 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
15 #if defined HAVE_CONFIG_H 16 #include "config.h" 17 #endif 18 19 #include "abi_common.h" 20 21 module m_gkk 22 23 use defs_basis 24 use m_abicore 25 use m_xmpi 26 use m_errors 27 use m_dtset 28 use m_ifc 29 use m_ebands 30 use m_ddb 31 use m_dvdb 32 use m_fft 33 use m_hamiltonian 34 use m_pawcprj 35 use m_wfk 36 use m_nctk 37 use m_dtfil 38 #ifdef HAVE_NETCDF 39 use netcdf 40 #endif 41 42 use defs_abitypes, only : MPI_type 43 use m_time, only : cwtime, sec2str 44 use m_io_tools, only : iomode_from_fname 45 use m_fstrings, only : itoa, sjoin, ktoa, ltoa, strcat 46 use m_symtk, only : littlegroup_q 47 use m_fftcore, only : get_kg 48 use defs_datatypes, only : ebands_t, pseudopotential_type 49 use m_crystal, only : crystal_t 50 use m_bz_mesh, only : findqg0 51 use m_cgtools, only : dotprod_g 52 use m_kg, only : getph 53 use m_pawang, only : pawang_type 54 use m_pawrad, only : pawrad_type 55 use m_pawtab, only : pawtab_type 56 use m_pawfgr, only : pawfgr_type 57 use m_eig2d, only : gkk_t, gkk_init, gkk_ncwrite, gkk_free 58 use m_wfd, only : wfd_init, wfd_t 59 use m_getgh1c, only : getgh1c, rf_transgrid_and_pack, getgh1c_setup 60 61 implicit none 62 63 private
m_gkk/eph_gkk [ Functions ]
[ Top ] [ m_gkk ] [ Functions ]
NAME
eph_gkk
FUNCTION
Compute electron-phonon coupling matrix elements.
INPUTS
wk0_path=String with the path to the GS unperturbed WFK file. ngfft(18),ngfftf(18)=Coarse and Fine FFT meshes. dtset<dataset_type>=All input variables for this dataset. ebands<ebands_t>=The GS KS band structure (energies, occupancies, k-weights...) dvdb<dbdb_type>=Database with the DFPT SCF potentials. ifc<ifc_type>=interatomic force constants and corresponding real space grid info. pawfgr <type(pawfgr_type)>=fine grid parameters and related data pawang<pawang_type)>=PAW angular mesh and related data. pawrad(ntypat*usepaw)<pawrad_type>=Paw radial mesh and related data. pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data. psps<pseudopotential_type>=Variables related to pseudopotentials. comm=MPI communicator.
OUTPUT
NOTES
SOURCE
99 subroutine eph_gkk(wfk0_path,wfq_path,dtfil,ngfft,ngfftf,dtset,cryst,ebands_k,ebands_kq,dvdb,ifc,& 100 pawfgr,pawang,pawrad,pawtab,psps,mpi_enreg,comm) 101 102 !Arguments ------------------------------------ 103 !scalars 104 character(len=*),intent(in) :: wfk0_path, wfq_path 105 integer,intent(in) :: comm 106 type(datafiles_type),intent(in) :: dtfil 107 type(dataset_type),intent(in) :: dtset 108 type(crystal_t),intent(in) :: cryst 109 type(ebands_t),intent(in) :: ebands_k, ebands_kq 110 type(dvdb_t),target,intent(inout) :: dvdb 111 type(pawang_type),intent(in) :: pawang 112 type(pseudopotential_type),intent(in) :: psps 113 type(pawfgr_type),intent(in) :: pawfgr 114 type(ifc_type),intent(in) :: ifc 115 type(mpi_type),intent(inout) :: mpi_enreg 116 !arrays 117 integer,intent(in) :: ngfft(18),ngfftf(18) 118 type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw) 119 type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw) 120 121 !Local variables ------------------------------ 122 !scalars 123 integer,parameter :: tim_getgh1c=1, berryopt0=0, useylmgr1=0, master=0, qptopt1 = 1 124 integer :: my_rank,nproc,mband,mband_kq,my_minb,my_maxb,nsppol,nkpt,nkpt_kq,idir,ipert 125 integer :: cplex,db_iqpt,natom,natom3,ipc,nspinor 126 integer :: ib1,ib2,band,ik,ikq,timerev_q 127 integer :: spin,istwf_k,istwf_kq,npw_k,npw_kq, comm_rpt 128 integer :: mpw,mpw_k,mpw_kq,ierr,my_kstart,my_kstop,ncid 129 integer :: n1,n2,n3,n4,n5,n6,nspden,ncerr 130 integer :: sij_opt,usecprj,usevnl,optlocal,optnl,opt_gvnlx1 131 integer :: nfft,nfftf,mgfft,mgfftf,nkpg,nkpg1, interpolated 132 real(dp) :: cpu,wall,gflops,ecut,eshift,eig0nk,dotr,doti 133 logical :: i_am_master, gen_eigenpb 134 type(wfd_t) :: wfd_k, wfd_kq 135 type(gs_hamiltonian_type) :: gs_hamkq 136 type(rf_hamiltonian_type) :: rf_hamkq 137 type(gkk_t) :: gkk2d 138 character(len=500) :: msg, what 139 character(len=fnlen) :: fname, gkkfilnam 140 !arrays 141 integer :: g0_k(3),symq(4,2,cryst%nsym) 142 integer,allocatable :: kg_k(:,:),kg_kq(:,:),nband(:,:),nband_kq(:,:),blkflg(:,:), wfd_istwfk(:) 143 real(dp) :: kk(3),kq(3),qpt(3),phfrq(3*cryst%natom) 144 real(dp) :: dvdb_qdamp(1) 145 real(dp),allocatable :: displ_cart(:,:,:),displ_red(:,:,:), eigens_kq(:,:,:) 146 real(dp),allocatable :: grad_berry(:,:),kinpw1(:),kpg1_k(:,:),kpg_k(:,:),dkinpw(:) 147 real(dp),allocatable :: ffnlk(:,:,:,:),ffnl1(:,:,:,:),ph3d(:,:,:),ph3d1(:,:,:) 148 real(dp),allocatable :: v1scf(:,:,:,:),gkk(:,:,:,:,:) 149 real(dp),allocatable :: bras(:,:,:),kets(:,:,:),h1_kets(:,:,:) 150 real(dp),allocatable :: ph1d(:,:),vlocal(:,:,:,:),vlocal1(:,:,:,:,:) 151 real(dp),allocatable :: ylm_kq(:,:),ylm_k(:,:),ylmgr_kq(:,:,:) 152 real(dp),allocatable :: dummy_vtrial(:,:),gvnlx1(:,:), gs1c(:,:), gkq_atm(:,:,:,:) 153 logical,allocatable :: bks_mask(:,:,:),bks_mask_kq(:,:,:),keep_ur(:,:,:),keep_ur_kq(:,:,:) 154 type(pawcprj_type),allocatable :: cwaveprj0(:,:) !natom,nspinor*usecprj) 155 156 !************************************************************************ 157 158 what = "(GKK files)"; if (dtset%eph_task == -2) what = "GKQ file" 159 write(msg, '(3a)') " Computation of electron-phonon coupling matrix elements ", trim(what), ch10 160 call wrtout([std_out, ab_out], msg, do_flush=.True.) 161 162 if (psps%usepaw == 1) then 163 ABI_ERROR("PAW not implemented") 164 ABI_UNUSED((/pawang%nsym, pawrad(1)%mesh_size/)) 165 end if 166 167 my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm); i_am_master = my_rank == master 168 169 ! Copy important dimensions 170 natom = cryst%natom 171 natom3 = 3 * natom 172 nsppol = ebands_k%nsppol 173 nspinor = ebands_k%nspinor 174 nspden = dtset%nspden 175 nkpt = ebands_k%nkpt 176 mband = ebands_k%mband 177 nkpt_kq = ebands_kq%nkpt 178 mband_kq = ebands_kq%mband 179 ecut = dtset%ecut 180 !write(std_out, *)"ebands dims (b, k, s): ", ebands_k%mband, ebands_k%nkpt, ebands_k%nsppol 181 !write(std_out, *)"ebands_kq dims (b, k, s): ", ebands_kq%mband, ebands_kq%nkpt, ebands_kq%nsppol 182 183 qpt = dtset%qptn(:) 184 185 nfftf = product(ngfftf(1:3)); mgfftf = maxval(ngfftf(1:3)) 186 nfft = product(ngfft(1:3)) ; mgfft = maxval(ngfft(1:3)) 187 n1=ngfft(1); n2=ngfft(2); n3=ngfft(3) 188 n4=ngfft(4); n5=ngfft(5); n6=ngfft(6) 189 190 ! Open the DVDB file 191 call dvdb%open_read(ngfftf, xmpi_comm_self) 192 193 ! Initialize the wave function descriptors. 194 ! For the time being, no memory distribution, each node has the full set of states. 195 my_minb = 1; my_maxb = mband 196 197 ABI_MALLOC(nband, (nkpt, nsppol)) 198 ABI_MALLOC(bks_mask,(mband, nkpt, nsppol)) 199 ABI_MALLOC(keep_ur,(mband, nkpt ,nsppol)) 200 nband=mband; bks_mask=.False.; keep_ur=.False. 201 202 ABI_MALLOC(nband_kq, (nkpt_kq, nsppol)) 203 ABI_MALLOC(bks_mask_kq,(mband_kq, nkpt_kq, nsppol)) 204 ABI_MALLOC(keep_ur_kq,(mband_kq, nkpt_kq ,nsppol)) 205 nband_kq=mband_kq; bks_mask_kq=.False.; keep_ur_kq=.False. 206 207 ! Distribute the k-points over the processors 208 call xmpi_split_work(nkpt,comm,my_kstart,my_kstop) 209 do ik=1,nkpt 210 if (.not. (ik >= my_kstart .and. ik <= my_kstop)) cycle 211 kk = ebands_k%kptns(:,ik) 212 kq = kk + qpt 213 ! Find the index of the k+q point 214 call findqg0(ikq,g0_k,kq,nkpt_kq,ebands_kq%kptns(:,:), [1,1,1]) 215 bks_mask(:,ik,:) = .True. 216 bks_mask_kq(:,ikq,:) = .True. 217 end do 218 219 ! Impose istwfk=1 for all k points. This is also done in respfn (see inkpts) 220 ! wfd_read_wfk will handle a possible conversion if WFK contains istwfk /= 1. 221 ABI_MALLOC(wfd_istwfk, (nkpt)) 222 wfd_istwfk = 1 223 224 ! Initialize the wavefunction descriptors 225 call wfd_init(wfd_k,cryst,pawtab,psps,keep_ur,mband,nband,nkpt,nsppol,bks_mask,& 226 nspden,nspinor,ecut,dtset%ecutsm,dtset%dilatmx,wfd_istwfk,ebands_k%kptns,ngfft,& 227 dtset%nloalg,dtset%prtvol,dtset%pawprtvol,comm) 228 ABI_FREE(wfd_istwfk) 229 230 call wfd_k%print(header="Wavefunctions on the k-points grid",mode_paral='PERS') 231 232 ABI_MALLOC(wfd_istwfk, (nkpt_kq)) 233 wfd_istwfk = 1 234 235 call wfd_init(wfd_kq,cryst,pawtab,psps,keep_ur_kq,mband_kq,nband_kq,nkpt_kq,nsppol,bks_mask_kq,& 236 nspden,nspinor,ecut,dtset%ecutsm,dtset%dilatmx,wfd_istwfk,ebands_kq%kptns,ngfft,& 237 dtset%nloalg,dtset%prtvol,dtset%pawprtvol,comm) 238 ABI_FREE(wfd_istwfk) 239 240 call wfd_kq%print(header="Wavefunctions on the q-shifted k-points grid",mode_paral='PERS') 241 242 ABI_FREE(nband) 243 ABI_FREE(bks_mask) 244 ABI_FREE(keep_ur) 245 ABI_FREE(nband_kq) 246 ABI_FREE(bks_mask_kq) 247 ABI_FREE(keep_ur_kq) 248 249 ! Read wavefunctions on the k-points grid and q-shifted k-points grid. 250 call wfd_k%read_wfk(wfk0_path, iomode_from_fname(wfk0_path)) 251 252 call wfd_kq%read_wfk(wfq_path, iomode_from_fname(wfq_path)) 253 254 ! ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information on the coarse grid. 255 ABI_MALLOC(ph1d, (2,3*(2*mgfft+1)*natom)) 256 call getph(cryst%atindx,natom,n1,n2,n3,ph1d,cryst%xred) 257 258 ! Find the appropriate value of mpw 259 call find_mpw(mpw_k, ebands_k%kptns(:,:), nsppol, nkpt, cryst%gmet,ecut,comm) 260 call find_mpw(mpw_kq, ebands_kq%kptns(:,:), nsppol, nkpt_kq, cryst%gmet,ecut,comm) 261 mpw = max(mpw_k, mpw_kq) 262 263 ! Allow PW-arrays dimensioned with mpw 264 ABI_MALLOC(kg_k, (3, mpw)) 265 ABI_MALLOC(kg_kq, (3, mpw)) 266 267 ! Spherical Harmonics for useylm==1. 268 ABI_MALLOC(ylm_k,(mpw, psps%mpsang*psps%mpsang*psps%useylm)) 269 ABI_MALLOC(ylm_kq,(mpw, psps%mpsang*psps%mpsang*psps%useylm)) 270 ABI_MALLOC(ylmgr_kq,(mpw, 3, psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)) 271 272 ! TODO FOR PAW 273 usecprj = 0 274 ABI_MALLOC(cwaveprj0, (natom, nspinor*usecprj)) 275 276 ! Prepare call to getgh1c 277 usevnl = 0 278 optlocal = 1 ! local part of H^(1) is computed in gh1c=<G|H^(1)|C> 279 optnl = 2 ! non-local part of H^(1) is totally computed in gh1c=<G|H^(1)|C> 280 opt_gvnlx1 = 0 ! gvnlx1 is output 281 ABI_MALLOC(gvnlx1, (2,usevnl)) 282 ABI_MALLOC(grad_berry, (2,nspinor*(berryopt0/4))) 283 284 ! This part is taken from dfpt_vtorho 285 !==== Initialize most of the Hamiltonian (and derivative) ==== 286 !1) Allocate all arrays and initialize quantities that do not depend on k and spin. 287 !2) Perform the setup needed for the non-local factors: 288 !* Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk. 289 !* PAW: Initialize the overlap coefficients and allocate the Dij coefficients. 290 291 call init_hamiltonian(gs_hamkq,psps,pawtab,nspinor,nsppol,nspden,natom,& 292 dtset%typat,cryst%xred,nfft,mgfft,ngfft,cryst%rprimd,dtset%nloalg,& 293 usecprj=usecprj,ph1d=ph1d,nucdipmom=dtset%nucdipmom,gpu_option=dtset%gpu_option,& 294 comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 295 296 ! Allocate vlocal. Note nvloc 297 ABI_MALLOC(vlocal,(n4,n5,n6,gs_hamkq%nvloc)) 298 ! Allocate work space arrays. 299 ABI_MALLOC(blkflg, (natom3,natom3)) 300 ABI_CALLOC(dummy_vtrial, (nfftf,nspden)) 301 302 call cwtime(cpu, wall, gflops, "start") 303 304 interpolated = 0 305 if (dtset%eph_use_ftinterp /= 0) then 306 ABI_WARNING(sjoin("Enforcing FT interpolation for q-point", ktoa(qpt))) 307 comm_rpt = xmpi_comm_self 308 call dvdb%ftinterp_setup(dtset%ddb_ngqpt, qptopt1, 1, dtset%ddb_shiftq, nfftf, ngfftf, comm_rpt) 309 cplex = 2 310 ABI_MALLOC(v1scf, (cplex, nfftf, nspden, dvdb%my_npert)) 311 call dvdb%ftinterp_qpt(qpt, nfftf, ngfftf, v1scf, dvdb%comm_rpt) 312 interpolated = 1 313 else 314 ! Find the index of the q-point in the DVDB. 315 db_iqpt = dvdb%findq(qpt) 316 if (db_iqpt /= -1) then 317 if (dtset%prtvol > 0) call wrtout(std_out, sjoin("Found: ",ktoa(qpt)," in DVDB with index ",itoa(db_iqpt))) 318 ! Read or reconstruct the dvscf potentials for all 3*natom perturbations. 319 ! This call allocates v1scf(cplex, nfftf, nspden, 3*natom)) 320 call dvdb%readsym_allv1(db_iqpt, cplex, nfftf, ngfftf, v1scf, comm) 321 else 322 ABI_WARNING(sjoin("Cannot find q-point:", ktoa(qpt), "in DVDB file")) 323 end if 324 end if 325 326 ! Examine the symmetries of the q wavevector 327 call littlegroup_q(cryst%nsym,qpt,symq,cryst%symrec,cryst%symafm,timerev_q,prtvol=dtset%prtvol) 328 329 ! Allocate vlocal1 with correct cplex. Note nvloc 330 ABI_MALLOC_OR_DIE(vlocal1,(cplex*n4,n5,n6,gs_hamkq%nvloc,natom3), ierr) 331 332 ABI_MALLOC(displ_cart, (2,3*cryst%natom,3*cryst%natom)) 333 ABI_MALLOC(displ_red, (2,3*cryst%natom,3*cryst%natom)) 334 335 if (dtset%eph_task == 2) then 336 ! Write GKK files (1 file for perturbation) 337 ABI_MALLOC(gkk, (2*mband*nsppol,nkpt,1,1,mband_kq)) 338 339 else if (dtset%eph_task == -2) then 340 ! Write GKQ file with all perturbations. gkq are given in the atom representation. 341 ! TODO: mband_kq == mband 342 ABI_MALLOC(gkq_atm, (2, mband_kq, mband, nkpt)) 343 if (i_am_master) then 344 call ifc%fourq(cryst, qpt, phfrq, displ_cart, out_displ_red=displ_red) 345 fname = strcat(dtfil%filnam_ds(4), "_GKQ.nc") 346 #ifdef HAVE_NETCDF 347 NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating GKQ file") 348 NCF_CHECK(cryst%ncwrite(ncid)) 349 ! Write bands on k mesh. 350 NCF_CHECK(ebands_ncwrite(ebands_k, ncid)) 351 ncerr = nctk_def_dims(ncid, [nctkdim_t('number_of_phonon_modes', natom3)], defmode=.True.) 352 NCF_CHECK(ncerr) 353 ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: & 354 "symdynmat", "symv1scf", "dvdb_add_lr", "interpolated"]) 355 NCF_CHECK(ncerr) 356 NCF_CHECK(nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "qdamp"])) 357 358 ! Define EPH arrays 359 ncerr = nctk_def_arrays(ncid, [ & 360 nctkarr_t('qpoint', "dp" , 'number_of_reduced_dimensions'), & 361 nctkarr_t('emacro_cart', "dp", 'number_of_cartesian_directions, number_of_cartesian_directions'), & 362 nctkarr_t('becs_cart', "dp", "number_of_cartesian_directions, number_of_cartesian_directions, number_of_atoms"), & 363 nctkarr_t("eigenvalues_kq", "dp", "max_number_of_states, number_of_kpoints, number_of_spins"), & 364 nctkarr_t('phfreqs', "dp", 'number_of_phonon_modes'), & 365 nctkarr_t('phdispl_cart', "dp", 'complex, number_of_phonon_modes, number_of_phonon_modes'), & 366 !nctkarr_t('phdispl_cart_qvers', "dp", 'complex, number_of_phonon_modes, number_of_phonon_modes'), & 367 nctkarr_t('phdispl_red', "dp", 'complex, number_of_phonon_modes, number_of_phonon_modes'), & 368 nctkarr_t("gkq_representation", "char", "character_string_length"), & 369 nctkarr_t('gkq', "dp", & 370 'complex, max_number_of_states, max_number_of_states, number_of_phonon_modes, number_of_kpoints, number_of_spins') & 371 ]) 372 NCF_CHECK(ncerr) 373 ! Write data. 374 NCF_CHECK(nctk_set_datamode(ncid)) 375 ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: & 376 "symdynmat", "symv1scf", "dvdb_add_lr", "interpolated"], & 377 [dtset%symdynmat, dtset%symv1scf, dtset%dvdb_add_lr, interpolated]) 378 NCF_CHECK(ncerr) 379 dvdb_qdamp = dvdb%qdamp 380 NCF_CHECK(nctk_write_dpscalars(ncid, [character(len=nctk_slen) :: "qdamp"], dvdb_qdamp)) 381 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "qpoint"), qpt)) 382 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "emacro_cart"), dvdb%dielt)) 383 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "becs_cart"), dvdb%zeff)) 384 ABI_MALLOC(eigens_kq, (ebands_kq%mband, nkpt, nsppol)) 385 do ik=1,nkpt 386 kk = ebands_k%kptns(:,ik) 387 kq = kk + qpt 388 ! Find the index of the k+q point 389 call findqg0(ikq, g0_k, kq, nkpt_kq, ebands_kq%kptns, [1,1,1]) 390 eigens_kq(:, ik, :) = ebands_kq%eig(:, ikq, :) 391 end do 392 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "eigenvalues_kq"), eigens_kq)) 393 ABI_FREE(eigens_kq) 394 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "phfreqs"), phfrq)) 395 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'phdispl_cart'), displ_cart)) 396 ! Add phonon displacement for qvers 397 !call ifc_fourq(ifc, cryst, qpt, phfrq, displ_cart, out_displ_red=displ_red, nanaqdir="reduced") 398 !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'phdispl_cart_qvers'), displ_cart)) 399 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'phdispl_red'), displ_red)) 400 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gkq_representation"), "atom")) 401 #endif 402 end if 403 else 404 ABI_ERROR(sjoin("Invalid value for eph_task:", itoa(dtset%eph_task))) 405 end if 406 407 ! Loop over all 3*natom perturbations. 408 do ipc=1,natom3 409 idir = mod(ipc-1, 3) + 1 410 ipert = (ipc - idir) / 3 + 1 411 write(msg, '(a,2i4)') " Treating ipert, idir = ", ipert, idir 412 call wrtout(std_out, msg, do_flush=.True.) 413 if (dtset%eph_task == 2) gkk = zero 414 415 do spin=1,nsppol 416 if (dtset%eph_task == -2) gkq_atm = zero 417 418 ! Set up local potential vlocal1 with proper dimensioning, from vtrial1 taking into account the spin. 419 call rf_transgrid_and_pack(spin,nspden,psps%usepaw,cplex,nfftf,nfft,ngfft,gs_hamkq%nvloc,& 420 pawfgr,mpi_enreg,dummy_vtrial,v1scf(:,:,:,ipc),vlocal,vlocal1(:,:,:,:,ipc)) 421 422 ! Continue to initialize the Hamiltonian 423 call gs_hamkq%load_spin(spin,vlocal=vlocal,with_nonlocal=.true.) 424 425 ! Allocate workspace for wavefunctions. Make npw larger than expected. 426 ABI_MALLOC(bras, (2, mpw*nspinor, mband)) 427 ABI_MALLOC(kets, (2, mpw*nspinor, mband)) 428 ABI_MALLOC(h1_kets, (2, mpw*nspinor, mband)) 429 430 ! GKA: This little block used to be right after the perturbation loop 431 ! Prepare application of the NL part. 432 call init_rf_hamiltonian(cplex,gs_hamkq,ipert,rf_hamkq,has_e1kbsc=.true.) 433 call rf_hamkq%load_spin(spin,vlocal1=vlocal1(:,:,:,:,ipc),with_nonlocal=.true.) 434 435 do ik=1,nkpt 436 ! Only do a subset a k-points 437 if (.not. (ik >= my_kstart .and. ik <= my_kstop)) cycle 438 439 kk = ebands_k%kptns(:,ik) 440 kq = kk + qpt 441 ! Find the index of the k+q point 442 call findqg0(ikq, g0_k, kq, nkpt_kq, ebands_kq%kptns, [1,1,1]) 443 444 ! Copy u_k(G) 445 istwf_k = wfd_k%istwfk(ik); npw_k = wfd_k%npwarr(ik) 446 ABI_CHECK(mpw >= npw_k, "mpw < npw_k") 447 kg_k(:,1:npw_k) = wfd_k%kdata(ik)%kg_k 448 do ib2=1,mband 449 call wfd_k%copy_cg(ib2, ik, spin, kets(1,1,ib2)) 450 end do 451 452 ! Copy u_kq(G) 453 istwf_kq = wfd_kq%istwfk(ikq); npw_kq = wfd_kq%npwarr(ikq) 454 ABI_CHECK(mpw >= npw_kq, "mpw < npw_kq") 455 kg_kq(:,1:npw_kq) = wfd_kq%kdata(ikq)%kg_k 456 do ib1=1,mband_kq 457 call wfd_kq%copy_cg(ib1, ikq, spin, bras(1,1,ib1)) 458 end do 459 460 ! if PAW, one has to solve a generalized eigenproblem 461 ! Be careful here because I will need sij_opt==-1 462 gen_eigenpb = (psps%usepaw==1) 463 sij_opt = 0; if (gen_eigenpb) sij_opt = 1 464 ABI_MALLOC(gs1c, (2,npw_kq*nspinor*((sij_opt+1)/2))) 465 466 ! GKA: Previous loop on 3*natom perturbations used to start here 467 ! This call is not optimal because there are quantities in out that do not depend on idir,ipert 468 call getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,kk,kq,idir,ipert,& ! In 469 cryst%natom,cryst%rmet,cryst%gprimd,cryst%gmet,istwf_k,& ! In 470 npw_k,npw_kq,useylmgr1,kg_k,ylm_k,kg_kq,ylm_kq,ylmgr_kq,& ! In 471 dkinpw,nkpg,nkpg1,kpg_k,kpg1_k,kinpw1,ffnlk,ffnl1,ph3d,ph3d1) ! Out 472 473 ! Calculate dvscf * psi_k, results stored in h1_kets on the k+q sphere. 474 ! Compute H(1) applied to GS wavefunction Psi(0) 475 do ib2=1,mband 476 eig0nk = ebands_k%eig(ib2,ik,spin) 477 ! Use scissor shift on 0-order eigenvalue 478 eshift = eig0nk - dtset%dfpt_sciss 479 480 call getgh1c(berryopt0,kets(:,:,ib2),cwaveprj0,h1_kets(:,:,ib2),& 481 grad_berry,gs1c,gs_hamkq,gvnlx1,idir,ipert,eshift,mpi_enreg,optlocal,& 482 optnl,opt_gvnlx1,rf_hamkq,sij_opt,tim_getgh1c,usevnl) 483 end do 484 485 ABI_FREE(kinpw1) 486 ABI_FREE(kpg1_k) 487 ABI_FREE(kpg_k) 488 ABI_FREE(dkinpw) 489 ABI_FREE(ffnlk) 490 ABI_FREE(ffnl1) 491 ABI_FREE(ph3d) 492 ABI_FREE(gs1c) 493 ABI_SFREE(ph3d1) 494 495 ! Calculate elphmat(j,i) = <psi_{k+q,j}|dvscf_q*psi_{k,i}> for this perturbation. 496 ! The array eig1_k contains: 497 ! 498 ! <u_(band,k+q)^(0)|H_(k+q,k)^(1)|u_(band,k)^(0)> (NC psps) 499 ! <u_(band,k+q)^(0)|H_(k+q,k)^(1)-(eig0_k+eig0_k+q)/2.S^(1)|u_(band,k)^(0)> (PAW) 500 do ib2=1,mband 501 do ib1=1,mband_kq 502 call dotprod_g(dotr,doti,istwf_kq,npw_kq*nspinor,2,bras(1,1,ib1),h1_kets(1,1,ib2),& 503 mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 504 band = 2*ib2-1 + (spin-1) * 2 * mband 505 506 if (dtset%eph_task == 2) then 507 gkk(band,ik,1,1,ib1) = dotr 508 gkk(band+1,ik,1,1,ib1) = doti 509 else 510 gkq_atm(:, ib1, ib2, ik) = [dotr, doti] 511 end if 512 513 end do 514 end do 515 516 end do ! ikpt 517 518 ABI_FREE(bras) 519 ABI_FREE(kets) 520 ABI_FREE(h1_kets) 521 call rf_hamkq%free() 522 523 if (dtset%eph_task == -2) then 524 ! Gather the k-points computed by all processes 525 call xmpi_sum_master(gkq_atm, master, comm, ierr) 526 if (i_am_master) then 527 ! Write the netCDF file. 528 #ifdef HAVE_NETCDF 529 ncerr = nf90_put_var(ncid, nctk_idname(ncid, "gkq"), gkq_atm, & 530 start=[1, 1, 1, ipc, 1, 1], count=[2, mband, mband, 1, nkpt, spin]) 531 NCF_CHECK(ncerr) 532 #endif 533 end if 534 end if 535 end do ! spin 536 537 if (dtset%eph_task == 2) then 538 ! Gather the k-points computed by all processes 539 call xmpi_sum_master(gkk,master,comm,ierr) 540 ! Init a gkk_t object 541 call gkk_init(gkk,gkk2d,mband,nsppol,nkpt,1,1) 542 ! Write the netCDF file. 543 call appdig(ipc,dtfil%fnameabo_gkk,gkkfilnam) 544 fname = strcat(gkkfilnam, ".nc") 545 #ifdef HAVE_NETCDF 546 if (i_am_master) then 547 NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating GKK file") 548 NCF_CHECK(cryst%ncwrite(ncid)) 549 NCF_CHECK(ebands_ncwrite(ebands_k, ncid)) 550 call gkk_ncwrite(gkk2d, qpt, 1.0_dp, ncid) 551 NCF_CHECK(nf90_close(ncid)) 552 end if 553 #endif 554 ! Free memory 555 call gkk_free(gkk2d) 556 end if 557 end do ! ipc (loop over 3*natom atomic perturbations) 558 559 call cwtime(cpu, wall, gflops, "stop") 560 write(msg, '(2a)') " Computation of gkq matrix elements with ", trim(what) 561 call wrtout([std_out, ab_out], msg, do_flush=.True.) 562 call wrtout(std_out, sjoin("cpu-time:", sec2str(cpu), ",wall-time:", sec2str(wall)), do_flush=.True.) 563 564 if (dtset%eph_task == -2 .and. i_am_master) then 565 #ifdef HAVE_NETCDF 566 NCF_CHECK(nf90_close(ncid)) 567 #endif 568 end if 569 570 ! =========== 571 ! Free memory 572 ! =========== 573 ABI_SFREE(gkk) 574 ABI_SFREE(gkq_atm) 575 ABI_FREE(displ_cart) 576 ABI_FREE(displ_red) 577 ABI_FREE(v1scf) 578 ABI_FREE(vlocal1) 579 ABI_FREE(gvnlx1) 580 ABI_FREE(grad_berry) 581 ABI_FREE(dummy_vtrial) 582 ABI_FREE(ph1d) 583 ABI_FREE(vlocal) 584 ABI_FREE(kg_k) 585 ABI_FREE(kg_kq) 586 ABI_FREE(ylm_k) 587 ABI_FREE(ylm_kq) 588 ABI_FREE(ylmgr_kq) 589 ABI_FREE(blkflg) 590 591 call gs_hamkq%free() 592 call wfd_k%free() 593 call wfd_kq%free() 594 call pawcprj_free(cwaveprj0) 595 ABI_FREE(cwaveprj0) 596 597 end subroutine eph_gkk
m_gkk/find_mpw [ Functions ]
[ Top ] [ m_gkk ] [ Functions ]
NAME
find_mpw
FUNCTION
Look at all k-points and spins to find the maximum number of plane waves.
INPUTS
OUTPUT
NOTES
SOURCE
910 subroutine find_mpw(mpw, kpts, nsppol, nkpt, gmet, ecut, comm) 911 912 !Arguments ------------------------------------ 913 !scalars 914 integer,intent(out) :: mpw 915 integer,intent(in) :: nsppol, nkpt 916 integer,intent(in) :: comm 917 real(dp),intent(in) :: ecut 918 !arrays 919 real(dp),intent(in) :: kpts(3,nkpt) 920 real(dp),intent(in) :: gmet(3,3) 921 922 !Local variables ------------------------------ 923 !scalars 924 integer :: my_rank, cnt, nproc, ierr 925 integer :: ispin, ikpt 926 integer :: my_mpw, onpw 927 integer,allocatable :: gtmp(:,:) 928 real(dp) :: kpt(3) 929 930 my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm) 931 932 mpw = 0; cnt=0 933 do ispin=1,nsppol 934 do ikpt=1,nkpt 935 cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle 936 kpt = kpts(:,ikpt) 937 call get_kg(kpt,1,ecut,gmet,onpw,gtmp) 938 ABI_FREE(gtmp) 939 mpw = max(mpw, onpw) 940 end do 941 end do 942 my_mpw = mpw; call xmpi_max(my_mpw, mpw, comm, ierr) 943 944 end subroutine find_mpw
m_gkk/ncwrite_v1qnu [ Functions ]
[ Top ] [ m_gkk ] [ Functions ]
NAME
ncwrite_v1qnu
FUNCTION
Compute \delta V_{q,nu)(r) and dump results to netcdf file.
This routine should be called by a single processor.
INPUT
dvdb<dbdb_type>=Database with the DFPT SCF potentials.
dtset<dataset_type>= Input variables.
ifc<ifc_type>=interatomic force constants and corresponding real space grid info.
out_ncpath=Name of the netcdf file.
OUTPUT
Only writing
SOURCE
621 subroutine ncwrite_v1qnu(dvdb, dtset, ifc, out_ncpath) 622 623 use m_bz_mesh, only : kpath_t, kpath_new 624 625 !Arguments ------------------------------------ 626 !scalars 627 type(dataset_type),target,intent(in) :: dtset 628 type(dvdb_t),intent(inout) :: dvdb 629 type(ifc_type),intent(in) :: ifc 630 character(len=*),intent(in) :: out_ncpath 631 632 !Local variables------------------------------- 633 !scalars 634 integer,parameter :: master = 0, qptopt1 = 1 635 integer :: db_iqpt, cplex, nfft, comm, ip, idir, ipert, my_rank, interpolated, comm_rpt 636 logical :: with_lr_model 637 #ifdef HAVE_NETCDF 638 integer :: ncid, ncerr 639 #endif 640 !arrays 641 integer :: ngfft(18) 642 real(dp) :: phfreqs(dvdb%natom3),qpt(3) 643 real(dp) :: displ_cart(2,3, dvdb%cryst%natom, dvdb%natom3), displ_red(2,dvdb%natom3,dvdb%natom3) 644 real(dp),allocatable :: v1scf(:,:,:,:), v1_qnu(:,:,:,:) 645 real(dp),allocatable :: v1lr_atm(:,:,:,:), v1lr_qnu(:,:,:,:) 646 #if 1 647 integer :: iq, nu, iatom, ii, jj, kk 648 real(dp) :: inv_qepsq, qtau, phre, phim, rtmp 649 type(kpath_t) :: qpath 650 real(dp) :: bounds(3,6), qpt_red(3), qpt_cart(3), glr(3) 651 real(dp) :: values(dvdb%natom3) 652 653 !************************************************************************ 654 655 ! +0.50000 +0.50000 +0.50000 # L 656 ! +0.00000 +0.00000 +0.00000 # $\Gamma$ 657 ! +0.50000 +0.00000 +0.50000 # X 658 ! +0.50000 +0.25000 +0.75000 # W 659 ! +0.37500 +0.37500 +0.75000 # K 660 ! +0.00000 +0.00000 +0.00000 # $\Gamma$ 661 ! +0.37500 +0.37500 +0.75000 # K 662 663 ! +0.62500 +0.25000 +0.62500 # U 664 ! +0.50000 +0.50000 +0.50000 # L 665 ! +0.37500 +0.37500 +0.75000 # K 666 ! +0.62500 +0.25000 +0.62500 # U 667 ! +0.50000 +0.00000 +0.50000 # X 668 669 bounds(:, 1) = tol3 * [+0.50000, +0.50000, +0.50000] ! # L 670 bounds(:, 2) = tol3 * [+0.00000, +0.00000, +0.00000] ! # $\Gamma$ 671 bounds(:, 3) = tol3 * [+0.50000, +0.00000, +0.50000] ! # X 672 bounds(:, 4) = tol3 * [+0.37500, +0.37500, +0.75000] ! # K 673 bounds(:, 5) = tol3 * [+0.00000, +0.00000, +0.00000] ! # $\Gamma$ 674 bounds(:, 6) = tol3 * [+0.50000, +0.25000, +0.75000] ! # W 675 676 qpath = kpath_new(bounds, dvdb%cryst%gprimd, dtset%ndivsm) 677 678 do iq=1,qpath%npts 679 qpt_red = qpath%points(:, iq) 680 qpt_cart = two_pi * matmul(dvdb%cryst%gprimd, qpt_red) 681 inv_qepsq = one / dot_product(qpt_cart, matmul(ifc%dielt, qpt_cart)) 682 call ifc%fourq(dvdb%cryst, qpt_red, phfreqs, displ_cart) 683 do nu=1, dvdb%natom3 684 glr = zero 685 do iatom=1, dvdb%cryst%natom 686 ! Phase factor exp(-i (q+G) . tau) 687 qtau = - two_pi * dot_product(qpt_red, dvdb%cryst%xred(:,iatom)) 688 phre = cos(qtau); phim = sin(qtau) 689 do jj=1,3 690 do ii=1,3 691 do kk=1,3 692 rtmp = dvdb%qstar(ii, jj, kk, iatom) * qpt_cart(ii) * qpt_cart(jj) 693 glr(1) = glr(1) + rtmp * (displ_cart(1, kk, iatom, nu) * phre - displ_cart(2, kk, iatom, nu) * phim) 694 glr(2) = glr(2) + rtmp * (displ_cart(2, kk, iatom, nu) * phre + displ_cart(1, kk, iatom, nu) * phre) 695 end do 696 end do 697 end do 698 end do 699 glr = half * (glr / inv_qepsq) * (four_pi / dvdb%cryst%ucvol) 700 values(nu) = (glr(1) ** 2 + glr(2) ** 2) / (two * phfreqs(nu)) 701 end do ! nu 702 write(std_out, "(i0, 4(f9.6), /, (es18.6, 1x))") iq, qpt_red, phfreqs(nu), (values(nu), nu=1, 3*dvdb%natom) 703 end do ! iqpt 704 705 call qpath%free() 706 return 707 #endif 708 709 my_rank = xmpi_comm_rank(dvdb%comm) 710 comm = dvdb%comm 711 qpt = dtset%qptn 712 713 call wrtout(std_out, sjoin(" Writing Delta V_{q,nu)(r) potentials to file:", out_ncpath), do_flush=.True.) 714 call wrtout([std_out, ab_out], sjoin(ch10, "- Results stored in: ", out_ncpath)) 715 call wrtout(std_out, sjoin(" Using qpt:", ktoa(qpt))) 716 !call wrtout([std_out, ab_out], " Use `abiopen.py out_V1QAVG.nc -e` to visualize results") 717 call dvdb%print(unit=std_out) 718 719 ! Define FFT mesh 720 ngfft = dvdb%ngfft 721 nfft = product(ngfft(1:3)) 722 723 if (dtset%eph_task == -16) then 724 call wrtout([std_out, ab_out], " Assuming q-point already in the DVDB file. No interpolation.") 725 interpolated = 0 726 727 else if (dtset%eph_task == +16) then 728 call wrtout([std_out, ab_out], " Using Fourier interpolation.") 729 comm_rpt = xmpi_comm_self 730 call dvdb%ftinterp_setup(dtset%ddb_ngqpt, qptopt1, 1, dtset%ddb_shiftq, nfft, ngfft, comm_rpt) 731 interpolated = 1 732 else 733 ABI_ERROR(sjoin("Invalid value for eph_task:", itoa(dtset%eph_task))) 734 end if 735 736 with_lr_model = .True. 737 738 ! Create netcdf file. 739 #ifdef HAVE_NETCDF 740 if (my_rank == master) then 741 NCF_CHECK(nctk_open_create(ncid, out_ncpath, comm)) 742 NCF_CHECK(dvdb%cryst%ncwrite(ncid)) 743 744 ! Add other dimensions. 745 ncerr = nctk_def_dims(ncid, [ & 746 nctkdim_t("nfft", nfft), nctkdim_t("nspden", dvdb%nspden), & 747 nctkdim_t("natom3", 3 * dvdb%cryst%natom)], defmode=.True.) 748 NCF_CHECK(ncerr) 749 750 ! Define arrays 751 ncerr = nctk_def_arrays(ncid, [ & 752 nctkarr_t("ngfft", "int", "three"), & 753 nctkarr_t("qpt", "dp", "three"), & 754 nctkarr_t("phfreqs", "dp", "natom3"), & 755 nctkarr_t("displ_cart", "dp", "two, natom3, natom3"), & 756 nctkarr_t("v1_qnu", "dp", "two, nfft, nspden, natom3")]) 757 NCF_CHECK(ncerr) 758 759 if (with_lr_model) then 760 NCF_CHECK(nctk_def_arrays(ncid, [nctkarr_t("v1lr_qnu", "dp", "two, nfft, nspden, natom3")])) 761 end if 762 763 NCF_CHECK(nctk_set_datamode(ncid)) 764 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ngfft"), ngfft(1:3))) 765 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "qpt"), qpt)) 766 end if 767 #endif 768 769 ABI_MALLOC(v1_qnu, (2, nfft, dvdb%nspden, dvdb%natom3)) 770 if (with_lr_model) then 771 ABI_MALLOC(v1lr_atm, (2, nfft, dvdb%nspden, dvdb%natom3)) 772 ABI_MALLOC(v1lr_qnu, (2, nfft, dvdb%nspden, dvdb%natom3)) 773 end if 774 775 ! Get phonon freqs and displacemented for this q-point. 776 call ifc%fourq(dvdb%cryst, qpt, phfreqs, displ_cart, out_displ_red=displ_red) 777 778 if (interpolated == 0) then 779 ! Find the index of the q-point in the DVDB. 780 db_iqpt = dvdb%findq(qpt) 781 if (db_iqpt /= -1) then 782 ! Read or reconstruct the dvscf potentials for all 3*natom perturbations. 783 ! This call allocates v1scf(cplex, nfft, nspden, 3*natom)) 784 call dvdb%readsym_allv1(db_iqpt, cplex, nfft, ngfft, v1scf, comm) 785 else 786 ABI_ERROR(sjoin("Cannot find q-point:", ktoa(qpt), "in DVDB file")) 787 end if 788 else 789 790 cplex = 2 791 ABI_MALLOC(v1scf, (cplex, nfft, dvdb%nspden, dvdb%my_npert)) 792 call dvdb%ftinterp_qpt(qpt, nfft, ngfft, v1scf, dvdb%comm_rpt) 793 end if 794 795 ! Compute scattering potential in phonon representations instead ot atomic one. 796 ! v1_qnu = \sum_{ka} phdispl{ka}(q,nu) D_{ka,q} V_scf(r) 797 ! NOTE: prefactor 1/sqrt(2 w(q,nu)) is not included in the potentials saved to file. 798 ! v1_qnu(2, nfft, nspden, natom3), v1scf(cplex, nfft, nspden, natom3) 799 call v1atm_to_vqnu(cplex, nfft, dvdb%nspden, dvdb%natom3, v1scf, displ_red, v1_qnu) 800 801 if (with_lr_model) then 802 ! Compute LR model in the atomic representation then compute phonon representation in v1lr_qnu. 803 v1lr_atm = zero 804 do idir=1,3 805 do ipert=1,dvdb%natom 806 ip = (ipert - 1) * 3 + idir 807 call dvdb%get_v1r_long_range(qpt, idir, ipert, nfft, ngfft, v1lr_atm(:,:,1,ip)) 808 if (dvdb%nspden == 2) v1lr_atm(:,:,2,ip) = v1lr_atm(:,:,1,ip) 809 end do 810 end do 811 call v1atm_to_vqnu(2, nfft, dvdb%nspden, dvdb%natom3, v1lr_atm, displ_red, v1lr_qnu) 812 end if 813 814 #ifdef HAVE_NETCDF 815 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "phfreqs"), phfreqs)) 816 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "displ_cart"), displ_cart)) 817 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "v1_qnu"), v1_qnu)) 818 if (with_lr_model) then 819 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "v1lr_qnu"), v1lr_qnu)) 820 end if 821 #endif 822 823 ABI_FREE(v1scf) 824 ABI_FREE(v1_qnu) 825 ABI_SFREE(v1lr_atm) 826 ABI_SFREE(v1lr_qnu) 827 828 #ifdef HAVE_NETCDF 829 NCF_CHECK(nf90_close(ncid)) 830 #endif 831 call dvdb%close() 832 833 call wrtout(std_out, "dvqnu file written", do_flush=.True.) 834 835 end subroutine ncwrite_v1qnu
m_gkk/v1atm_to_vqnu [ Functions ]
[ Top ] [ m_gkk ] [ Functions ]
NAME
v1atm_to_vqnu
FUNCTION
Receive potentials in atomic representation and return potential in phonon representation
INPUTS
OUTPUT
SOURCE
853 subroutine v1atm_to_vqnu(cplex, nfft, nspden, natom3, v1_atm, displ_red, v1_qnu) 854 855 !Arguments ------------------------------------ 856 !scalars 857 integer,intent(in) :: cplex, nfft, nspden, natom3 858 !arrays 859 real(dp),intent(in) :: v1_atm(cplex, nfft, nspden, natom3) 860 real(dp),intent(out) :: v1_qnu(2, nfft, nspden, natom3) 861 real(dp),intent(in) :: displ_red(2, natom3, natom3) 862 863 !Local variables------------------------------- 864 !scalars 865 integer :: nu, ip, ispden 866 867 !************************************************************************ 868 869 do nu=1,natom3 870 ! v1_qnu = \sum_{ka} phdispl{ka}(q,nu) D_{ka,q} V_scf(r) 871 ! NOTE: prefactor 1/sqrt(2 w(q,nu)) is not included in the potentials saved to file. 872 ! v1_qnu(2, nfft, nspden, natom3), v1_atm(cplex, nfft, nspden, natom3) 873 v1_qnu(:, :, :, nu) = zero 874 do ip=1,natom3 875 do ispden=1,nspden 876 if (cplex == 2) then 877 v1_qnu(1, :, ispden, nu) = v1_qnu(1, :, ispden, nu) + & 878 displ_red(1,ip,nu) * v1_atm(1,:,ispden,ip) - displ_red(2,ip,nu) * v1_atm(2,:,ispden,ip) 879 v1_qnu(2, :, ispden, nu) = v1_qnu(2, :, ispden, nu) + & 880 displ_red(2,ip,nu) * v1_atm(1,:,ispden,ip) + displ_red(1,ip,nu) * v1_atm(2,:,ispden,ip) 881 else 882 ! Gamma point. d(q) = d(-q)* --> d is real. 883 v1_qnu(1, :, ispden, nu) = v1_qnu(1, :, ispden, nu) + displ_red(1,ip,nu) * v1_atm(1,:,ispden,ip) 884 end if 885 end do 886 end do 887 end do 888 889 end subroutine v1atm_to_vqnu