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