! Copyright (C) 2001-2009 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!-----------------------------------------------------------------------
PROGRAM gw
  !-----------------------------------------------------------------------
  ! ... This is the main driver of the GW code.
  ! ... It reads all the quantities calculated by pwscf, it
  ! ... checks if some recover file is present and determines
  ! ... which calculation needs to be done. Finally, it makes
  ! ... a loop over the q points. At a generic q, if necessary, it
  ! ... recalculates the band structure by calling pwscf again.
  ! NC = norm conserving pseudopotentials
  ! US = ultrasoft pseudopotentials
  ! PAW = projector augmented-wave
  ! [1] LDA, [2] [1]+GGA, [3] [2]+LSDA/sGGA, [4] [3]+Spin-orbit/nonmagnetic,
  ! [5] [4]+Spin-orbit/magnetic

  USE io_global,        ONLY : stdout
  USE wvfct,            ONLY : nbnd

 !HL epsil, trans, elph,
  USE disp,             ONLY : nqs
  USE output,           ONLY : fildrho
  USE check_stop,       ONLY : check_stop_init
  USE gw_restart,       ONLY : gw_writefile, destroy_status_run
  USE save_gw,          ONLY : clean_input_variables
  USE mp_global,        ONLY: mp_startup, nimage
  USE path_io_routines, ONLY : io_path_start

  !USE control_gw,       ONLY : bands_computed, nbnd_occ
  USE control_gw,      ONLY : done_bands, reduce_io, recover, tmp_dir_gw, &
                              ext_restart, bands_computed, bands_computed, nbnd_occ

  USE input_parameters,ONLY : pseudo_dir
  USE io_files,        ONLY : prefix, tmp_dir
  USE control_flags,   ONLY : restart
  USE qpoint,           ONLY : xq

  USE save_gw,         ONLY : tmp_dir_save
  USE environment,      ONLY: environment_start
  USE freq_gw,          ONLY : nfs, nwsigma
  USE units_gw,         ONLY : iuncoul, iungreen, lrgrn, lrcoul, iunsigma, lrsigma
  USE basis,            ONLY : starting_wfc, starting_pot, startingconfig
  USE gwsigma,          ONLY : nr1sex, nr2sex, nr3sex, nrsex, nlsex, ecutsex, &
                               nr1sco, nr2sco, nr3sco, nrsco, nlsco, ecutsco, &
                               ngmsig
   
  USE gvect,           ONLY : nl
  
  IMPLICIT NONE
  INTEGER :: iq, ik
  INTEGER :: ios
  LOGICAL :: do_band, do_iq, setup_pw, exst
  CHARACTER (LEN=9)   :: code = 'GW'
  CHARACTER (LEN=256) :: auxdyn

   ! Initialize MPI, clocks, print initial messages
   ! /Modules/mp.f90

#ifdef __PARA
   CALL mp_startup ( )
   IF (nimage>1) CALL io_path_start()
#endif
   
   !/Modules/environment.f90 prints out all parallel information and opening message.
    CALL environment_start ( code )
    CALL gwq_readin()

    WRITE(stdout, '(/5x, "Finished reading variables")') 

   !HL
   !Check stop init Modules/check_stop.f90
   !This module contains functions to check if the code should
   !be smoothly stopped.

     CALL check_stop_init()

   !This routine checks the initial status of the GW run, initializes the qmesh, and prepares
   !the control of the dispersion calculation. 

     CALL check_initial_status(auxdyn)

!    Generate frequency grid for GW convolution. 
     CALL freqbins()

!    CALL refold() 
     CALL clean_pw( .FALSE. )
     CALL allocate_fft()
     CALL ggen()

!   Generating Exchange and Correlation grid.
    CALL sig_fft_g(nr1sco, nr2sco, nr3sco, nrsco, ecutsco, 1)
    CALL sig_fft_g(nr1sex, nr2sex, nr3sex, nrsex, ecutsex, 2)

    CALL clean_pw( .FALSE. )

!   Only writing correlation energy to disk so ngmsigma = ngmsco.   

!   Coulomb file
    iuncoul = 28
    lrcoul = 2 * ngmsig * ngmsig * nfs
    CALL diropn (iuncoul, 'coul', lrcoul, exst)

!   Green's function file
    iungreen = 31
    lrgrn  = 2 * ngmsig * ngmsig
    CALL diropn (iungreen, 'green', lrgrn, exst)

!   Sigma file
    iunsigma = 32
    lrsigma = 2 * ngmsig * ngmsig * nwsigma
    CALL diropn(iunsigma, 'sigma', lrsigma, exst)

!   Should sigma_ex need to be written to file:
!   iunsex = 33
!   lrsex = 2 * ngmsex * ngmsex
!   CALL diropn(iunsex, 'sigma_ex', lrsex, exst)
!   CALCULATE W(r,r';iw)
GOTO 123
    DO iq = 1, nqs
        CALL prepare_q(do_band, do_iq, setup_pw, iq)
        CALL run_pwscf(do_band)
        CALL initialize_gw()

       !CALCULATE W(r,r';iw):
        CALL coulomb(iq)
        CALL clean_pw_gw(iq)
    END DO
    WRITE(stdout, '("Finished Calculating Screened Coulomb")') 

   WRITE(stdout, '(/5x, "GREEN LINEAR SYSTEM SOLVER")')

   DO ik = 1, 1
!     For Green's Fxn, need to generate a new list of 
!     Symmetry operations for k-points not at gamma...
!     CALL setup_green(ik)
      DO iq = 1, nqs
        CALL prepare_kmq(do_band, do_iq, setup_pw, iq, ik)
        CALL run_pwscf(do_band)
        CALL initialize_gw()
       !CALCULATE G(r,r'; w) 
        CALL green_linsys(ik, iq)
        CALL clean_pw_gw(iq)
      ENDDO
   ENDDO

   WRITE(stdout,'("Finished Calculating Greens Function")') 

!  For sigma_c and sigma_exchange should be able to just generate 
!  the eigenvalues/eigenvectors once...  need to get the sign right.

123 CONTINUE
   DO ik = 1, 1
     !CALCULATE Sigma_corr(r,r';w) = iG(r,r';w + w')(W(r,r';w') - v(r,r'))
     !CALL sigma_c(ik)
     !CALCULATE Sigma_ex(r,r') = iG(r,r')v(r,r')
      CALL sigma_exch(ik)
   ENDDO

   WRITE(6, '("Finished CALCULATING SIGMA")') 

   DO ik = 1, 1
       CALL sigma_matel(ik) 
   ENDDO
    
   STOP
   CALL gw_writefile('init',0)
   CALL clean_input_variables()
   CALL collect_grid_files()
   CALL destroy_status_run()
   IF (bands_computed) CALL print_clock_pw()
   CALL stop_gw( .TRUE. )
END PROGRAM gw