! 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, ionode_id, ionode USE wvfct, ONLY : nbnd,npwx USE disp, ONLY : nqs, num_k_pts, xk_kpoints, w_of_q_start 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, npool, intra_image_comm, inter_image_comm, & nproc_pool, mpime, nproc, my_pool_id, me_pool, & mp_global_end, inter_pool_comm USE mp, ONLY: mp_barrier, mp_bcast, mp_sum, mp_end USE parallel_include, ONLY: mpi_comm_world USE path_io_routines, ONLY : io_path_start USE control_gw, ONLY : done_bands, reduce_io, recover, tmp_dir_gw,& ext_restart, bands_computed, bands_computed, nbnd_occ, lgamma,& do_coulomb, do_sigma_c, do_sigma_exx,do_sigma_exxG, do_green,& do_sigma_matel, do_q0_only, multishift, do_sigma_extra,& solve_direct, tinvert, lrpa, do_serial, do_epsil, do_imag, do_pade_coul 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, nwgreen USE units_gw, ONLY : iuncoul, iungreen, lrgrn, lrcoul, iunsigma, lrsigma, lrsex, iunsex,& iunresid, lrresid, iunalphabeta, lralphabeta, iunsigext, lrsigext USE basis, ONLY : starting_wfc, starting_pot, startingconfig USE gwsigma, ONLY : nr1sex, nr2sex, nr3sex, nrsex, nlsex, ecutsex, & nr1sco, nr2sco, nr3sco, nrsco, nlsco, ecutsco, & ngmsig, ngmsex, ecutsig, ngmsco, ngmgrn, ngmpol, nbnd_sig USE gvect, ONLY : nl, g USE kinds, ONLY : DP USE gwsymm, ONLY : ngmunique, ig_unique, use_symm, sym_friend, sym_ig USE, INTRINSIC :: ieee_arithmetic IMPLICIT NONE INTEGER :: iq, ik INTEGER :: ios LOGICAL :: do_band, do_iq, setup_pw, exst CHARACTER (LEN=9) :: code = 'GW' CHARACTER (LEN=256) :: auxdyn REAL(DP) :: Gb !Variables to split Planewaves perturbations across processors. INTEGER :: igstart, igstop, ngpool, ngr, igs, ig COMPLEX(DP), ALLOCATABLE :: scrcoul_g(:,:,:,:) !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 ) WRITE(stdout, '(/5x, "Reading variables")') 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 clean_pw( .FALSE. ) CALL allocate_fft() CALL ggen() ! Generating Exchange and Correlation grid. CALL sig_fft_g(nr1sco, nr2sco, nr3sco, nrsco, ecutsig, 1) write(6,'(4x,"Exchange cutoffs: ")') CALL sig_fft_g(nr1sex, nr2sex, nr3sex, nrsex, ecutsex, 2) Gb = (DBLE(2)**26) WRITE( stdout, '(8x,"Greens fxn: ", f10.2," Gb", 5x,"(",i7,",",i7,",",i7,")")') & ((DBLE(ngmgrn)**2)*float(nwgreen))/Gb, ngmgrn, ngmgrn, nwgreen WRITE( stdout, '(8x,"Sigma(r,rp) fxn: ", f10.2," Gb", 5x,"(",i7,",",i7,")")') & ((DBLE(nrsco)**2))/Gb, nrsco, nrsco write(6,*) write(6,'(4x, "Screening:")') if(lrpa) then write(6,'(4x,"RPA Screening")') else write(6,'(4x,"Using Approx. Vertex Correction")') endif write(6,*) ALLOCATE ( scrcoul_g( ngmpol, ngmpol, nfs, 1)) ALLOCATE ( ig_unique( ngmpol) ) ALLOCATE ( sym_ig(ngmpol)) ALLOCATE ( sym_friend(ngmpol)) iuncoul = 28 lrcoul = 2 * ngmpol * ngmpol * nfs iungreen = 31 lrgrn = 2 * ngmgrn * ngmgrn iunsigma = 32 lrsigma = 2*ngmsco*ngmsco if(multishift) then iunresid = 34 lrresid = 2*npwx !Could probably keep the alphabeta coefficients in memory. iunalphabeta = 35 lralphabeta = 4 endif IF (ionode) THEN iuncoul = 28 lrcoul = 2 * ngmpol * ngmpol * nfs CALL diropn (iuncoul, 'coul', lrcoul, exst) ! Green's function file iungreen = 31 lrgrn = 2 * ngmgrn * ngmgrn CALL diropn (iungreen, 'green', lrgrn, exst) ! Sigma file iunsigma = 32 if(do_serial) then lrsigma = 2 * ngmsco * ngmsco else lrsigma = 2 * ngmsco * ngmsco * nwsigma endif 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) ! Should sigma_extra need to be written to file: iunsigext = 36 lrsigext = 2 * ngmsco * ngmsco CALL diropn(iunsigext, 'sig_ext', lrsigext, exst) ENDIF IF(do_pade_coul) THEN write(6, *) write(6,'(4x,"Analytic Continuation of Epsilon")') if(ionode) CALL pade_eps() call mp_barrier(inter_pool_comm) GOTO 126 ENDIF CALL clean_pw( .FALSE. ) IF(do_coulomb) THEN DO iq = w_of_q_start, nqs scrcoul_g(:,:,:,:) = (0.0d0, 0.0d0) !Prepare k, k+q grids, run nscf calculation, find small group of q. CALL prepare_q(do_band, do_iq, setup_pw, iq) CALL run_pwscf(do_band) CALL initialize_gw() !Determine the unique G vectors in the small group of q if symmetry is being used. !If not then all the vectors up to ngmpol (correlation cutoff) are treated as unique and calculated. if(use_symm) then WRITE(6,'("")') WRITE(6,'("SYMMETRIZING COULOMB Perturbations")') WRITE(6,'("")') CALL stern_symm() else ngmunique = ngmpol do ig = 1, ngmpol ig_unique(ig) = ig enddo endif !Distribute unique G-vectors between processors: #ifdef __PARA npool = nproc / nproc_pool if (npool.gt.1) then ! number of g-vec per pool and reminder ngpool = ngmunique / npool ngr = ngmunique - ngpool * npool ! the remainder goes to the first ngr pools if ( my_pool_id < ngr ) ngpool = ngpool + 1 igs = ngpool * my_pool_id + 1 if ( my_pool_id >= ngr ) igs = igs + ngr ! the index of the first and the last g vec in this pool igstart = igs igstop = igs - 1 + ngpool write (stdout,'(/4x,"Max n. of PW perturbations per pool = ",i5)') igstop-igstart+1 else #endif igstart = 1 igstop = ngmunique #ifdef __PARA endif #endif !CALCULATE W(G,G';iw): if((igstop-igstart+1).ne.0) then CALL coulomb(iq, igstart, igstop, scrcoul_g) endif !COLLECT G-VECTORS: CALL mp_barrier(inter_pool_comm) CALL mp_sum ( scrcoul_g, inter_pool_comm ) CALL mp_barrier(inter_pool_comm) !Write W_{q}(G,G';iw) to file: IF (ionode) THEN CALL unfold_w(scrcoul_g,iq) !If solve direct now need to invert epsilon: IF(solve_direct.and.tinvert) CALL invert_epsilon(scrcoul_g, iq) CALL davcio(scrcoul_g, lrcoul, iuncoul, iq, +1, ios) ENDIF CALL mp_barrier(inter_pool_comm) CALL clean_pw_gw(iq) CALL mp_barrier(inter_pool_comm) if(do_q0_only) GOTO 124 if(do_epsil.and.(iq.eq.nqs)) GOTO 126 END DO WRITE(stdout, '("Finished Calculating Screened Coulomb")') ENDIF 124 CONTINUE !Free up memory from scrcoul_g call mp_barrier(inter_pool_comm) DEALLOCATE( scrcoul_g ) DEALLOCATE( ig_unique ) DO ik = 1, 1 !DO ik = 1, num_k_pts if(do_green.and.multishift) CALL diropn(iunresid, 'resid', lrresid, exst) if(do_green.and.multishift) CALL diropn(iunalphabeta, 'alphbet', lralphabeta, exst) xq(:) = xk_kpoints(:, ik) WRITE(6,'(4x,"Sigma_k", 3f12.7)') xk_kpoints(:,ik) do_iq=.TRUE. lgamma = ( xq(1) == 0.D0 .AND. xq(2) == 0.D0 .AND. xq(3) == 0.D0 ) setup_pw = .TRUE. do_band = .TRUE. ! Generates small group of k and then forms IBZ_{k}. CALL run_pwscf_green(do_band) CALL initialize_gw() ! CALCULATE G(r,r'; w) ! WRITE(stdout, '(/5x, "GREEN LINEAR SYSTEM SOLVER")') if(do_green) write(6,'("Do green_linsys")') if(do_green.and.(.not.multishift)) then write(6,'("Green Not-MultiShift, Serial sigma_C", i4)') ik CALL green_linsys(ik) call mp_barrier() if(do_sigma_c) CALL sigma_c(ik) endif ! CALCULATE Sigma_corr(r,r';w) = i\int G(r,r'; w + w')(W(r,r';w') - v(r,r')) dw' if(do_green.and.multishift.and.do_serial) then write(6,'("Green Linsys_Shift Serial", i4)') ik if(.not.do_imag) then CALL green_linsys_serial(ik) else if(do_imag) then CALL green_linsys_serial_im(ik) endif call mp_barrier() else if (multishift.and.(.not.do_serial))then if(do_green) CALL green_linsys_shift(ik) if(do_sigma_c) CALL sigma_c(ik) call mp_barrier() endif if(do_green.and.multishift) then CLOSE(UNIT = iunresid, STATUS = 'DELETE') CLOSE(UNIT = iunalphabeta, STATUS = 'DELETE') endif if(ionode) then !CALCULATE Sigma_ex(r,r') = iG(r,r')v(r,r') if(do_sigma_exx.and..not.(do_sigma_exxG)) CALL sigma_exch(ik) endif if(ionode) WRITE(6, '("Finished CALCULATING SIGMA")') CALL mp_barrier(inter_pool_comm) CALL clean_pw_gw(ik) CALL mp_barrier(inter_pool_comm) ENDDO DO ik = 1, 1 if (do_sigma_matel) then nbnd = nbnd_sig if(nbnd_sig.lt.nbnd) write(6, '("nbnd_sig is less than the number of occupied states. Setting default nbnd=nbnd+4.")') if(nbnd_sig.lt.nbnd) nbnd = nbnd + 4 xq(:) = xk_kpoints(:, ik) do_iq=.TRUE. lgamma = ( xq(1) == 0.D0 .AND. xq(2) == 0.D0 .AND. xq(3) == 0.D0 ) setup_pw = .TRUE. do_band = .TRUE. CALL run_pwscf_green(do_band) CALL initialize_gw() WRITE(6,'(4x,"Sigma exchange")') if(ionode.and.do_sigma_exxG) CALL sigma_exchg(ik) CALL mp_barrier(inter_pool_comm) CALL sigma_matel(ik) CALL mp_barrier(inter_pool_comm) CALL clean_pw_gw(ik) endif ENDDO 126 CONTINUE call mp_barrier(inter_pool_comm) 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