Line data Source code
1 : !!****m* ABINIT/m_outwant
2 : !! NAME
3 : !! m_outwant
4 : !!
5 : !! FUNCTION
6 : !! Interface with want code.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2005-2026 ABINIT group (CMorari)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
14 : !!
15 : !! SOURCE
16 :
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_matrix, only : matr3inv
33 :
34 : implicit none
35 :
36 : private
37 : public :: outwant
38 :
39 : contains
40 : !!***
41 :
42 : !!****f* m_outwant/outwant
43 : !! NAME
44 : !! outwant
45 : !!
46 : !! FUNCTION
47 : !! This routine creates an output file containing all the
48 : !! information needed to run WanT as a post-processing program
49 : !! The resulting file is 'launch.dat'.
50 : !!
51 : !! The routine writes to the disk (unformatted file unitwnt) the following information:
52 : !!
53 : !! alat - lattice parameter
54 : !! rprim - primitive translation vectors
55 : !! ntypat - nr of atom types in elementary cell
56 : !! tpat - nr of types of atoms in the elementary cell
57 : !! xcart - cartesian coordinates of the atoms in the elem. cell
58 : !! ecut - energy cut-off
59 : !! mband - # of bands taken in calculation (same for each K-pt)
60 : !! nk(3) - # of k-pts for each direction (uniform grid in the WHOLE BZ)
61 : !! s0(3) - the origin of the K-space
62 : !! kg_tmp(3,mpw*mkmem ) - reduced planewave coordinates
63 : !! imax - Maximum index of a G vector among all k points (see explanation bellow)
64 : !! nkpt - total no of K-pts
65 : !! nsppol - nr of spin polarisations (1 or 2)
66 : !! eig(mband, nkp_tot) - eigenvalues/band/K_point
67 : !! ngfft(3) - nr of points used for FFT in each direction
68 : !! wfc(i)- cmplx(cg(1,i),cg(2,i)) - wavefunction
69 : !!
70 : !! INPUTS
71 : !! dtset <type(dataset_type)>=all input variables for this dataset
72 : !! eig(mband*nkpt*nsppol) = array for holding eigenvalues (Hartree)
73 : !! cg(2,mcg) = planewave coefficients of wavefunction
74 : !! kg(3, mpw*mkmem) = reduced planewave coordinates
75 : !! npwarr(nkpt) = number of planewaves in basis at this k-point
76 : !! mband = maximum number of bands
77 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
78 : !! nkpt = number of k - points
79 : !! nsppol = 1 for unpolarized, 2 for spin polarized
80 : !! nspinor = number of spinorial components of the wavefunction (on current proc)
81 : !! mkmem = number of k points treated by this node.
82 : !! mpw = maximum dimensioned size of npw
83 : !! prtwant = if set to 1, print 0 in S0 output
84 : !!
85 : !! OUTPUT
86 : !! (only writing)
87 : !!
88 : !! SOURCE
89 :
90 1 : 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 1 : integer,allocatable :: iwfi(:,:),kg_tmp(:,:),tpat(:)
117 : real(dp) :: drprim(3,3),gmat(3,3),gmod(3),s0(3),t1(3),t2(3)
118 1 : real(dp),allocatable :: xcoord(:,:,:)
119 1 : complex,allocatable :: wfc(:)
120 :
121 : ! ***************************************************************************
122 :
123 : !WARNING: not tested for nsppol,nspinor >1
124 : !
125 : !Initialisations
126 1 : nameat = ' '
127 1 : bandtot=mband*nkpt*nsppol*dtset%nspinor
128 1 : filewnt='launch.dat'
129 :
130 1 : write(message,'(3a)')ch10,' Opening file for WanT input: ',trim(filewnt)
131 1 : call wrtout(std_out,message,'COLL')
132 :
133 : !Open the file
134 1 : if (open_file(filewnt,message,newunit=unitwnt,form='unformatted', status='unknown') /=0) then
135 0 : ABI_ERROR(message)
136 : end if
137 :
138 : !Comments
139 1 : if(prtwant>1) then
140 0 : write(std_out,*) 'Wrong value for prtwant. Reseting to 1'
141 0 : prtwant=1
142 1 : elseif(prtwant==1) then
143 4 : do i=1,3
144 4 : s0(i)=0._dp
145 : end do
146 : end if
147 :
148 : !Discussion of 'alat' ABINIT/ WanT
149 1 : 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 1 : alat=dtset%acell_orig(1,1)
152 4 : do i=1,3
153 13 : do j=1,3
154 12 : drprim( i, j) = dtset%rprim_orig( i, j, 1 )
155 : end do
156 : end do
157 : else
158 : ! Redefining the drprim( i, j)
159 0 : alat=dtset%acell_orig(1,1)
160 0 : do i=1,3
161 0 : do j=1,3
162 0 : 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 4 : nk(:)=1
171 1 : ifind=0
172 1 : icount=2
173 5 : do i=1,nkpt
174 : index1=0
175 16 : do j=1,3
176 16 : if(dtset%kptns(j,i)<tol8) index1=index1+1
177 : end do
178 5 : if(index1==3) then
179 1 : index=i
180 1 : ifind=1
181 1 : cycle
182 : end if
183 : end do
184 1 : if(ifind==0) then
185 0 : write(std_out,*) 'GAMMA POINT NOT IN THE LIST OF KPTS?'
186 0 : do ii=1,nkpt
187 0 : write(std_out,*) (dtset%kptns(j,ii),j=1,3)
188 : end do
189 0 : ABI_ERROR("fatal error")
190 : end if
191 :
192 1 : call matr3inv(drprim,gmat)
193 :
194 : !Modules for each vector in recip. space; nb: g(index coord, index point)
195 4 : do j=1,3
196 3 : gmod(j)=0.D0
197 12 : do i=1,3
198 12 : gmod(j)=gmod(j)+gmat(i,j)**2
199 : end do
200 4 : gmod(j)=sqrt(gmod(j))
201 : end do
202 1 : if(nkpt==2) then
203 0 : do j=1,3
204 0 : do ii=1,3
205 0 : t1(ii)=dtset%kptns(ii,1)-dtset%kptns(ii,2)
206 : end do
207 : tt=0._dp
208 0 : do iij=1,3
209 0 : t2(iij)=0._dp
210 0 : do ii=1,3
211 0 : t2(iij)=t2(iij)+t1(ii)*gmat(ii,iij)
212 : end do
213 0 : tt=tt + t2(iij)**2
214 : end do
215 0 : tt=sqrt(tt)
216 0 : scal=0._dp
217 0 : do ii=1,3
218 0 : scal=scal+t2(ii)*gmat(j,ii)
219 : end do
220 0 : scal=abs(scal)
221 : ! Compare scal(tt,gmat) with simple product of modules -> paralel or not
222 0 : if(abs(scal-tt*gmod(j))<tol8) nk(j)=2
223 : end do
224 :
225 1 : elseif(nkpt>2) then
226 :
227 5 : do i=1,nkpt
228 5 : if(i.ne.index) then
229 12 : do ii=1,3
230 12 : t1(ii)=dtset%kptns(ii,index)-dtset%kptns(ii,i)
231 : end do
232 : tt=0._dp
233 12 : do iij=1,3
234 9 : t2(iij)=0._dp
235 36 : do ii=1,3
236 36 : t2(iij)=t2(iij)+t1(ii)*gmat(ii,iij)
237 : end do
238 12 : tt=tt + t2(iij)**2
239 : end do
240 3 : tt=sqrt(tt)
241 : ! check for each direction in the BZ
242 12 : do j=1,3
243 : scal=0._dp
244 36 : do ii=1,3
245 36 : scal=scal+t2(ii)*gmat(j,ii)
246 : end do
247 9 : scal=abs(scal)
248 : ! Compare scal(t1,gmat) with simple product of modules -> paralel or not
249 12 : 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 1 : index=1
255 4 : do i=1,3
256 4 : index=index*nk(i)
257 : end do
258 :
259 1 : if(index.ne.nkpt) then
260 0 : write(message,'(a,2i0)')' OutwanT: Wrong assignemt of kpts', index,nkpt
261 0 : 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 3 : ABI_MALLOC(tpat,(dtset%ntypat))
267 2 : tpat(:)=zero
268 2 : do i=1,dtset%natom
269 3 : do j=1,dtset%ntypat
270 2 : if(dtset%typat(i)==j) tpat(j)=tpat(j)+1
271 : end do
272 : end do
273 2 : maxat=maxval(tpat(:))
274 4 : ABI_MALLOC(xcoord,(3,maxat,dtset%ntypat))
275 1 : index=1
276 2 : do i=1, dtset%ntypat
277 3 : do k=1,tpat(i)
278 4 : do j=1,3
279 4 : xcoord(j,k,i)=dtset%xred_orig(j,index,1)
280 : end do
281 2 : 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 3 : ABI_MALLOC(kg_tmp,(3,mpw*nkpt))
292 4 : ABI_MALLOC(iwfi,(nkpt,mpw))
293 7873 : kg_tmp(:,:)=zero
294 2461 : iwfi(:,:)=zero
295 1 : imax=npwarr(1)
296 1 : index=0
297 :
298 5 : do i=1, nkpt
299 4 : if(i>1) then
300 3 : index=index+npwarr(i-1)
301 : end if
302 1892 : do j=1, npwarr(i)
303 1891 : if(i.eq.1) then
304 459 : iwfi(i,j)=j
305 1836 : if(mkmem>0) kg_tmp(:,j)=kg(:,j)
306 : else
307 1428 : ifind=0
308 5712 : if(mkmem>0) ikg(:)=kg(:,index+j)
309 :
310 754939 : do k=1,imax
311 754939 : if(ikg(1)==kg_tmp(1,k)) then
312 57817 : if(ikg(2)==kg_tmp(2,k)) then
313 4361 : if(ikg(3)==kg_tmp(3,k)) then
314 1279 : ifind=1
315 1279 : iwfi(i,j)=k
316 : end if
317 : end if
318 : end if
319 : end do
320 :
321 1428 : if(ifind==0) then
322 149 : imax=imax+1
323 596 : kg_tmp(:,imax)=ikg(:)
324 149 : iwfi(i,j)=imax
325 : end if
326 : end if
327 : end do
328 : end do
329 1 : ngm=imax
330 :
331 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
332 : !PART ONE: writing the header
333 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
334 1 : write(unitwnt) alat
335 1 : write(unitwnt) ( drprim( i, 1 ), i = 1, 3 ) ! save A1
336 1 : write(unitwnt) ( drprim( i, 2 ), i = 1, 3 ) ! save A2
337 1 : 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 1 : write(unitwnt) dtset%ntypat
343 : !write(std_out,*) dtset%ntypat, 'NTYPAT', tpat
344 :
345 2 : do i = 1, dtset%ntypat
346 1 : write(unitwnt) tpat(i), nameat
347 2 : 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 1 : ABI_FREE(tpat)
352 1 : ABI_FREE(xcoord)
353 :
354 : !energy cut-off in Rydberg (WANT option)
355 1 : write (unitwnt) 2._dp*dtset%ecut, mband
356 : !write(std_out,*) 2._dp*dtset%ecut, mband
357 1 : 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 609 : write (unitwnt) ( kg_tmp( 1, i ), kg_tmp( 2, i ), kg_tmp( 3, i ), i = 1, ngm )
360 1 : write (unitwnt) mpw, mband, dtset%nkpt/dtset%nsppol
361 : !write(std_out,*) mpw, mband, dtset%nkpt/dtset%nsppol
362 :
363 5 : do i=1, nkpt
364 5 : write(unitwnt) (iwfi(i,j), j=1,mpw)
365 : end do
366 1 : ABI_FREE(kg_tmp)
367 :
368 : !Eigenvalues in HARTREE
369 1 : write (unitwnt) ( eig( i ), i = 1, bandtot)
370 1 : write (unitwnt) ( npwarr( ik ), ik = 1, nkpt )
371 5 : write (unitwnt) ( mband, ik = 1, nkpt )
372 4 : 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 1 : if (open_file('band.gpl',message,newunit=iunit,status='unknown') /=0) then
381 0 : ABI_ERROR(message)
382 : end if
383 :
384 1 : index=1
385 11 : do i=1,mband
386 10 : index=1
387 50 : do j=1,nkpt
388 40 : write(iunit,*) index, Ha_eV*eig(i+(j-1)*mband), eig(i+(j-1)*mband)
389 50 : index=index+1
390 : end do
391 11 : write(iunit,*)
392 : end do
393 :
394 1 : close(iunit)
395 :
396 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
397 : !PART TWO: Writing the wavefunction
398 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
399 : !not used
400 1 : ngw_=0
401 1 : ik_=0
402 1 : nk_=0
403 1 : ispin_=0
404 1 : nspin_=0
405 1 : scal_=1._dp
406 : !!!!!!!!!!!!!!!!!!!!!!!!!!
407 1 : iwf = 1
408 1 : iwf_k=1
409 3 : ABI_MALLOC(wfc,(imax))
410 :
411 : !Loop over k-pt
412 5 : do nkp=1,nkpt
413 : ! Not relevant
414 4 : write(unitwnt) twrite, ik_, section_name
415 : ! Only 'mband' is relevant here
416 4 : write(unitwnt) ngw_, mband, ik_, nk_, nk_,ispin_, nspin_, scal_
417 4 : write(unitwnt) imax
418 : ! Not relevant
419 4 : write(unitwnt) twrite
420 : ! Loop over bands
421 :
422 : ! Preparing WF
423 44 : do k=1,mband
424 40 : if(mkmem >0) then
425 24360 : wfc(:)=zero
426 : ! From cg to wf:
427 18910 : do i=iwf, iwf+npwarr(nkp)-1
428 18870 : index=i-iwf+1
429 18910 : 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 0 : message = 'Wrong mkmem in outwant'
434 0 : ABI_ERROR(message)
435 : end if
436 44 : write(unitwnt) (wfc(ig), ig=1,imax)
437 : ! End loop over bands
438 : end do
439 :
440 : ! Not relevant
441 4 : write(unitwnt) twrite
442 : ! Not relevant
443 45 : do i=1,mband
444 44 : write(unitwnt) i
445 : end do
446 :
447 : ! End loop over k-pts
448 : end do
449 :
450 1 : ABI_FREE(iwfi)
451 1 : ABI_FREE(wfc)
452 :
453 1 : call wrtout(std_out,'Closing file','COLL')
454 1 : close(unit=unitwnt)
455 :
456 1 : end subroutine outwant
457 : !!***
458 :
459 : end module m_outwant
460 : !!***
|
