1!
2! Copyright (C) 2001-2013 Quantum ESPRESSO group
3! This file is distributed under the terms of the
4! GNU General Public License. See the file `License'
5! in the root directory of the present distribution,
6! or http://www.gnu.org/copyleft/gpl.txt .
7!
8!
9
10
11!this is the main part of the GWW program
12
13   PROGRAM gww
14
15   USE  input_gw,            ONLY : input_options, read_input_gww
16   USE  io_global,           ONLY : stdout, ionode
17   USE  self_energy_storage
18   USE  expansion
19   USE  energies_gww
20   USE  start_end
21   USE  mp_world,           ONLY : mpime, world_comm
22   USE  para_gww
23   USE  times_gw
24   USE  w_divergence
25   USE  mp,                  ONLY : mp_barrier
26   USE  contour
27      USE io_files,  ONLY : prefix, tmp_dir
28
29   implicit none
30
31   TYPE(input_options)   :: options
32   TYPE(self_storage)    :: ss
33   TYPE(self_expansion)  :: se
34   TYPE(self_on_real)    :: sr
35   TYPE(quasi_particles) :: qp
36   TYPE(times_freqs)     :: tf
37   TYPE(gv_time)         :: gt
38   TYPE(w_expectation)   :: we
39   TYPE(w_poles)         :: wp
40
41   INTEGER :: ispin
42   CHARACTER(5) :: name_proc
43   INTEGER :: ie
44   REAL(kind=DP) :: energy
45
46#if defined(_OPENMP)
47   INTEGER :: omp_get_num_threads, omp_get_max_threads
48   EXTERNAL omp_set_num_threads, omp_get_num_threads, omp_get_max_threads
49#endif
50
51   tmp_dir=' '
52
53!setup MPI environment
54
55   call startup
56
57   !CALL remove_stack_limit ( )
58
59#if defined(_OPENMP)
60     ntids=omp_get_max_threads()
61     ! call omp_set_num_threads(1)
62#endif
63
64
65
66#if defined(_OPENMP)
67   write(stdout,*)  'ntids = ', ntids
68#endif
69
70
71
72!initialize arrays
73   call initialize_quasi_particle(qp)
74
75
76!  read in input structure
77
78   call read_input_gww(options)
79#if defined(__MPI)
80   if(options%l_verbose) then
81      write(name_proc,'(5i1)') &
82           & (mpime+1)/10000,mod(mpime+1,10000)/1000,mod(mpime+1,1000)/100,mod(mpime+1,100)/10,mod(mpime+1,10)
83      OPEN( UNIT = stdout, FILE = trim(tmp_dir)//trim(prefix)//'-out_'//name_proc, STATUS = 'UNKNOWN' )
84   endif
85#endif
86
87   FLUSH(stdout)
88   if(options%grid_freq/=5.and.options%grid_freq/=6) then
89      call setup_para_gww(options%n, options%max_i, options%i_min, options%i_max)
90   else
91      call setup_para_gww(options%n+(1+2*options%second_grid_i)*options%second_grid_n, options%max_i, options%i_min, options%i_max)
92   endif
93   FLUSH(stdout)
94! setup time/frequency grid if required
95   call setup_timefreq(tf,options)
96
97!Step 0
98!calculates the exchange energies
99
100
101
102   if(options%starting_point <=1) then
103      call  go_exchange_main( options, qp)
104      call write_quasi_particles(qp, options,.false.)
105   else
106      call read_quasi_particles(qp,options,.false.)
107   endif
108
109
110!Step 1
111!create the Green function G_0 in imaginary time and save on file
112!it also calculates here the exchage energies
113
114
115   if(options%starting_point <= 1 .and. options%ending_point >= 1) then
116      if(.not.options%lpola_file .and. .not. options%lvcprim_file) then
117         call go_green(tf,options, qp)
118      endif
119   endif
120
121
122
123!Step 2
124!create the polarization in imaginary time and save on  file
125
126
127!loop on spin
128
129   do ispin=1,options%nspin
130
131      if(options%starting_point <= 2 .and. options%ending_point >=2 ) then
132         if(options%l_t_wannier) then
133            call calculate_compact_pola_lanczos(options,ispin)
134         endif
135
136      endif
137
138      if(options%starting_point <= 3 .and. options%ending_point >= 3 ) then
139         write(stdout,*) "*******************************"
140         write(stdout,*) "     RESTART FROM POINT 3"
141         write(stdout,*) "*******************************"
142!Step 3
143!FFT of polarization to imaginary frequency and save on file
144
145
146         call do_polarization_lanczos(tf,options,ispin)
147
148      endif
149   enddo
150
151  if(options%starting_point<=4  .and. options%ending_point >= 4) then
152
153
154!Step 3.1
155!calculate dresses interaction W, and save on file
156
157
158
159     write(stdout,*) 'Call go_dressed_w'
160
161
162     call go_dressed_w(options)
163
164
165
166     call write_quasi_particles(qp,options,.false.)
167
168
169
170  endif
171
172
173!Step 3.2
174!FFT of W to imaginary time and save on file
175  if(options%starting_point<=5  .and. options%ending_point >= 5) then
176
177
178
179     call read_quasi_particles(qp,options,.false.)
180
181     if(.not.  options%l_self_lanczos) then
182        write(stdout,*) 'Call FFT'
183        call go_fft_para2(tf, options)
184!if required do fft of gt structure
185        if(options%w_divergence==2) then
186           write(stdout,*) 'Go fft gt'
187           call initialize_gv_time(gt)
188           write(stdout,*) 'Go fft gt 1'
189           call read_gv_time(gt)
190           write(stdout,*) 'Go fft gt 1.5'
191           call fft_gv_time(gt,tf)
192           write(stdout,*) 'Go fft gt2'
193           call write_gv_time(gt)
194           write(stdout,*) 'Go fft gt3'
195           call free_memory_gv_time(gt)
196        endif
197     else
198        call do_reducible_pola(tf ,options)
199
200     endif
201  endif
202
203
204   if(options%starting_point <= 6  .and. options%ending_point >= 6) then
205!Step 4
206      write(stdout,*) '*******************************'
207      write(stdout,*) '     RESTART FROM POINT 6'
208      write(stdout,*) '*******************************'
209
210
211      if(options%n_real_axis>0) then
212         call initialize_w_expectation(we)
213         call create_w_expectation(we, tf, options)
214         call write_w_expectation(we)
215         call free_memory_w_expectation(we)
216      endif
217
218
219      if(.not.  options%l_self_lanczos) then
220!calculate the expectation value of Sigma in imaginary time and save on file
221         call create_self_ontime(tf, ss,options,qp)
222
223         if(options%lconduction.and. .not.options%lvcprim_file .and. .not.options%l_self_beta) then
224            if(.not.options%lcprim_file) then
225               call addconduction_self_ontime(ss, options)
226            else
227               call addconduction_self_ontime_file(ss, tf, options)
228            endif
229         endif
230         if(options%l_self_upper) then
231            call selfenergy_ontime_upper(ss, tf ,options)
232         endif
233
234         if(options%debug) call write_storage(tf,ss)
235
236         if(options%l_fft_timefreq) then
237            call fft_storage(ss)
238         else
239            if(tf%grid_fit==0) then
240               call fft_storage_grid(tf,ss)
241            else
242               call fft_storage_grid_fit(tf, ss)
243            endif
244         endif
245         if(options%debug) call write_storage(tf,ss)
246         call write_self_storage_ondisk(ss, options)
247      else
248!lanczos calculation of self-energy
249         if(options%n_real_axis==0) then
250            if(.not.options%l_self_time) then
251               call do_self_lanczos(ss, tf ,options)
252            else
253               if(.not.options%l_full) then
254                  call do_self_lanczos_time(ss, tf ,options,.false.,0.d0)
255               else
256                  call do_self_lanczos_full(ss, tf ,options,.false.,0.d0)
257               endif
258               call fft_storage_grid_fit(tf, ss)
259            endif
260            call write_self_storage_ondisk(ss, options)
261         else
262            call do_self_on_real(options,tf,ss,sr)
263            call write_self_on_real(sr,0)
264            call free_memory_self_on_real(sr)
265         endif
266         call write_self_storage_ondisk(ss, options)
267      endif
268
269
270   endif
271
272
273
274
275
276   if(options%starting_point <= 7 .and. options%ending_point >= 7) then
277!Step 7
278! fit self_energy with a multipole expansion
279      call read_self_storage_ondisk(ss, options)
280      call create_self_energy_fit( tf, se, ss, options,sr,.false.)
281      call mp_barrier( world_comm )
282      call print_fit_onfile(tf, se,ss)
283      call mp_barrier( world_comm )
284      call free_memory_self_storage(ss)
285      call mp_barrier( world_comm )
286      call create_quasi_particles(options,qp,se)
287      call mp_barrier( world_comm )
288      call write_self_expansion(se)
289      call free_memory_self_expansion(se)
290      call mp_barrier( world_comm )
291      call printout_quasi_particles(qp)
292   endif
293
294
295   if(options%starting_point <= 8 .and. options%ending_point >= 8) then
296!if the whole self_energy matrix has been calculate do use it for obtaining QPEs and QPAs
297      if(options%whole_s) then
298         call initialize_self_expansion(se)
299         call read_self_expansion(se)
300
301         call create_quasi_particles_off(options,qp,se)
302
303         call printout_quasi_particles_off(qp)
304         call free_memory_self_expansion(se)
305      endif
306   endif
307
308    if(options%starting_point <= 9 .and. options%ending_point >= 9) then
309!here does analytic continuation for contour integration
310       call initialize_w_expectation(we)
311       call initialize_w_poles(wp)
312       call read_w_expectation(we)
313       call create_w_poles(we,wp,options)
314       call write_w_poles(wp)
315       call free_memory_w_expectation(we)
316       call free_memory_w_poles(wp)
317    endif
318
319    if(options%starting_point <= 10 .and. options%ending_point >= 10) then
320!adds poles
321       call initialize_w_poles(wp)
322       call initialize_self_on_real(sr)
323       call read_w_poles(wp)
324       call read_self_on_real(sr,0)
325       !call self_on_real_print(sr)
326!NOT_TO_BE_INCLUDED_START
327       call do_contour(sr,wp,options)
328!NOT_TO_BE_INCLUDED_END
329       call write_self_on_real(sr,1)
330       call self_on_real_print(sr)
331       call free_memory_w_poles(wp)
332       call free_memory_self_on_real(sr)
333    endif
334    if(options%starting_point <= 11 .and. options%ending_point >= 11) then
335       call initialize_self_on_real(sr)
336       call read_self_on_real(sr,1)
337       call create_quasi_particle_on_real(options,qp,sr)
338       call printout_quasi_particles(qp)
339       call free_memory_self_on_real(sr)
340    endif
341!stops MPI
342    call free_memory_times_freqs(tf)
343    call free_memory_para_gww
344
345    call stop_run
346
347   stop
348   END PROGRAM gww
349
350