TABLE OF CONTENTS
ABINIT/m_outwant [ Modules ]
NAME
m_outwant
FUNCTION
Interface with want code.
COPYRIGHT
Copyright (C) 2005-2024 ABINIT group (CMorari) 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.
SOURCE
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 module m_outwant 24 25 use defs_basis 26 use m_errors 27 use m_abicore 28 use m_hdr 29 use m_dtset 30 31 use m_io_tools, only : open_file 32 use m_symtk, only : matr3inv 33 34 implicit none 35 36 private 37 public :: outwant 38 39 contains
m_outwant/outwant [ Functions ]
[ Top ] [ m_outwant ] [ Functions ]
NAME
outwant
FUNCTION
This routine creates an output file containing all the information needed to run WanT as a post-processing program The resulting file is 'launch.dat'. The routine writes to the disk (unformatted file unitwnt) the following information: alat - lattice parameter rprim - primitive translation vectors ntypat - nr of atom types in elementary cell tpat - nr of types of atoms in the elementary cell xcart - cartesian coordinates of the atoms in the elem. cell ecut - energy cut-off mband - # of bands taken in calculation (same for each K-pt) nk(3) - # of k-pts for each direction (uniform grid in the WHOLE BZ) s0(3) - the origin of the K-space kg_tmp(3,mpw*mkmem ) - reduced planewave coordinates imax - Maximum index of a G vector among all k points (see explanation bellow) nkpt - total no of K-pts nsppol - nr of spin polarisations (1 or 2) eig(mband, nkp_tot) - eigenvalues/band/K_point ngfft(3) - nr of points used for FFT in each direction wfc(i)- cmplx(cg(1,i),cg(2,i)) - wavefunction
INPUTS
dtset <type(dataset_type)>=all input variables for this dataset eig(mband*nkpt*nsppol) = array for holding eigenvalues (Hartree) cg(2,mcg) = planewave coefficients of wavefunction kg(3, mpw*mkmem) = reduced planewave coordinates npwarr(nkpt) = number of planewaves in basis at this k-point mband = maximum number of bands mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol nkpt = number of k - points nsppol = 1 for unpolarized, 2 for spin polarized nspinor = number of spinorial components of the wavefunction (on current proc) mkmem = number of k points treated by this node. mpw = maximum dimensioned size of npw prtwant = if set to 1, print 0 in S0 output
OUTPUT
(only writing)
SOURCE
90 subroutine outwant(dtset,eig,cg,kg,npwarr,mband,mcg,nkpt,nsppol,mkmem,mpw,prtwant) 91 92 !Arguments ------------------------------------ 93 !scalars 94 integer :: mband,mcg,mkmem,mpw,nkpt,nsppol,prtwant 95 type(dataset_type),intent(in) :: dtset 96 !arrays 97 integer :: kg(3,mpw*mkmem),npwarr(nkpt) 98 real(dp) :: cg(2,mcg),eig(mband*nkpt*nsppol) 99 100 !Local variables------------------------------- 101 ! the following variables are not used; they are written to 'launch.dat' 102 ! in order to be compatible with the WANT format 103 !scalars 104 integer :: bandtot,i,icount,ifind,ig,ii,iij,ik,ik_,imax 105 integer :: index,index1,ispin_,iunit,iwf,iwf_k,j,k 106 integer :: maxat,ngm,ngw_,nk_,nkp,nspin_ 107 integer :: unitwnt 108 real(dp) :: alat,scal,scal_,tt 109 logical :: twrite=.true. 110 character(len=20) :: section_name 111 character(len=3) :: nameat 112 character(len=500) :: message 113 character(len=fnlen) :: filewnt 114 !arrays 115 integer :: ikg(3),nk(3) 116 integer,allocatable :: iwfi(:,:),kg_tmp(:,:),tpat(:) 117 real(dp) :: drprim(3,3),gmat(3,3),gmod(3),s0(3),t1(3),t2(3) 118 real(dp),allocatable :: xcoord(:,:,:) 119 complex,allocatable :: wfc(:) 120 121 ! *************************************************************************** 122 123 !WARNING: not tested for nsppol,nspinor >1 124 ! 125 !Initialisations 126 nameat = ' ' 127 bandtot=mband*nkpt*nsppol*dtset%nspinor 128 filewnt='launch.dat' 129 130 write(message,'(3a)')ch10,' Opening file for WanT input: ',trim(filewnt) 131 call wrtout(std_out,message,'COLL') 132 133 !Open the file 134 if (open_file(filewnt,message,newunit=unitwnt,form='unformatted', status='unknown') /=0) then 135 ABI_ERROR(message) 136 end if 137 138 !Comments 139 if(prtwant>1) then 140 write(std_out,*) 'Wrong value for prtwant. Reseting to 1' 141 prtwant=1 142 elseif(prtwant==1) then 143 do i=1,3 144 s0(i)=0._dp 145 end do 146 end if 147 148 !Discussion of 'alat' ABINIT/ WanT 149 if(dtset%acell_orig(1,1)==dtset%acell_orig(2,1).and.& 150 dtset%acell_orig(1,1)==dtset%acell_orig(3,1)) then 151 alat=dtset%acell_orig(1,1) 152 do i=1,3 153 do j=1,3 154 drprim( i, j) = dtset%rprim_orig( i, j, 1 ) 155 end do 156 end do 157 else 158 ! Redefining the drprim( i, j) 159 alat=dtset%acell_orig(1,1) 160 do i=1,3 161 do j=1,3 162 drprim( i, j) = dtset%rprim_orig( i, j, 1 )*dtset%acell_orig(j, 1)/alat 163 end do 164 end do 165 end if 166 167 !Now finding the no of k-pt for each direction PARALEL with the 168 !generators of the first B.Z. 169 !First decide if we have the Gamma point in the list; its index in the list is ... index 170 nk(:)=1 171 ifind=0 172 icount=2 173 do i=1,nkpt 174 index1=0 175 do j=1,3 176 if(dtset%kptns(j,i)<tol8) index1=index1+1 177 end do 178 if(index1==3) then 179 index=i 180 ifind=1 181 cycle 182 end if 183 end do 184 if(ifind==0) then 185 write(std_out,*) 'GAMMA POINT NOT IN THE LIST OF KPTS?' 186 do ii=1,nkpt 187 write(std_out,*) (dtset%kptns(j,ii),j=1,3) 188 end do 189 ABI_ERROR("fatal error") 190 end if 191 192 call matr3inv(drprim,gmat) 193 194 !Modules for each vector in recip. space; nb: g(index coord, index point) 195 do j=1,3 196 gmod(j)=0.D0 197 do i=1,3 198 gmod(j)=gmod(j)+gmat(i,j)**2 199 end do 200 gmod(j)=sqrt(gmod(j)) 201 end do 202 if(nkpt==2) then 203 do j=1,3 204 do ii=1,3 205 t1(ii)=dtset%kptns(ii,1)-dtset%kptns(ii,2) 206 end do 207 tt=0._dp 208 do iij=1,3 209 t2(iij)=0._dp 210 do ii=1,3 211 t2(iij)=t2(iij)+t1(ii)*gmat(ii,iij) 212 end do 213 tt=tt + t2(iij)**2 214 end do 215 tt=sqrt(tt) 216 scal=0._dp 217 do ii=1,3 218 scal=scal+t2(ii)*gmat(j,ii) 219 end do 220 scal=abs(scal) 221 ! Compare scal(tt,gmat) with simple product of modules -> paralel or not 222 if(abs(scal-tt*gmod(j))<tol8) nk(j)=2 223 end do 224 225 elseif(nkpt>2) then 226 227 do i=1,nkpt 228 if(i.ne.index) then 229 do ii=1,3 230 t1(ii)=dtset%kptns(ii,index)-dtset%kptns(ii,i) 231 end do 232 tt=0._dp 233 do iij=1,3 234 t2(iij)=0._dp 235 do ii=1,3 236 t2(iij)=t2(iij)+t1(ii)*gmat(ii,iij) 237 end do 238 tt=tt + t2(iij)**2 239 end do 240 tt=sqrt(tt) 241 ! check for each direction in the BZ 242 do j=1,3 243 scal=0._dp 244 do ii=1,3 245 scal=scal+t2(ii)*gmat(j,ii) 246 end do 247 scal=abs(scal) 248 ! Compare scal(t1,gmat) with simple product of modules -> paralel or not 249 if(abs(scal-tt*gmod(j))<tol8) nk(j)=nk(j)+1 250 end do 251 end if 252 end do 253 end if 254 index=1 255 do i=1,3 256 index=index*nk(i) 257 end do 258 259 if(index.ne.nkpt) then 260 write(message,'(a,2i0)')' OutwanT: Wrong assignemt of kpts', index,nkpt 261 ABI_ERROR(message) 262 end if 263 264 !End counting/assigning no of kpts/direction 265 !Reordering the coordinates of all atoms - xcoord array 266 ABI_MALLOC(tpat,(dtset%ntypat)) 267 tpat(:)=zero 268 do i=1,dtset%natom 269 do j=1,dtset%ntypat 270 if(dtset%typat(i)==j) tpat(j)=tpat(j)+1 271 end do 272 end do 273 maxat=maxval(tpat(:)) 274 ABI_MALLOC(xcoord,(3,maxat,dtset%ntypat)) 275 index=1 276 do i=1, dtset%ntypat 277 do k=1,tpat(i) 278 do j=1,3 279 xcoord(j,k,i)=dtset%xred_orig(j,index,1) 280 end do 281 index=index+1 282 end do 283 end do 284 ! 285 !Defining the kg_tmp list 286 !Preparing the output of reduced coords., in a single list (kg_tmp(3,imax)) 287 !We start with kg_tmp(:,i)=kg(:,i=1,npwarr(1)) then the new coordinates are added 288 !ONLY if they are not allready in the list. An index is associated 289 !for each kg_tmp which allow us to recover kg(3,mpw*nkpt) from 290 !the smaller list kg_tmp(3, imax) 291 ABI_MALLOC(kg_tmp,(3,mpw*nkpt)) 292 ABI_MALLOC(iwfi,(nkpt,mpw)) 293 kg_tmp(:,:)=zero 294 iwfi(:,:)=zero 295 imax=npwarr(1) 296 index=0 297 298 do i=1, nkpt 299 if(i>1) then 300 index=index+npwarr(i-1) 301 end if 302 do j=1, npwarr(i) 303 if(i.eq.1) then 304 iwfi(i,j)=j 305 if(mkmem>0) kg_tmp(:,j)=kg(:,j) 306 else 307 ifind=0 308 if(mkmem>0) ikg(:)=kg(:,index+j) 309 310 do k=1,imax 311 if(ikg(1)==kg_tmp(1,k)) then 312 if(ikg(2)==kg_tmp(2,k)) then 313 if(ikg(3)==kg_tmp(3,k)) then 314 ifind=1 315 iwfi(i,j)=k 316 end if 317 end if 318 end if 319 end do 320 321 if(ifind==0) then 322 imax=imax+1 323 kg_tmp(:,imax)=ikg(:) 324 iwfi(i,j)=imax 325 end if 326 end if 327 end do 328 end do 329 ngm=imax 330 331 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 332 !PART ONE: writing the header 333 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 334 write(unitwnt) alat 335 write(unitwnt) ( drprim( i, 1 ), i = 1, 3 ) ! save A1 336 write(unitwnt) ( drprim( i, 2 ), i = 1, 3 ) ! save A2 337 write(unitwnt) ( drprim( i, 3 ), i = 1, 3 ) ! save A3 338 !write(std_out,* ) ( drprim( i, 1 ), i = 1, 3 ) ! save A1 339 !write(std_out,* ) ( drprim( i, 2 ), i = 1, 3 ) ! save A2 340 !write(std_out,* ) ( drprim( i, 3 ), i = 1, 3 ) ! save A3 341 342 write(unitwnt) dtset%ntypat 343 !write(std_out,*) dtset%ntypat, 'NTYPAT', tpat 344 345 do i = 1, dtset%ntypat 346 write(unitwnt) tpat(i), nameat 347 write(unitwnt) ((xcoord(j,k,i),j=1,3), k=1, tpat(i)) 348 ! write(std_out,*) tpat(i), nameat 349 ! write(std_out,*) ((xcoord(j,k,i),j=1,3),k=1,tpat(i)), 'XCART' 350 end do 351 ABI_FREE(tpat) 352 ABI_FREE(xcoord) 353 354 !energy cut-off in Rydberg (WANT option) 355 write (unitwnt) 2._dp*dtset%ecut, mband 356 !write(std_out,*) 2._dp*dtset%ecut, mband 357 write (unitwnt) ( nk(i), i = 1, 3 ), ( s0(j), j = 1, 3 ),ngm 358 !write(std_out,*) ( nk(i), i = 1, 3 ), ( s0(j), j = 1, 3 ),imax 359 write (unitwnt) ( kg_tmp( 1, i ), kg_tmp( 2, i ), kg_tmp( 3, i ), i = 1, ngm ) 360 write (unitwnt) mpw, mband, dtset%nkpt/dtset%nsppol 361 !write(std_out,*) mpw, mband, dtset%nkpt/dtset%nsppol 362 363 do i=1, nkpt 364 write(unitwnt) (iwfi(i,j), j=1,mpw) 365 end do 366 ABI_FREE(kg_tmp) 367 368 !Eigenvalues in HARTREE 369 write (unitwnt) ( eig( i ), i = 1, bandtot) 370 write (unitwnt) ( npwarr( ik ), ik = 1, nkpt ) 371 write (unitwnt) ( mband, ik = 1, nkpt ) 372 write (unitwnt) (dtset%ngfft(i),i=1,3), imax, imax 373 !write(std_out,*) ( eig( i ), i = 1, bandtot ) 374 !write(std_out,*) ( npwarr( ik ), ik = 1, nkpt ) 375 !write(std_out,*) ( mband, ik = 1, nkpt ) 376 !write(std_out,*) (dtset%ngfft(i),i=1,3), imax ,imax 377 !a list with the band structure; usefull for 'windows' and 'disentangle' programs 378 !from WanT distribution 379 380 if (open_file('band.gpl',message,newunit=iunit,status='unknown') /=0) then 381 ABI_ERROR(message) 382 end if 383 384 index=1 385 do i=1,mband 386 index=1 387 do j=1,nkpt 388 write(iunit,*) index, Ha_eV*eig(i+(j-1)*mband), eig(i+(j-1)*mband) 389 index=index+1 390 end do 391 write(iunit,*) 392 end do 393 394 close(iunit) 395 396 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 397 !PART TWO: Writing the wavefunction 398 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 399 !not used 400 ngw_=0 401 ik_=0 402 nk_=0 403 ispin_=0 404 nspin_=0 405 scal_=1._dp 406 !!!!!!!!!!!!!!!!!!!!!!!!!! 407 iwf = 1 408 iwf_k=1 409 ABI_MALLOC(wfc,(imax)) 410 411 !Loop over k-pt 412 do nkp=1,nkpt 413 ! Not relevant 414 write(unitwnt) twrite, ik_, section_name 415 ! Only 'mband' is relevant here 416 write(unitwnt) ngw_, mband, ik_, nk_, nk_,ispin_, nspin_, scal_ 417 write(unitwnt) imax 418 ! Not relevant 419 write(unitwnt) twrite 420 ! Loop over bands 421 422 ! Preparing WF 423 do k=1,mband 424 if(mkmem >0) then 425 wfc(:)=zero 426 ! From cg to wf: 427 do i=iwf, iwf+npwarr(nkp)-1 428 index=i-iwf+1 429 wfc(iwfi(nkp,index))=cmplx(cg(1,i), cg(2,i), kind(0._dp)) 430 end do 431 iwf=iwf+npwarr(nkp) 432 else 433 message = 'Wrong mkmem in outwant' 434 ABI_ERROR(message) 435 end if 436 write(unitwnt) (wfc(ig), ig=1,imax) 437 ! End loop over bands 438 end do 439 440 ! Not relevant 441 write(unitwnt) twrite 442 ! Not relevant 443 do i=1,mband 444 write(unitwnt) i 445 end do 446 447 ! End loop over k-pts 448 end do 449 450 ABI_FREE(iwfi) 451 ABI_FREE(wfc) 452 453 call wrtout(std_out,'Closing file','COLL') 454 close(unit=unitwnt) 455 456 end subroutine outwant