Branch data 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 : : !!***
|
