1!
2! Copyright (C) 2001-2016 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!-----------------------------------------------------------------------
9PROGRAM lr_main
10  !---------------------------------------------------------------------
11  !
12  ! This is the main driver of the TDDFPT code
13  ! for Absorption Spectroscopy.
14  ! It applys the Lanczos algorithm to the matrix
15  ! of equations coming from TDDFPT.
16  !
17  ! Brent Walker, ICTP, 2004
18  ! Dario Rocca, SISSA, 2006
19  ! Osman Baris Malcioglu, SISSA, 2008
20  ! Simone Binnie, SISSA, 2011
21  ! Xiaochuan Ge, SISSA, 2013
22  ! Iurii Timrov, SISSA, 2015
23  !
24  USE lr_lanczos,            ONLY : one_lanczos_step
25  USE io_global,             ONLY : stdout
26  USE kinds,                 ONLY : dp
27  USE lr_variables,          ONLY : restart, restart_step, itermax, lr_verbosity,  &
28                                  & evc1, evc1_old,norm0, charge_response, n_ipol, &
29                                  & d0psi, LR_iteration, LR_polarization, &
30                                  & plot_type, no_hxc, nbnd_total, project, F,R, &
31                                  & itermax_int, revc0, lr_io_level, code1
32  USE charg_resp,            ONLY : lr_calc_w_T, read_wT_beta_gamma_z, lr_project_init,&
33                                  & lr_dump_rho_tot_compat1, lr_dump_rho_tot_cube,&
34                                  & lr_dump_rho_tot_xyzd,lr_dump_rho_tot_xcrys,&
35                                  & lr_dump_rho_tot_pxyd,chi,lr_calc_R,w_t_norm0_store,&
36                                  & resonance_condition, lr_dump_rho, lr_calc_project
37  USE ions_base,             ONLY : tau,nat,atm,ityp
38  USE environment,           ONLY : environment_start
39  USE mp_global,             ONLY : nimage, mp_startup, inter_bgrp_comm, &
40                                    ibnd_start, ibnd_end
41  USE wvfct,                 ONLY : nbnd
42  USE wavefunctions,  ONLY : psic
43  USE check_stop,            ONLY : check_stop_now, check_stop_init
44  USE funct,                 ONLY : dft_is_hybrid
45  USE fft_base,              ONLY : dffts
46  USE uspp,                  ONLY : okvan
47  USE mp_bands,              ONLY : ntask_groups
48  !
49  IMPLICIT NONE
50  !
51  ! Local variables
52  !
53  INTEGER            :: ip,na,pol_index,ibnd
54  INTEGER            :: iter_restart,iteration
55  LOGICAL            :: rflag
56  COMPLEX(kind=dp)   :: temp
57  LOGICAL, EXTERNAL  :: test_restart
58  !
59  pol_index = 1
60  !
61#if defined(__MPI)
62  CALL mp_startup ( )
63#endif
64  !
65  CALL environment_start ( code1 )
66  !
67  CALL start_clock('lr_main')
68  !
69  IF (lr_verbosity > 5) THEN
70     WRITE(stdout,'("<lr_main>")')
71  ENDIF
72  !
73  ! Reading input file and PWSCF xml, some initialisation;
74  ! Read the input variables for TDDFPT;
75  ! Allocate space for all quantities already computed
76  ! by PWscf, and read them from the data file;
77  ! Define the tmp_dir directory.
78  !
79  CALL lr_readin ( )
80  !
81  ! Writing a summary of plugin variables
82  !
83  CALL plugin_summary()
84  !
85  CALL check_stop_init()
86  !
87  ! Initialisation
88  !
89  CALL lr_init_nfo()
90  !
91  ! Allocate the arrays
92  !
93  CALL lr_alloc_init()
94  !
95  ! Print a preamble info about the run
96  !
97  CALL lr_print_preamble()
98  !
99  IF ( ntask_groups > 1 ) WRITE(stdout,'(5X,"Task groups is activated...")' )
100  !
101  ! Read ground state wavefunctions from PWscf.
102  !
103  CALL lr_read_wf()
104  !
105  ! Band groups parallelization (if activated)
106  !
107  CALL divide(inter_bgrp_comm,nbnd,ibnd_start,ibnd_end)
108  !
109  ! Set up initial response orbitals (starting Lanczos vectors)
110  !
111  IF ( test_restart(1) ) THEN
112     CALL lr_read_d0psi()
113  ELSE
114     CALL lr_solve_e()
115  ENDIF
116  !
117  !do ip = 1, n_ipol
118  !  temp = wfc_dot(ibnd)
119  !enddo
120  !
121  DEALLOCATE( psic )
122  !
123  ! Projection analysis
124  !
125  IF (project) CALL lr_project_init()
126  !
127  ! Calculate a derivative of the XC potential
128  !
129  CALL lr_dv_setup()
130  !
131  ! S. Binnie: Write coordinates of the read atom (just in case).
132  !
133  IF (lr_verbosity > 1) THEN
134     !
135     WRITE(stdout,'(/,5X,"Positions of atoms in internal coordinates")')
136     !
137     DO na = 1, nat
138        WRITE(stdout,'(5X,A3,2X,3F15.8)') atm(ityp(na)), tau(1:3,na)
139     ENDDO
140     !
141  ENDIF
142  !
143  WRITE(stdout,'(/,5X,"LANCZOS LINEAR-RESPONSE SPECTRUM CALCULATION")')
144  WRITE(stdout,'(5X," ")')
145  WRITE(stdout,'(5x,"Number of Lanczos iterations = ",i6)') itermax
146  !
147  ! Lanczos loop where the real work happens
148  !
149  DO ip = 1, n_ipol
150     !
151     IF (n_ipol/=1) THEN
152        LR_polarization = ip
153        pol_index = LR_polarization
154     ENDIF
155     !
156     ! This is for the response of the charge density
157     !
158     IF (charge_response == 1) THEN
159        !
160        ! Read precalculated beta, gamma, and z.
161        !
162        CALL read_wT_beta_gamma_z()
163        CALL lr_calc_w_T()
164        !
165     ENDIF
166     !
167     ! Read the starting Lanczos vectors d0psi from the file which
168     ! was written above by lr_solve_e.
169     !
170     CALL lr_read_d0psi()
171     !
172     ! Normalization of the starting Lanczos vectors,
173     ! or reading of the data from the restart file.
174     !
175     IF (test_restart(2)) THEN
176        !
177        CALL lr_restart(iter_restart,rflag)
178        !
179        WRITE(stdout,'(/5x,"Restarting Lanczos loop",1x,i8)') LR_polarization
180        !
181     ELSE
182        !
183        ! The two starting Lanczos vectors are equal.
184        !
185        evc1(:,:,:,1) = d0psi(:,:,:,pol_index)
186        !
187        ! Xiaochuan Ge: The new structure of the Lanczos algorithm
188        ! does not need the normalisation of the starting Lanczos
189        ! vectors here.
190        !
191        !CALL lr_normalise( evc1(:,:,:,1), norm0(pol_index) )
192        !
193        evc1(:,:,:,2) = evc1(:,:,:,1)
194        !
195        evc1_old(:,:,:,1) = (0.0d0,0.0d0)
196        evc1_old(:,:,:,2) = (0.0d0,0.0d0)
197        !
198        iter_restart = 1
199        !
200        WRITE(stdout,'(/5x,"Starting Lanczos loop",1x,i8)') LR_polarization
201        !
202     ENDIF
203     !
204     ! d0psi = S * d0psi
205     ! This is needed in lr_lanczos for the dot product
206     ! in the calculation of the zeta-coefficients.
207     !
208     IF (okvan) CALL sd0psi()
209     !
210     ! Loop on the Lanczos iterations
211     !
212     lancz_loop1 : DO iteration = iter_restart, itermax
213        !
214        LR_iteration = iteration
215        !
216        WRITE(stdout,'(/5x,"Lanczos iteration:",1x,i6,3x,"Pol:",i1,i8)') LR_iteration, ip
217        !
218        CALL one_lanczos_step()
219        !
220        IF ( lr_io_level > 0 .and. (mod(LR_iteration,restart_step)==0 .or. &
221                           & LR_iteration==itermax .or. LR_iteration==1) ) &
222                           CALL lr_write_restart()
223        !
224        ! Check to see if the wall time limit has been exceeded.
225        ! If it has exit gracefully saving the last set of Lanczos
226        ! iterations.
227        !
228        IF ( check_stop_now() ) THEN
229           !
230           CALL lr_write_restart()
231           !
232           ! Deallocate PW variables.
233           !
234           CALL clean_pw( .FALSE. )
235           CALL stop_clock('lr_main')
236           CALL print_clock_lr()
237           CALL stop_lr( .FALSE. )
238           !
239        ENDIF
240        !
241     ENDDO lancz_loop1
242     !
243     IF (charge_response == 1) THEN
244        !
245        ! Write the apropriate charge density plot.
246        !
247        CALL lr_dump_rho (plot_type)
248        !
249     ENDIF
250     !
251     IF (project) THEN
252        !
253        ! Calculate projections onto virtual states if required.
254        !
255        CALL lr_calc_project(ip)
256        !
257     ENDIF
258     !
259  ENDDO
260  !
261  WRITE(stdout,'(5x,"End of Lanczos iterations")')
262  !
263  IF (project .AND. n_ipol == 3) THEN
264     !
265     !  Final projection and report (if required).
266     !
267     CALL lr_calc_project(4)
268     !
269  ENDIF
270  !
271  ! Deallocate PW variables
272  !
273  CALL clean_pw( .FALSE. )
274  !
275  WRITE(stdout,'(5x,"Finished linear response calculation...")')
276  !
277  CALL stop_clock('lr_main')
278  !
279  CALL print_clock_lr()
280  !
281  IF (lr_verbosity > 5) THEN
282     WRITE(stdout,'("<end of lr_main>")')
283  ENDIF
284  !
285  CALL stop_lr( .TRUE. )
286
287CONTAINS
288
289SUBROUTINE lr_print_preamble()
290
291    USE lr_variables,        ONLY : no_hxc, d0psi_rs
292    USE uspp,                ONLY : okvan
293    USE funct,               ONLY : dft_is_hybrid
294    USE martyna_tuckerman,   ONLY : do_comp_mt
295    USE control_flags,       ONLY : do_makov_payne
296
297    IMPLICIT NONE
298    !
299    WRITE( stdout, '(/5x,"=-----------------------------------------------------------------=")')
300    WRITE( stdout, '(/5x,"Please cite the TDDFPT project as:")')
301    WRITE( stdout, '(7x,"O. B. Malcioglu, R. Gebauer, D. Rocca, and S. Baroni,")')
302    WRITE( stdout, '(7x,"Comput. Phys. Commun. 182, 1744 (2011)")')
303    WRITE( stdout, '(5x,"and")' )
304    WRITE( stdout, '(7x,"X. Ge, S. J. Binnie, D. Rocca, R. Gebauer, and S. Baroni,")')
305    WRITE( stdout, '(7x,"Comput. Phys. Commun. 185, 2080 (2014)")')
306    WRITE( stdout, '(5x,"in publications and presentations arising from this work.")' )
307    WRITE( stdout, '(/5x,"=-----------------------------------------------------------------=")')
308    !
309    IF (okvan) WRITE( stdout, '(/5x,"Ultrasoft (Vanderbilt) Pseudopotentials")' )
310    !
311    IF (do_comp_mt) THEN
312       WRITE( stdout, '(/5x,"Martyna-Tuckerman periodic-boundary correction is used")' )
313    ELSEIF (do_makov_payne) THEN
314       WRITE( stdout, '(/5x,"WARNING! Makov-Payne periodic-boundary correction was activated in PWscf,",  &
315                      & /5x,"but it is of no use for TDDFPT. It just corrects the total energy in PWscf", &
316                      & /5x,"(post-processing correction) and nothing more, thus no effect on TDDFPT.", &
317                      & /5x,"You can try to use the Martyna-Tuckerman correction scheme instead.")' )
318    ENDIF
319    !
320    IF (no_hxc)  THEN
321       WRITE(stdout,'(5x,"No Hartree/Exchange/Correlation")')
322    ELSEIF (dft_is_hybrid() .AND. .NOT.d0psi_rs) THEN
323       WRITE(stdout, '(/5x,"Use of exact-exchange enabled. Note the EXX correction to the [H,X]", &
324                     & /5x,"commutator is NOT included hence the f-sum rule will be violated.",   &
325                     & /5x,"You can try to use the variable d0psi_rs=.true. (see the documentation).")' )
326    ENDIF
327    !
328    RETURN
329    !
330END SUBROUTINE lr_print_preamble
331
332END PROGRAM lr_main
333!-----------------------------------------------------------------------
334