TABLE OF CONTENTS


ABINIT/m_outwant [ Modules ]

[ Top ] [ Modules ]

NAME

 m_outwant

FUNCTION

 Interface with want code.

COPYRIGHT

 Copyright (C) 2005-2022 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