Line data Source code
1 : !!****p* ABINIT/macroave
2 : !! NAME
3 : !! macroave
4 : !!
5 : !! FUNCTION
6 : !! **********************************************************************
7 : !! The MACROAVE program implements the macroscopic average technique,
8 : !! introduced by A. Baldereschi and coworkers
9 : !! (A. Baldereschi, S. Baroni, and R. Resta, Phys. Rev. Lett. 61, 734 (1988) [[cite:Baldereschi1988]]).
10 : !! This is an extremely powerful method that relates
11 : !! microscopic quantities, typical outputs of first-principles codes,
12 : !! with macroscopic magnitudes, needed to perform electrostatic analysis.
13 : !! Within this methodology, we will be able to wash out all the
14 : !! wiggles of the rapidly-varying functions of position (resembling
15 : !! the underlying atomic structure) of the microscopic quantities,
16 : !! blowing up only the macroscopic features.
17 : !! It can be used to compute band offsets, work functions, effective
18 : !! charges, and high frequency dielectric constants, among others
19 : !! interesting physical properties.
20 : !! Ref: L. Colombo, R. Resta and S. Baroni Phys Rev B 44, 5572 (1991) [[cite:Colombo1991]].
21 : !! Coded by P. Ordejon and J. Junquera, April 1999.
22 : !! Modified by J. Junquera, November 2001.
23 : !! **********************************************************************
24 : !!
25 : !! COPYRIGHT
26 : !! Copyright (C) 1999-2008 (P. Ordejon, J. Junquera, J. Soler, A. Garcia)
27 : !! This file is distributed under the terms of the
28 : !! GNU General Public License, see ~abinit/COPYING
29 : !! or http://www.gnu.org/copyleft/gpl.txt .
30 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
31 : !!
32 : !! INPUTS
33 : !! (main routine)
34 : !!
35 : !! OUTPUT
36 : !! (main routine)
37 : !!
38 : !! SOURCE
39 :
40 : #if defined HAVE_CONFIG_H
41 : #include "config.h"
42 : #endif
43 :
44 : #include "abi_common.h"
45 :
46 8 : program macroave
47 :
48 8 : use defs_basis
49 : use m_xmpi
50 : use m_abicore
51 : use m_errors
52 : use m_nctk
53 : use netcdf
54 : use m_hdr
55 : use m_macroave
56 :
57 : use m_fstrings, only : sjoin, strcat, endswith
58 : use m_io_tools, only : file_exists, open_file
59 :
60 : implicit none
61 :
62 : !Arguments -----------------------------------
63 :
64 : !Local variables-------------------------------
65 : !no_abirules
66 : ! --------- PARAMETERS -------------------------------------------------
67 : ! INTEGER NP : Parameter needed to define maximum number of points
68 : ! for FFT grid.
69 : ! Number of points = (2**NP)
70 : ! INTEGER N : Maximum number of complex data point for FFT.
71 : ! MUST be a power of 2
72 : ! REAL*8 HARTREE: Conversion factor from Hartrees to Rydbergs
73 : ! 1 hartree = 2 Ry
74 : ! REAL*8 RYDBERG: Conversion factor from Rydbergs to eV
75 : ! ----------------------------------------------------------------------
76 : integer, parameter :: np=12
77 : integer, parameter :: n=2**np
78 : real(dp), parameter :: hartree=two
79 : real(dp), parameter :: rydberg=Ha_eV *half
80 : ! --------- VARIABLES --------------------------------------------------
81 : integer :: i,ii,ij,ip,is,j,iomode
82 : integer :: nconv,npoints,npt,nsm,nspin,nz,nspden
83 : integer :: unit1,unit2,unit3,unit4,varid
84 : integer :: mesh(3)
85 : character(len=10) :: code,interp
86 : character(len=15) :: SNAME,inpdata
87 : character(len=fnlen) :: fnamerho,fnamedelv,fnameplave
88 : character(len=500) :: msg
89 : character(len=nctk_slen) :: varname
90 : logical :: siesta,abinit,potential,charge,totalcharge
91 : logical :: found,linear,splin
92 8 : real,allocatable :: rhos(:,:)
93 8 : real(dp),allocatable :: rho(:,:)
94 : real(dp) :: cell(3,3),dcell(3,3)
95 : real(dp) :: l,sur,ds,length,convfac,qren,qtot,lav1,lav2,vol
96 8 : real(dp),allocatable :: z(:),rhoz(:),d2rhoz(:),drhodz(:)
97 : real(dp) :: data(2*n),th(2*n),re(n),im(n)!,v(2*n)
98 : real(dp) :: x,delta,yp1,ypn,phi
99 : complex(dp) :: a,b,c
100 : ! ABINIT variables
101 : integer :: fform
102 : integer :: comm,nproc,my_rank
103 8 : type(hdr_type) :: hdr
104 : ! end ABINIT variables
105 :
106 : !************************************************************************
107 : !CHARACTER CODE : First principles-code used to generate the
108 : !electrostatic potential at the points of a grid
109 : !in real space. It is read from file 'macroave.in'
110 : !CHARACTER SNAME : System Label
111 : !(If code = ABINIT, then SNAME = FNAMERHO)
112 : !CHARACTER INPDATA : Calculate the band offset from the charge
113 : !density or from the electrostatic potential?
114 : !CHARACTER FNAMERHO : Name of the file where the electrostatic
115 : !potential at the mesh is stored
116 : !CHARACTER FNAMEPLAVE : Name of the file where the planar average of the
117 : !electronic potential will be stored
118 : !CHARACTER FNAMEDELV : Name of the file where the profile
119 : !of the electrostatic potential will be stored
120 : !LOGICAL SIESTA : Have you used SIESTA to get the electrostatic potential
121 : !or the charge density?
122 : !LOGICAL ABINIT : Have you used ABINIT to get the electrostatic potential
123 : !or the charge density?
124 : !LOGICAL LINEAR : Linear interpolation to get the charge
125 : !density/potential in the FFT grid.
126 : !LOGICAL SPLIN : Cubic spline interpolation to get the charge
127 : !density/potential in the FFT grid.
128 : !LOGICAL POTENTIAL : We are going to compute the band offset from
129 : !the electrostatic potential
130 : !LOGICAL CHARGE : We are going to compute the band offset from
131 : !the charge density
132 : !LOGICAL FOUND : Were data found? (only when task in iorho ='read')
133 : !INTEGER NATOMS : Number of atoms in unit cell
134 : !INTEGER NSPIN : Number of spin polarizations (1 or 2)
135 : !INTEGER MESH(3) : Number of mesh divisions of each lattice vectors,
136 : !INCLUDING subgrid
137 : !INTEGER NSM : Number of sub-mesh points per mesh point
138 : !(not used in this version)
139 : !INTEGER NPOINTS : Number of mesh subdivisions in the normal direction
140 : !to the interface
141 : !INTEGER NCONV : Number of convolutions required to calculate the
142 : !macroscopic average
143 : !INTEGER NPT : Total number of mesh points (included subpoints)
144 : !REAL*8 CELL(3,3) : Unit cell lattice vectors (a.u.) CELL(IXYZ,IVECT)
145 : !REAL*8 DS : Differential area per point of the mesh
146 : !REAL*8 SUR : Area of a plane parallel to the interface
147 : !REAL*8 LENGTH : Distance between two planes parallel to the interface
148 : !REAL*8 L : Length of the cell in the direction nomal
149 : !to the interface (a.u.)
150 : !REAL*8 LAV1 : Linear period of the electrostatic potential
151 : !in the bulklike region for the first material
152 : !REAL*8 LAV2 : Linear period of the electrostatic potential
153 : !in the bulklike region for the second material
154 : !REAL*8 CONVFAC : Conversion factor for the output units
155 : !REAL*8 QTOT : Total electronic charge in the unit cell
156 : !REAL*8 QREN : Total electronic charge calculated from
157 : !the input data file
158 : !REAL*4 Z(MESH(3)) : Z-coordinate of the planes where the elec. density
159 : !is averaged
160 : !REAL*4 RHO : Electron density
161 : !Notice single precision in this version
162 : !REAL*8 RHOZ : Planar average of the electron density
163 : !REAL*8 DATA(2*N) : Fourier coefficients of the planar average density
164 : !REAL*8 TH(2*N) : Fourier coefficients of the step functions
165 : !REAL*8 V(2*N) : Fourier coefficients of the potential
166 : !REAL*4 VREEC(N) : Real part of the electronic potential in real space
167 : !REAL*4 VIMEC(N) : Imag. part of the electronic potential in real space
168 : !REAL*4 VRENC(N) : Real part of the nuclear potential in real space
169 : !REAL*4 VIMNC(N) : Imag. part of the nuclear potential in real space
170 : !*********************************************************************
171 :
172 : !Change communicator for I/O (mandatory!)
173 8 : call abi_io_redirect(new_io_comm=xmpi_world)
174 :
175 : !Initialize MPI
176 8 : call xmpi_init()
177 8 : comm = xmpi_world
178 8 : nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
179 :
180 : !Initialize memory profiling if it is activated
181 : !if a full abimem.mocc report is desired, set the argument of abimem_init to "2" instead of "0"
182 : !note that abimem.mocc files can easily be multiple GB in size so don't use this option normally
183 : #ifdef HAVE_MEM_PROFILING
184 : call abimem_init(0)
185 : #endif
186 :
187 : !Reading input data from a file ---------------------------------------
188 8 : if (open_file('macroave.in', msg, newunit=UNIT1, STATUS='OLD') /= 0) then
189 0 : ABI_ERROR(msg)
190 : end if
191 :
192 8 : READ(UNIT1,'(A)')CODE
193 8 : READ(UNIT1,'(A)')INPDATA
194 8 : READ(UNIT1,'(A)')SNAME
195 8 : READ(UNIT1,*)NCONV
196 8 : READ(UNIT1,*)LAV1
197 8 : READ(UNIT1,*)LAV2
198 8 : READ(UNIT1,*)QTOT
199 8 : READ(UNIT1,*)INTERP
200 8 : close(UNIT1)
201 :
202 : !Which code has been used to get the electrostatic potential? ---------
203 8 : if ( CODE == 'siesta' .OR. CODE == 'SIESTA' .OR.&
204 : & CODE == 'Siesta' ) then
205 0 : SIESTA = .TRUE.
206 0 : ABINIT = .FALSE.
207 : else if ( CODE == 'abinit' .OR. CODE == 'ab-init' .OR.&
208 : & CODE == 'ABINIT' .OR. CODE == 'AB-INIT' .OR.&
209 8 : & CODE == 'Abinit' .OR. CODE == 'Ab-init') then
210 8 : SIESTA = .FALSE.
211 8 : ABINIT = .TRUE.
212 : else
213 0 : ABI_ERROR(sjoin('macroave: Unknown code: ', CODE))
214 : end if
215 :
216 : !Are we going to compute the band offset from the charge density or
217 : !from the electrostatic potential? ------------------------------------
218 8 : if ( INPDATA == 'potential' .OR. INPDATA == 'POTENTIAL' .OR.&
219 : & INPDATA == 'Potential' ) then
220 4 : POTENTIAL = .TRUE.
221 4 : CHARGE = .FALSE.
222 4 : TOTALCHARGE = .FALSE.
223 4 : else if ( INPDATA == 'charge' .OR. INPDATA == 'CHARGE' .OR.&
224 : & INPDATA == 'Charge' ) then
225 4 : POTENTIAL = .FALSE.
226 4 : CHARGE = .TRUE.
227 4 : TOTALCHARGE = .FALSE.
228 : else if ( INPDATA == 'totalcharge' .OR. &
229 : & INPDATA == 'Totalcharge' .OR.&
230 0 : & INPDATA == 'TotalCharge' .OR.&
231 : & INPDATA == 'TOTALCHARGE' ) then
232 0 : POTENTIAL = .FALSE.
233 0 : CHARGE = .FALSE.
234 0 : TOTALCHARGE = .TRUE.
235 : else
236 0 : ABI_ERROR(sjoin('macroave: Unknown input data ', INPDATA))
237 : end if
238 :
239 : !What kind of interpolation will we use to get the charge density/
240 : !potential in a FFT grid? ---------------------------------------------
241 8 : if ( INTERP == 'linear' .OR. INTERP == 'Linear' .OR.&
242 : & INTERP == 'LINEAR' ) then
243 2 : LINEAR = .TRUE.
244 2 : SPLIN = .FALSE.
245 6 : else if ( INTERP == 'spline' .OR. INTERP == 'Spline' .OR.&
246 : & INTERP == 'SPLINE' ) then
247 6 : LINEAR = .FALSE.
248 6 : SPLIN = .TRUE.
249 : end if
250 :
251 : !Reading charge density from a file -----------------------------------
252 8 : if ( SIESTA ) then
253 0 : if (POTENTIAL) then
254 0 : FNAMERHO = strcat(trim(SNAME),'.VH')
255 0 : elseif (CHARGE) then
256 0 : FNAMERHO = strcat(trim(SNAME),'.RHO')
257 0 : elseif (TOTALCHARGE) then
258 0 : FNAMERHO = strcat(trim(SNAME),'.TOCH')
259 : end if
260 8 : else if ( ABINIT ) then
261 8 : FNAMERHO = trim(SNAME)
262 : end if
263 :
264 8 : if (SIESTA) then
265 0 : NSM = 1
266 0 : NPT = 0
267 0 : NSPIN = 0
268 : CALL IORHO( 'READ', trim(FNAMERHO), DCELL, MESH, NSM, NPT, NSPIN,&
269 0 : & RHOS, FOUND )
270 0 : if (FOUND) then
271 0 : ABI_MALLOC( RHOS,(NPT,NSPIN))
272 0 : ABI_MALLOC( RHO,(NPT,NSPIN))
273 : CALL IORHO( 'READ', trim(FNAMERHO), DCELL, MESH, NSM, NPT, NSPIN,&
274 0 : & RHOS, FOUND )
275 0 : do I = 1, 3
276 0 : do J = 1, 3
277 0 : CELL(J,I) = DCELL(J,I)
278 : end do
279 : end do
280 : ! Transform the density or the potential read from SIESTA
281 : ! from a single precision variable to a double precision variable
282 0 : do IS = 1, NSPIN
283 0 : do IP = 1, NPT
284 0 : RHO(IP,IS) = RHOS(IP,IS) * 1.0D0
285 : end do
286 : end do
287 :
288 : else
289 0 : ABI_ERROR(sjoin('macroave: file not found: ', FNAMERHO))
290 : end if
291 :
292 8 : else if (ABINIT) then
293 :
294 8 : if (nctk_try_fort_or_ncfile(FNAMERHO, msg) /= 0) then
295 0 : ABI_ERROR(msg)
296 : end if
297 8 : iomode = IO_MODE_FORTRAN; if (endswith(FNAMERHO, ".nc")) iomode = IO_MODE_ETSF
298 :
299 8 : if (iomode == IO_MODE_FORTRAN) then
300 0 : if (open_file(FNAMERHO, msg, newunit=unit2, form="unformatted", status="old") /= 0) then
301 0 : ABI_ERROR(msg)
302 : end if
303 0 : call hdr%fort_read(unit2, fform)
304 0 : ABI_CHECK(FFORM /= 0, "fform == 0")
305 : else
306 8 : NCF_CHECK(nctk_open_read(unit2, fnamerho, xmpi_comm_self))
307 8 : call hdr%ncread(unit2, fform)
308 : end if
309 :
310 : ! For debugging
311 : ! call hdr%echo(fform,4,std_out)
312 :
313 32 : do I = 1, 3
314 24 : MESH(I) = HDR%NGFFT(I)
315 104 : do J = 1, 3
316 96 : CELL(J,I) = HDR%RPRIMD(J,I)
317 : end do
318 : end do
319 8 : NSPIN = HDR%NSPPOL
320 8 : nspden = hdr%nspden
321 8 : ABI_CHECK(hdr%nspinor == 1, "nspinor == 2 not coded")
322 8 : call hdr%free()
323 :
324 8 : NPT = MESH(1) * MESH(2) * MESH(3)
325 8 : ABI_MALLOC( RHO,(NPT,NSPIN))
326 :
327 8 : if (iomode == IO_MODE_FORTRAN) then
328 0 : do IS = 1, NSPIN
329 0 : READ(UNIT2) (RHO(IP,IS),IP=1,NPT)
330 : end do
331 0 : close(UNIT2)
332 : else
333 8 : varname = varname_from_fname(fnamerho)
334 8 : NCF_CHECK(nf90_inq_varid(unit2, varname, varid))
335 : ! [cplex, n1, n2, n3, nspden]
336 48 : NCF_CHECK(nf90_get_var(unit2, varid, rho, start=[1,1,1,1,1], count=[1, mesh(1), mesh(2), mesh(3), nspden]))
337 : end if
338 :
339 : ! Units for the potential in Ab-init are in Hartrees,
340 : ! so we transform them into Ry. No transformation is
341 : ! needed for the charge density
342 : ! (it is directly read in electrons/bohr**3).
343 :
344 8 : if (POTENTIAL) then
345 8 : do IS = 1, NSPIN
346 28808 : do IP = 1, NPT
347 28804 : RHO(IP,IS) = RHO(IP,IS) * HARTREE
348 : end do
349 : end do
350 : end if
351 : end if
352 :
353 : !Initialize some variables (we suppose cells with a c axis orthogonal to a and b) -------------
354 :
355 8 : L = CELL(3,3)
356 8 : SUR = SURPLA( CELL ) ! surface of unit cell in xy plane, perpendicular to z
357 8 : VOL = VOLCEL( CELL )
358 8 : DS = SUR/( MESH(1) * MESH(2) ) ! this seems adapted to arbitrary in-plane cells
359 8 : LENGTH = L/DBLE(N)
360 8 : NPOINTS = MESH(3)
361 :
362 8 : ABI_MALLOC(Z,(NPOINTS+1))
363 8 : ABI_MALLOC(RHOZ,(NPOINTS+1))
364 8 : ABI_MALLOC(D2RHOZ,(NPOINTS+1))
365 8 : ABI_MALLOC(DRHODZ,(N))
366 :
367 480 : RHOZ(1:NPOINTS+1) = 0.D0
368 480 : D2RHOZ(1:NPOINTS+1) = 0.D0
369 32776 : DRHODZ(1:N) = 0.D0
370 :
371 8 : if (POTENTIAL) then
372 4 : CONVFAC = RYDBERG
373 4 : else if (CHARGE) then
374 4 : CONVFAC = 1.0D0
375 0 : else if (TOTALCHARGE) then
376 0 : CONVFAC = 1.0D0
377 : end if
378 :
379 :
380 : !Loop over all points and calculate the planar average ----------------
381 : !Warning: The planar average is only done for the first component of
382 : !RHO. Spin polarization is not implemented yet ---------------
383 51400 : do IP = 1, NPT
384 51392 : NZ = (IP-1) / (MESH(1)*MESH(2)) + 1
385 51400 : RHOZ(NZ) = RHOZ(NZ) + RHO(IP,1)*DS
386 : end do
387 :
388 472 : do IP = 1, NPOINTS
389 472 : RHOZ(IP) = RHOZ(IP) / SUR
390 : end do
391 :
392 472 : do IP = 1, NPOINTS
393 472 : Z(IP) = (IP-1)*CELL(3,3)/DBLE(NPOINTS)
394 : end do
395 :
396 : !Calculate electrostatic potential or electronic charge density -------
397 : !in fft grid, interpolating the planar average calculated before ------
398 8 : if (SPLIN) then
399 6 : Z(NPOINTS+1) = L
400 6 : RHOZ(NPOINTS+1) = RHOZ(1)
401 6 : DELTA = L/DBLE(NPOINTS)
402 6 : YP1 = ( RHOZ(2) - RHOZ(NPOINTS) ) / (2.0D0*DELTA)
403 6 : YPN = YP1
404 6 : CALL MACROAV_SPLINE(DELTA, RHOZ, NPOINTS+1, YP1, YPN, D2RHOZ)
405 6 : I = 0
406 6 : do II = 1, 2*N-1, 2
407 24576 : I = I + 1
408 24576 : X = (I-1)*L/DBLE(N)
409 : CALL MACROAV_SPLINT( DELTA, RHOZ, D2RHOZ, NPOINTS+1, X, DATA(II), &
410 24576 : & DRHODZ(I) )
411 24576 : DATA(II+1) = 0.D0
412 : end do
413 2 : else if (LINEAR) then
414 2 : I = 0
415 8194 : do II = 1,2*N-1,2
416 8192 : I = I + 1
417 8192 : X = (I-1)*L/DBLE(N)
418 307122 : do IJ = 1, NPOINTS
419 307010 : if (X == Z(IJ)) then
420 14 : DATA(II) = RHOZ(IJ)
421 14 : DATA(II+1) = 0.D0
422 14 : GOTO 20
423 : end if
424 307108 : if (Z(IJ) > X) then
425 : DATA(II) = RHOZ(IJ-1) +&
426 : & (X-Z(IJ-1))*(RHOZ(IJ)-RHOZ(IJ-1))/&
427 8066 : & (Z(IJ)-Z(IJ-1))
428 8066 : DATA(II+1) = 0.D0
429 8066 : GOTO 20
430 : end if
431 : end do
432 : DATA(II)=RHOZ(NPOINTS) +&
433 : & (X-Z(NPOINTS))*(RHOZ(1)-RHOZ(NPOINTS))/&
434 112 : & (Z(NPOINTS)-Z(NPOINTS-1))
435 112 : DATA(II+1) = 0.D0
436 2 : 20 CONTINUE
437 : end do
438 : end if
439 :
440 : !Renormalize the charge density ---------------------------------------
441 8 : if (CHARGE .OR. TOTALCHARGE) then
442 4 : QREN = 0.D0
443 4 : do IP = 1, 2*N-1, 2
444 16384 : QREN = QREN + DATA(IP)*LENGTH*SUR
445 : end do
446 16388 : do IP = 1, 2*N-1, 2
447 16388 : if (CHARGE) then
448 16384 : DATA(IP) = DATA(IP) * QTOT/QREN
449 0 : elseif(TOTALCHARGE) then
450 0 : DATA(IP) = DATA(IP) - QREN/VOL
451 : end if
452 : end do
453 4 : QREN = 0.D0
454 4 : do IP = 1, 2*N-1, 2
455 16384 : QREN = QREN + DATA(IP)*LENGTH*SUR
456 : end do
457 : ! For debugging
458 : ! write(std_out,*)' QREN = ', QREN
459 : end if
460 : !...
461 :
462 : !Print planar average of the electrostatic potential or ---------------
463 : !the electronic charge density ----------------------------------------
464 8 : FNAMEPLAVE = strcat(SNAME,'.PAV')
465 8 : if (open_file(FNAMEPLAVE,msg,newunit=UNIT3,STATUS='UNKNOWN') /= 0) then
466 0 : ABI_ERROR(msg)
467 : end if
468 8 : I = 0
469 32776 : do II = 1, 2*N-1, 2
470 32768 : I = I+1
471 32768 : X=(I-1)*L/DBLE(N)
472 : ! WRITE(UNIT3,'(3F20.12)')X,
473 : ! . DATA(II)*CONVFAC,DATA(II+1)*CONVFAC
474 32768 : WRITE(UNIT3,'(2F20.12)')X,&
475 65544 : & DATA(II)*CONVFAC
476 : end do
477 8 : close(UNIT3)
478 : !...
479 :
480 :
481 : !Calculate Fourier transform of the electrostatic potential or
482 : !the electronic density -----------------------------------------------
483 8 : CALL FOUR1(DATA,N,1)
484 : !...
485 :
486 : !Calculate macroscopic average of the electrostatic potential or the
487 : !electronic charge density taking the convolution with two step functions.
488 : !In Fourier space, it is a product of the Fourier transform components -
489 : !The decompositions in the sum over II is due to the special way in which
490 : !the data are stored in subroutine four1( see Fig. 12.2.2, in
491 : !'Numerical Recipes, The Art of Scientific Computing'
492 : !by W.H. Press, S.A. Teukolsky, W.T. Veterling and B.P. Flannery,
493 : !Cambridge U.P. 1987 and 1992.
494 :
495 :
496 8 : CALL THETAFT(N,L,LAV1,TH)
497 :
498 8 : do II = 1, N+1, 2
499 16392 : A = DATA(II)*(1.D0,0.D0) + DATA(II+1)*(0.D0,1.D0)
500 16392 : B = TH(II)*(1.D0,0.D0) + TH(II+1)*(0.D0,1.D0)
501 16392 : C = A*B
502 16392 : DATA(II) = REAL(C)*L/DBLE(N)
503 16392 : DATA(II+1) = AIMAG(C)*L/DBLE(N)
504 : end do
505 :
506 8 : do II = N+3, 2*N-1, 2
507 16376 : A = DATA(II)*(1.D0,0.D0) + DATA(II+1)*(0.D0,1.D0)
508 16376 : B = TH(II)*(1.D0,0.D0) + TH(II+1)*(0.D0,1.D0)
509 16376 : C = A*B
510 16376 : DATA(II) = REAL(C)*L/DBLE(N)
511 16376 : DATA(II+1) = AIMAG(C)*L/DBLE(N)
512 : end do
513 :
514 :
515 8 : if (NCONV == 2) then
516 2 : CALL THETAFT(N,L,LAV2,TH)
517 :
518 2 : do II = 1, N+1, 2
519 4098 : A = DATA(II)*(1.D0,0.D0) + DATA(II+1)*(0.D0,1.D0)
520 4098 : B = TH(II)*(1.D0,0.D0) + TH(II+1)*(0.D0,1.D0)
521 4098 : C = A*B
522 4098 : DATA(II) = REAL(C)*L/DBLE(N)
523 4098 : DATA(II+1) = AIMAG(C)*L/DBLE(N)
524 : ! if ( POISON ) then
525 : ! IG = (II-1) / 2
526 : ! GSQ= (2.D0*PI*IG/L)**2
527 : ! if(GSQ > 0.D0) then
528 : ! V(II) = DATA(II) * (4.D0*PI/GSQ) * HARTREE * RYDBERG
529 : ! V(II+1) = DATA(II+1) * (4.D0*PI/GSQ) * HARTREE * RYDBERG
530 : ! else
531 : ! V(II) = 0.D0
532 : ! V(II+1) = 0.D0
533 : ! endif
534 : ! endif
535 : end do
536 :
537 2 : do II = N+3, 2*N-1, 2
538 4094 : A = DATA(II)*(1.D0,0.D0) + DATA(II+1)*(0.D0,1.D0)
539 4094 : B = TH(II)*(1.D0,0.D0) + TH(II+1)*(0.D0,1.D0)
540 4094 : C = A*B
541 4094 : DATA(II) = REAL(C)*L/DBLE(N)
542 4094 : DATA(II+1) = AIMAG(C)*L/DBLE(N)
543 : ! if ( POISON ) then
544 : ! IG = (-2*N+II-1) / 2
545 : ! GSQ= (2.D0*PI*IG/L)**2
546 : ! if(GSQ > 0.D0) then
547 : ! V(II) = DATA(II) * (4.D0*PI/GSQ) * HARTREE * RYDBERG
548 : ! V(II+1) = DATA(II+1) * (4.D0*PI/GSQ) * HARTREE * RYDBERG
549 : ! else
550 : ! V(II) = 0.D0
551 : ! V(II+1) = 0.D0
552 : ! endif
553 : ! endif
554 : end do
555 :
556 : end if
557 : !...
558 :
559 : !Transform average electronic density and potential to real space -----
560 : !The decompositions in the sum over J is due to the special way in which
561 : !the data are stored in subroutine four1( see Fig. 12.2.2, in
562 : !'Numerical Recipes, The Art of Scientific Computing'
563 : !by W.H. Press, S.A. Teukolsky, W.T. Veterling and B.P. Flannery,
564 : !Cambridge U.P. 1987 and 1992.
565 :
566 32776 : do II = 1, N
567 32768 : RE(II) = 0.D0
568 32768 : IM(II) = 0.D0
569 : ! if ( POISON ) then
570 : ! VREEC(II) = 0.D0
571 : ! VIMEC(II) = 0.D0
572 : ! endif
573 32768 : do J = 1, N+1, 2
574 67141632 : PHI = -2.D0 * PI * (II-1) * ( (J-1)/2 ) / DBLE(N)
575 : RE(II)=RE(II)+(1.D0/DBLE(N))*(DATA(J)*COS(PHI)&
576 67141632 : & -DATA(J+1)*SIN(PHI))
577 : IM(II)=IM(II)+(1.D0/DBLE(N))*(DATA(J)*SIN(PHI)&
578 67141632 : & +DATA(J+1)*COS(PHI))
579 : ! if ( POISON ) then
580 : ! VREEC(II)=VREEC(II)+(1.D0/DBLE(N))*(V(J)*COS(PHI)
581 : ! . -V(J+1)*SIN(PHI))
582 : ! VIMEC(II)=VIMEC(II)+(1.D0/DBLE(N))*(V(J)*SIN(PHI)
583 : ! . +V(J+1)*COS(PHI))
584 : ! endif
585 : end do
586 :
587 32776 : do J = N+3, 2*N-1, 2
588 67076096 : PHI = -2.0D0 * PI * (II-1) * ((-2*N+J-1)/2) / DBLE(N)
589 : RE(II)=RE(II)+(1.D0/DBLE(N))*(DATA(J)*COS(PHI)&
590 67076096 : & -DATA(J+1)*SIN(PHI))
591 : IM(II)=IM(II)+(1.D0/DBLE(N))*(DATA(J)*SIN(PHI)&
592 67076096 : & +DATA(J+1)*COS(PHI))
593 : ! if ( POISON ) then
594 : ! VREEC(II)=VREEC(II)+(1.D0/DBLE(N))*(V(J)*COS(PHI)
595 : ! . -V(J+1)*SIN(PHI))
596 : ! VIMEC(II)=VIMEC(II)+(1.D0/DBLE(N))*(V(J)*SIN(PHI)
597 : ! . +V(J+1)*COS(PHI))
598 : ! endif
599 : end do
600 : end do
601 : !...
602 :
603 : !Print averaged electronic charge density and potential ---------------
604 8 : FNAMEDELV = strcat( trim(SNAME),'.MAV')
605 8 : if (open_file(FNAMEDELV, msg, newunit=UNIT4, STATUS='UNKNOWN') /= 0) then
606 0 : ABI_ERROR(msg)
607 : end if
608 32776 : do I = 1, N
609 32768 : X=(I-1)*L/DBLE(N)
610 : ! WRITE(UNIT4,'(3F20.5)')X,
611 : ! . RE(I)*CONVFAC ,IM(I)*CONVFAC
612 32768 : WRITE(UNIT4,'(2F20.12)')X,&
613 65544 : & RE(I)*CONVFAC
614 : end do
615 8 : close(UNIT4)
616 : !...
617 :
618 : !Print electrostatic potential ----------------------------------------
619 : !if (POISON) then
620 : !FNAMEVEC = strcat( SNAME,'.VEC')
621 : !if (open_file(FNAMEVEC, msg newunit=UNIT5, STATUS='UNKNOWN') /= 0) then
622 : ! ABI_ERROR(msg)
623 : !end if
624 : !do I = 1, N
625 : !X=(I-1)*L/DBLE(N)
626 : !WRITE(UNIT5,'(3F20.12)')X, VREEC(I), VIMEC(I)
627 : !enddo
628 : !close(unit5)
629 : !endif
630 : !...
631 :
632 8 : ABI_FREE(Z)
633 8 : ABI_FREE(RHOZ)
634 8 : ABI_FREE(DRHODZ)
635 8 : ABI_FREE(D2RHOZ)
636 8 : if (allocated(rho)) then
637 8 : ABI_FREE(rho)
638 : end if
639 :
640 : !Write information on file about the memory before ending mpi module, if memory profiling is enabled
641 8 : call abinit_doctor("__macroave")
642 :
643 8 : call xmpi_end()
644 :
645 0 : end program macroave
646 : !!***
|
