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