TABLE OF CONTENTS


ABINIT/m_outwant [ Modules ]

[ Top ] [ Modules ]

NAME

 m_outwant

FUNCTION

 Interface with want code.

COPYRIGHT

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

PARENTS

CHILDREN

SOURCE

22 #if defined HAVE_CONFIG_H
23 #include "config.h"
24 #endif
25 
26 #include "abi_common.h"
27 
28 module m_outwant
29 
30  use defs_basis
31  use defs_abitypes
32  use m_errors
33  use m_abicore
34  use m_hdr
35 
36  use m_io_tools,   only : open_file
37  use m_symtk,      only : matr3inv
38 
39  implicit none
40 
41  private
42  public ::  outwant
43 
44 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)

PARENTS

      outscfcv

CHILDREN

      matr3inv,wrtout

SOURCE

101 subroutine outwant(dtset,eig,cg,kg,npwarr,mband,mcg,nkpt,nsppol,mkmem,mpw,prtwant)
102 
103 
104 !This section has been created automatically by the script Abilint (TD).
105 !Do not modify the following lines by hand.
106 #undef ABI_FUNC
107 #define ABI_FUNC 'outwant'
108 !End of the abilint section
109 
110  implicit none
111 
112 !Arguments ------------------------------------
113 !scalars
114  integer :: mband,mcg,mkmem,mpw,nkpt,nsppol,prtwant
115  type(dataset_type),intent(in) :: dtset
116 !arrays
117  integer :: kg(3,mpw*mkmem),npwarr(nkpt)
118  real(dp) :: cg(2,mcg),eig(mband*nkpt*nsppol)
119 
120 !Local variables-------------------------------
121 ! the following variables are not used; they are written to 'launch.dat'
122 ! in order to be compatible with the WANT format
123 !scalars
124  integer :: bandtot,i,icount,ifind,ig,ii,iij,ik,ik_,imax
125  integer :: index,index1,ispin_,iunit,iwf,iwf_k,j,k
126  integer :: maxat,ngm,ngw_,nk_,nkp,nspin_
127  integer :: unitwnt
128  real(dp) :: alat,scal,scal_,tt
129  logical :: twrite=.true.
130  character(len=20) :: section_name
131  character(len=3) :: nameat
132  character(len=500) :: message
133  character(len=fnlen) :: filewnt
134 !arrays
135  integer :: ikg(3),nk(3)
136  integer,allocatable :: iwfi(:,:),kg_tmp(:,:),tpat(:)
137  real(dp) :: drprim(3,3),gmat(3,3),gmod(3),s0(3),t1(3),t2(3)
138  real(dp),allocatable :: xcoord(:,:,:)
139  complex,allocatable :: wfc(:)
140 
141 ! ***************************************************************************
142 
143 !WARNING: not tested for nsppol,nspinor >1
144 !
145 !Initialisations
146  nameat = ' '
147  bandtot=mband*nkpt*nsppol*dtset%nspinor
148  filewnt='launch.dat'
149 
150  write(message,'(3a)')ch10,' Opening file for WanT input: ',trim(filewnt)
151  call wrtout(std_out,message,'COLL')
152 
153 !Open the file
154  if (open_file(filewnt,message,newunit=unitwnt,form='unformatted', status='unknown') /=0) then
155    MSG_ERROR(message)
156  end if
157 
158 !Comments
159  if(prtwant>1) then
160    write(std_out,*) 'Wrong value for prtwant. Reseting to 1'
161    prtwant=1
162  elseif(prtwant==1) then
163    do i=1,3
164      s0(i)=0._dp
165    end do
166  end if
167 
168 !Discussion of 'alat' ABINIT/ WanT
169  if(dtset%acell_orig(1,1)==dtset%acell_orig(2,1).and.&
170  dtset%acell_orig(1,1)==dtset%acell_orig(3,1)) then
171    alat=dtset%acell_orig(1,1)
172    do i=1,3
173      do j=1,3
174        drprim( i, j) = dtset%rprim_orig( i, j, 1 )
175      end do
176    end do
177  else
178 !  Redefining the drprim( i, j)
179    alat=dtset%acell_orig(1,1)
180    do i=1,3
181      do j=1,3
182        drprim( i, j) = dtset%rprim_orig( i, j, 1 )*dtset%acell_orig(j, 1)/alat
183      end do
184    end do
185  end if
186 
187 !Now finding the no of k-pt for each direction PARALEL with the
188 !generators of the first B.Z.
189 !First decide if we have the Gamma point in the list; its index in the list is ... index
190  nk(:)=1
191  ifind=0
192  icount=2
193  do i=1,nkpt
194    index1=0
195    do j=1,3
196      if(dtset%kptns(j,i)<tol8) index1=index1+1
197    end do
198    if(index1==3) then
199      index=i
200      ifind=1
201      cycle
202    end if
203  end do
204  if(ifind==0) then
205    write(std_out,*) 'GAMMA POINT NOT IN THE LIST OF KPTS?'
206    do ii=1,nkpt
207      write(std_out,*) (dtset%kptns(j,ii),j=1,3)
208    end do
209    MSG_ERROR("fatal error")
210  end if
211 
212  call matr3inv(drprim,gmat)
213 
214 !Modules for each vector in recip. space; nb: g(index coord, index point)
215  do j=1,3
216    gmod(j)=0.D0
217    do i=1,3
218      gmod(j)=gmod(j)+gmat(i,j)**2
219    end do
220    gmod(j)=sqrt(gmod(j))
221  end do
222  if(nkpt==2) then
223    do j=1,3
224      do ii=1,3
225        t1(ii)=dtset%kptns(ii,1)-dtset%kptns(ii,2)
226      end do
227      tt=0._dp
228      do iij=1,3
229        t2(iij)=0._dp
230        do ii=1,3
231          t2(iij)=t2(iij)+t1(ii)*gmat(ii,iij)
232        end do
233        tt=tt + t2(iij)**2
234      end do
235      tt=sqrt(tt)
236      scal=0._dp
237      do ii=1,3
238        scal=scal+t2(ii)*gmat(j,ii)
239      end do
240      scal=abs(scal)
241 !    Compare scal(tt,gmat) with simple product of modules -> paralel or not
242      if(abs(scal-tt*gmod(j))<tol8) nk(j)=2
243    end do
244 
245  elseif(nkpt>2) then
246 
247    do i=1,nkpt
248      if(i.ne.index) then
249        do ii=1,3
250          t1(ii)=dtset%kptns(ii,index)-dtset%kptns(ii,i)
251        end do
252        tt=0._dp
253        do iij=1,3
254          t2(iij)=0._dp
255          do ii=1,3
256            t2(iij)=t2(iij)+t1(ii)*gmat(ii,iij)
257          end do
258          tt=tt + t2(iij)**2
259        end do
260        tt=sqrt(tt)
261 !      check for each direction in the BZ
262        do j=1,3
263          scal=0._dp
264          do ii=1,3
265            scal=scal+t2(ii)*gmat(j,ii)
266          end do
267          scal=abs(scal)
268 !        Compare scal(t1,gmat) with simple product of modules -> paralel or not
269          if(abs(scal-tt*gmod(j))<tol8) nk(j)=nk(j)+1
270        end do
271      end if
272    end do
273  end if
274  index=1
275  do i=1,3
276    index=index*nk(i)
277  end do
278 
279  if(index.ne.nkpt) then
280    write(message,'(a,2i0)')' OutwanT: Wrong assignemt of kpts', index,nkpt
281    MSG_ERROR(message)
282  end if
283 
284 !End counting/assigning no of kpts/direction
285 !Reordering the coordinates of all atoms - xcoord array
286  ABI_ALLOCATE(tpat,(dtset%ntypat))
287  tpat(:)=zero
288  do i=1,dtset%natom
289    do j=1,dtset%ntypat
290      if(dtset%typat(i)==j) tpat(j)=tpat(j)+1
291    end do
292  end do
293  maxat=maxval(tpat(:))
294  ABI_ALLOCATE(xcoord,(3,maxat,dtset%ntypat))
295  index=1
296  do i=1, dtset%ntypat
297    do k=1,tpat(i)
298      do j=1,3
299        xcoord(j,k,i)=dtset%xred_orig(j,index,1)
300      end do
301      index=index+1
302    end do
303  end do
304 !
305 !Defining the kg_tmp list
306 !Preparing the output of reduced coords., in a single list (kg_tmp(3,imax))
307 !We start with kg_tmp(:,i)=kg(:,i=1,npwarr(1)) then the new coordinates are added
308 !ONLY if they are not allready in the list. An index is associated
309 !for each kg_tmp which allow us to recover kg(3,mpw*nkpt) from
310 !the smaller list kg_tmp(3, imax)
311  ABI_ALLOCATE(kg_tmp,(3,mpw*nkpt))
312  ABI_ALLOCATE(iwfi,(nkpt,mpw))
313  kg_tmp(:,:)=zero
314  iwfi(:,:)=zero
315  imax=npwarr(1)
316  index=0
317 
318  do i=1,  nkpt
319    if(i>1) then
320      index=index+npwarr(i-1)
321    end if
322    do j=1, npwarr(i)
323      if(i.eq.1) then
324        iwfi(i,j)=j
325        if(mkmem>0) kg_tmp(:,j)=kg(:,j)
326      else
327        ifind=0
328        if(mkmem>0) ikg(:)=kg(:,index+j)
329 
330        do k=1,imax
331          if(ikg(1)==kg_tmp(1,k)) then
332            if(ikg(2)==kg_tmp(2,k)) then
333              if(ikg(3)==kg_tmp(3,k)) then
334                ifind=1
335                iwfi(i,j)=k
336              end if
337            end if
338          end if
339        end do
340 
341        if(ifind==0) then
342          imax=imax+1
343          kg_tmp(:,imax)=ikg(:)
344          iwfi(i,j)=imax
345        end if
346      end if
347    end do
348  end do
349  ngm=imax
350 
351 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
352 !PART ONE: writing the header
353 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
354  write(unitwnt) alat
355  write(unitwnt) ( drprim( i, 1 ), i = 1, 3 )  ! save A1
356  write(unitwnt) ( drprim( i, 2 ), i = 1, 3 )  ! save A2
357  write(unitwnt) ( drprim( i, 3 ), i = 1, 3 )  ! save A3
358 !write(std_out,* ) ( drprim( i, 1 ), i = 1, 3 )  ! save A1
359 !write(std_out,* ) ( drprim( i, 2 ), i = 1, 3 )  ! save A2
360 !write(std_out,* ) ( drprim( i, 3 ), i = 1, 3 )  ! save A3
361 
362  write(unitwnt) dtset%ntypat
363 !write(std_out,*) dtset%ntypat, 'NTYPAT', tpat
364 
365  do i = 1, dtset%ntypat
366    write(unitwnt) tpat(i), nameat
367    write(unitwnt) ((xcoord(j,k,i),j=1,3), k=1, tpat(i))
368 !  write(std_out,*) tpat(i), nameat
369 !  write(std_out,*) ((xcoord(j,k,i),j=1,3),k=1,tpat(i)), 'XCART'
370  end do
371  ABI_DEALLOCATE(tpat)
372  ABI_DEALLOCATE(xcoord)
373 
374 !energy cut-off in Rydberg (WANT option)
375  write (unitwnt) 2._dp*dtset%ecut, mband
376 !write(std_out,*)   2._dp*dtset%ecut, mband
377  write (unitwnt) ( nk(i), i = 1, 3 ), ( s0(j), j = 1, 3 ),ngm
378 !write(std_out,*) ( nk(i), i = 1, 3 ), ( s0(j), j = 1, 3 ),imax
379  write (unitwnt) ( kg_tmp( 1, i ), kg_tmp( 2, i ), kg_tmp( 3, i ), i = 1, ngm )
380  write (unitwnt) mpw, mband, dtset%nkpt/dtset%nsppol
381 !write(std_out,*) mpw, mband,  dtset%nkpt/dtset%nsppol
382 
383  do i=1, nkpt
384    write(unitwnt) (iwfi(i,j), j=1,mpw)
385  end do
386  ABI_DEALLOCATE(kg_tmp)
387 
388 !Eigenvalues in HARTREE
389  write (unitwnt)  ( eig( i ), i = 1, bandtot)
390  write (unitwnt) ( npwarr( ik ), ik = 1, nkpt )
391  write (unitwnt) ( mband, ik = 1, nkpt )
392  write (unitwnt) (dtset%ngfft(i),i=1,3), imax, imax
393 !write(std_out,*)  ( eig( i ), i = 1, bandtot )
394 !write(std_out,*) ( npwarr( ik ), ik = 1, nkpt )
395 !write(std_out,*) ( mband, ik = 1, nkpt )
396 !write(std_out,*) (dtset%ngfft(i),i=1,3), imax ,imax
397 !a list with the band structure; usefull for 'windows' and 'disentangle' programs
398 !from WanT distribution
399 
400  if (open_file('band.gpl',message,newunit=iunit,status='unknown') /=0) then
401    MSG_ERROR(message)
402  end if
403 
404  index=1
405  do i=1,mband
406    index=1
407    do j=1,nkpt
408      write(iunit,*) index, Ha_eV*eig(i+(j-1)*mband), eig(i+(j-1)*mband)
409      index=index+1
410    end do
411    write(iunit,*)
412  end do
413 
414  close(iunit)
415 
416 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
417 !PART TWO: Writing the wavefunction
418 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
419 !not used
420  ngw_=0
421  ik_=0
422  nk_=0
423  ispin_=0
424  nspin_=0
425  scal_=1._dp
426 !!!!!!!!!!!!!!!!!!!!!!!!!!
427  iwf = 1
428  iwf_k=1
429  ABI_ALLOCATE(wfc,(imax))
430 
431 !Loop over k-pt
432  do nkp=1,nkpt
433 !  Not relevant
434    write(unitwnt) twrite, ik_, section_name
435 !  Only 'mband' is relevant here
436    write(unitwnt) ngw_, mband, ik_, nk_, nk_,ispin_, nspin_, scal_
437    write(unitwnt) imax
438 !  Not relevant
439    write(unitwnt) twrite
440 !  Loop over bands
441 
442 !  Preparing WF
443    do k=1,mband
444      if(mkmem >0) then
445        wfc(:)=zero
446 !      From cg to wf:
447        do i=iwf, iwf+npwarr(nkp)-1
448          index=i-iwf+1
449          wfc(iwfi(nkp,index))=cmplx(cg(1,i), cg(2,i), kind(0._dp))
450        end do
451        iwf=iwf+npwarr(nkp)
452      else
453        message = 'Wrong mkmem in outwant'
454        MSG_ERROR(message)
455      end if
456      write(unitwnt) (wfc(ig), ig=1,imax)
457 !    End loop over bands
458    end do
459 
460 !  Not relevant
461    write(unitwnt) twrite
462 !  Not relevant
463    do i=1,mband
464      write(unitwnt) i
465    end do
466 
467 !  End loop over k-pts
468  end do
469 
470  ABI_DEALLOCATE(iwfi)
471  ABI_DEALLOCATE(wfc)
472 
473  call wrtout(std_out,'Closing file','COLL')
474  close(unit=unitwnt)
475 
476 end subroutine outwant