1!
2! Copyright (C) 2002-2011 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!----------------------------------------------------------------------------
10SUBROUTINE init_run()
11  !----------------------------------------------------------------------------
12  !
13  ! ... this routine initialise the cp code and allocates (calling the
14  ! ... appropriate routines) the memory
15  !
16  USE kinds,                    ONLY : DP
17  USE control_flags,            ONLY : nbeg, nomore, lwf, iverbosity, iprint, &
18                                       ndr, ndw, tfor, tprnfor, tpre, ts_vdw, &
19                                       force_pairing
20  USE cp_electronic_mass,       ONLY : emass, emass_cutoff
21  USE ions_base,                ONLY : na, nax, nat, nsp, iforce, amass, cdms, ityp
22  USE ions_positions,           ONLY : tau0, taum, taup, taus, tausm, tausp, &
23                                       vels, velsm, velsp, fion, fionm
24  USE gvecw,                    ONLY : ngw, ngw_g, g2kin, g2kin_init
25  USE smallbox_gvec,            ONLY : ngb
26  USE gvect,                    ONLY : gstart, gg
27  USE fft_base,                 ONLY : dfftp, dffts
28  USE electrons_base,           ONLY : nspin, nbsp, nbspx, nupdwn, f
29  USE uspp,                     ONLY : nkb, vkb, deeq, becsum,nkbus
30  USE core,                     ONLY : rhoc
31  USE wavefunctions,     ONLY : c0_bgrp, cm_bgrp, phi_bgrp
32  USE ensemble_dft,             ONLY : tens, z0t
33  USE cg_module,                ONLY : tcg
34  USE electrons_base,           ONLY : nudx, nbnd
35  USE efield_module,            ONLY : tefield, tefield2
36  USE uspp_param,               ONLY : nhm
37  USE ions_nose,                ONLY : xnhp0, xnhpm, vnhp, nhpcl, nhpdim
38  USE cell_base,                ONLY : h, hold, hnew, velh, tpiba2, ibrav, &
39                                       alat, celldm, at, bg
40  USE cp_main_variables,        ONLY : lambda, lambdam, lambdap, ema0bg, &
41                                       sfac, eigr, taub, &
42                                       irb, eigrb, rhog, rhos, rhor,     &
43                                       acc, acc_this_run, wfill, &
44                                       edft, nfi, vpot, ht0, htm, iprint_stdout
45  USE cp_main_variables,        ONLY : allocate_mainvar, idesc
46  USE energies,                 ONLY : eself, enl, ekin, etot, enthal, ekincm
47  USE dener,                    ONLY : detot
48  USE time_step,                ONLY : dt2, delt, tps
49  USE electrons_nose,           ONLY : xnhe0, xnhem, vnhe
50  USE electrons_base,           ONLY : nbspx_bgrp
51  USE cell_nose,                ONLY : xnhh0, xnhhm, vnhh
52  USE funct,                    ONLY : dft_is_meta, dft_is_hybrid
53  USE metagga_cp,               ONLY : crosstaus, dkedtaus, gradwfc
54  !
55  USE efcalc,                   ONLY : clear_nbeg
56  USE local_pseudo,             ONLY : allocate_local_pseudo
57  USE cp_electronic_mass,       ONLY : emass_precond
58  USE wannier_subroutines,      ONLY : wannier_startup
59  USE cp_interfaces,            ONLY : readfile
60  USE ions_base,                ONLY : ions_cofmass
61  USE ensemble_dft,             ONLY : id_matrix_init, allocate_ensemble_dft, h_matrix_init
62  USE efield_module,            ONLY : allocate_efield, allocate_efield2
63  USE cg_module,                ONLY : allocate_cg
64  USE wannier_module,           ONLY : allocate_wannier
65  USE io_files,                 ONLY : tmp_dir, create_directory, restart_dir
66  USE io_global,                ONLY : ionode, stdout
67  USE printout_base,            ONLY : printout_base_init
68  USE wave_types,               ONLY : wave_descriptor_info
69  USE orthogonalize_base,       ONLY : mesure_diag_perf, mesure_mmul_perf
70  USE ions_base,                ONLY : ions_reference_positions, cdmi
71  USE mp_bands,                 ONLY : nbgrp
72  USE mp,                       ONLY : mp_barrier
73  USE wrappers
74  USE ldaU_cp
75  USE control_flags,            ONLY : lwfpbe0nscf         ! exx_wf related
76  USE wavefunctions,     ONLY : cv0                 ! exx_wf related
77  USE wannier_base,             ONLY : vnbsp               ! exx_wf related
78  !!!USE cp_restart,               ONLY : cp_read_wfc_Kong    ! exx_wf related
79  USE input_parameters,         ONLY : ref_cell
80  USE cell_base,                ONLY : ref_tpiba2, init_tpiba2
81  USE tsvdw_module,             ONLY : tsvdw_initialize
82  USE exx_module,               ONLY : exx_initialize
83  !
84  IMPLICIT NONE
85  !
86  INTEGER            :: i
87  CHARACTER(LEN=256) :: dirname
88  REAL(DP)           :: a1(3), a2(3), a3(3)
89  LOGICAL            :: ftest
90  !
91  !
92  CALL start_clock( 'initialize' )
93  !
94  ! ... initialize directories
95  !
96
97  IF( nbeg < 0 ) THEN
98     CALL create_directory( tmp_dir )
99  END IF
100  !
101  CALL plugin_initialization()
102
103  IF( nbgrp > 1 .AND. force_pairing ) &
104     CALL errore( ' init_run ', ' force_pairing with parallelization over bands not implemented yet ', 1 )
105  !
106  ! ... Open files containing MD information
107  !
108  CALL printout_base_init( )
109  !
110  ! ... Create main restart directory
111  !
112  dirname = restart_dir( ndw )
113  CALL create_directory( dirname )
114  !
115  ! ... initialize g-vectors, fft grids
116  ! ... The number of g-vectors are based on the input celldm!
117  !
118  CALL init_dimensions()
119  !
120  ! ... initialization of plugin variables and arrays
121  !
122  CALL plugin_init_base()
123  !
124  ! ... initialize atomic positions and cell
125  !
126  CALL init_geometry()
127  !
128  ! ... mesure performances of parallel routines
129  !
130  CALL mesure_mmul_perf( nudx )
131  !
132  CALL mesure_diag_perf( nudx )
133  !
134  IF ( lwf ) CALL clear_nbeg( nbeg )
135  !
136  !=======================================================================
137  !     allocate and initialize local and nonlocal potentials
138  !=======================================================================
139  !
140  CALL allocate_local_pseudo( dffts%ngm, nsp )
141  !
142  CALL nlinit()
143  !
144  !=======================================================================
145  !     allocation of all arrays not already allocated in init and nlinit
146  !=======================================================================
147  !
148  CALL allocate_mainvar( ngw, ngw_g, ngb, dffts%ngm, dfftp%ngm, dfftp%nr1,dfftp%nr2,dfftp%nr3, dfftp%nr1x, &
149                         dfftp%nr2x, dfftp%my_nr3p, dfftp%nnr, dffts%nnr, nat, nax, nsp,   &
150                         nspin, nbsp, nbspx, nupdwn, nkb, gstart, nudx, &
151                         tpre, nbspx_bgrp )
152  !
153  !=======================================================================
154  !     Initialization of the TS-vdW code (RAD)
155  !=======================================================================
156  !
157  IF (ts_vdw) CALL tsvdw_initialize()
158  !
159  !=======================================================================
160  !     Initialization of the exact exchange code (exx_module)
161  !=======================================================================
162  !exx_wf related
163  IF ( dft_is_hybrid() .AND. lwf ) THEN
164    !
165    CALL exx_initialize()
166    !
167  END IF
168  !
169  !  initialize wave functions descriptors and allocate wf
170  !
171  IF(lwfpbe0nscf) ALLOCATE(cv0( ngw, vnbsp ) )   ! Lingzhu Kong
172  ALLOCATE( c0_bgrp( ngw, nbspx ) )
173  ALLOCATE( cm_bgrp( ngw, nbspx ) )
174  ALLOCATE( phi_bgrp( ngw, nbspx ) )
175  !
176  IF ( iverbosity > 1 ) THEN
177     !
178     CALL wave_descriptor_info( wfill, 'wfill', stdout )
179     !
180  END IF
181  !
182  ! Depending on the verbosity set the frequency of
183  ! verbose information to stdout
184  !
185  IF( iverbosity < 0 ) iprint_stdout = 100 * iprint
186  IF( iverbosity ==0 .OR. iverbosity == 1 ) iprint_stdout = 10 * iprint
187  IF( iverbosity > 1 ) iprint_stdout = iprint
188  !
189  acc          = 0.D0
190  acc_this_run = 0.D0
191  !
192  edft%ent  = 0.D0
193  edft%esr  = 0.D0
194  edft%evdw = 0.D0
195  edft%ekin = 0.D0
196  edft%enl  = 0.D0
197  edft%etot = 0.D0
198  !
199  ALLOCATE( becsum(  nhm*(nhm+1)/2, nat, nspin ) )
200  ALLOCATE( deeq( nhm, nhm, nat, nspin ) )
201  !
202  ALLOCATE( vkb( ngw, nkb ) )
203  !
204  IF ( dft_is_meta() .AND. tens ) &
205     CALL errore( ' init_run ', 'ensemble_dft not implemented for metaGGA', 1 )
206  !
207  IF ( dft_is_meta() .AND. nbgrp > 1 ) &
208     CALL errore( ' init_run ', 'band parallelization not implemented for metaGGA', 1 )
209  !
210  IF ( dft_is_meta() .AND. tpre ) THEN
211     !
212     ALLOCATE( crosstaus( dffts%nnr, 6, nspin ) )
213     ALLOCATE( dkedtaus(  dffts%nnr, 3, 3, nspin ) )
214     ALLOCATE( gradwfc(   dffts%nnr, 3 ) )
215     !
216     if (nspin.ne.1) &
217       CALL errore( ' init_run ', 'spin-polarized stress not implemented for metaGGA', 1 )
218     !
219  END IF
220  !
221  IF ( lwf ) THEN
222     IF( nbgrp > 1 ) &
223        CALL errore( ' init_run ', ' wannier with band parallelization not implemented ', 1 )
224     CALL allocate_wannier( nbsp, dffts%nnr, nspin, dfftp%ngm )
225  END IF
226  !
227  IF ( tens .OR. tcg ) THEN
228     IF( nbgrp > 1 ) &
229        CALL errore( ' init_run ', ' ensemble_dft with band parallelization not implemented ', 1 )
230     CALL allocate_ensemble_dft( nkb, nbsp, ngw, nudx, nspin, nbspx, &
231                                 dffts%nnr, nat, idesc )
232  END IF
233  !
234  IF ( tcg ) THEN
235     CALL allocate_cg( ngw, nbspx,nkbus )
236  END IF
237  !
238  IF ( tefield ) THEN
239     IF( nbgrp > 1 ) &
240        CALL errore( ' init_run ', ' efield with band paralleliztion not implemented ', 1 )
241     CALL allocate_efield( ngw, ngw_g, nbspx, nhm, nax, nsp )
242  END IF
243  IF ( tefield2 ) THEN
244     IF( nbgrp > 1 ) &
245        CALL errore( ' init_run ', ' efield with band paralleliztion not implemented ', 1 )
246     CALL allocate_efield2( ngw, nbspx, nhm, nax, nsp )
247  END IF
248  !
249  IF ( ALLOCATED( deeq ) ) deeq(:,:,:,:) = 0.D0
250  !
251  IF ( ALLOCATED( lambda  ) ) lambda  = 0.D0
252  IF ( ALLOCATED( lambdam ) ) lambdam = 0.D0
253  !
254  taum  = tau0
255  taup  = 0.D0
256  tausm = taus
257  tausp = 0.D0
258  vels  = 0.D0
259  velsm = 0.D0
260  velsp = 0.D0
261  !
262  hnew = h
263  !
264  IF(lwfpbe0nscf) cv0=( 0.D0, 0.D0 )    ! Lingzhu Kong
265  cm_bgrp  = ( 0.D0, 0.D0 )
266  c0_bgrp  = ( 0.D0, 0.D0 )
267  phi_bgrp = ( 0.D0, 0.D0 )
268  !
269  IF ( tens ) then
270     CALL id_matrix_init( idesc, nspin )
271     CALL h_matrix_init( idesc, nspin )
272  ENDIF
273  !
274  a1(:)=at(:,1)*alat; a2(:)=at(:,2)*alat; a3(:)=at(:,3)*alat
275  IF ( lwf ) CALL wannier_startup( ibrav, alat, a1, a2, a3, &
276                                   bg(:,1), bg(:,2), bg(:,3) )
277  !
278  ! ... Calculate: ema0bg = ecutmass /  MAX( 1.0d0, (2pi/alat)^2 * |G|^2 )
279  !
280  IF ( ref_cell ) THEN
281    WRITE( stdout,'(/,3X,"Reference cell parameters are used in electron mass preconditioning")' )
282    WRITE( stdout,'(3X,"ref_tpiba2=",F14.8)' ) ref_tpiba2
283    CALL g2kin_init( gg, ref_tpiba2 )
284    CALL emass_precond( ema0bg, g2kin, ngw, ref_tpiba2, emass_cutoff )
285    WRITE( stdout,'(3X,"current_tpiba2=",F14.8)' ) tpiba2
286    CALL g2kin_init( gg, tpiba2 )
287  ELSE
288    WRITE( stdout,'(/,3X,"Cell parameters from input file are used in electron mass preconditioning")' )
289    WRITE( stdout,'(3X,"init_tpiba2=",F14.8)' ) init_tpiba2
290    CALL g2kin_init( gg, init_tpiba2 )
291    CALL emass_precond( ema0bg, g2kin, ngw, init_tpiba2, emass_cutoff )
292    !WRITE( stdout,'(3X,"current_tpiba2=",F14.8)' ) tpiba2 !BS : DEBUG
293    CALL g2kin_init( gg, tpiba2 )
294  END IF
295  !
296  CALL print_legend( )
297  !
298  CALL ldaU_init()
299  !
300  IF ( nbeg < 0 ) THEN
301     !
302     !======================================================================
303     !     Initialize from scratch nbeg = -1
304     !======================================================================
305     !
306     nfi = 0
307     !
308     CALL from_scratch( )
309     !
310  ELSE
311     !
312     !======================================================================
313     !     nbeg = 0, nbeg = 1
314     !======================================================================
315     !
316     !======================================================================
317     ! Kong, read the valence orbitals
318     IF(lwfpbe0nscf) THEN
319       !!! CALL cp_read_wfc_Kong( 36, tmp_dir, 1, 1, 1, 1, cv0, 'v' )
320       CALL errore( 'init_run', 'cp_read_wfc_Kong no longer available', 1)
321     ENDIF
322     !======================================================================
323     i = 1
324     CALL start_clock( 'init_readfile' )
325     CALL readfile( i, h, hold, nfi, c0_bgrp, cm_bgrp, taus,   &
326                    tausm, vels, velsm, acc, lambda, lambdam, xnhe0, xnhem, &
327                    vnhe, xnhp0, xnhpm, vnhp,nhpcl,nhpdim,ekincm, xnhh0, xnhhm,&
328                    vnhh, velh, fion, tps, z0t, f )
329     !
330     CALL from_restart( )
331     !
332     CALL stop_clock( 'init_readfile' )
333     !
334  END IF
335  !
336  !=======================================================================
337  !     restart with new averages and nfi=0
338  !=======================================================================
339  !
340  ! ... reset some variables if nbeg < 0
341  ! ... ( new simulation or step counter reset to 0 )
342  !
343  IF ( nbeg <= 0 ) THEN
344     !
345     acc = 0.D0
346     nfi = 0
347     !
348  END IF
349  !
350  IF ( .NOT. tfor .AND. .NOT. tprnfor ) fion(:,:) = 0.D0
351  !
352  nomore = nomore + nfi
353  !
354  !  Set center of mass for scaled coordinates
355  !
356  CALL ions_cofmass( taus, amass, nat, ityp, cdms )
357  !
358  IF ( nbeg <= 0 .OR. lwf ) THEN
359     !
360     CALL ions_reference_positions( tau0 ) ! BS: screws up msd calculation for lwf ...
361     !
362  END IF
363  !
364  CALL stop_clock( 'initialize' )
365  !
366  RETURN
367  !
368END SUBROUTINE init_run
369