1!************************************************************************* 2 3Program EPM2BGW 4!************************************************************************* 5! EPM calculates the band structure along various directions in 6! k-space using a plane-wave basis and a fixed effective potential. 7! Several choices of empirical pseudopotentials are provided. 8!************************************************************************* 9! Written by William Mattson and Richard M. Martin of the University 10! of Illinois, based upon program written by K. Glassford and I. Souza. 11! Modified by G. Samsonidze for use with BerkeleyGW (August 2008). 12!************************************************************************* 13! Note: atomic units (a.u.) are used throughout the program. 14!************************************************************************* 15 16#include "f_defs.h" 17 18 use sysParams, only : double, zero, one 19 use TagHandlerMod, only : tagHandlerT, TagHandlerInit, TagHandlerDestroy, & 20 & FindTag 21 use StructureMod, only : StructureT, StructureInit, StructureDestroy 22 use kPointsMod, only : kPointsT, KPointsInit, KPointsDestroy 23 use graphMod, only : graphT, PlotBands 24 use hamiltonianMod, only : hamiltonianArrayT, HInit, HDestroy, HDiagonalize, & 25 & HPrint, HCheck 26 use eigenStatesMod, only : eigenStatesT, EigenStatesInit, EigenStatesDestroy, & 27 & PrintEigenvalues 28 use pwHamMod, only : hamInfoT, HamInfoInit, HamInfoDestroy, & 29 & HGenerate, FindHMaxSize, GridInit, & 30 & PotentialRGenerate, KineticGenerate, & 31 & CalculateDensity, GridDestroy 32 use cgParamsMod, only : cgParamsT, cgParamsInit 33 use typeMod, only : GridT 34 use ConjGradMod, only : cgEigenSystem,cgFFTDiagonalise 35 use densityArrayMod, only : densityArrayT, DensityArrayInit, & 36 & DensityArrayDestroy 37 38 use nrtype_m 39 use message_m 40 use wfn_rho_vxc_io_m, only : write_binary_header, write_binary_gvectors, & 41 & write_binary_real_data, write_binary_complex_data, write_matrix_elements 42 use fftw_m, only : check_FFT_size 43 use check_inversion_m, only : check_inversion 44 use symmetries_m, only : get_symmetries 45 46 implicit none 47 48 type(tagHandlerT), pointer :: tagHandler 49 type(hamInfoT), pointer :: hamInfo 50 type(StructureT), pointer :: structure1 51 type(kPointsT), pointer :: kPoints 52 type(eigenStatesT), pointer :: eigenStates 53 type(hamiltonianArrayT), pointer :: hamiltonian 54 type(graphT), pointer :: graphinfo 55 type(cgParamsT) :: cgParams 56 type(GridT), pointer :: Grid 57 type(densityArrayT), pointer :: densityArray 58 59 integer, allocatable :: gVectorsK(:,:) 60 real(double), allocatable :: eigenVectorsReal(:,:,:,:) 61 complex(double), allocatable :: eigenVectorsComplex(:,:,:,:) 62 real(double), allocatable :: dataReal(:,:) 63 complex(double), allocatable :: dataComplex(:,:), vxc_mtxel(:,:) 64 65 integer, pointer :: atyp(:) 66 real(double), pointer :: apos(:,:) 67 integer, pointer :: ngk(:) 68 real(double), pointer :: kw(:) 69 real(double), pointer :: kpt(:,:) 70 integer, pointer :: ifmin(:,:) 71 integer, pointer :: ifmax(:,:) 72 real(double), pointer :: energies(:,:,:) 73 real(double), pointer :: occupations(:,:,:) 74 logical :: wfng_flag, rhog_flag, vxcg_flag, vxc_flag, v_of_q, disable_symmetries 75 integer :: real_or_complex, ReportNo, error, i, j, k, & 76 & ik, ispin, ib, ig, cell_symmetry, nat, nsym, nspin, & 77 & vxc_diag_nmin, vxc_diag_nmax, vxc_offdiag_nmin, & 78 & vxc_offdiag_nmax, ndiag, noffdiag, FFTgrid(3), & 79 & rotation(3, 3, 48), spin_index(2), idiag, ioffdiag, spacegroup, & 80 & identity(3,3) 81 integer, allocatable :: diag(:), offdiag(:) 82 real(double) :: celvol, recvol, al, bl, ecutwfn, ecutrho, & 83 & gcutm, abstol, temperature, efermi, & 84 & a(3, 3), b(3, 3), adot(3, 3), bdot(3, 3), translation(3, 48), ff(6) 85 character(len=3) :: sheader 86 character(len=256) :: wfng_file, rhog_file, vxcg_file, vxc_file, & 87 & filename, tmpstr 88 character*21 :: symbol 89 90 integer, parameter :: Nfac = 3 91 92! Initialization 93 94 nullify(tagHandler) 95 nullify(hamInfo) 96 nullify(structure1) 97 nullify(kPoints) 98 nullify(eigenStates) 99 nullify(hamiltonian) 100 nullify(graphinfo) 101 nullify(Grid) 102 nullify(densityArray) 103 104 nullify(atyp) 105 nullify(apos) 106 nullify(ngk) 107 nullify(kw) 108 nullify(kpt) 109 nullify(ifmin) 110 nullify(ifmax) 111 nullify(energies) 112 nullify(occupations) 113 114 ReportNo = 110 115 116 filename = 'input.tmp' 117 call open_file(unit = 10, file = trim(filename), status = 'replace') 118 i = 0 119 do while (i .eq. 0) 120 read(unit = 5, fmt = 999, iostat = i) tmpstr 121 write(unit = 10, fmt = 999) tmpstr 122 enddo 123 call close_file(unit = 10) 124 125 ! Necessary if using conjugate gradient method 126 call Random_Seed 127 128 ! Initialize "tag handler" to read the input from file 'filename'. 129 call TagHandlerInit( tagHandler, 10, filename) 130 131 ! Input parameter real_or_complex for BerkeleyGW files 132 call FindTag( tagHandler, 'real_or_complex', error) 133 if(error .eq. 0) then 134 read(unit = tagHandler%fileno, fmt = *, iostat = error) real_or_complex 135 if(error .eq. 0) then 136 if(real_or_complex .ne. 1 .and. real_or_complex .ne. 2) error = 1 137 endif 138 endif 139 if(error .ne. 0) real_or_complex = 2 140 write(*,'(" real_or_complex ",i1)')real_or_complex 141 142 ! Input parameter disable_symmetries for BerkeleyGW wavefunction file 143 call FindTag( tagHandler, 'disable_symmetries', error) 144 if(error .eq. 0) then 145 read(unit = tagHandler%fileno, fmt = *, iostat = error) disable_symmetries 146 endif 147 if(error .ne. 0) disable_symmetries = .false. 148 write(*,'(" disable_symmetries ",l1)')disable_symmetries 149 150 ! Input parameter wfng_flag for BerkeleyGW wavefunction file 151 call FindTag( tagHandler, 'wfng_flag', error) 152 if(error .eq. 0) then 153 read(unit = tagHandler%fileno, fmt = *, iostat = error) wfng_flag 154 endif 155 if(error .ne. 0) wfng_flag = .false. 156 write(*,'(" wfng_flag ",l1)')wfng_flag 157 158 ! Input parameter wfng_file for BerkeleyGW wavefunction file 159 call FindTag( tagHandler, 'wfng_file', error) 160 if(error .eq. 0) then 161 read(unit = tagHandler%fileno, fmt = *, iostat = error) wfng_file 162 endif 163 if(error .ne. 0) write(wfng_file,'("WFN")') 164 write(*,'(" wfng_file ",a)')trim(wfng_file) 165 166 ! Input parameter rhog_flag for BerkeleyGW charge density file 167 call FindTag( tagHandler, 'rhog_flag', error) 168 if(error .eq. 0) then 169 read(unit = tagHandler%fileno, fmt = *, iostat = error) rhog_flag 170 endif 171 if(error .ne. 0) rhog_flag = .false. 172 write(*,'(" rhog_flag ",l1)')rhog_flag 173 174 ! Input parameter rhog_file for BerkeleyGW charge density file 175 call FindTag( tagHandler, 'rhog_file', error) 176 if(error .eq. 0) then 177 read(unit = tagHandler%fileno, fmt = *, iostat = error) rhog_file 178 endif 179 if(error .ne. 0) write(rhog_file,'("RHO")') 180 write(*,'(" rhog_file ",a)')trim(rhog_file) 181 182 ! Input parameter vxcg_flag for BerkeleyGW exchange-correlation potential file 183 call FindTag( tagHandler, 'vxcg_flag', error) 184 if(error .eq. 0) then 185 read(unit = tagHandler%fileno, fmt = *, iostat = error) vxcg_flag 186 endif 187 if(error .ne. 0) vxcg_flag = .false. 188 write(*,'(" vxcg_flag ",l1)')vxcg_flag 189 190 ! Input parameter vxcg_file for BerkeleyGW exchange-correlation potential file 191 call FindTag( tagHandler, 'vxcg_file', error) 192 if(error .eq. 0) then 193 read(unit = tagHandler%fileno, fmt = *, iostat = error) vxcg_file 194 endif 195 if(error .ne. 0) write(vxcg_file,'("VXC")') 196 write(*,'(" vxcg_file ",a)')trim(vxcg_file) 197 198 ! Input parameter vxc_flag for BerkeleyGW exchange-correlation matrix element file 199 call FindTag( tagHandler, 'vxc_flag', error) 200 if(error .eq. 0) then 201 read(unit = tagHandler%fileno, fmt = *, iostat = error) vxc_flag 202 endif 203 if(error .ne. 0) vxc_flag = .false. 204 write(*,'(" vxc_flag ",l1)')vxc_flag 205 206 ! Input parameter vxc_file for BerkeleyGW exchange-correlation matrix element file 207 call FindTag( tagHandler, 'vxc_file', error) 208 if(error .eq. 0) then 209 read(unit = tagHandler%fileno, fmt = *, iostat = error) vxc_file 210 endif 211 if(error .ne. 0) write(vxc_file,'("vxc.dat")') 212 write(*,'(" vxc_file ",a)')trim(vxc_file) 213 214 ! Input parameter vxc_diag_nmin for BerkeleyGW exchange-correlation matrix element file 215 call FindTag( tagHandler, 'vxc_diag_nmin', error) 216 if(error .eq. 0) then 217 read(unit = tagHandler%fileno, fmt = *, iostat = error) vxc_diag_nmin 218 endif 219 if(error .ne. 0) vxc_diag_nmin = 0 220 write(*,'(" vxc_diag_nmin ",i2)')vxc_diag_nmin 221 222 ! Input parameter vxc_diag_nmax for BerkeleyGW exchange-correlation matrix element file 223 call FindTag( tagHandler, 'vxc_diag_nmax', error) 224 if(error .eq. 0) then 225 read(unit = tagHandler%fileno, fmt = *, iostat = error) vxc_diag_nmax 226 endif 227 if(error .ne. 0) vxc_diag_nmax = 0 228 write(*,'(" vxc_diag_nmax ",i2)')vxc_diag_nmax 229 230 ! Input parameter vxc_offdiag_nmin for BerkeleyGW exchange-correlation matrix element file 231 call FindTag( tagHandler, 'vxc_offdiag_nmin', error) 232 if(error .eq. 0) then 233 read(unit = tagHandler%fileno, fmt = *, iostat = error) vxc_offdiag_nmin 234 endif 235 if(error .ne. 0) vxc_offdiag_nmin = 0 236 write(*,'(" vxc_offdiag_nmin ",i2)')vxc_offdiag_nmin 237 238 ! Input parameter vxc_offdiag_nmax for BerkeleyGW exchange-correlation matrix element file 239 call FindTag( tagHandler, 'vxc_offdiag_nmax', error) 240 if(error .eq. 0) then 241 read(unit = tagHandler%fileno, fmt = *, iostat = error) vxc_offdiag_nmax 242 endif 243 if(error .ne. 0) vxc_offdiag_nmax = 0 244 write(*,'(" vxc_offdiag_nmax ",i2)')vxc_offdiag_nmax 245 246 ! for nspin = 2, a second copy of the wavefunctions and eigenvalues will be written 247 call FindTag( tagHandler, 'nspin', error) 248 if(error .eq. 0) then 249 read(unit = tagHandler%fileno, fmt = *, iostat = error) nspin 250 endif 251 if(error .ne. 0) nspin = 1 252 write(*,'(" nspin ",i2)') nspin 253 254 ! electronic temperature 255 call FindTag( tagHandler, 'temperature', error) 256 if(error .eq. 0) then 257 read(unit = tagHandler%fileno, fmt = *, iostat = error) temperature 258 endif 259 if(error .ne. 0) temperature = 0d0 260 write(*,'(" temperature ",f13.6," Ry")') temperature 261 temperature = temperature * 0.5d0 ! convert to Ha 262 263 ! Use form-factors from input file or hard-coded V(q) potentials 264 call FindTag( tagHandler,'FormFactors', error) 265 if(error .eq. 0) then 266 read(unit = tagHandler%fileno, fmt = *, iostat = error) ff 267 if(error /= 0) stop 'Error, could not read FormFactors.' 268 v_of_q = .false. 269 write(*,'("Using Form Factors V_3,8,11^S & V_3,4,11^A")') 270 else 271 ff(:) = 0.0d0 272 v_of_q = .true. 273 write(*,'("Using Potentials V(q) for Species")') 274 endif 275 276 ! Read LAPACK Absolute Tolerance 277 call FindTag( tagHandler, 'AbsoluteTolerance', error) 278 if(error .eq. 0) then 279 read(unit = tagHandler%fileno, fmt = *, iostat = error) abstol 280 endif 281 if(error .ne. 0) abstol = -one 282 write(*,'(" AbsoluteTolerance ",e13.6)')abstol 283 284 ! Read the Conjugate Gradient Parameter 285 call cgParamsInit(cgParams,tagHandler) 286 287 ! This will call DimensionsInit, LatticeInit, and AtomsInit. 288 call StructureInit( structure1, tagHandler ) 289 290 call KPointsInit( kPoints, graphInfo, structure1%ndim, structure1%lattice, tagHandler) 291 292 ! Gvector list and structurefactor 293 call HamInfoInit( hamInfo, structure1, tagHandler, v_of_q, ff ) 294 295 ! Find hMaxSize = max size of H at any k point 296 call FindHMaxSize( hamInfo, structure1, kPoints ) 297 298 ! Allocates array for H and work arrays 299 call HInit( hamiltonian, hamInfo%hMaxSize ) 300 hamiltonian%abstol = abstol 301 302 ! Allocates array for eigenStates array 303 call EigenStatesInit( eigenStates, kPoints%numBands, kPoints%numKPoints, & 304 & hamInfo%hMaxSize, -1.0d0, 1) 305 eigenStates%eigenVectors = zero 306 307 ! Initialise the FFT Grid and Plans 308 Call GridInit(Grid, hamInfo, eigenStates, structure1%ndim) 309 310! DAS -- do not write this useless file if not using CG 311 if(cgParams%Switch.eq.1) then 312 call open_file(unit=ReportNo,file='Report.cg',status='replace',position='rewind') 313 write(ReportNo,*) 'Eigenvalues:' 314 write(ReportNo,*) 'The FFT Grid Size +/-',Grid%Size(1) 315 write(ReportNo,*) 'The number of GVectors:',Grid%NumgVec 316 write(ReportNo,*) 'Eigenvalues:' 317 endif 318 319 ! Calculate the Potential and FFT to Vr 320 Call PotentialRGenerate(hamInfo, Grid) 321 322 if(cgParams%Switch .eq. 1) then 323 do k = 1, kPoints%numKPoints 324 call HGenerate( hamInfo, structure1, kPoints, eigenStates, & 325 & k, hamiltonian%h, hamiltonian%sqH, hamiltonian%hSize ) 326 327 write(ReportNo, *) 'k point ', k 328 329 call cgEigenSystem( hamiltonian, eigenStates, k, cgParams%tol, & 330 & cgParams%period ,ReportNo ) 331 332 write(*, '(" k point ", i4, " out of ", i4)') k, kPoints%numKPoints 333 334 end do ! k 335 else 336 do k = 1, kPoints%numKPoints 337 338 call HGenerate( hamInfo, structure1, kPoints, eigenStates, & 339 & k, hamiltonian%h, hamiltonian%sqH, hamiltonian%hSize ) 340 341! DAS -- do not write this useless file 342! write(ReportNo, *) 'k point ', k 343 344 call HDiagonalize( hamiltonian, eigenStates, k ) 345 346! DAS -- make sure stupid LAPACK routines have genuinely solved the Hamiltonian, else die 347#ifdef DEBUG 348 call HCheck( hamiltonian, eigenStates, k ) 349#endif 350 351 write(*, '(" k point ", i4, " out of ", i4)') k, kPoints%numKPoints 352 353 end do ! k 354 end if 355 356 if(cgParams%Switch.eq.1) then 357 call close_file(unit=ReportNo) 358 endif 359 360! Initialize data structures for BerkeleyGW output 361 362 a(:,:) = structure1%lattice%aLatVec(:,:) 363 al = sqrt(a(1,1)**2 + a(2,1)**2 + a(3,1)**2) 364 a(:,:) = a(:,:) / al 365 celvol = a(1,1) * (a(2,2) * a(3,3) - a(2,3) * a(3,2)) - & 366 a(2,1) * (a(1,2) * a(3,3) - a(1,3) * a(3,2)) + & 367 a(3,1) * (a(1,2) * a(2,3) - a(1,3) * a(2,2)) 368 bl = 2.0d0 * PI_D / al 369 b(1,1) = (a(2,2) * a(3,3) - a(3,2) * a(2,3)) / celvol 370 b(2,1) = (a(3,2) * a(1,3) - a(1,2) * a(3,3)) / celvol 371 b(3,1) = (a(1,2) * a(2,3) - a(2,2) * a(1,3)) / celvol 372 b(1,2) = (a(2,3) * a(3,1) - a(3,3) * a(2,1)) / celvol 373 b(2,2) = (a(3,3) * a(1,1) - a(1,3) * a(3,1)) / celvol 374 b(3,2) = (a(1,3) * a(2,1) - a(2,3) * a(1,1)) / celvol 375 b(1,3) = (a(2,1) * a(3,2) - a(3,1) * a(2,2)) / celvol 376 b(2,3) = (a(3,1) * a(1,2) - a(1,1) * a(3,2)) / celvol 377 b(3,3) = (a(1,1) * a(2,2) - a(2,1) * a(1,2)) / celvol 378 celvol = abs(celvol) * al**3 379 recvol = (2.0d0 * PI_D)**3 / celvol 380 adot(:,:) = 0.0d0 381 do i=1,3 382 do j=1,3 383 do k=1,3 384 adot(j,i) = adot(j,i) + a(k,j) * a(k,i) 385 enddo 386 enddo 387 enddo 388 adot(:,:) = adot(:,:) * al**2 389 bdot(:,:) = 0.0d0 390 do i=1,3 391 do j=1,3 392 do k=1,3 393 bdot(j,i) = bdot(j,i) + b(k,j) * b(k,i) 394 enddo 395 enddo 396 enddo 397 bdot(:,:) = bdot(:,:) * bl**2 398 399 ! kinetic energy cutoff for wave functions 400 ! convert from Hartree to Rydberg 401 ecutwfn=2.0d0*hamInfo%gVectors%energyCutOff 402 ! kinetic energy cutoff for charge density and potential 403 ecutrho=4.0d0*ecutwfn 404 gcutm=ecutrho/bl**2 405 406 call FindTag( tagHandler, 'FFTGrid', error) 407 if(error .eq. 0) then 408 read(unit = tagHandler%fileno, fmt = *, iostat = error) (FFTgrid(j),j=1,3) 409 if(error /= 0) stop 'Error, could not read FFTgrid.' 410 if(any(FFTgrid(:) <= 0)) stop 'Error, FFTGrid values must be positive.' 411 else 412 do i=1,3 413 FFTgrid(i)=int(2.0d0*sqrt(gcutm)*sqrt(a(1,i)**2+a(2,i)**2+a(3,i)**2))+1 414 do while (.not.check_FFT_size(FFTgrid(i), Nfac)) 415 FFTgrid(i)=FFTgrid(i)+1 416 enddo 417 enddo 418 endif 419 write(*,'(a,3i6)') 'FFT grid: ', FFTgrid(:) 420 421 nat=structure1%atomBasis%totalNumAtoms 422 SAFE_ALLOCATE(atyp, (nat)) 423 SAFE_ALLOCATE(apos, (3, nat)) 424 atyp(:) = structure1%atomBasis%atomicNumber(structure1%atomBasis%speciesOfAtom(:)) 425 apos(:,:) = structure1%atomBasis%atomCart(:,:)/al 426 427! gsm: if using form-factors make sgam_at think we have different species 428 if (.not. v_of_q .and. abs(ff(4)) + abs(ff(5)) + abs(ff(6)) .gt. TOL_Small) & 429 atyp(1) = atyp(1) + 1000 430 431 if(disable_symmetries) then 432 cell_symmetry = 0 433 nsym = 1 434 write(*,'(a)') 'Symmetries have been disabled.' 435 identity = reshape((/1, 0, 0, 0, 1, 0, 0, 0, 1/), shape(identity)) 436 rotation(:,:,nsym) = identity 437 else 438 call get_symmetries(nat, atyp, structure1%atomBasis%atomLat, a, FFTgrid, cell_symmetry, nsym, & 439 rotation, translation, spacegroup, symbol) 440 write(*,'(a,i3,a,a)') 'Space group ', spacegroup, ', symbol ', trim(symbol) 441 endif 442 443! gsm: restore atomic species back to normal 444 if (.not. v_of_q .and. abs(ff(4)) + abs(ff(5)) + abs(ff(6)) .gt. TOL_Small) & 445 atyp(1) = atyp(1) - 1000 446 447 call check_inversion(real_or_complex, nsym, rotation(:,:,:), nspin, .true., .true., tnp = translation) 448 449 if (wfng_flag) then 450 451 write(*,*) 452 453 SAFE_ALLOCATE(ngk, (kPoints%numKPoints)) 454 SAFE_ALLOCATE(kw, (kPoints%numKPoints)) 455 SAFE_ALLOCATE(kpt, (3, kPoints%numKPoints)) 456 ngk(1:kPoints%numKPoints)=eigenStates%numBasisVectors(1:kPoints%numKPoints) 457 kw(1:kPoints%numKPoints)=kPoints%kPointWeights(1:kPoints%numKPoints) 458 kpt(1:3,1:kPoints%numKPoints)=kPoints%kPointsLat(1:3,1:kPoints%numKPoints) 459 460 SAFE_ALLOCATE(ifmin, (kPoints%numKPoints,nspin)) 461 SAFE_ALLOCATE(ifmax, (kPoints%numKPoints,nspin)) 462 ifmin(1:kPoints%numKPoints,1:nspin)=1 463 ifmax(1:kPoints%numKPoints,1:nspin)=kPoints%numOccupiedBands 464 465 466 ! Convert eigenvalues from Hartree to Rydberg for BerkeleyGW 467 SAFE_ALLOCATE(energies, (kPoints%numBands, kPoints%numKPoints, nspin)) 468 do ispin=1,nspin 469 do ik=1,kPoints%numKPoints 470 do ib=1,kPoints%numBands 471 energies(ib,ik,ispin)=eigenStates%eigenValues(ib,ik)*2.0d0 472 enddo 473 enddo 474 enddo 475 476 if(temperature > TOL_Small) then 477 efermi = (minval(energies(kPoints%numOccupiedBands+1:,:,:)) + maxval(energies(1:kPoints%numOccupiedBands,:,:)))*0.5d0 478 write(*,'(a,f13.6,a)') 'Fermi energy (mid-gap): ', efermi*2d0, ' Ry' 479 endif 480 481 SAFE_ALLOCATE(occupations, (kPoints%numBands, kPoints%numKPoints, nspin)) 482 do ispin=1,nspin 483 do ik=1,kPoints%numKPoints 484 do ib=1,kPoints%numBands 485 if(temperature > TOL_Small) then 486 occupations(ib, ik, ispin) = 1 / (exp((energies(ib,ik,ispin) - efermi)/temperature) + 1) 487 else 488 if(ib < ifmin(ik, ispin) .or. ib > ifmax(ik, ispin)) then 489 occupations(ib, ik, ispin) = 0.0d0 490 else 491 occupations(ib, ik, ispin) = 1.0d0 492 endif 493 endif 494 enddo 495 enddo 496 enddo 497 498 ! Convert eigenvalues from Hartree to Rydberg for stdout 499 do ik=1,kPoints%numKPoints 500 do ib=1,kPoints%numBands 501! eigenStates%eigenValues(ib,ik)=eigenStates%eigenValues(ib,ik)*2.0d0*RYD 502 eigenStates%eigenValues(ib,ik)=eigenStates%eigenValues(ib,ik)*2.0d0 503 enddo 504 call PrintEigenvalues(eigenStates,ik) 505 write(*,*) 506 enddo 507 508 SAFE_ALLOCATE(gVectorsK, (3, eigenStates%maxBasisVectors)) 509 510 if(real_or_complex .eq. 1) then 511 write(*,'(" Applying Gram-Schmidt process...")') 512 SAFE_ALLOCATE(eigenVectorsReal, (eigenStates%maxBasisVectors, nspin, eigenStates%numStates, eigenStates%numKPoints)) 513 call real_wfng(eigenStates%numKPoints, eigenStates%numStates, & 514 & nspin, eigenStates%maxBasisVectors, eigenStates%numBasisVectors, & 515 & eigenStates%eigenValues, eigenStates%eigenVectors, eigenVectorsReal) 516 write(*,'(" ...done!")') 517 write(*,*) 518 else 519 SAFE_ALLOCATE(eigenVectorsComplex, (eigenStates%maxBasisVectors, nspin, eigenStates%numStates, eigenStates%numKPoints)) 520 do ik=1,kPoints%numKPoints 521 do ib=1,kPoints%numBands 522 do ispin=1,nspin 523 do ig=1,eigenStates%numBasisVectors(ik) 524 eigenVectorsComplex(ig, ispin, ib, ik) = & 525 eigenStates%eigenVectors(ig, ib, ik) 526 enddo 527 enddo 528 enddo 529 enddo 530 endif 531 532 write(*,'(" Writing BerkeleyGW wavefunction file...")') 533 534 call open_file(20,file=trim(wfng_file),status='replace',form='unformatted') 535 536 sheader = 'WFN' 537 call write_binary_header(20, sheader, real_or_complex, nspin, & 538 hamInfo%gVectors%numGVectors, nsym, cell_symmetry, nat, & 539 kPoints%numKPoints, kPoints%numBands, eigenStates%maxBasisVectors, & 540 ecutrho, ecutwfn, FFTgrid, kPoints%kPointNumbers, kPoints%kPointOffsets, & 541 celvol, al, a, adot, recvol, bl, b, bdot, rotation, translation, & 542 atyp, apos, ngk, kw, kpt, ifmin, ifmax, energies, occupations) 543 544 call write_binary_gvectors(20, hamInfo%gVectors%numGVectors, & 545 hamInfo%gVectors%numGVectors, hamInfo%gVectors%gVectors(:,:)) 546 547 do ik=1,kPoints%numKPoints 548 do ig=1,eigenStates%numBasisVectors(ik) 549 gVectorsK(:,ig)= & 550 hamInfo%gVectors%gVectors(:,eigenStates%basisIndex(ig,ik)) 551 enddo 552 call write_binary_gvectors(20, eigenStates%numBasisVectors(ik), & 553 eigenStates%maxBasisVectors, gVectorsK(:,:)) 554 555 do ib=1,kPoints%numBands 556 if(real_or_complex .eq. 1) then 557 call write_binary_real_data(20, eigenStates%numBasisVectors(ik), & 558 eigenStates%maxBasisVectors, nspin, eigenVectorsReal(:, :, ib, ik)) 559 else 560 call write_binary_complex_data(20, eigenStates%numBasisVectors(ik), & 561 eigenStates%maxBasisVectors, nspin, eigenVectorsComplex(:, :, ib, ik)) 562 endif 563 enddo 564 565 enddo 566 567 call close_file(20) 568 569 write(*,'(" ...done!")') 570 571 SAFE_DEALLOCATE(gVectorsK) 572 if (real_or_complex .eq. 1) then 573 SAFE_DEALLOCATE(eigenVectorsReal) 574 else 575 SAFE_DEALLOCATE(eigenVectorsComplex) 576 endif 577 578 endif 579 580 if (rhog_flag) then 581 582 write(*,*) 583 584 write(*,'(" Calculating charge density...")') 585 586 call DensityArrayInit( DensityArray, structure1%ndim, & 587 & FFTgrid, hamInfo%gVectors%numGVectors) 588 589 call CalculateDensity( hamInfo, kPoints, eigenStates, densityArray ) 590 591 write(*,'(" ...done!")') 592 593 write(*,*) 594 595 if(real_or_complex .eq. 1) then 596 SAFE_ALLOCATE(dataReal, (hamInfo%gVectors%numGVectors, nspin)) 597 do ispin=1,nspin 598 do ig=1,hamInfo%gVectors%numGVectors 599 dataReal(ig,ispin)= & 600 dble(densityArray%rhog(ig)*cmplx(2.0d0/dble(nspin),0.0d0,kind=double)) 601 enddo 602 enddo 603 else 604 SAFE_ALLOCATE(dataComplex, (hamInfo%gVectors%numGVectors, nspin)) 605 do ispin=1,nspin 606 do ig=1,hamInfo%gVectors%numGVectors 607 dataComplex(ig,ispin)= & 608 densityArray%rhog(ig)*cmplx(2.0d0/dble(nspin),0.0d0,kind=double) 609 enddo 610 enddo 611 endif 612 613 write(*,'(" Writing BerkeleyGW charge density file...")') 614 615 call open_file(20,file=trim(rhog_file),status='replace',form='unformatted') 616 617 sheader = 'RHO' 618 call write_binary_header(20, sheader, real_or_complex, nspin, & 619 hamInfo%gVectors%numGVectors, nsym, cell_symmetry, nat, & 620 kPoints%numKPoints, kPoints%numBands, eigenStates%maxBasisVectors, & 621 ecutrho, ecutwfn, FFTgrid, kPoints%kPointNumbers, kPoints%kPointOffsets, & 622 celvol, al, a, adot, recvol, bl, b, bdot, rotation, translation, & 623 atyp, apos, ngk, kw, kpt, ifmin, ifmax, energies, occupations) 624 625 call write_binary_gvectors(20, hamInfo%gVectors%numGVectors, & 626 hamInfo%gVectors%numGVectors, hamInfo%gVectors%gVectors(:,:)) 627 628 if(real_or_complex .eq. 1) then 629 call write_binary_real_data(20, hamInfo%gVectors%numGVectors, & 630 haminfo%gVectors%numGVectors, nspin, dataReal(:, :)) 631 else 632 call write_binary_complex_data(20, hamInfo%gVectors%numGVectors, & 633 haminfo%gVectors%numGVectors, nspin, dataComplex(:, :)) 634 endif 635 636 call close_file(20) 637 638 write(*,'(" ...done!")') 639 640 call DensityArrayDestroy( densityArray ) 641 642 if(real_or_complex .eq. 1) then 643 SAFE_DEALLOCATE(dataReal) 644 else 645 SAFE_DEALLOCATE(dataComplex) 646 endif 647 648 endif 649 650 if (vxcg_flag) then 651 652 write(*,*) 653 654 if(real_or_complex .eq. 1) then 655 SAFE_ALLOCATE(dataReal, (hamInfo%gVectors%numGVectors, nspin)) 656 do ispin=1,nspin 657 do ig=1,hamInfo%gVectors%numGVectors 658 dataReal(ig,ispin)=0.0d0 659 enddo 660 enddo 661 else 662 SAFE_ALLOCATE(dataComplex, (hamInfo%gVectors%numGVectors, nspin)) 663 do ispin=1,nspin 664 do ig=1,hamInfo%gVectors%numGVectors 665 dataComplex(ig,ispin)=(0.0d0,0.0d0) 666 enddo 667 enddo 668 endif 669 670 write(*,'(" Writing BerkeleyGW exchange-correlation potential file...")') 671 672 call open_file(20,file=trim(vxcg_file),status='replace',form='unformatted') 673 674 sheader = 'VXC' 675 call write_binary_header(20, sheader, real_or_complex, nspin, & 676 hamInfo%gVectors%numGVectors, nsym, cell_symmetry, nat, & 677 kPoints%numKPoints, kPoints%numBands, eigenStates%maxBasisVectors, & 678 ecutrho, ecutwfn, FFTgrid, kPoints%kPointNumbers, kPoints%kPointOffsets, & 679 celvol, al, a, adot, recvol, bl, b, bdot, rotation, translation, & 680 atyp, apos, ngk, kw, kpt, ifmin, ifmax, energies, occupations) 681 682 call write_binary_gvectors(20, hamInfo%gVectors%numGVectors, & 683 hamInfo%gVectors%numGVectors, hamInfo%gVectors%gVectors(:,:)) 684 685 if(real_or_complex .eq. 1) then 686 call write_binary_real_data(20, hamInfo%gVectors%numGVectors, & 687 hamInfo%gVectors%numGVectors, nspin, dataReal(:, :)) 688 else 689 call write_binary_complex_data(20, hamInfo%gVectors%numGVectors, & 690 hamInfo%gVectors%numGVectors, nspin, dataComplex(:, :)) 691 endif 692 693 call close_file(20) 694 695 write(*,'(" ...done!")') 696 697 if(real_or_complex .eq. 1) then 698 SAFE_DEALLOCATE(dataReal) 699 else 700 SAFE_DEALLOCATE(dataComplex) 701 endif 702 703 endif 704 705 SAFE_DEALLOCATE_P(atyp) 706 SAFE_DEALLOCATE_P(apos) 707 SAFE_DEALLOCATE_P(ngk) 708 SAFE_DEALLOCATE_P(kw) 709 SAFE_DEALLOCATE_P(kpt) 710 SAFE_DEALLOCATE_P(ifmin) 711 SAFE_DEALLOCATE_P(ifmax) 712 SAFE_DEALLOCATE_P(energies) 713 SAFE_DEALLOCATE_P(occupations) 714 715 if (vxc_flag) then 716 717 write(*,*) 718 719 if (vxc_diag_nmin .lt. 1) vxc_diag_nmin = 1 720 if (vxc_diag_nmax .gt. kPoints%numBands) vxc_diag_nmax = kPoints%numBands 721 ndiag = MAX (vxc_diag_nmax - vxc_diag_nmin + 1, 0) 722 if (vxc_offdiag_nmin .lt. 1) vxc_offdiag_nmin = 1 723 if (vxc_offdiag_nmax .gt. kPoints%numBands) vxc_offdiag_nmax = kPoints%numBands 724 noffdiag = MAX (vxc_offdiag_nmax - vxc_offdiag_nmin + 1, 0) 725 noffdiag = noffdiag**2 726 727 do ispin = 1, nspin 728 spin_index(ispin) = ispin 729 enddo 730 SAFE_ALLOCATE(diag, (ndiag)) 731 do idiag = 1, ndiag 732 diag(idiag) = vxc_diag_nmin + idiag - 1 733 enddo 734 SAFE_ALLOCATE(offdiag, (noffdiag)) 735 do ioffdiag = 1, noffdiag 736 offdiag(ioffdiag) = vxc_offdiag_nmin + ioffdiag - 1 737 enddo 738 SAFE_ALLOCATE(vxc_mtxel, (ndiag + noffdiag, nspin)) 739 vxc_mtxel(:,:) = CMPLX(0d0, 0d0) 740 741 write(*,'(" Writing BerkeleyGW exchange-correlation matrix element file...")') 742 743 call open_file(20,file=trim(vxc_file),status='unknown',form='formatted') 744 do ik=1,kPoints%numKPoints 745 call write_matrix_elements(20, kPoints%kPointsLat(:,ik), nspin, ndiag, noffdiag, & 746 spin_index, diag, offdiag, offdiag, vxc_mtxel(:,:)) 747 enddo 748 call close_file(20) 749 750 SAFE_DEALLOCATE(diag) 751 SAFE_DEALLOCATE(offdiag) 752 SAFE_DEALLOCATE(vxc_mtxel) 753 754 write(*,'(" ...done!")') 755 756 endif 757 758 write(*,*) 759 760 ! Clean up memory 761 call EigenStatesDestroy( eigenStates ) 762 763 call HDestroy( hamiltonian ) 764 765 call GridDestroy( Grid) 766 767 call HamInfoDestroy( hamInfo ) 768 769 call KPointsDestroy( kPoints ) 770 771 call StructureDestroy( structure1 ) 772 773 ! Clean up and close the input file 774 call TagHandlerDestroy( tagHandler ) 775 776 call open_file(unit = 10, file = filename, status = 'old') 777 call close_file(unit = 10, delete = .true.) 778 779 999 format(a) 780 781contains 782 783!************************************************************************* 784subroutine real_wfng (nk, nb, ns, ngkmax, ng, en, wfc, wfr) 785!************************************************************************* 786! Constructs real wavefunctions in G-space for systems 787! with inversion symmetry by applying Gram-Schmidt process. 788! Based on paratecSGL/src/para/gwreal.f90 789!************************************************************************* 790 791 use SysParams, only : double 792 use message_m 793 794 implicit none 795 796 integer, intent(in) :: nk 797 integer, intent(in) :: nb 798 integer, intent(in) :: ns 799 integer, intent(in) :: ngkmax 800 integer, intent(in) :: ng(nk) 801 real(double), intent(in) :: en(:, :) !< (nb, nk) 802 complex(double), intent(in) :: wfc(:, :, :) !< (ngkmax, nb, nk) 803 real(double), intent(out) :: wfr(:, :, :, :) !< (ngkmax, ns, nb, nk) 804 805 real(double), parameter :: eps2 = 1.0d-2 806 real(double), parameter :: eps5 = 1.0d-5 807 real(double), parameter :: eps6 = 1.0d-6 808 809 integer :: i, j, k, ik, ib, jb, is, ig, deg, mdeg, inc 810 integer :: dimension_span, reduced_span 811 real(double) :: x 812 integer, allocatable :: inum (:) 813 integer, allocatable :: imap (:, :) 814 integer, allocatable :: inull(:) 815 integer, allocatable :: null_map(:, :) 816 real(double), allocatable :: psi(:, :) 817 real(double), allocatable :: phi(:, :) 818 real(double), allocatable :: vec(:) 819 820 ! determine size of degenerate subspace 821 mdeg = 1 822 do ik = 1, nk 823 do ib = 1, nb 824 deg = 1 825 do jb = ib + 1, nb 826 if (abs(en(ib, ik) - en(jb, ik)) .lt. & 827 eps5 * dble(jb - ib + 1)) deg = deg + 1 828 enddo 829 if (deg .gt. mdeg) mdeg = deg 830 enddo 831 enddo 832 mdeg = mdeg * 2 833 834 SAFE_ALLOCATE(imap, (nb, nk)) 835 SAFE_ALLOCATE(inum, (nk)) 836 SAFE_ALLOCATE(inull, (nb)) 837 SAFE_ALLOCATE(null_map, (mdeg, nb)) 838 839 do ik = 1, nk 840 inum(ik) = 1 841 do ib = 1, nb 842 if (ib .eq. nb) then 843 imap(inum(ik), ik) = ib 844 inum(ik) = inum(ik) + 1 845 elseif (abs(en(ib, ik) - en(ib + 1, ik)) .gt. eps5) then 846 imap(inum(ik), ik) = ib 847 inum(ik) = inum(ik) + 1 848 endif 849 enddo 850 inum(ik) = inum(ik) - 1 851 enddo 852 853 SAFE_ALLOCATE(psi, (ngkmax, mdeg)) 854 SAFE_ALLOCATE(phi, (ngkmax, mdeg)) 855 SAFE_ALLOCATE(vec, (ngkmax)) 856 857 do ik = 1, nk 858 inc = 1 859 do i = 1, inum(ik) 860 inull(i) = 1 861 do ib = inc, imap(i, ik) 862 x = 0.0d0 863 do ig = 1, ng(ik) 864 x = x + dble(wfc(ig, ib, ik))**2 865 enddo 866 if (x .lt. eps2) null_map(inull(i), i) = 0 867 if (x .gt. eps2) null_map(inull(i), i) = 1 868 inull(i) = inull(i) + 1 869 x = 0.0d0 870 do ig = 1, ng(ik) 871 x = x + IMAG(wfc(ig, ib, ik))**2 872 enddo 873 if (x .lt. eps2) null_map(inull(i), i) = 0 874 if (x .gt. eps2) null_map(inull(i), i) = 1 875 inull(i) = inull(i) + 1 876 enddo 877 inull(i) = inull(i) - 1 878 inc = imap(i, ik) + 1 879 enddo 880 inc = 1 881 ib = 1 882 do i = 1, inum(ik) 883 k = 1 884 do j = 1, 2 * (imap(i, ik) - inc) + 1, 2 885 if (null_map(j, i) .eq. 1 .or. null_map(j + 1, i) .eq. 1) then 886 if (null_map(j, i) .eq. 1) then 887 do ig = 1, ng(ik) 888 phi(ig, k) = dble(wfc(ig, ib, ik)) 889 enddo 890 k = k + 1 891 endif 892 if (null_map(j + 1, i) .eq. 1) then 893 do ig = 1, ng(ik) 894 phi(ig, k) = IMAG(wfc(ig, ib, ik)) 895 enddo 896 k = k + 1 897 endif 898 ib = ib + 1 899 endif 900 enddo 901 dimension_span = k - 1 902 if (dimension_span .eq. 0) then 903 write(0,201)ik,inc 904 stop 905 endif 906 do j = 1, dimension_span 907 x = 0.0d0 908 do ig = 1, ng(ik) 909 x = x + phi(ig, j)**2 910 enddo 911 x = sqrt(x) 912 do ig = 1, ng(ik) 913 phi(ig, j) = phi(ig, j) / x 914 enddo 915 enddo 916! 917! the Gram-Schmidt process begins 918! 919 reduced_span = 1 920 do ig = 1, ng(ik) 921 psi(ig, 1) = phi(ig, 1) 922 enddo 923 do j = 1, dimension_span - 1 924 do ig = 1, ng(ik) 925 vec(ig) = phi(ig, j + 1) 926 enddo 927 do k = 1, reduced_span 928 x = 0.0d0 929 do ig = 1, ng(ik) 930 x = x + phi(ig, j + 1) * psi(ig, k) 931 enddo 932 do ig = 1, ng(ik) 933 vec(ig) = vec(ig) - psi(ig, k) * x 934 enddo 935 enddo 936 x = 0.0d0 937 do ig = 1, ng(ik) 938 x = x + vec(ig)**2 939 enddo 940 x = sqrt(x) 941 if (x .gt. eps6) then 942 reduced_span = reduced_span + 1 943 do ig = 1, ng(ik) 944 psi(ig, reduced_span) = vec(ig) / x 945 enddo 946 endif 947 enddo 948! 949! the Gram-Schmidt process ends 950! 951 if (reduced_span .lt. imap(i, ik) - inc + 1) then 952 write(0,202) ik,inc 953 stop 954 endif 955 do ib = inc, imap(i, ik) 956 do is = 1, ns 957 do ig = 1, ng(ik) 958 wfr(ig, is, ib, ik) = psi(ig, ib - inc + 1) 959 enddo 960 enddo 961 enddo 962 inc = imap(i, ik) + 1 963 enddo 964 enddo 965 966 SAFE_DEALLOCATE(inum) 967 SAFE_DEALLOCATE(imap) 968 SAFE_DEALLOCATE(inull) 969 SAFE_DEALLOCATE(null_map) 970 SAFE_DEALLOCATE(psi) 971 SAFE_DEALLOCATE(phi) 972 SAFE_DEALLOCATE(vec) 973 974 return 975 976 201 format(1x,"ERROR: failed Gram-Schmidt dimension span",/,14x, & 977 & "k-point =",1x,i4.4,1x,"band =",1x,i4.4,/) 978 202 format(1x,"ERROR: failed Gram-Schmidt reduced span",/,14x, & 979 & "k-point =",1x,i4.4,1x,"band =",1x,i4.4,/) 980 981end subroutine real_wfng 982 983end Program EPM2BGW 984 985