!
! 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 control_gw,      ONLY : bands_computed 
  !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 environment,     ONLY: environment_start
  USE units_gw,    ONLY : iuncoul, iungreen
  !
  IMPLICIT NONE
  !
  INTEGER :: iq, ik
  INTEGER :: ios
  LOGICAL :: do_band, do_iq, setup_pw
  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()


  ! Checking the status of the calculation and if necessary initialize the q mesh
  ! This routine checks the initial status of the GW run and prepares
  ! the control of the dispersion calculation. 
  
  CALL check_initial_status(auxdyn)

 !Generate frequency grid for GW convolution and G-vector mapping.
  CALL freqbins()
  CALL refold()

  DO iq = 1, nqs
     !
     CALL prepare_q(do_band, do_iq, setup_pw, iq)
     !
     !  If this q is not done in this run, cycle
     !
     IF (.NOT.do_iq) CYCLE
     !
     !  If necessary the bands are recalculated
     !
     IF (setup_pw) WRITE (stdout, '("RUNNING NSCF CALCULATION")')
     IF (setup_pw) CALL run_pwscf(do_band)

! Initialize the quantities which do not depend on
! the linear response of the system

     CALL initialize_gw()

! CALCULATE W(r,r';iw)
! We are going to do one NSCF step to generate a k grid and then the k and k+q for the linear system
! will be obtained by shuffling the {k} grid with the grid_shuffle routine. 

     CALL coulomb(iq)

! if(fildrho /= ' ' )call punch_plot_gw() routine for plotting charge density which could be useful.
!     CALL clean_pw_gw(iq)

  END DO

!CALCULATE G(r, r';iw) 
  WRITE(stdout, '("GREEN LINEAR SYSTEM SOLVER")')
  DO ik = 1, 1
     DO iq = 1, nqs
       CALL green_linsys(ik, iq)
       !rec0 = (ik0-1) * nq + (iq-1) + 1
       !greenf_g = greenf(1:ngms,1:ngms,:)
       !write ( iungreen, rec = rec0, iostat = ios) greenf_g
     ENDDO
  ENDDO

! Form GW product. 

  DO ik = 1, 1
    CALL gw_product(ik, iq)
  ENDDO

! Take Matrix Elements
! do ik0 = 1, nk0
!   CALL sigma_matel() 
! enddo 

  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. )
  !
  STOP
  !
END PROGRAM gw