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