TABLE OF CONTENTS
ABINIT/get_npert_rbz [ Functions ]
NAME
get_npert_rbz
FUNCTION
Get the number of effective pertubation done in looper3, nkpt_rbz, nband_rbz
COPYRIGHT
Copyright (C) 1998-2017 ABINIT group (FJ) 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
dtset <type(dataset_type)>=all input variables for this dataset
OUTPUT
npert=number of effective pertubation done in looper3 nkpt_rbz= nkpt in the reduced brillouin zone nband_rbz= nband in the reduced brillouin zone
PARENTS
finddistrproc,initmpi_pert,mpi_setup
CHILDREN
irreducible_set_pert,littlegroup_pert,littlegroup_q,mati3inv,metric mkrdim,symatm,symkpt,wrtout
SOURCE
33 #if defined HAVE_CONFIG_H 34 #include "config.h" 35 #endif 36 37 #include "abi_common.h" 38 39 subroutine get_npert_rbz(dtset,nband_rbz,nkpt_rbz,npert) 40 41 use defs_basis 42 use defs_abitypes 43 use m_errors 44 use m_profiling_abi 45 46 !This section has been created automatically by the script Abilint (TD). 47 !Do not modify the following lines by hand. 48 #undef ABI_FUNC 49 #define ABI_FUNC 'get_npert_rbz' 50 use interfaces_14_hidewrite 51 use interfaces_32_util 52 use interfaces_41_geometry 53 !End of the abilint section 54 55 implicit none 56 57 !Arguments ------------------------------------ 58 !scalars 59 integer,intent(out) :: npert 60 !arrays 61 integer,pointer :: nkpt_rbz(:) 62 real(dp),pointer :: nband_rbz(:,:) 63 type(dataset_type),intent(in) :: dtset 64 65 !Local variables------------------------------- 66 !scalars 67 integer :: icase,idir,ikpt,ikpt1,ipert,isppol,isym,maxidir,mpert,nband_k,nsym1,timrev,timrev_pert 68 integer :: to_compute_this_pert 69 real(dp) :: tolsym8,ucvol 70 character(len=500) :: message 71 !arrays 72 integer :: rfdir(9),rf2dir(9),rf2_dir1(3),rf2_dir2(3) 73 integer,allocatable :: indkpt1(:,:),indsym(:,:,:),pertsy(:,:),rfpert(:),symq(:,:,:),symrec(:,:,:) 74 integer, allocatable :: pert_tmp(:,:), pert_calc(:,:) 75 integer,allocatable :: symaf1(:),symrc1(:,:,:),symrl1(:,:,:),symrl1_tmp(:,:,:) 76 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3) 77 real(dp),allocatable :: tnons1_tmp(:,:),wtk_folded(:) 78 79 ! ************************************************************************* 80 81 !Define the set of admitted perturbations 82 mpert=dtset%natom+6 83 if(dtset%natom+10/=0.or.dtset%natom+11/=0) mpert=dtset%natom+11 84 85 ABI_ALLOCATE(symrec,(3,3,dtset%nsym)) 86 !Get the symmetry matrices in terms of reciprocal basis 87 do isym=1,dtset%nsym 88 call mati3inv(dtset%symrel(:,:,isym),symrec(:,:,isym)) 89 end do 90 91 ABI_ALLOCATE(indsym,(4,dtset%nsym,dtset%natom)) 92 !Obtain a list of rotated atom labels: 93 tolsym8=tol8 94 call symatm(indsym,dtset%natom,dtset%nsym,symrec,dtset%tnons,tolsym8,dtset%typat,dtset%xred_orig) 95 96 ABI_ALLOCATE(symq,(4,2,dtset%nsym)) 97 timrev=1 98 call littlegroup_q(dtset%nsym,dtset%qptn,symq,symrec,dtset%symafm,timrev) 99 100 !Initialize the list of perturbations rfpert 101 ABI_ALLOCATE(rfpert,(mpert)) 102 rfpert(:)=0 103 if(dtset%rfphon==1)rfpert(dtset%rfatpol(1):dtset%rfatpol(2))=1 104 105 if(dtset%rfddk==1)rfpert(dtset%natom+1)=1 106 if(dtset%rfddk==2)rfpert(dtset%natom+6)=1 107 108 if(dtset%rf2_dkdk/=0) rfpert(dtset%natom+10)=1 109 if(dtset%rf2_dkde/=0) rfpert(dtset%natom+11)=1 110 111 if(dtset%rfelfd==1.or.dtset%rfelfd==2)rfpert(dtset%natom+1)=1 112 if(dtset%rfelfd==1.or.dtset%rfelfd==3)rfpert(dtset%natom+2)=1 113 114 if(dtset%rfstrs==1.or.dtset%rfstrs==3)rfpert(dtset%natom+3)=1 115 if(dtset%rfstrs==2.or.dtset%rfstrs==3)rfpert(dtset%natom+4)=1 116 117 if(dtset%rfuser==1.or.dtset%rfuser==3)rfpert(dtset%natom+6)=1 118 if(dtset%rfuser==2.or.dtset%rfuser==3)rfpert(dtset%natom+7)=1 119 120 if(dtset%rfmagn==1) rfpert(dtset%natom+5)=1 121 122 ABI_ALLOCATE(pertsy,(3,mpert)) 123 call irreducible_set_pert(indsym,mpert,dtset%natom,dtset%nsym,pertsy,dtset%rfdir,rfpert,symq,symrec,dtset%symrel) 124 npert=0 125 ! ABI_ALLOCATE(pert_tmp,(3*mpert)) 126 127 ! do ipert=1,mpert 128 ! do idir=1,3 129 ! if( rfpert(ipert)==1 .and. dtset%rfdir(idir) == 1 )then 130 ! if (pertsy(idir,ipert)==1.or.& 131 !& (dtset%prepanl == 1.and.ipert == dtset%natom+2).or.& 132 !& (dtset%prepgkk == 1.and.ipert <= dtset%natom) ) then 133 ! npert = npert+1; 134 ! pert_tmp(npert) = idir+(ipert-1)*3; 135 ! else 136 ! write(message, '(a,a,i0,a,i0,a,a,a,a,a,a)' )ch10,& 137 !& 'The perturbation idir=',idir,' ipert=',ipert,' is',ch10,& 138 !& 'symmetric of a previously calculated perturbation.',ch10,& 139 !& 'So, its SCF calculation is not needed.',ch10 140 ! call wrtout(std_out,message,'COLL') 141 ! end if ! Test of existence of symmetry of perturbation 142 ! end if ! Test of existence of perturbation 143 ! end do 144 ! end do 145 146 !Initialize rf2dir : 147 rf2_dir1(1:3)=dtset%rf2_pert1_dir(1:3) 148 rf2_dir2(1:3)=dtset%rf2_pert2_dir(1:3) 149 !Diagonal terms : 150 rf2dir(1) = rf2_dir1(1)*rf2_dir2(1) 151 rf2dir(2) = rf2_dir1(2)*rf2_dir2(2) 152 rf2dir(3) = rf2_dir1(3)*rf2_dir2(3) 153 !Upper triangular terms : 154 rf2dir(4) = rf2_dir1(2)*rf2_dir2(3) 155 rf2dir(5) = rf2_dir1(1)*rf2_dir2(3) 156 rf2dir(6) = rf2_dir1(1)*rf2_dir2(2) 157 !Lower triangular terms : 158 rf2dir(7) = rf2_dir1(3)*rf2_dir2(2) 159 rf2dir(8) = rf2_dir1(3)*rf2_dir2(1) 160 rf2dir(9) = rf2_dir1(2)*rf2_dir2(1) 161 162 !Determine existence of pertubations and of pertubation symmetries 163 !Create array with pertubations which have to be calculated 164 ABI_ALLOCATE(pert_tmp,(2,3*(dtset%natom+6)+18)) 165 166 do ipert=1,mpert 167 if (ipert<dtset%natom+10) then 168 maxidir = 3 169 rfdir(1:3) = dtset%rfdir(:) 170 rfdir(4:9) = 0 171 else 172 maxidir = 9 173 rfdir(1:9) = rf2dir(:) 174 end if 175 do idir=1,maxidir 176 if( rfpert(ipert)==1 .and. rfdir(idir) == 1 )then 177 to_compute_this_pert = 0 178 if (ipert>=dtset%natom+10) then 179 to_compute_this_pert = 1 180 else if ((pertsy(idir,ipert)==1).or.& 181 & ((dtset%prepanl == 1).and.(ipert == dtset%natom+2)).or.& 182 & ((dtset%prepgkk == 1).and.(ipert <= dtset%natom)) ) then 183 to_compute_this_pert = 1 184 end if 185 if (to_compute_this_pert /= 0) then 186 npert = npert+1; 187 pert_tmp(1,npert) = ipert 188 pert_tmp(2,npert) = idir 189 else 190 write(message, '(a,a,i4,a,i4,a,a,a,a,a,a)' )ch10,& 191 & ' The perturbation idir=',idir,' ipert=',ipert,' is',ch10,& 192 & ' symmetric of a previously calculated perturbation.',ch10,& 193 & ' So, its SCF calculation is not needed.',ch10 194 call wrtout(std_out,message,'COLL') 195 end if ! Test of existence of symmetry of perturbation 196 end if ! Test of existence of perturbation 197 end do 198 end do 199 ABI_ALLOCATE(pert_calc,(2,npert)) 200 do icase=1,npert 201 pert_calc(:,icase)=pert_tmp(:,icase) 202 end do 203 ABI_DEALLOCATE(pert_tmp) 204 ABI_DEALLOCATE(pertsy) 205 ABI_DEALLOCATE(rfpert) 206 207 ! Write YAML doc with the list of irreducible perturbations. Example. 208 ! 209 !--- !IrredPerts 210 !# List of irreducible perturbations 211 !irred_perts: 212 ! - qpt: [ 0.0000000000000000, 0.0000000000000000, 0.0000000000000000] 213 ! ipert : 1 214 ! idir : 1 215 ! - qpt: [ 0.0000000000000000, 0.0000000000000000, 0.0000000000000000] 216 ! ipert : 2 217 ! idir : 1 218 !.. 219 write(std_out,'(a)')"--- !IrredPerts" 220 write(std_out,'(a)')'# List of irreducible perturbations' 221 write(std_out,'(a)')'irred_perts:' 222 223 do icase=1,npert 224 ! pert = pert_tmp(icase) 225 226 ! if (pert <= dtset%natom*3) then 227 ! idir = mod(pert, 3) 228 ! if (idir==0) idir=3 229 ! ipert=((pert-idir) / 3 + 1) 230 ! else 231 ! idir = mod(pert, 3) 232 ! if (idir==0) idir=3 233 ! ipert = dtset%natom + ((pert - 3*dtset%natom - 1) / 3) + 1 234 ! end if 235 ipert = pert_calc(1,icase) 236 idir = pert_calc(2,icase) 237 238 write(std_out,'(a,3(f20.16,a))')" - qpt: [ ",dtset%qptn(1),", ", dtset%qptn(2),", ", dtset%qptn(3),"]" 239 write(std_out,'(a,i0)')" ipert: ",ipert 240 write(std_out,'(a,i0)')" idir: ",idir 241 end do 242 243 write(std_out,'(a)')"..." 244 245 ! ABI_ALLOCATE(pert_calc,(npert)) 246 ! do icase=1,npert 247 ! pert_calc(icase) = pert_tmp(icase) 248 ! end do 249 250 call mkrdim(dtset%acell_orig(1:3,1),dtset%rprim_orig(1:3,1:3,1),rprimd) 251 call metric(gmet,gprimd,std_out,rmet,rprimd,ucvol) 252 253 ABI_ALLOCATE(nkpt_rbz,(npert)) 254 ABI_ALLOCATE(indkpt1,(dtset%nkpt,npert)) 255 indkpt1=0 256 257 do icase=1,npert 258 ! if (pert_calc(icase) <= dtset%natom*3) then 259 ! idir = mod(pert_calc(icase),3) 260 ! if (idir==0) idir=3 261 ! ipert=( (pert_calc(icase)-idir) / 3 + 1) 262 ! else 263 ! ipert = dtset%natom + ((pert_calc(icase) - 3*dtset%natom - 1) / 3) + 1 264 ! idir = mod(pert_calc(icase),3) 265 ! if (idir==0) idir=3 266 ! end if 267 ipert = pert_calc(1,icase) 268 idir = pert_calc(2,icase) 269 270 ABI_ALLOCATE(symrl1_tmp,(3,3,dtset%nsym)) 271 ABI_ALLOCATE(symaf1,(dtset%nsym)) 272 ABI_ALLOCATE(tnons1_tmp,(3,dtset%nsym)) 273 ! MJV TODO: check whether prepgkk should be used here 274 if (dtset%prepanl /= 1 .and. dtset%berryopt /=4 .and. dtset%berryopt /=6 .and. dtset%berryopt /=7 .and. & 275 & dtset%berryopt /=14 .and. dtset%berryopt /=16 .and. dtset%berryopt /=17) then !!HONG 276 call littlegroup_pert(gprimd,idir,indsym,std_out,ipert,dtset%natom,dtset%nsym,nsym1,2,& 277 & dtset%symafm,symaf1,symq,symrec,& 278 & dtset%symrel,symrl1_tmp,0,dtset%tnons,tnons1_tmp) 279 else 280 nsym1 = 1 281 end if 282 ABI_DEALLOCATE(tnons1_tmp) 283 ABI_DEALLOCATE(symaf1) 284 285 ABI_ALLOCATE(symrc1,(3,3,nsym1)) 286 ABI_ALLOCATE(symrl1,(3,3,nsym1)) 287 symrl1(:,:,1:nsym1)=symrl1_tmp(:,:,1:nsym1) 288 ABI_DEALLOCATE(symrl1_tmp) 289 do isym=1,nsym1 290 call mati3inv(symrl1(:,:,isym),symrc1(:,:,isym)) 291 end do 292 ABI_DEALLOCATE(symrl1) 293 294 ABI_ALLOCATE(wtk_folded,(dtset%nkpt)) 295 timrev_pert=timrev 296 if(dtset%ieig2rf>0) then 297 call symkpt(0,gmet,indkpt1(:,icase),std_out,dtset%kptns,dtset%nkpt,nkpt_rbz(icase),& 298 & 1,symrc1,0,dtset%wtk,wtk_folded) 299 else 300 ! For the time being, the time reversal symmetry is not used 301 ! for ddk, elfd, mgfd perturbations. 302 if(ipert==dtset%natom+1 .or. ipert==dtset%natom+10 .or. ipert==dtset%natom+11 .or. & 303 & ipert==dtset%natom+2 .or. dtset%berryopt==4 .or. dtset%berryopt==6 .or. dtset%berryopt==7 & 304 & .or. dtset%berryopt==14 .or. dtset%berryopt==16 .or. dtset%berryopt==17 )timrev_pert=0 !!HONG 305 call symkpt(0,gmet,indkpt1(:,icase),std_out,dtset%kptns,dtset%nkpt,nkpt_rbz(icase),& 306 nsym1,symrc1,timrev_pert,dtset%wtk,wtk_folded) 307 end if 308 ABI_DEALLOCATE(wtk_folded) 309 ABI_DEALLOCATE(symrc1) 310 end do 311 312 ABI_ALLOCATE(nband_rbz,(maxval(nkpt_rbz)*dtset%nsppol,npert)) 313 nband_rbz=zero 314 do icase=1,npert 315 do isppol=1,dtset%nsppol 316 ikpt1=1 317 do ikpt=1,dtset%nkpt 318 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 319 ! Must test against ikpt1/=nkpt_rbz+1, before evaluate indkpt1(ikpt1) 320 if(ikpt1/=nkpt_rbz(icase)+1)then 321 if(ikpt==indkpt1(ikpt1,icase))then 322 nband_rbz(ikpt1+(isppol-1)*nkpt_rbz(icase),icase)=nband_k 323 end if 324 end if 325 end do 326 end do 327 328 end do 329 330 ABI_DEALLOCATE(indkpt1) 331 ABI_DEALLOCATE(symq) 332 ABI_DEALLOCATE(symrec) 333 ABI_DEALLOCATE(indsym) 334 ABI_DEALLOCATE(pert_calc) 335 336 end subroutine get_npert_rbz