1! 2! Copyright (C) 2001-2020 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 wfcinit() 11 !---------------------------------------------------------------------------- 12 ! 13 ! ... This routine computes an estimate of the starting wavefunctions 14 ! ... from superposition of atomic wavefunctions and/or random wavefunctions. 15 ! ... It also open needed files or memory buffers 16 ! 17 USE io_global, ONLY : stdout, ionode, ionode_id 18 USE basis, ONLY : natomwfc, starting_wfc 19 USE bp, ONLY : lelfield 20 USE klist, ONLY : xk, nks, ngk, igk_k 21 USE control_flags, ONLY : io_level, lscf 22 USE fixed_occ, ONLY : one_atom_occupations 23 USE ldaU, ONLY : lda_plus_u, U_projection, wfcU, lda_plus_u_kind 24 USE lsda_mod, ONLY : lsda, current_spin, isk 25 USE io_files, ONLY : nwordwfc, nwordwfcU, iunhub, iunwfc,& 26 diropn, xmlfile, restart_dir 27 USE buffers, ONLY : open_buffer, close_buffer, get_buffer, save_buffer 28 USE uspp, ONLY : nkb, vkb 29 USE wavefunctions, ONLY : evc 30 USE wvfct, ONLY : nbnd, current_k 31 USE wannier_new, ONLY : use_wannier 32 USE pw_restart_new, ONLY : read_collected_wfc 33 USE mp, ONLY : mp_bcast, mp_sum 34 USE mp_images, ONLY : intra_image_comm 35 USE qexsd_module, ONLY : qexsd_readschema 36 USE qes_types_module, ONLY : output_type 37 USE qes_libs_module, ONLY : qes_reset 38 ! 39 IMPLICIT NONE 40 ! 41 INTEGER :: ik, ierr, exst_sum 42 LOGICAL :: exst, exst_mem, exst_file, opnd_file, twfcollect_file 43 CHARACTER (LEN=256) :: dirname 44 TYPE ( output_type ) :: output_obj 45 ! 46 ! 47 CALL start_clock( 'wfcinit' ) 48 ! 49 ! ... Orthogonalized atomic functions needed for DFT+U and other cases 50 ! 51 IF ( use_wannier .OR. one_atom_occupations ) CALL orthoatwfc ( use_wannier ) 52 IF ( lda_plus_u ) CALL orthoUwfc() 53 ! 54 ! ... open files/buffer for wavefunctions (nwordwfc set in openfil) 55 ! ... io_level > 1 : open file, otherwise: open buffer 56 ! 57 CALL open_buffer( iunwfc, 'wfc', nwordwfc, io_level, exst_mem, exst_file ) 58 ! 59 IF ( TRIM(starting_wfc) == 'file') THEN 60 ! Check whether all processors have found a file when opening a buffer 61 IF (exst_file) THEN 62 exst_sum = 0 63 ELSE 64 exst_sum = 1 65 END IF 66 CALL mp_sum (exst_sum, intra_image_comm) 67 ! 68 ! Check whether wavefunctions are collected (info in xml file) 69 dirname = restart_dir ( ) 70 IF (ionode) CALL qexsd_readschema ( xmlfile(), ierr, output_obj ) 71 CALL mp_bcast(ierr, ionode_id, intra_image_comm) 72 IF ( ierr <= 0 ) THEN 73 ! xml file is valid 74 IF (ionode) twfcollect_file = output_obj%band_structure%wf_collected 75 CALL mp_bcast(twfcollect_file, ionode_id, intra_image_comm) 76 CALL qes_reset ( output_obj ) 77 ELSE 78 ! xml file not found or not valid 79 twfcollect_file = .FALSE. 80 END IF 81 ! 82 IF ( twfcollect_file ) THEN 83 ! 84 DO ik = 1, nks 85 CALL read_collected_wfc ( dirname, ik, evc ) 86 CALL save_buffer ( evc, nwordwfc, iunwfc, ik ) 87 END DO 88 ! 89 ELSE IF ( exst_sum /= 0 ) THEN 90 ! 91 WRITE( stdout, '(5X,"Cannot read wfcs: file not found")' ) 92 IF (exst_file) THEN 93 CALL close_buffer(iunwfc, 'delete') 94 CALL open_buffer(iunwfc,'wfc', nwordwfc, io_level, exst_mem, exst_file) 95 END IF 96 starting_wfc = 'atomic+random' 97 ! 98 ELSE 99 ! 100 ! ... wavefunctions are read from file (or buffer) not here but 101 ! ...in routine c_bands. If however there is a single k-point, 102 ! ... c_bands doesn't read wavefunctions, so we read them here 103 ! ... (directly from file to avoid a useless buffer allocation) 104 ! 105 IF ( nks == 1 ) THEN 106 INQUIRE (unit = iunwfc, opened = opnd_file) 107 IF ( .NOT.opnd_file ) CALL diropn( iunwfc, 'wfc', 2*nwordwfc, exst ) 108 CALL davcio ( evc, 2*nwordwfc, iunwfc, nks, -1 ) 109 IF ( .NOT.opnd_file ) CLOSE ( UNIT=iunwfc, STATUS='keep' ) 110 END IF 111 END IF 112 END IF 113 ! 114 ! ... state what will happen 115 ! 116 IF ( TRIM(starting_wfc) == 'file' ) THEN 117 ! 118 WRITE( stdout, '(5X,"Starting wfcs from file")' ) 119 ! 120 ELSE IF ( starting_wfc == 'atomic' ) THEN 121 ! 122 IF ( natomwfc >= nbnd ) THEN 123 WRITE( stdout, '(5X,"Starting wfcs are ",I4," atomic wfcs")' ) natomwfc 124 ELSE 125 WRITE( stdout, '(5X,"Starting wfcs are ",I4," atomic + ", & 126 & I4," random wfcs")' ) natomwfc, nbnd-natomwfc 127 END IF 128 ! 129 ELSE IF ( TRIM(starting_wfc) == 'atomic+random' .AND. natomwfc > 0) THEN 130 ! 131 IF ( natomwfc >= nbnd ) THEN 132 WRITE( stdout, '(5X,"Starting wfcs are ",I4," randomized atomic wfcs")')& 133 natomwfc 134 ELSE 135 WRITE( stdout, '(5X,"Starting wfcs are ",I4," randomized atomic wfcs + "& 136 & ,I4," random wfcs")' ) natomwfc, nbnd-natomwfc 137 END IF 138 ! 139 ELSE 140 ! 141 WRITE( stdout, '(5X,"Starting wfcs are random")' ) 142 ! 143 END IF 144 ! 145 ! ... exit here if starting from file or for non-scf calculations. 146 ! ... In the latter case the starting wavefunctions are not 147 ! ... calculated here but just before diagonalization (to reduce I/O) 148 ! 149 IF ( ( .NOT. lscf .AND. .NOT. lelfield ) .OR. TRIM(starting_wfc) == 'file' ) THEN 150 ! 151 CALL stop_clock( 'wfcinit' ) 152 RETURN 153 ! 154 END IF 155 ! 156 ! ... calculate and write all starting wavefunctions to buffer 157 ! 158 DO ik = 1, nks 159 ! 160 ! ... Hpsi initialization: k-point index, spin, kinetic energy 161 ! 162 current_k = ik 163 IF ( lsda ) current_spin = isk(ik) 164 call g2_kin (ik) 165 ! 166 ! ... More Hpsi initialization: nonlocal pseudopotential projectors |beta> 167 ! 168 IF ( nkb > 0 ) CALL init_us_2( ngk(ik), igk_k(1,ik), xk(1,ik), vkb ) 169 ! 170 ! ... Needed for DFT+U 171 ! 172 IF ( nks > 1 .AND. lda_plus_u .AND. (U_projection .NE. 'pseudo') ) & 173 CALL get_buffer( wfcU, nwordwfcU, iunhub, ik ) 174 ! 175 ! DFT+U+V: calculate the phase factor at a given k point 176 ! 177 IF (lda_plus_u .AND. lda_plus_u_kind.EQ.2) CALL phase_factor(ik) 178 ! 179 ! ... calculate starting wavefunctions (calls Hpsi) 180 ! 181 CALL init_wfc ( ik ) 182 ! 183 ! ... write starting wavefunctions to file 184 ! 185 IF ( nks > 1 .OR. (io_level > 1) .OR. lelfield ) & 186 CALL save_buffer ( evc, nwordwfc, iunwfc, ik ) 187 ! 188 END DO 189 ! 190 CALL stop_clock( 'wfcinit' ) 191 RETURN 192 ! 193END SUBROUTINE wfcinit 194! 195!---------------------------------------------------------------------------- 196SUBROUTINE init_wfc ( ik ) 197 !---------------------------------------------------------------------------- 198 ! 199 ! ... This routine computes starting wavefunctions for k-point ik 200 ! 201 USE kinds, ONLY : DP 202 USE bp, ONLY : lelfield 203 USE becmod, ONLY : allocate_bec_type, deallocate_bec_type, & 204 bec_type, becp 205 USE constants, ONLY : tpi 206 USE basis, ONLY : natomwfc, starting_wfc 207 USE gvect, ONLY : g, gstart 208 USE klist, ONLY : xk, ngk, igk_k 209 USE wvfct, ONLY : nbnd, npwx, et 210 USE uspp, ONLY : nkb, okvan 211 USE noncollin_module, ONLY : npol 212 USE wavefunctions, ONLY : evc 213 USE random_numbers, ONLY : randy 214 USE mp_bands, ONLY : intra_bgrp_comm, inter_bgrp_comm, & 215 nbgrp, root_bgrp_id 216 USE mp, ONLY : mp_bcast 217 USE funct, ONLY : dft_is_hybrid, stop_exx 218 ! 219 IMPLICIT NONE 220 ! 221 INTEGER, INTENT(in) :: ik 222 ! 223 INTEGER :: ibnd, ig, ipol, n_starting_wfc, n_starting_atomic_wfc 224 LOGICAL :: lelfield_save 225 ! 226 REAL(DP) :: rr, arg 227 REAL(DP), ALLOCATABLE :: etatom(:) ! atomic eigenvalues 228 ! 229 COMPLEX(DP), ALLOCATABLE :: wfcatom(:,:,:) ! atomic wfcs for initialization 230 ! 231 ! 232 IF ( starting_wfc(1:6) == 'atomic' ) THEN 233 ! 234 n_starting_wfc = MAX( natomwfc, nbnd ) 235 n_starting_atomic_wfc = natomwfc 236 ! 237 ELSE IF ( starting_wfc == 'random' ) THEN 238 ! 239 n_starting_wfc = nbnd 240 n_starting_atomic_wfc = 0 241 ! 242 ELSE 243 ! 244 ! ...case 'file' should not be done here 245 ! 246 CALL errore ( 'init_wfc', & 247 'invalid value for startingwfc: ' // TRIM ( starting_wfc ) , 1 ) 248 ! 249 END IF 250 ! 251 ALLOCATE( wfcatom( npwx, npol, n_starting_wfc ) ) 252 ! 253 IF ( starting_wfc(1:6) == 'atomic' ) THEN 254 ! 255 CALL start_clock( 'wfcinit:atomic' ); !write(*,*) 'start wfcinit:atomic' ; FLUSH(6) 256 CALL atomic_wfc( ik, wfcatom ) 257 CALL stop_clock( 'wfcinit:atomic' ); !write(*,*) 'stop wfcinit:atomic' ; FLUSH(6) 258 ! 259 IF ( starting_wfc == 'atomic+random' .AND. & 260 n_starting_wfc == n_starting_atomic_wfc ) THEN 261 ! 262 ! ... in this case, introduce a small randomization of wavefunctions 263 ! ... to prevent possible "loss of states" 264 ! 265 DO ibnd = 1, n_starting_atomic_wfc 266 ! 267 DO ipol = 1, npol 268 ! 269 DO ig = 1, ngk(ik) 270 ! 271 rr = randy() 272 arg = tpi * randy() 273 ! 274 wfcatom(ig,ipol,ibnd) = wfcatom(ig,ipol,ibnd) * & 275 ( 1.0_DP + 0.05_DP * CMPLX( rr*COS(arg), rr*SIN(arg) ,kind=DP) ) 276 ! 277 END DO 278 ! 279 END DO 280 ! 281 END DO 282 ! 283 END IF 284 ! 285 END IF 286 ! 287 ! ... if not enough atomic wfc are available, 288 ! ... fill missing wfcs with random numbers 289 ! 290 DO ibnd = n_starting_atomic_wfc + 1, n_starting_wfc 291 ! 292 DO ipol = 1, npol 293 ! 294 wfcatom(:,ipol,ibnd) = (0.0_dp, 0.0_dp) 295 ! 296 DO ig = 1, ngk(ik) 297 ! 298 rr = randy() 299 arg = tpi * randy() 300 ! 301 wfcatom(ig,ipol,ibnd) = & 302 CMPLX( rr*COS( arg ), rr*SIN( arg ) ,kind=DP) / & 303 ( ( xk(1,ik) + g(1,igk_k(ig,ik)) )**2 + & 304 ( xk(2,ik) + g(2,igk_k(ig,ik)) )**2 + & 305 ( xk(3,ik) + g(3,igk_k(ig,ik)) )**2 + 1.0_DP ) 306 END DO 307 ! 308 END DO 309 ! 310 END DO 311 312 ! when band parallelization is active, the first band group distributes 313 ! the wfcs to the others making sure all bgrp have the same starting wfc 314 ! FIXME: maybe this should be done once evc are computed, not here? 315 ! 316 IF( nbgrp > 1 ) CALL mp_bcast( wfcatom, root_bgrp_id, inter_bgrp_comm ) 317 ! 318 ! ... Diagonalize the Hamiltonian on the basis of atomic wfcs 319 ! 320 ALLOCATE( etatom( n_starting_wfc ) ) 321 ! 322 ! ... Allocate space for <beta|psi> 323 ! 324 CALL allocate_bec_type ( nkb, n_starting_wfc, becp, intra_bgrp_comm ) 325 ! 326 ! ... the following trick is for electric fields with Berry's phase: 327 ! ... by setting lelfield = .false. one prevents the calculation of 328 ! ... electric enthalpy in the Hamiltonian (cannot be calculated 329 ! ... at this stage: wavefunctions at previous step are missing) 330 ! 331 lelfield_save = lelfield 332 lelfield = .FALSE. 333 ! 334 ! ... subspace diagonalization (calls Hpsi) 335 ! 336 IF ( dft_is_hybrid() ) CALL stop_exx() 337 CALL start_clock( 'wfcinit:wfcrot' ); !write(*,*) 'start wfcinit:wfcrot' ; FLUSH(6) 338 CALL rotate_wfc ( npwx, ngk(ik), n_starting_wfc, gstart, nbnd, wfcatom, npol, okvan, evc, etatom ) 339 CALL stop_clock( 'wfcinit:wfcrot' ); !write(*,*) 'stop wfcinit:wfcrot' ; FLUSH(6) 340 ! 341 lelfield = lelfield_save 342 ! 343 ! ... copy the first nbnd eigenvalues 344 ! ... eigenvectors are already copied inside routine rotate_wfc 345 ! 346 et(1:nbnd,ik) = etatom(1:nbnd) 347 ! 348 CALL deallocate_bec_type ( becp ) 349 DEALLOCATE( etatom ) 350 DEALLOCATE( wfcatom ) 351 ! 352 RETURN 353 ! 354END SUBROUTINE init_wfc 355