TABLE OF CONTENTS
ABINIT/getkgrid [ Functions ]
NAME
getkgrid
FUNCTION
Compute the grid of k points in the irreducible Brillouin zone. Note that nkpt (and nkpthf) can be computed by calling this routine with nkpt=0, provided that kptopt/=0. If downsampling is present, also compute a downsampled k grid.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, MM) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
chksymbreak= if 1, will check whether the k point grid is symmetric (for kptopt=1,2 and 4), and stop if not. iout=unit number for echoed output . 0 if no output is wished. iscf= ( <= 0 =>non-SCF), >0 => SCF) MG: FIXME I don't understand why we have to pass the value iscf. kptopt=option for the generation of k points (defines whether spatial symmetries and/or time-reversal can be used) msym=default maximal number of symmetries nsym=number of symmetries rprimd(3,3)=dimensional real space primitive translations (bohr) symafm(nsym)=(anti)ferromagnetic part of symmetry operations symrel(3,3,nsym)=symmetry operations in real space in terms of primitive translations vacuum(3)=for each direction, 0 if no vacuum, 1 if vacuum [downsampling(3) = input variable that governs the downsampling]
OUTPUT
kptrlen=length of the smallest real space supercell vector associated with the lattice of k points. nkpt_computed=number of k-points in the IBZ computed in the present routine If nkpt/=0 the following are also output : kpt(3,nkpt)=reduced coordinates of k points. wtk(nkpt)=weight assigned to each k point. [fullbz(3,nkpt_fullbz)]=k-points generated in the full Brillouin zone. In output: allocated array with the list of k-points in the BZ. [kpthf(3,nkpthf)]=k-points generated in the full Brillouin zone, possibly downsampled (for Fock).
NOTES
msym not needed since nsym is the last index.
SIDE EFFECTS
Input/Output nkpt=number of k points (might be zero, see output description) kptrlatt(3,3)=k-point lattice specification nshiftk=actual number of k-point shifts in shiftk shiftk(3,210)=shift vectors for k point generation [nkpthf = number of k points in the full BZ, for the Fock operator]
PARENTS
ep_setupqpt,getshell,inkpts,inqpt,m_ab7_kpoints,m_bz_mesh,m_kpts nonlinear,testkgrid,thmeig
CHILDREN
mati3inv,matr3inv,metric,smallprim,smpbz,symkpt
SOURCE
62 #if defined HAVE_CONFIG_H 63 #include "config.h" 64 #endif 65 66 #include "abi_common.h" 67 68 subroutine getkgrid(chksymbreak,iout,iscf,kpt,kptopt,kptrlatt,kptrlen,& 69 & msym,nkpt,nkpt_computed,nshiftk,nsym,rprimd,shiftk,symafm,symrel,vacuum,wtk,& 70 & fullbz,nkpthf,kpthf,downsampling) ! optional 71 72 use defs_basis 73 use m_profiling_abi 74 use m_errors 75 76 use m_geometry, only : metric 77 78 !This section has been created automatically by the script Abilint (TD). 79 !Do not modify the following lines by hand. 80 #undef ABI_FUNC 81 #define ABI_FUNC 'getkgrid' 82 use interfaces_32_util 83 use interfaces_41_geometry 84 use interfaces_56_recipspace, except_this_one => getkgrid 85 !End of the abilint section 86 87 implicit none 88 89 !Arguments ------------------------------------ 90 !scalars 91 integer,intent(in) :: chksymbreak,iout,iscf,kptopt,msym,nkpt,nsym 92 integer,intent(inout),optional :: nkpthf 93 integer,intent(inout) :: nshiftk 94 integer,intent(inout) :: nkpt_computed !vz_i 95 real(dp),intent(out) :: kptrlen 96 !arrays 97 integer,intent(in) :: symafm(msym),symrel(3,3,msym),vacuum(3) 98 integer,optional,intent(in) :: downsampling(3) 99 integer,intent(inout) :: kptrlatt(3,3) 100 real(dp),intent(in) :: rprimd(3,3) 101 real(dp),intent(inout) :: shiftk(3,210) 102 real(dp),intent(inout) :: kpt(3,nkpt) !vz_i 103 real(dp),intent(inout) :: wtk(nkpt) 104 real(dp),optional,allocatable,intent(out) :: fullbz(:,:) 105 real(dp),optional,intent(out) :: kpthf(:,:) 106 107 !Local variables------------------------------- 108 !scalars 109 integer, parameter :: max_number_of_prime=47,mshiftk=210 110 integer :: brav,decreased,found,ii,ikpt,iprime,ishiftk,isym,jshiftk,kshiftk,mkpt,mult 111 integer :: nkpthf_computed,nkpt_fullbz,nkptlatt,nshiftk2,nsym_used,option 112 integer :: test_prime,timrev 113 real(dp) :: length2,ucvol,ucvol_super 114 character(len=500) :: message 115 !arrays 116 integer, parameter :: prime_factor(max_number_of_prime)=(/2,3,5,7,9, 11,13,17,19,23,& 117 & 29,31,37,41,43, 47,53,59,61,67,& 118 & 71,73,79,83,89, 97,101,103,107,109,& 119 & 113,127,131,137,139, 149,151,157,163,167,& 120 & 173,179,181,191,193, 197,199/) 121 integer :: kptrlatt2(3,3) 122 integer,allocatable :: belong_chain(:),generator(:),indkpt(:),number_in_chain(:) 123 integer,allocatable :: repetition_factor(:),symrec(:,:,:) 124 ! real(dp) :: cart(3,3) 125 real(dp) :: dijk(3),delta_dmult(3),dmult(3),fact_vacuum(3),gmet(3,3) 126 real(dp) :: gmet_super(3,3),gprimd(3,3),gprimd_super(3,3),klatt2(3,3) 127 real(dp) :: klatt3(3,3),kptrlattr(3,3),ktransf(3,3),ktransf_invt(3,3) 128 real(dp) :: metmin(3,3),minim(3,3),rmet(3,3),rmet_super(3,3),rprimd_super(3,3) 129 real(dp),allocatable :: deltak(:,:),kpt_fullbz(:,:),shiftk2(:,:),shiftk3(:,:),spkpt(:,:),wtk_folded(:),wtk_fullbz(:) 130 131 ! ************************************************************************* 132 133 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 134 135 if (kptopt==1.or.kptopt==4) then 136 ! Cannot use antiferromagnetic symmetry operations to decrease the number of k points 137 nsym_used=0 138 do isym=1,nsym 139 if(symafm(isym)==1)nsym_used=nsym_used+1 140 end do 141 ABI_ALLOCATE(symrec,(3,3,nsym_used)) 142 nsym_used=0 143 do isym=1,nsym ! Get the symmetry matrices in terms of reciprocal basis 144 if(symafm(isym)==1)then 145 nsym_used=nsym_used+1 146 call mati3inv(symrel(:,:,isym),symrec(:,:,nsym_used)) 147 end if 148 end do 149 else if (kptopt==2) then 150 !Use only the time-reversal 151 nsym_used=1 152 ABI_ALLOCATE(symrec,(3,3,1)) 153 symrec(1:3,1:3,1)=0 154 do ii=1,3 155 symrec(ii,ii,1)=1 156 end do 157 end if 158 159 kptrlatt2(:,:)=kptrlatt(:,:) 160 nshiftk2=nshiftk 161 ABI_ALLOCATE(shiftk2,(3,mshiftk)) 162 ABI_ALLOCATE(shiftk3,(3,mshiftk)) 163 shiftk2(:,:)=shiftk(:,:) 164 165 !Find a primitive k point lattice, if possible, by decreasing the number of shifts. 166 if(nshiftk2/=1)then 167 168 do ! Loop to be repeated if there has been a successful reduction of nshiftk2 169 170 ABI_ALLOCATE(deltak,(3,nshiftk2)) 171 ABI_ALLOCATE(repetition_factor,(nshiftk2)) 172 ABI_ALLOCATE(generator,(nshiftk2)) 173 ABI_ALLOCATE(belong_chain,(nshiftk2)) 174 ABI_ALLOCATE(number_in_chain,(nshiftk2)) 175 176 decreased=0 177 deltak(1,1:nshiftk2)=shiftk2(1,1:nshiftk2)-shiftk2(1,1) 178 deltak(2,1:nshiftk2)=shiftk2(2,1:nshiftk2)-shiftk2(2,1) 179 deltak(3,1:nshiftk2)=shiftk2(3,1:nshiftk2)-shiftk2(3,1) 180 deltak(:,:)=deltak(:,:)-floor(deltak(:,:)+tol8) 181 182 ! Identify for each shift, the smallest repetition prime factor that yields a reciprocal lattice vector. 183 repetition_factor(:)=0 184 repetition_factor(1)=1 185 do ishiftk=2,nshiftk2 186 do iprime=1,max_number_of_prime 187 test_prime=prime_factor(iprime) 188 dmult(:)=test_prime*deltak(:,ishiftk) 189 if(sum(abs( dmult(:)-nint(dmult(:)) ))<tol8)then 190 repetition_factor(ishiftk)=test_prime 191 exit 192 end if 193 end do 194 end do 195 196 ! Initialize the selection of tentative generators 197 generator(:)=1 198 do ishiftk=1,nshiftk2 199 if(repetition_factor(ishiftk)==0 .or. repetition_factor(ishiftk)==1)generator(ishiftk)=0 200 end do 201 202 ! Try different shifts as generators, by order of increasing repetition factor, 203 ! provided they are equal or bigger than 2 204 do iprime=1,max_number_of_prime 205 do ishiftk=2,nshiftk2 206 ! Note that ishiftk=1 is never a generator. It is the reference starting point. 207 if(generator(ishiftk)==1 .and. repetition_factor(ishiftk)==prime_factor(iprime))then 208 ! Test the generator : is it indeed closed ? 209 if(prime_factor(iprime)/=2)then 210 do mult=2,prime_factor(iprime)-1 211 dmult(:)=mult*deltak(:,ishiftk) 212 found=0 213 do jshiftk=1,nshiftk2 214 delta_dmult(:)=deltak(:,jshiftk)-dmult(:) 215 if(sum(abs(delta_dmult(:)-nint(delta_dmult(:)) ))<tol8)then 216 found=1 217 exit 218 end if 219 end do 220 if(found==0)exit 221 end do 222 if(found==0)generator(ishiftk)=0 223 end if 224 if(generator(ishiftk)==0)cycle 225 else 226 cycle 227 end if 228 ! Now, test whether all k points can be found in all possible chains 229 belong_chain(:)=0 230 do jshiftk=1,nshiftk2 231 ! Initialize a chain starting from a k point not yet in a chain 232 if(belong_chain(jshiftk)==0)then 233 number_in_chain(:)=0 ! Not a member of the chain (yet) 234 number_in_chain(jshiftk)=1 ! The first point in chain 235 do mult=1,prime_factor(iprime)-1 236 dmult(:)=mult*deltak(:,ishiftk) 237 found=0 238 do kshiftk=jshiftk+1,nshiftk2 239 delta_dmult(:)=deltak(:,kshiftk)-deltak(:,jshiftk)-dmult(:) 240 if(sum(abs(delta_dmult(:)-nint(delta_dmult(:)) ))<tol8)then 241 found=1 242 number_in_chain(kshiftk)=mult+1 243 exit 244 end if 245 end do 246 if(found==0)then 247 generator(ishiftk)=0 248 exit 249 end if 250 end do 251 if(generator(ishiftk)==1)then 252 ! Store the chain 253 do kshiftk=1,nshiftk2 254 if(number_in_chain(kshiftk)/=0)belong_chain(kshiftk)=number_in_chain(kshiftk) 255 end do 256 else 257 exit 258 end if 259 end if 260 end do 261 262 if(generator(ishiftk)==0)cycle 263 264 ! For the generator based on ishiftk, all the k points have been found to belong to one chain. 265 ! All the initializing k points in the different chains have belong_chain(:)=1 . 266 ! They must be kept, and the others thrown away. 267 ktransf(:,:)=0.0_dp 268 ktransf(1,1)=1.0_dp 269 ktransf(2,2)=1.0_dp 270 ktransf(3,3)=1.0_dp 271 ! Replace one of the unit vectors by the shift vector deltak(:,ishiftk). 272 ! However, must pay attention not to make linear combinations. 273 ! Also, choose positive sign for first-non-zero value. 274 if(abs(deltak(1,ishiftk)-nint(deltak(1,ishiftk)))>tol8)then 275 if(deltak(1,ishiftk)>0)ktransf(:,1)= deltak(:,ishiftk) 276 if(deltak(1,ishiftk)<0)ktransf(:,1)=-deltak(:,ishiftk) 277 else if(abs(deltak(2,ishiftk)-nint(deltak(2,ishiftk)))>tol8)then 278 if(deltak(2,ishiftk)>0)ktransf(:,2)= deltak(:,ishiftk) 279 if(deltak(2,ishiftk)<0)ktransf(:,2)=-deltak(:,ishiftk) 280 else if(abs(deltak(3,ishiftk)-nint(deltak(3,ishiftk)))>tol8)then 281 if(deltak(3,ishiftk)>0)ktransf(:,3)= deltak(:,ishiftk) 282 if(deltak(3,ishiftk)<0)ktransf(:,3)=-deltak(:,ishiftk) 283 end if 284 ! Copy the integers to real(dp) 285 kptrlattr(:,:)=kptrlatt2(:,:) 286 ! Go to reciprocal space 287 call matr3inv(kptrlattr,klatt2) 288 ! Make the transformation 289 do ii=1,3 290 klatt3(:,ii)=ktransf(1,ii)*klatt2(:,1)+ktransf(2,ii)*klatt2(:,2)+ktransf(3,ii)*klatt2(:,3) 291 end do 292 ! Back to real space 293 call matr3inv(klatt3,kptrlattr) 294 ! real(dp) to integer 295 kptrlatt2(:,:)=nint(kptrlattr(:,:)) 296 ! Prepare the transformation of the shifts 297 call matr3inv(ktransf,ktransf_invt) 298 decreased=1 299 kshiftk=0 300 do jshiftk=1,nshiftk2 301 if(belong_chain(jshiftk)==1)then 302 kshiftk=kshiftk+1 303 ! Place the shift with index jshiftk in place of the one in kshiftk, 304 ! also transform the shift from the old to the new coordinate system 305 shiftk3(:,kshiftk)=ktransf_invt(1,:)*shiftk2(1,jshiftk)+& 306 & ktransf_invt(2,:)*shiftk2(2,jshiftk)+& 307 & ktransf_invt(3,:)*shiftk2(3,jshiftk) 308 end if 309 end do 310 nshiftk2=nshiftk2/prime_factor(iprime) 311 shiftk2(:,1:nshiftk2)=shiftk3(:,1:nshiftk2)-floor(shiftk3(:,1:nshiftk2)+tol8) 312 if(kshiftk/=nshiftk2)then 313 MSG_BUG('The search for a primitive k point lattice contains a bug.') 314 end if 315 316 ! If this trial shift was successful, must exit the loop on trial ishiftk, 317 ! and reinitialize the global loop 318 if(decreased==1)exit 319 end do ! ishiftk 320 if(decreased==1)exit 321 end do ! iprime 322 323 ABI_DEALLOCATE(belong_chain) 324 ABI_DEALLOCATE(deltak) 325 ABI_DEALLOCATE(number_in_chain) 326 ABI_DEALLOCATE(repetition_factor) 327 ABI_DEALLOCATE(generator) 328 329 if(decreased==0 .or. nshiftk2==1)exit 330 331 end do ! Infinite loop 332 333 end if ! End nshiftk being 1 or larger 334 335 !Impose shiftk coordinates to be in [0,1[ 336 do ishiftk=1,nshiftk2 337 do ii=1,3 338 if(shiftk2(ii,ishiftk)>one-tol8) shiftk2(ii,ishiftk)=shiftk2(ii,ishiftk)-1.0_dp 339 if(shiftk2(ii,ishiftk)<-tol8) shiftk2(ii,ishiftk)=shiftk2(ii,ishiftk)+1.0_dp 340 end do 341 end do 342 343 !Compute the number of k points in the G-space unit cell 344 nkptlatt=kptrlatt2(1,1)*kptrlatt2(2,2)*kptrlatt2(3,3) & 345 & +kptrlatt2(1,2)*kptrlatt2(2,3)*kptrlatt2(3,1) & 346 & +kptrlatt2(1,3)*kptrlatt2(2,1)*kptrlatt2(3,2) & 347 & -kptrlatt2(1,2)*kptrlatt2(2,1)*kptrlatt2(3,3) & 348 & -kptrlatt2(1,3)*kptrlatt2(2,2)*kptrlatt2(3,1) & 349 & -kptrlatt2(1,1)*kptrlatt2(2,3)*kptrlatt2(3,2) 350 351 !Check whether the number of k points is positive, 352 !otherwise, change the handedness of kptrlatt2 353 if(nkptlatt<=0)then 354 ! write(std_out,*)' getkgrid : nkptlatt is negative !' 355 kptrlatt2(:,3)=-kptrlatt2(:,3) 356 nkptlatt=-nkptlatt 357 do ishiftk=1,nshiftk2 358 shiftk2(3,ishiftk)=-shiftk2(3,ishiftk) 359 end do 360 end if 361 362 !Determine the smallest supercell R-vector whose contribution 363 !is not taken correctly into account in the k point integration. 364 !Increase enormously the size of the cell when vacuum is present. 365 fact_vacuum(:)=1 366 if(vacuum(1)==1)fact_vacuum(1)=1000.0_dp 367 if(vacuum(2)==1)fact_vacuum(2)=1000.0_dp 368 if(vacuum(3)==1)fact_vacuum(3)=1000.0_dp 369 do ii=1,3 370 rprimd_super(:,ii)=fact_vacuum(1)*rprimd(:,1)*kptrlatt2(1,ii)+& 371 & fact_vacuum(2)*rprimd(:,2)*kptrlatt2(2,ii)+& 372 & fact_vacuum(3)*rprimd(:,3)*kptrlatt2(3,ii) 373 end do 374 375 call metric(gmet_super,gprimd_super,-1,rmet_super,rprimd_super,ucvol_super) 376 call smallprim(metmin,minim,rprimd_super) 377 length2=min(metmin(1,1),metmin(2,2),metmin(3,3)) 378 kptrlen=sqrt(length2) 379 380 !write(message,'(a,es16.6)' )' getkgrid : length of smallest supercell vector (bohr)=',kptrlen 381 !call wrtout(std_out,message,'COLL') 382 ! If the number of shifts has been decreased, determine the set of kptrlatt2 vectors 383 ! with minimal length (without using fact_vacuum) 384 ! It is worth to determine the minimal set of vectors so that the kptrlatt that is output 385 ! does not seem screwy, although correct but surprising. 386 if(nshiftk/=nshiftk2)then 387 do ii=1,3 388 rprimd_super(:,ii)=rprimd(:,1)*kptrlatt2(1,ii)+rprimd(:,2)*kptrlatt2(2,ii)+rprimd(:,3)*kptrlatt2(3,ii) 389 end do 390 call metric(gmet_super,gprimd_super,-1,rmet_super,rprimd_super,ucvol_super) 391 ! Shift vectors in cartesian coordinates (reciprocal space) 392 do ishiftk=1,nshiftk2 393 shiftk3(:,ishiftk)=gprimd_super(:,1)*shiftk2(1,ishiftk)+& 394 & gprimd_super(:,2)*shiftk2(2,ishiftk)+& 395 & gprimd_super(:,3)*shiftk2(3,ishiftk) 396 end do 397 call smallprim(metmin,minim,rprimd_super) 398 call metric(gmet_super,gprimd_super,-1,rmet_super,minim,ucvol_super) 399 ! This is the new kptrlatt2 400 do ii=1,3 401 dijk(:)=gprimd(1,:)*minim(1,ii)+& 402 & gprimd(2,:)*minim(2,ii)+& 403 & gprimd(3,:)*minim(3,ii) 404 kptrlatt2(:,ii)=nint(dijk(:)) 405 end do 406 ! Shifts in the new set of kptrlatt vectors 407 do ishiftk=1,nshiftk2 408 shiftk2(:,ishiftk)=minim(1,:)*shiftk3(1,ishiftk)+& 409 & minim(2,:)*shiftk3(2,ishiftk)+& 410 & minim(3,:)*shiftk3(3,ishiftk) 411 end do 412 end if 413 414 !brav=1 is able to treat all bravais lattices. 415 brav=1 416 mkpt=nkptlatt*nshiftk2 417 418 ABI_ALLOCATE(spkpt,(3,mkpt)) 419 option=0 420 if(iout/=0)option=1 421 422 if (PRESENT(downsampling))then 423 call smpbz(brav,iout,kptrlatt2,mkpt,nkpthf_computed,nshiftk2,option,shiftk2,spkpt,downsampling=downsampling) 424 if (PRESENT(kpthf) .and. nkpthf/=0) then ! Returns list of k-points in the Full BZ, possibly downsampled for Fock 425 kpthf = spkpt(:,1:nkpthf) 426 end if 427 nkpthf=nkpthf_computed 428 429 end if 430 431 call smpbz(brav,iout,kptrlatt2,mkpt,nkpt_fullbz,nshiftk2,option,shiftk2,spkpt) 432 433 if (PRESENT(fullbz)) then ! Returns list of k-points in the Full BZ. 434 ABI_ALLOCATE(fullbz,(3,nkpt_fullbz)) 435 fullbz = spkpt(:,1:nkpt_fullbz) 436 end if 437 438 if(kptopt==1 .or. kptopt==2 .or. kptopt==4)then 439 440 ABI_ALLOCATE(indkpt,(nkpt_fullbz)) 441 ABI_ALLOCATE(kpt_fullbz,(3,nkpt_fullbz)) 442 ABI_ALLOCATE(wtk_fullbz,(nkpt_fullbz)) 443 ABI_ALLOCATE(wtk_folded,(nkpt_fullbz)) 444 445 kpt_fullbz(:,:)=spkpt(:,1:nkpt_fullbz) 446 wtk_fullbz(1:nkpt_fullbz)=1.0_dp/dble(nkpt_fullbz) 447 448 timrev=1;if (kptopt==4) timrev=0 449 450 call symkpt(chksymbreak,gmet,indkpt,iout,kpt_fullbz,nkpt_fullbz,& 451 & nkpt_computed,nsym_used,symrec,timrev,wtk_fullbz,wtk_folded) 452 453 ABI_DEALLOCATE(symrec) 454 ABI_DEALLOCATE(wtk_fullbz) 455 456 else if(kptopt==3)then 457 nkpt_computed=nkpt_fullbz 458 end if 459 460 !The number of k points has been computed from kptopt, kptrlatt, nshiftk, shiftk, 461 !and the eventual symmetries, it is presently called nkpt_computed. 462 463 !Check that the argument nkpt is coherent with nkpt_computed, if nkpt/=0. 464 if(nkpt/=nkpt_computed .and. nkpt/=0)then 465 write(message, '(a,i6,5a,i6,7a)') & 466 & 'The argument nkpt=',nkpt,', does not match',ch10,& 467 & 'the number of k points generated by kptopt, kptrlatt, shiftk,',ch10,& 468 & 'and the eventual symmetries, that is, nkpt=',nkpt_computed,'.',ch10,& 469 & 'However, note that it might be due to the user,',ch10,& 470 & 'if nkpt is explicitely defined in the input file.',ch10,& 471 & 'In this case, please check your input file.' 472 MSG_BUG(message) 473 end if 474 475 if(kptopt==1 .or. kptopt==2 .or. kptopt==4)then 476 477 if(nkpt/=0)then 478 do ikpt=1,nkpt 479 kpt(:,ikpt)=kpt_fullbz(:,indkpt(ikpt)) 480 if(iscf>=0 .or. iscf==-3 .or. iscf==-1.or.iscf==-2)wtk(ikpt)=wtk_folded(indkpt(ikpt)) 481 end do 482 end if 483 484 ABI_DEALLOCATE(indkpt) 485 ABI_DEALLOCATE(kpt_fullbz) 486 ABI_DEALLOCATE(spkpt) 487 ABI_DEALLOCATE(wtk_folded) 488 489 else if(kptopt==3)then 490 491 if(nkpt/=0)then 492 kpt(:,1:nkpt)=spkpt(:,1:nkpt) 493 if(iscf>1 .or. iscf==-3 .or. iscf==-1.or.iscf==-2)wtk(1:nkpt)=1.0_dp/dble(nkpt) 494 end if 495 ABI_DEALLOCATE(spkpt) 496 497 end if 498 499 kptrlatt(:,:)=kptrlatt2(:,:) 500 nshiftk=nshiftk2 501 shiftk(:,1:nshiftk)=shiftk2(:,1:nshiftk) 502 ABI_DEALLOCATE(shiftk2) 503 ABI_DEALLOCATE(shiftk3) 504 505 end subroutine getkgrid