1! 2! Copyright (C) 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!---------------------------------------------------------------------------- 9MODULE pw_restart_new 10!---------------------------------------------------------------------------- 11 ! 12 ! ... New PWscf I/O using xml schema and (optionally) hdf5 binaries 13 ! ... Parallel execution: the xml file is written by one processor only 14 ! ... ("ionode_id"), read by all processors ; 15 ! ... the wavefunction files are written / read by one processor per pool, 16 ! ... collected on / distributed to all other processors in pool 17 ! 18 USE kinds, ONLY: dp 19 USE qes_types_module, ONLY : output_type, parallel_info_type, & 20 general_info_type, input_type, gateInfo_type, dipoleOutput_type, & 21 BerryPhaseOutput_type, hybrid_type, vdw_type, dftU_type, smearing_type 22 USE qes_write_module, ONLY: qes_write 23 USE qes_reset_module, ONLY: qes_reset 24 USE qes_bcast_module,ONLY : qes_bcast 25 USE qexsd_module, ONLY: qexsd_xf, qexsd_openschema, qexsd_closeschema, & 26 qexsd_readschema 27 USE qexsd_input, ONLY : qexsd_input_obj, qexsd_init_k_points_ibz, & 28 qexsd_init_occupations, qexsd_init_smearing 29 USE qexsd_init, ONLY: qexsd_init_convergence_info, qexsd_init_algorithmic_info, & 30 qexsd_init_atomic_species, qexsd_init_atomic_structure, & 31 qexsd_init_symmetries, qexsd_init_basis_set, qexsd_init_dft, & 32 qexsd_init_magnetization,qexsd_init_band_structure, & 33 qexsd_init_dipole_info, qexsd_init_total_energy, & 34 qexsd_init_vdw, qexsd_init_forces, qexsd_init_stress, & 35 qexsd_init_outputElectricField, qexsd_init_outputPBC, & 36 qexsd_init_gate_info, qexsd_init_hybrid, qexsd_init_dftU, & 37 qexsd_occ_obj, qexsd_bp_obj, qexsd_start_k_obj 38 USE qexsd_copy, ONLY : qexsd_copy_parallel_info, & 39 qexsd_copy_algorithmic_info, qexsd_copy_atomic_species, & 40 qexsd_copy_atomic_structure, qexsd_copy_symmetry, & 41 qexsd_copy_basis_set, qexsd_copy_dft, qexsd_copy_efield, & 42 qexsd_copy_band_structure, qexsd_copy_magnetization, & 43 qexsd_copy_kpoints 44 USE io_global, ONLY : ionode, ionode_id 45 USE io_files, ONLY : iunpun, xmlfile 46 ! 47 IMPLICIT NONE 48 ! 49 CHARACTER(LEN=6), EXTERNAL :: int_to_char 50 PRIVATE 51 PUBLIC :: pw_write_schema, pw_write_binaries 52 PUBLIC :: read_xml_file, read_collected_wfc 53 ! 54 CONTAINS 55 !------------------------------------------------------------------------ 56 SUBROUTINE pw_write_schema( only_init, wf_collect ) 57 !------------------------------------------------------------------------ 58 ! 59 ! only_init = T write only variables that are known after the 60 ! initial steps of initialization (e.g. structure) 61 ! = F write the complete xml file 62 ! wf_collect = T if final wavefunctions in portable format are written, 63 ! F if wavefunctions are either not written or are written 64 ! in binary non-portable form (for checkpointing) 65 ! NB: wavefunctions are not written here in any case 66 ! 67 USE control_flags, ONLY : istep, conv_ions, & 68 lscf, gamma_only, & 69 tqr, tq_smoothing, tbeta_smoothing, & 70 noinv, smallmem, & 71 llondon, lxdm, ts_vdw, scf_error, n_scf_steps 72 USE constants, ONLY : e2 73 USE realus, ONLY : real_space 74 USE uspp, ONLY : okvan 75 USE paw_variables, ONLY : okpaw 76 USE uspp_param, ONLY : upf 77 USE cell_base, ONLY : at, bg, alat, ibrav 78 USE ions_base, ONLY : nsp, ityp, atm, nat, tau, zv, amass 79 USE noncollin_module, ONLY : noncolin, npol 80 USE io_files, ONLY : psfile, pseudo_dir 81 USE klist, ONLY : nks, nkstot, xk, ngk, wk, & 82 lgauss, ngauss, smearing, degauss, nelec, & 83 two_fermi_energies, nelup, neldw, tot_charge, ltetra 84 USE start_k, ONLY : nk1, nk2, nk3, k1, k2, k3, & 85 nks_start, xk_start, wk_start 86 USE gvect, ONLY : ngm, ngm_g, g 87 USE fft_base, ONLY : dfftp 88 USE basis, ONLY : natomwfc 89 USE gvecs, ONLY : ngms_g, dual 90 USE fft_base, ONLY : dffts 91 USE wvfct, ONLY : npwx, et, wg, nbnd 92 USE ener, ONLY : ef, ef_up, ef_dw, vtxc, etxc, ewld, etot, & 93 ehart, eband, demet, edftd3, elondon, exdm 94 USE tsvdw_module, ONLY : EtsvdW 95 USE gvecw, ONLY : ecutwfc 96 USE fixed_occ, ONLY : tfixed_occ, f_inp 97 USE ldaU, ONLY : lda_plus_u, lda_plus_u_kind, U_projection, & 98 Hubbard_lmax, Hubbard_l, Hubbard_U, Hubbard_J, & 99 Hubbard_l_back, Hubbard_l1_back, Hubbard_V, & 100 Hubbard_alpha, Hubbard_alpha_back, nsg, & 101 Hubbard_J0, Hubbard_beta, Hubbard_U_back, & 102 is_hubbard, is_hubbard_back, backall 103 USE spin_orb, ONLY : lspinorb, domag 104 USE symm_base, ONLY : nrot, nsym, invsym, s, ft, irt, & 105 t_rev, sname, time_reversal, no_t_rev,& 106 spacegroup 107 USE lsda_mod, ONLY : nspin, isk, lsda, starting_magnetization, magtot, absmag 108 USE noncollin_module, ONLY : angle1, angle2, i_cons, mcons, bfield, magtot_nc, & 109 lambda 110 USE funct, ONLY : get_dft_short, get_nonlocc_name, dft_is_nonlocc 111 USE scf, ONLY : rho 112 USE force_mod, ONLY : lforce, sumfor, force, sigma, lstres 113 USE extfield, ONLY : tefield, dipfield, edir, etotefield, & 114 emaxpos, eopreg, eamp, el_dipole, ion_dipole,& 115 gate, zgate, relaxz, block, block_1,& 116 block_2, block_height, etotgatefield ! TB 117 USE mp, ONLY : mp_sum 118 USE mp_bands, ONLY : intra_bgrp_comm 119 USE funct, ONLY : get_exx_fraction, dft_is_hybrid, & 120 get_gau_parameter, & 121 get_screening_parameter, exx_is_active 122 USE exx_base, ONLY : x_gamma_extrapolation, nq1, nq2, nq3, & 123 exxdiv_treatment, yukawa, ecutvcut 124 USE exx, ONLY : ecutfock, local_thr 125 USE london_module, ONLY : scal6, lon_rcut, c6_i 126 USE xdm_module, ONLY : xdm_a1=>a1i, xdm_a2=>a2i 127 USE tsvdw_module, ONLY : vdw_isolated, vdw_econv_thr 128 USE input_parameters, ONLY : verbosity, calculation, ion_dynamics, starting_ns_eigenvalue, & 129 vdw_corr, london, k_points, assume_isolated, & 130 input_parameters_occupations => occupations, dftd3_threebody, & 131 dftd3_version 132 USE bp, ONLY : lelfield, lberry, el_pol, ion_pol 133 ! 134 USE rap_point_group, ONLY : elem, nelem, name_class 135 USE rap_point_group_so, ONLY : elem_so, nelem_so, name_class_so 136 USE bfgs_module, ONLY : bfgs_get_n_iter 137 USE fcp_variables, ONLY : lfcpopt, lfcpdyn, fcp_mu 138 USE control_flags, ONLY : conv_elec, conv_ions, ldftd3, do_makov_payne 139 USE Coul_cut_2D, ONLY : do_cutoff_2D 140 USE esm, ONLY : do_comp_esm 141 USE martyna_tuckerman, ONLY : do_comp_mt 142 USE run_info, ONLY : title 143 ! 144 IMPLICIT NONE 145 ! 146 LOGICAL, INTENT(IN) :: only_init, wf_collect 147 ! 148 CHARACTER(LEN=26) :: dft_name 149 CHARACTER(LEN=8) :: smearing_loc 150 CHARACTER(LEN=8), EXTERNAL :: schema_smearing 151 INTEGER :: i, ig, ngg, ipol 152 INTEGER :: npwx_g, ispin 153 INTEGER, ALLOCATABLE :: ngk_g(:) 154 INTEGER :: iclass, isym, ielem 155 CHARACTER(LEN=15) :: symop_2_class(48) 156 LOGICAL :: opt_conv_ispresent, dft_is_vdw, empirical_vdw 157 INTEGER :: n_opt_steps, n_scf_steps_, h_band 158 REAL(DP),TARGET :: h_energy 159 TYPE(gateInfo_type),TARGET :: gate_info_temp 160 TYPE(gateInfo_type),POINTER :: gate_info_ptr 161 TYPE(dipoleOutput_type),TARGET :: dipol_obj 162 TYPE(dipoleOutput_type),POINTER :: dipol_ptr 163 TYPE(BerryPhaseOutput_type), POINTER :: bp_obj_ptr 164 TYPE(hybrid_type), POINTER :: hybrid_obj 165 TYPE(vdW_type), POINTER :: vdw_obj 166 TYPE(dftU_type), POINTER :: dftU_obj 167 REAL(DP), TARGET :: lumo_tmp, ef_targ, dispersion_energy_term 168 REAL(DP), POINTER :: lumo_energy, ef_point 169 REAL(DP), ALLOCATABLE :: ef_updw(:) 170 ! 171 ! 172 ! 173 TYPE(output_type) :: output_obj 174 REAL(DP),POINTER :: degauss_, demet_, efield_corr, potstat_corr, gatefield_corr 175 LOGICAL, POINTER :: optimization_has_converged 176 LOGICAL, TARGET :: conv_opt 177 LOGICAL :: scf_has_converged 178 INTEGER :: itemp = 1 179 REAL(DP),ALLOCATABLE :: london_c6_(:), bp_el_pol(:), bp_ion_pol(:), U_opt(:), J0_opt(:), alpha_opt(:), & 180 J_opt(:,:), beta_opt(:), U_back_opt(:), alpha_back_opt(:) 181 INTEGER,ALLOCATABLE :: Hubbard_l_back_opt(:), Hubbard_l1_back_opt(:) 182 LOGICAL, ALLOCATABLE :: backall_opt(:) 183 CHARACTER(LEN=3),ALLOCATABLE :: species_(:) 184 CHARACTER(LEN=20),TARGET :: dft_nonlocc_ 185 INTEGER,TARGET :: dftd3_version_ 186 CHARACTER(LEN=20),TARGET :: vdw_corr_, pbc_label 187 CHARACTER(LEN=20),POINTER :: non_local_term_pt, vdw_corr_pt 188 REAL(DP),TARGET :: temp(20), lond_rcut_, lond_s6_, ts_vdw_econv_thr_, xdm_a1_, xdm_a2_, ectuvcut_,& 189 scr_par_, loc_thr_ 190 REAL(DP),POINTER :: vdw_term_pt, ts_thr_pt, london_s6_pt, london_rcut_pt, xdm_a1_pt, xdm_a2_pt, & 191 ts_vdw_econv_thr_pt, ectuvcut_opt, scr_par_opt, loc_thr_p, h_energy_ptr 192 LOGICAL,TARGET :: dftd3_threebody_, ts_vdw_isolated_, domag_ 193 LOGICAL,POINTER :: ts_isol_pt, dftd3_threebody_pt, ts_vdw_isolated_pt, domag_opt 194 INTEGER,POINTER :: dftd3_version_pt 195 TYPE(smearing_type),TARGET :: smear_obj 196 TYPE(smearing_type),POINTER:: smear_obj_ptr 197 198 NULLIFY( degauss_, demet_, efield_corr, potstat_corr, gatefield_corr) 199 NULLIFY( gate_info_ptr, dipol_ptr, bp_obj_ptr, hybrid_obj, vdw_obj, dftU_obj, lumo_energy, ef_point) 200 NULLIFY ( optimization_has_converged, non_local_term_pt, & 201 vdw_corr_pt, vdw_term_pt, ts_thr_pt, london_s6_pt, london_rcut_pt, & 202 xdm_a1_pt, xdm_a2_pt, ts_vdw_econv_thr_pt, ts_isol_pt, & 203 dftd3_threebody_pt, ts_vdw_isolated_pt, dftd3_version_pt ) 204 NULLIFY ( ectuvcut_opt, scr_par_opt, loc_thr_p, h_energy_ptr, smear_obj_ptr, domag_opt) 205 206 ! 207 ! Global PW dimensions need to be properly computed, reducing across MPI tasks 208 ! If local PW dimensions are not available, set to 0 209 ! 210 ALLOCATE( ngk_g( nkstot ) ) 211 ngk_g(:) = 0 212 IF ( ALLOCATED (ngk) ) THEN 213 ngk_g(1:nks) = ngk(:) 214 CALL mp_sum( ngk_g(1:nks), intra_bgrp_comm ) 215 CALL ipoolrecover( ngk_g, 1, nkstot, nks ) 216 END IF 217 ! BEWARE: only the first pool has ngk_g for all k-points 218 ! 219 ! ... compute the maximum number of G vector among all k points 220 ! 221 npwx_g = MAXVAL( ngk_g(1:nkstot) ) 222 ! 223 ! XML descriptor 224 ! 225 IF ( ionode ) THEN 226 ! 227 ! ... here we init the variables and finally write them to file 228 ! 229!------------------------------------------------------------------------------- 230! ... HEADER 231!------------------------------------------------------------------------------- 232 ! 233 output_obj%tagname="output" 234 output_obj%lwrite = .TRUE. 235 output_obj%lread = .TRUE. 236 ! 237!------------------------------------------------------------------------------- 238! ... CONVERGENCE_INFO 239!------------------------------------------------------------------------------- 240 SELECT CASE (TRIM( calculation )) 241 CASE ( "relax","vc-relax" ) 242 conv_opt = conv_ions 243 optimization_has_converged => conv_opt 244 IF (TRIM( ion_dynamics) == 'bfgs' ) THEN 245 n_opt_steps = bfgs_get_n_iter('bfgs_iter ') 246 ELSE 247 n_opt_steps = istep 248 END IF 249 scf_has_converged = conv_elec 250 n_scf_steps_ = n_scf_steps 251 CASE ("nscf", "bands" ) 252 n_opt_steps = 0 253 scf_has_converged = .FALSE. 254 n_scf_steps_ = 1 255 CASE default 256 n_opt_steps = 0 257 scf_has_converged = conv_elec 258 n_scf_steps_ = n_scf_steps 259 END SELECT 260 ! 261 call qexsd_init_convergence_info(output_obj%convergence_info, & 262 SCf_HAS_CONVERGED = scf_has_converged, & 263 OPTIMIZATION_HAS_CONVERGED = optimization_has_converged,& 264 N_SCF_STEPS = n_scf_steps_, SCF_ERROR=scf_error/e2,& 265 N_OPT_STEPS = n_opt_steps, GRAD_NORM = sumfor) 266 output_obj%convergence_info_ispresent = .TRUE. 267 ! 268 269!------------------------------------------------------------------------------- 270! ... ALGORITHMIC_INFO 271!------------------------------------------------------------------------------- 272 ! 273 CALL qexsd_init_algorithmic_info(output_obj%algorithmic_info, & 274 REAL_SPACE_BETA = real_space, REAL_SPACE_Q=tqr , USPP=okvan, PAW=okpaw) 275 ! 276!------------------------------------------------------------------------------- 277! ... ATOMIC_SPECIES 278!------------------------------------------------------------------------------- 279 ! 280 ! while amass's are always present, starting_mag should not be passed 281 ! for nspin==1 or contrained magnetization calculations 282 ! 283 IF (noncolin) THEN 284 CALL qexsd_init_atomic_species(output_obj%atomic_species, nsp, atm, psfile, & 285 amass, STARTING_MAGNETIZATION = starting_magnetization, & 286 ANGLE1=angle1, ANGLE2=angle2) 287 ELSE IF (nspin==2) THEN 288 CALL qexsd_init_atomic_species(output_obj%atomic_species, nsp, atm, psfile, & 289 amass, STARTING_MAGNETIZATION=starting_magnetization) 290 ELSE 291 CALL qexsd_init_atomic_species(output_obj%atomic_species, nsp, atm,psfile, & 292 amass) 293 END IF 294 output_obj%atomic_species%pseudo_dir = TRIM(pseudo_dir) 295 output_obj%atomic_species%pseudo_dir_ispresent = .TRUE. 296 ! 297!------------------------------------------------------------------------------- 298! ... ATOMIC_STRUCTURE 299!------------------------------------------------------------------------------- 300 ! 301 CALL qexsd_init_atomic_structure(output_obj%atomic_structure, nsp, atm, ityp, & 302 nat, alat*tau, alat, alat*at(:,1), alat*at(:,2), alat*at(:,3), ibrav) 303 ! 304!------------------------------------------------------------------------------- 305! ... SYMMETRIES 306!------------------------------------------------------------------------------- 307 ! 308 symop_2_class="not found" 309 IF (TRIM (verbosity) == 'medium' .OR. TRIM(verbosity) == 'high') THEN 310 IF ( noncolin ) THEN 311 symmetries_so_loop:DO isym = 1, nrot 312 classes_so_loop:DO iclass = 1, 24 313 elements_so_loop:DO ielem=1, nelem_so(iclass) 314 IF ( elem_so(ielem,iclass) == isym) THEN 315 symop_2_class(isym) = name_class_so(iclass) 316 EXIT symmetries_so_loop 317 END IF 318 END DO elements_so_loop 319 END DO classes_so_loop 320 END DO symmetries_so_loop 321 ! 322 ELSE 323 symmetries_loop:DO isym = 1, nrot 324 classes_loop:DO iclass = 1, 12 325 elements_loop:DO ielem=1, nelem (iclass) 326 IF ( elem(ielem,iclass) == isym) THEN 327 symop_2_class(isym) = name_class(iclass) 328 EXIT classes_loop 329 END IF 330 END DO elements_loop 331 END DO classes_loop 332 END DO symmetries_loop 333 END IF 334 END IF 335 CALL qexsd_init_symmetries(output_obj%symmetries, nsym, nrot, spacegroup,& 336 s, ft, sname, t_rev, nat, irt,symop_2_class(1:nrot), verbosity, & 337 noncolin) 338 output_obj%symmetries_ispresent=.TRUE. 339 ! 340!------------------------------------------------------------------------------- 341! ... BASIS SET 342!------------------------------------------------------------------------------- 343 ! 344 CALL qexsd_init_basis_set(output_obj%basis_set, gamma_only, ecutwfc/e2, ecutwfc*dual/e2, & 345 dfftp%nr1, dfftp%nr2, dfftp%nr3, dffts%nr1, dffts%nr2, dffts%nr3, & 346 .FALSE., dfftp%nr1, dfftp%nr2, dfftp%nr3, ngm_g, ngms_g, npwx_g, & 347 bg(:,1), bg(:,2), bg(:,3) ) 348 ! 349!------------------------------------------------------------------------------- 350! ... DFT 351!------------------------------------------------------------------------------- 352 ! 353 IF (dft_is_hybrid() ) THEN 354 ALLOCATE ( hybrid_obj) 355 IF (get_screening_parameter() > 0.0_DP) THEN 356 scr_par_ = get_screening_parameter() 357 scr_par_opt=> scr_par_ 358 END IF 359 IF (ecutvcut > 0.0_DP) THEN 360 ectuvcut_ = ecutvcut/e2 361 ectuvcut_opt => ectuvcut_ 362 END IF 363 IF ( local_thr > 0._DP) THEN 364 loc_thr_ = local_thr 365 loc_thr_p => loc_thr_ 366 END IF 367 CALL qexsd_init_hybrid(hybrid_obj, DFT_IS_HYBRID = .TRUE., NQ1 = nq1 , NQ2 = nq2, NQ3 =nq3, ECUTFOCK = ecutfock/e2, & 368 EXX_FRACTION = get_exx_fraction(), SCREENING_PARAMETER = scr_par_opt, & 369 EXXDIV_TREATMENT = exxdiv_treatment, X_GAMMA_EXTRAPOLATION = x_gamma_extrapolation,& 370 ECUTVCUT = ectuvcut_opt, LOCAL_THR = loc_thr_p ) 371 END IF 372 373 empirical_vdw = (llondon .OR. ldftd3 .OR. lxdm .OR. ts_vdw ) 374 dft_is_vdw = dft_is_nonlocc() 375 IF ( dft_is_vdw .OR. empirical_vdw ) THEN 376 ALLOCATE (vdw_obj) 377 IF ( empirical_vdw) THEN 378 vdw_term_pt => dispersion_energy_term 379 vdw_corr_ = TRIM(vdw_corr) 380 vdw_corr_pt => vdw_corr_ 381 IF (llondon ) THEN 382 dispersion_energy_term = elondon/e2 383 lond_s6_ = scal6 384 london_s6_pt => lond_s6_ 385 lond_rcut_ = lon_rcut 386 london_rcut_pt => lond_rcut_ 387 IF (ANY( c6_i(1:nsp) .NE. -1._DP )) THEN 388 ALLOCATE (london_c6_(nsp), species_(nsp)) 389 london_c6_(1:nsp) = c6_i(1:nsp) 390 species_(1:nsp) = atm(1:nsp) 391 END IF 392 ! 393 ELSE IF ( lxdm ) THEN 394 dispersion_energy_term = exdm/e2 395 xdm_a1_ = xdm_a1 396 xdm_a1_pt => xdm_a1_ 397 xdm_a2_ = xdm_a2 398 xdm_a2_pt => xdm_a2_ 399 ! 400 ELSE IF ( ldftd3) THEN 401 dispersion_energy_term = edftd3/e2 402 dftd3_version_ = dftd3_version 403 dftd3_version_pt => dftd3_version_ 404 dftd3_threebody_ = dftd3_threebody 405 dftd3_threebody_pt => dftd3_threebody_ 406 ELSE IF ( ts_vdw ) THEN 407 dispersion_energy_term = 2._DP * EtsvdW/e2 408 ts_vdw_isolated_ = vdw_isolated 409 ts_vdw_isolated_pt => ts_vdw_isolated_ 410 ts_vdw_econv_thr_ = vdw_econv_thr 411 ts_vdw_econv_thr_pt => ts_vdw_econv_thr_ 412 END IF 413 ELSE 414 vdw_corr_ = 'none' 415 vdw_corr_pt => vdw_corr_ 416 END IF 417 IF (dft_is_vdw) THEN 418 dft_nonlocc_ = TRIM(get_nonlocc_name()) 419 non_local_term_pt => dft_nonlocc_ 420 END IF 421 CALL qexsd_init_vdw(vdw_obj, non_local_term_pt, vdw_corr_pt, vdw_term_pt, & 422 ts_thr_pt, ts_isol_pt, london_s6_pt, LONDON_C6 = london_c6_, & 423 LONDON_RCUT = london_rcut_pt, XDM_A1 = xdm_a1_pt, XDM_A2 = xdm_a2_pt,& 424 DFTD3_VERSION = dftd3_version_pt, DFTD3_THREEBODY = dftd3_threebody_pt) 425 END IF 426 IF ( lda_plus_u) THEN 427 ALLOCATE (dftU_obj) 428 CALL check_and_allocate_real(U_opt, Hubbard_U) 429 CALL check_and_allocate_real(J0_opt, Hubbard_J0) 430 CALL check_and_allocate_real(alpha_opt, Hubbard_alpha) 431 CALL check_and_allocate_real(beta_opt, Hubbard_beta) 432 CALL check_and_allocate_real(U_back_opt, Hubbard_U_back) 433 CALL check_and_allocate_real(alpha_back_opt, Hubbard_alpha_back) 434 CALL check_and_allocate_integer(Hubbard_l_back_opt, Hubbard_l_back) 435 CALL check_and_allocate_integer(Hubbard_l1_back_opt, Hubbard_l1_back) 436 CALL check_and_allocate_logical(backall_opt, backall) 437 IF ( ANY(Hubbard_J(:,1:nsp) /= 0.0_DP)) THEN 438 ALLOCATE (J_opt(3,nsp)) 439 J_opt(:, 1:nsp) = Hubbard_J(:, 1:nsp) 440 END IF 441 ! 442 ! Currently Hubbard_V, rho%nsb, and nsg are not written (read) to (from) XML 443 ! 444 CALL qexsd_init_dftU (dftU_obj, NSP = nsp, PSD = upf(1:nsp)%psd, SPECIES = atm(1:nsp), ITYP = ityp(1:nat), & 445 IS_HUBBARD = is_hubbard, IS_HUBBARD_BACK = is_hubbard_back, & 446 BACKALL = backall, HUBB_L_BACK = Hubbard_l_back_opt, HUBB_L1_BACK = Hubbard_l1_back_opt, & 447 NONCOLIN = noncolin, LDA_PLUS_U_KIND = lda_plus_u_kind, U_PROJECTION_TYPE = U_projection, & 448 U =U_opt, U_back = U_back_opt, J0 = J0_opt, J = J_opt, & 449 alpha = alpha_opt, beta = beta_opt, alpha_back = alpha_back_opt, & 450 starting_ns = starting_ns_eigenvalue, Hub_ns = rho%ns, Hub_ns_nc = rho%ns_nc) 451 END IF 452 dft_name = get_dft_short() 453 ! 454 CALL qexsd_init_dft (output_obj%dft, dft_name, hybrid_obj, vdw_obj, dftU_obj) 455 IF (ASSOCIATED (hybrid_obj)) THEN 456 CALL qes_reset(hybrid_obj) 457 DEALLOCATE (hybrid_obj) 458 END IF 459 IF (ASSOCIATED (vdw_obj)) THEN 460 CALL qes_reset(vdw_obj) 461 DEALLOCATE (vdw_obj) 462 END IF 463 IF (ASSOCIATED (dftU_obj)) THEN 464 CALL qes_reset( dftU_obj) 465 DEALLOCATE (dftU_obj) 466 END IF 467 ! 468!------------------------------------------------------------------------------- 469! ... PERIODIC BOUNDARY CONDITIONS 470!------------------------------------------------------------------------------- 471 ! 472 IF (ANY([do_makov_payne, do_comp_mt, do_comp_esm, do_cutoff_2D])) THEN 473 output_obj%boundary_conditions_ispresent=.TRUE. 474 IF (do_makov_payne) THEN 475 pbc_label = 'makov_payne' 476 ELSE IF ( do_comp_mt) THEN 477 pbc_label = 'martyna_tuckerman' 478 ELSE IF ( do_comp_esm) THEN 479 pbc_label = 'esm' 480 ELSE IF ( do_cutoff_2D) THEN 481 pbc_label = '2D' 482 ELSE 483 CALL errore ('pw_restart_new.f90: ', 'internal error line 470', 1) 484 END IF 485 CALL qexsd_init_outputPBC(output_obj%boundary_conditions, TRIM(pbc_label) ) 486 ENDIF 487 ! 488!------------------------------------------------------------------------------- 489! ... MAGNETIZATION 490!------------------------------------------------------------------------------- 491 ! 492 IF (noncolin) THEN 493 domag_ = domag 494 domag_opt=> domag_ 495 END IF 496 CALL qexsd_init_magnetization(output_obj%magnetization, lsda, noncolin, lspinorb, & 497 magtot, magtot_nc, absmag, domag_opt ) 498 ! 499 500!-------------------------------------------------------------------------------------- 501! ... BAND STRUCTURE 502!------------------------------------------------------------------------------------- 503 ! 504 ! skip if not yet computed 505 ! 506 IF ( only_init ) GO TO 10 507 ! 508 IF ( .NOT. ( lgauss .OR. ltetra )) THEN 509 CALL get_homo_lumo( h_energy, lumo_tmp) 510 h_energy = h_energy/e2 511 h_energy_ptr => h_energy 512 IF ( lumo_tmp .LT. 1.d+6 ) THEN 513 lumo_tmp = lumo_tmp/e2 514 lumo_energy => lumo_tmp 515 END IF 516 END IF 517 IF (nks_start == 0 .AND. nk1*nk2*nk3 > 0 ) THEN 518 CALL qexsd_init_k_points_ibz(qexsd_start_k_obj, "automatic", calculation, & 519 nk1, nk2, nk3, k1, k2, k3, nks_start, xk_start, wk_start, alat, at(:,1), .TRUE.) 520 ELSE 521 CALL qexsd_init_k_points_ibz(qexsd_start_k_obj, k_points, calculation, & 522 nk1, nk2, nk3, k1, k2, k3, nks_start, xk_start, wk_start, alat, at(:,1), .TRUE.) 523 END IF 524 qexsd_start_k_obj%tagname = 'starting_kpoints' 525 IF ( TRIM (qexsd_input_obj%tagname) == 'input') THEN 526 qexsd_occ_obj = qexsd_input_obj%bands%occupations 527 ELSE 528 CALL qexsd_init_occupations ( qexsd_occ_obj, input_parameters_occupations, nspin) 529 END IF 530 qexsd_occ_obj%tagname = 'occupations_kind' 531 IF ( two_fermi_energies ) THEN 532 ALLOCATE ( ef_updw (2) ) 533 IF (TRIM(input_parameters_occupations) == 'fixed') THEN 534 ef_updw(1) = MAXVAL(et(INT(nelup),1:nkstot/2))/e2 535 ef_updw(2) = MAXVAL(et(INT(neldw),nkstot/2+1:nkstot))/e2 536 ELSE 537 ef_updw = [ef_up/e2, ef_dw/e2] 538 END IF 539 ELSE 540 ! The Fermi energy is written also for insulators because it can 541 ! be useful for further postprocessing, especially of bands 542 ! (for an insulator the Fermi energy is equal to the HOMO/VBMax) 543 ef_targ = ef/e2 544 ef_point => ef_targ 545 END IF 546 547 548 IF ( lgauss ) THEN 549 IF (TRIM(qexsd_input_obj%tagname) == 'input') THEN 550 smear_obj = qexsd_input_obj%bands%smearing 551 ELSE 552 smearing_loc = schema_smearing( smearing ) 553 CALL qexsd_init_smearing(smear_obj, smearing_loc, degauss/e2) 554 END IF 555 smear_obj_ptr => smear_obj 556 END IF 557 ! 558 559 CALL qexsd_init_band_structure( output_obj%band_structure,lsda,noncolin,lspinorb, nelec, natomwfc, & 560 et, wg, nkstot, xk, ngk_g, wk, SMEARING = smear_obj_ptr, & 561 STARTING_KPOINTS = qexsd_start_k_obj, OCCUPATIONS_KIND = qexsd_occ_obj, & 562 WF_COLLECTED = wf_collect, NBND = nbnd, FERMI_ENERGY = ef_point, EF_UPDW = ef_updw,& 563 HOMO = h_energy_ptr, LUMO = lumo_energy ) 564 ! 565 IF (lgauss) CALL qes_reset (smear_obj) 566 CALL qes_reset (qexsd_start_k_obj) 567 CALL qes_reset (qexsd_occ_obj) 568 ! 569!------------------------------------------------------------------------------------------- 570! ... TOTAL ENERGY 571!------------------------------------------------------------------------------------------- 572 ! 573 IF ( degauss > 0.0d0 ) THEN 574 ! 575 itemp = itemp + 1 576 temp(itemp) = degauss/e2 577 degauss_ => temp(itemp) 578 ! 579 itemp = itemp+1 580 temp(itemp) = demet/e2 581 demet_ => temp(itemp) 582 END IF 583 IF ( tefield ) THEN 584 itemp = itemp+1 585 temp(itemp) = etotefield/e2 586 efield_corr => temp(itemp) 587 END IF 588 IF (lfcpopt .OR. lfcpdyn ) THEN 589 itemp = itemp +1 590 temp(itemp) = ef * tot_charge/e2 591 potstat_corr => temp(itemp) 592 output_obj%FCP_tot_charge_ispresent = .TRUE. 593 output_obj%FCP_tot_charge = tot_charge 594 output_obj%FCP_force_ispresent = .TRUE. 595 !FIXME ( decide what units to use here ) 596 output_obj%FCP_force = fcp_mu - ef 597 END IF 598 IF ( gate) THEN 599 itemp = itemp + 1 600 temp(itemp) = etotgatefield/e2 601 gatefield_corr => temp(itemp) 602 END IF 603 604 CALL qexsd_init_total_energy(output_obj%total_energy, etot/e2, eband/e2, ehart/e2, vtxc/e2, & 605 etxc/e2, ewld/e2, degauss_, demet_, efield_corr, potstat_corr,& 606 gatefield_corr, DISPERSION_CONTRIBUTION = vdw_term_pt) 607 ! 608 NULLIFY(degauss_, demet_, efield_corr, potstat_corr, gatefield_corr) 609 itemp = 0 610 ! 611!--------------------------------------------------------------------------------------------- 612! ... FORCES 613!---------------------------------------------------------------------------------------------- 614 ! 615 IF ( lforce .and. conv_elec ) THEN 616 output_obj%forces_ispresent = .TRUE. 617 CALL qexsd_init_forces(output_obj%forces,nat,force,lforce) 618 ELSE 619 output_obj%forces_ispresent = .FALSE. 620 output_obj%forces%lwrite = .FALSE. 621 END IF 622 ! 623!------------------------------------------------------------------------------------------------ 624! ... STRESS 625!------------------------------------------------------------------------------------------------ 626 IF ( lstres .and. conv_elec ) THEN 627 output_obj%stress_ispresent=.TRUE. 628 CALL qexsd_init_stress(output_obj%stress, sigma, lstres ) 629 ELSE 630 output_obj%stress_ispresent=.FALSE. 631 output_obj%stress%lwrite=.FALSE. 632 END IF 633!------------------------------------------------------------------------------------------------- 634! ... ELECTRIC FIELD 635!------------------------------------------------------------------------------------------------- 636 output_obj%electric_field_ispresent = ( gate .OR. lelfield .OR. lberry .OR. tefield ) 637 638 IF ( gate ) THEN 639 CALL qexsd_init_gate_info(gate_info_temp,"gateInfo", etotgatefield/e2, zgate, nelec, & 640 alat, at, bg, zv, ityp) 641 gate_info_ptr => gate_info_temp 642 END IF 643 IF ( lelfield ) THEN 644 ALLOCATE (bp_el_pol(2), bp_ion_pol(3) ) 645 bp_el_pol = el_pol 646 bp_ion_pol(1:3) = ion_pol(1:3) 647 END IF 648 IF ( tefield .AND. dipfield) THEN 649 CALL qexsd_init_dipole_info(dipol_obj, el_dipole, ion_dipole, edir, eamp, & 650 emaxpos, eopreg ) 651 dipol_ptr => dipol_obj 652 END IF 653 IF ( lberry ) bp_obj_ptr => qexsd_bp_obj 654 IF (output_obj%electric_field_ispresent) & 655 CALL qexsd_init_outputElectricField(output_obj%electric_field, lelfield, tefield, dipfield, & 656 lberry, BP_OBJ = bp_obj_ptr, EL_POL = bp_el_pol, ION_POL = bp_ion_pol, & 657 GATEINFO = gate_info_ptr, DIPOLE_OBJ = dipol_ptr) 658 ! 659 IF (ASSOCIATED(gate_info_ptr)) THEN 660 CALL qes_reset (gate_info_ptr) 661 NULLIFY(gate_info_ptr) 662 ENDIF 663 IF (ASSOCIATED (dipol_ptr) ) THEN 664 CALL qes_reset (dipol_ptr) 665 NULLIFY(dipol_ptr) 666 ENDIF 667 NULLIFY ( bp_obj_ptr) 668!------------------------------------------------------------------------------- 669! ... ACTUAL WRITING 670!------------------------------------------------------------------------------- 671 10 CONTINUE 672 ! 673 CALL qexsd_openschema( xmlfile(), iunpun, 'PWSCF', title ) 674 CALL qes_write (qexsd_xf,output_obj) 675 CALL qes_reset (output_obj) 676 CALL qexsd_closeschema() 677 ! 678!------------------------------------------------------------------------------- 679 ! 680 END IF 681 DEALLOCATE (ngk_g) 682 ! 683 RETURN 684 ! 685 CONTAINS 686 SUBROUTINE check_and_allocate_real(alloc, mydata) 687 IMPLICIT NONE 688 REAL(DP),ALLOCATABLE :: alloc(:) 689 REAL(DP) :: mydata(:) 690 IF ( ANY(mydata(1:nsp) /= 0.0_DP)) THEN 691 ALLOCATE(alloc(nsp)) 692 alloc(1:nsp) = mydata(1:nsp) 693 END IF 694 RETURN 695 END SUBROUTINE check_and_allocate_real 696 ! 697 SUBROUTINE check_and_allocate_integer(alloc, mydata) 698 IMPLICIT NONE 699 INTEGER,ALLOCATABLE :: alloc(:) 700 INTEGER :: mydata(:) 701 IF ( ANY(mydata(1:nsp) /= -1)) THEN 702 ALLOCATE(alloc(nsp)) 703 alloc(1:nsp) = mydata(1:nsp) 704 END IF 705 RETURN 706 END SUBROUTINE check_and_allocate_integer 707 ! 708 SUBROUTINE check_and_allocate_logical(alloc, mydata) 709 IMPLICIT NONE 710 LOGICAL,ALLOCATABLE :: alloc(:) 711 LOGICAL :: mydata(:) 712 IF ( ANY(mydata(1:nsp))) THEN 713 ALLOCATE(alloc(nsp)) 714 alloc(1:nsp) = mydata(1:nsp) 715 END IF 716 RETURN 717 END SUBROUTINE check_and_allocate_logical 718 ! 719 END SUBROUTINE pw_write_schema 720 ! 721 !------------------------------------------------------------------------ 722 SUBROUTINE pw_write_binaries( ) 723 !------------------------------------------------------------------------ 724 ! 725 USE mp, ONLY : mp_sum, mp_max 726 USE io_base, ONLY : write_wfc 727 USE io_files, ONLY : restart_dir, iunwfc, nwordwfc 728 USE cell_base, ONLY : tpiba, alat, bg 729 USE control_flags, ONLY : gamma_only, smallmem 730 USE gvect, ONLY : ig_l2g 731 USE noncollin_module, ONLY : noncolin, npol 732 USE buffers, ONLY : get_buffer 733 USE wavefunctions, ONLY : evc 734 USE klist, ONLY : nks, nkstot, xk, ngk, igk_k 735 USE gvect, ONLY : ngm, g, mill 736 USE fft_base, ONLY : dfftp 737 USE basis, ONLY : natomwfc 738 USE wvfct, ONLY : npwx, et, wg, nbnd 739 USE lsda_mod, ONLY : nspin, isk, lsda 740 USE mp_pools, ONLY : intra_pool_comm, inter_pool_comm 741 USE mp_bands, ONLY : me_bgrp, root_bgrp, intra_bgrp_comm, & 742 root_bgrp_id, my_bgrp_id 743 USE wrappers, ONLY : f_mkdir_safe 744 ! 745 IMPLICIT NONE 746 ! 747 INTEGER :: ios, ig, ngg, ipol, ispin 748 INTEGER :: ik, ik_g, ike, iks, npw_g 749 INTEGER, EXTERNAL :: global_kpoint_index 750 INTEGER, ALLOCATABLE :: ngk_g(:), mill_k(:,:) 751 INTEGER, ALLOCATABLE :: igk_l2g(:), igk_l2g_kdip(:) 752 CHARACTER(LEN=2), DIMENSION(2) :: updw = (/ 'up', 'dw' /) 753 CHARACTER(LEN=256) :: dirname 754 CHARACTER(LEN=320) :: filename 755 ! 756 dirname = restart_dir () 757 ! 758 ! ... check that restart_dir exists on all processors that write 759 ! ... wavefunctions; create one if restart_dir is not found. This 760 ! ... is needed for k-point parallelization, in the case of non-parallel 761 ! ... scratch file systems, that are not visible to all processors 762 ! 763 IF ( my_bgrp_id == root_bgrp_id .AND. me_bgrp == root_bgrp ) THEN 764 ios = f_mkdir_safe( TRIM(dirname) ) 765 END IF 766 ! 767 ! ... write wavefunctions and k+G vectors 768 ! 769 iks = global_kpoint_index (nkstot, 1) 770 ike = iks + nks - 1 771 ! 772 ! ... ngk_g: global number of k+G vectors 773 ! 774 ALLOCATE( ngk_g( nks ) ) 775 ngk_g(1:nks) = ngk(1:nks) 776 CALL mp_sum( ngk_g, intra_bgrp_comm) 777 ! 778 ! ... The igk_l2g array yields the correspondence between the 779 ! ... local k+G index and the global G index 780 ! 781 ALLOCATE ( igk_l2g( npwx ) ) 782 ! 783 ! ... the igk_l2g_kdip local-to-global map yields the correspondence 784 ! ... between the global order of k+G and the local index for k+G. 785 ! 786 ALLOCATE ( igk_l2g_kdip( npwx ) ) 787 ! 788 ALLOCATE ( mill_k( 3, npwx ) ) 789 ! 790 k_points_loop: DO ik = 1, nks 791 ! 792 ! ik_g is the index of k-point ik in the global list 793 ! 794 ik_g = ik + iks - 1 795 ! 796 ! ... Compute the igk_l2g array from previously computed arrays 797 ! ... igk_k (k+G indices) and ig_l2g (local to global G index map) 798 ! 799 igk_l2g = 0 800 DO ig = 1, ngk (ik) 801 igk_l2g(ig) = ig_l2g(igk_k(ig,ik)) 802 END DO 803 ! 804 ! ... npw_g is the maximum G vector index among all processors 805 ! 806 npw_g = MAXVAL( igk_l2g(1:ngk(ik)) ) 807 CALL mp_max( npw_g, intra_pool_comm ) 808 ! 809 igk_l2g_kdip = 0 810 CALL gk_l2gmap_kdip( npw_g, ngk_g(ik), ngk(ik), igk_l2g, & 811 igk_l2g_kdip ) 812 ! 813 ! ... mill_k(:,i) contains Miller indices for (k+G)_i 814 ! 815 DO ig = 1, ngk (ik) 816 mill_k(:,ig) = mill(:,igk_k(ig,ik)) 817 END DO 818 ! 819 ! ... read wavefunctions - do not read if already in memory (nsk==1) 820 ! 821 IF ( nks > 1 ) CALL get_buffer ( evc, nwordwfc, iunwfc, ik ) 822 ! 823 IF ( nspin == 2 ) THEN 824 ! 825 ! ... LSDA: spin mapped to k-points, isk(ik) tracks up and down spin 826 ! 827 ik_g = MOD ( ik_g-1, nkstot/2 ) + 1 828 ispin = isk(ik) 829 filename = TRIM(dirname) // 'wfc' // updw(ispin) // & 830 & TRIM(int_to_char(ik_g)) 831 ! 832 ELSE 833 ! 834 ispin = 1 835 filename = TRIM(dirname) // 'wfc' // TRIM(int_to_char(ik_g)) 836 ! 837 ENDIF 838 ! 839 ! ... Only the first band group of each pool writes 840 ! ... No warranty it works for more than one band group 841 ! 842 IF ( my_bgrp_id == root_bgrp_id ) CALL write_wfc( iunpun, & 843 filename, root_bgrp, intra_bgrp_comm, ik_g, tpiba*xk(:,ik), & 844 ispin, nspin, evc, npw_g, gamma_only, nbnd, & 845 igk_l2g_kdip(:), ngk(ik), tpiba*bg(:,1), tpiba*bg(:,2), & 846 tpiba*bg(:,3), mill_k, 1.D0 ) 847 ! 848 END DO k_points_loop 849 ! 850 DEALLOCATE ( mill_k ) 851 DEALLOCATE ( igk_l2g_kdip ) 852 DEALLOCATE ( igk_l2g ) 853 DEALLOCATE ( ngk_g ) 854 ! 855 RETURN 856 ! 857 END SUBROUTINE pw_write_binaries 858 ! 859 !----------------------------------------------------------------------- 860 SUBROUTINE gk_l2gmap_kdip( npw_g, ngk_g, ngk, igk_l2g, igk_l2g_kdip, igwk ) 861 !----------------------------------------------------------------------- 862 ! 863 ! ... This subroutine maps local G+k index to the global G vector index 864 ! ... the mapping is used to collect wavefunctions subsets distributed 865 ! ... across processors. 866 ! ... This map is used to obtained the G+k grids related to each kpt 867 ! 868 USE mp_bands, ONLY : intra_bgrp_comm 869 USE mp, ONLY : mp_sum 870 ! 871 IMPLICIT NONE 872 ! 873 ! ... Here the dummy variables 874 ! 875 INTEGER, INTENT(IN) :: npw_g, ngk_g, ngk 876 INTEGER, INTENT(IN) :: igk_l2g(ngk) 877 INTEGER, INTENT(OUT) :: igk_l2g_kdip(ngk) 878 INTEGER, OPTIONAL, INTENT(OUT) :: igwk(ngk_g) 879 ! 880 INTEGER, ALLOCATABLE :: igwk_(:), itmp(:), igwk_lup(:) 881 INTEGER :: ig, ig_, ngg 882 ! 883 ! 884 ALLOCATE( itmp( npw_g ) ) 885 ALLOCATE( igwk_( ngk_g ) ) 886 ! 887 itmp(:) = 0 888 igwk_(:) = 0 889 ! 890 DO ig = 1, ngk 891 itmp(igk_l2g(ig)) = igk_l2g(ig) 892 END DO 893 ! 894 CALL mp_sum( itmp, intra_bgrp_comm ) 895 ! 896 ngg = 0 897 DO ig = 1, npw_g 898 ! 899 IF ( itmp(ig) == ig ) THEN 900 ! 901 ngg = ngg + 1 902 igwk_(ngg) = ig 903 ! 904 END IF 905 ! 906 END DO 907 ! 908 IF ( ngg /= ngk_g ) & 909 CALL errore( 'gk_l2gmap_kdip', 'unexpected dimension in ngg', 1 ) 910 ! 911 IF ( PRESENT( igwk ) ) THEN 912 ! 913 igwk(1:ngk_g) = igwk_(1:ngk_g) 914 ! 915 END IF 916 ! 917 ALLOCATE( igwk_lup( npw_g ) ) 918 ! 919!$omp parallel private(ig_, ig) 920!$omp workshare 921 igwk_lup = 0 922!$omp end workshare 923!$omp do 924 DO ig_ = 1, ngk_g 925 igwk_lup(igwk_(ig_)) = ig_ 926 END DO 927!$omp end do 928!$omp do 929 DO ig = 1, ngk 930 igk_l2g_kdip(ig) = igwk_lup(igk_l2g(ig)) 931 END DO 932!$omp end do 933!$omp end parallel 934 ! 935 DEALLOCATE( igwk_lup ) 936 ! 937 DEALLOCATE( itmp, igwk_ ) 938 ! 939 RETURN 940 ! 941 END SUBROUTINE gk_l2gmap_kdip 942 ! 943 !-------------------------------------------------------------------------- 944 SUBROUTINE read_xml_file ( wfc_is_collected ) 945 !------------------------------------------------------------------------ 946 ! 947 ! ... This routine allocates space for all quantities already computed 948 ! ... in the pwscf program and reads them from the data file. 949 ! ... All quantities that are initialized in subroutine "setup" when 950 ! ... starting from scratch should be initialized here when restarting 951 ! 952 USE kinds, ONLY : dp 953 USE constants, ONLY : e2 954 USE gvect, ONLY : ngm_g, ecutrho 955 USE gvecs, ONLY : ngms_g, dual 956 USE gvecw, ONLY : ecutwfc 957 USE fft_base, ONLY : dfftp, dffts 958 USE io_global, ONLY : stdout 959 USE io_files, ONLY : psfile, pseudo_dir, pseudo_dir_cur, & 960 restart_dir 961 USE mp_global, ONLY : nproc_file, nproc_pool_file, & 962 nproc_image_file, ntask_groups_file, & 963 nproc_bgrp_file, nproc_ortho_file 964 USE ions_base, ONLY : nat, nsp, ityp, amass, atm, tau, extfor 965 USE cell_base, ONLY : alat, at, bg, ibrav, celldm, omega 966 USE force_mod, ONLY : force 967 USE klist, ONLY : nks, nkstot, xk, wk, tot_magnetization, & 968 nelec, nelup, neldw, smearing, degauss, ngauss, lgauss, ltetra,& 969 two_fermi_energies 970 USE ktetra, ONLY : ntetra, tetra_type 971 USE start_k, ONLY : nks_start, xk_start, wk_start, & 972 nk1, nk2, nk3, k1, k2, k3 973 USE ener, ONLY : ef, ef_up, ef_dw 974 USE electrons_base, ONLY : nupdwn, set_nelup_neldw 975 USE wvfct, ONLY : npwx, nbnd, et, wg 976 USE extfield, ONLY : forcefield, forcegate, tefield, dipfield, & 977 edir, emaxpos, eopreg, eamp, el_dipole, ion_dipole, gate, zgate, & 978 relaxz, block, block_1, block_2, block_height 979 USE symm_base, ONLY : nrot, nsym, invsym, s, ft, irt, t_rev, & 980 sname, inverse_s, s_axis_to_cart, & 981 time_reversal, no_t_rev, nosym, checkallsym 982 USE ldaU, ONLY : lda_plus_u, lda_plus_u_kind, Hubbard_lmax, Hubbard_lmax_back, & 983 Hubbard_l, Hubbard_l_back, Hubbard_l1_back, backall, & 984 Hubbard_U, Hubbard_U_back, Hubbard_J, Hubbard_V, Hubbard_alpha, & 985 Hubbard_alpha_back, Hubbard_J0, Hubbard_beta, U_projection 986 USE funct, ONLY : set_exx_fraction, set_screening_parameter, & 987 set_gau_parameter, enforce_input_dft, & 988 start_exx, dft_is_hybrid 989 USE london_module, ONLY : scal6, lon_rcut, in_C6 990 USE tsvdw_module, ONLY : vdw_isolated 991 USE exx_base, ONLY : x_gamma_extrapolation, nq1, nq2, nq3, & 992 exxdiv_treatment, yukawa, ecutvcut 993 USE exx, ONLY : ecutfock, local_thr 994 USE control_flags, ONLY : noinv, gamma_only, tqr, llondon, ldftd3, & 995 lxdm, ts_vdw 996 USE Coul_cut_2D, ONLY : do_cutoff_2D 997 USE noncollin_module,ONLY : noncolin, npol, angle1, angle2, bfield, & 998 nspin_lsda, nspin_gga, nspin_mag 999 USE spin_orb, ONLY : domag, lspinorb 1000 USE lsda_mod, ONLY : nspin, isk, lsda, starting_magnetization,& 1001 current_spin 1002 USE realus, ONLY : real_space 1003 USE basis, ONLY : natomwfc 1004 USE uspp, ONLY : okvan 1005 USE paw_variables, ONLY : okpaw 1006 ! 1007 USE mp_images, ONLY : intra_image_comm 1008 USE mp, ONLY : mp_bcast 1009 ! 1010 IMPLICIT NONE 1011 LOGICAL, INTENT(OUT) :: wfc_is_collected 1012 ! 1013 INTEGER :: i, is, ik, ierr, dum1,dum2,dum3 1014 LOGICAL :: magnetic_sym, lvalid_input, lfixed 1015 CHARACTER(LEN=26) :: dft_name 1016 CHARACTER(LEN=20) :: vdw_corr, occupations 1017 CHARACTER(LEN=320):: filename 1018 REAL(dp) :: exx_fraction, screening_parameter 1019 TYPE (output_type) :: output_obj 1020 TYPE (parallel_info_type) :: parinfo_obj 1021 TYPE (general_info_type ) :: geninfo_obj 1022 TYPE (input_type) :: input_obj 1023 ! 1024 ! 1025 filename = xmlfile ( ) 1026 ! 1027 IF (ionode) CALL qexsd_readschema ( filename, & 1028 ierr, output_obj, parinfo_obj, geninfo_obj, input_obj) 1029 CALL mp_bcast(ierr, ionode_id, intra_image_comm) 1030 IF ( ierr > 0 ) THEN 1031 CALL errore ( 'read_xml_file', 'fatal error reading xml file', ierr ) 1032 ELSE IF ( ierr < 0 ) THEN 1033 input_obj%tagname = "not_read" 1034 ! ierr = -1 means that input_obj was not read: do not broadcast it 1035 ELSE 1036 CALL qes_bcast(input_obj, ionode_id, intra_image_comm) 1037 END IF 1038 CALL qes_bcast(output_obj, ionode_id, intra_image_comm) 1039 CALL qes_bcast(parinfo_obj, ionode_id, intra_image_comm) 1040 CALL qes_bcast(geninfo_obj, ionode_id, intra_image_comm) 1041 ! 1042 ! ... Now read all needed variables from xml objects 1043 ! 1044 wfc_is_collected = output_obj%band_structure%wf_collected 1045 lvalid_input = (TRIM(input_obj%tagname) == "input") 1046 ! 1047 CALL qexsd_copy_parallel_info (parinfo_obj, nproc_file, & 1048 nproc_pool_file, nproc_image_file, ntask_groups_file, & 1049 nproc_bgrp_file, nproc_ortho_file) 1050 ! 1051 pseudo_dir_cur = restart_dir ( ) 1052 CALL qexsd_copy_atomic_species ( output_obj%atomic_species, & 1053 nsp, atm, amass, starting_magnetization, angle1, angle2, & 1054 psfile, pseudo_dir ) 1055 IF ( pseudo_dir == ' ' ) pseudo_dir=pseudo_dir_cur 1056 !! Atomic structure section 1057 !! tau and ityp are allocated inside qexsd_copy_atomic_structure 1058 ! 1059 CALL qexsd_copy_atomic_structure (output_obj%atomic_structure, nsp, & 1060 atm, nat, tau, ityp, alat, at(:,1), at(:,2), at(:,3), ibrav ) 1061 ! 1062 !! More initializations needed for atomic structure: 1063 !! bring atomic positions and crystal axis into "alat" units; 1064 !! recalculate celldm; compute cell volume, reciprocal lattice vectors 1065 ! 1066 at = at / alat 1067 tau(:,1:nat) = tau(:,1:nat)/alat 1068 CALL at2celldm (ibrav,alat,at(:,1),at(:,2),at(:,3),celldm) 1069 CALL volume (alat,at(:,1),at(:,2),at(:,3),omega) 1070 !! 1071 !! Basis set section 1072 CALL qexsd_copy_basis_set ( output_obj%basis_set, gamma_only, ecutwfc,& 1073 ecutrho, dffts%nr1,dffts%nr2,dffts%nr3, dfftp%nr1,dfftp%nr2,dfftp%nr3, & 1074 dum1,dum2,dum3, ngm_g, ngms_g, npwx, bg(:,1), bg(:,2), bg(:,3) ) 1075 ecutwfc = ecutwfc*e2 1076 ecutrho = ecutrho*e2 1077 dual = ecutrho/ecutwfc 1078 ! FIXME: next line ensures exact consistency between reciprocal and 1079 ! direct lattice vectors, preventing weird phonon symmetry errors 1080 ! (due to lousy algorithms, extraordinarily sensitive to tiny errors) 1081 CALL recips ( at(1,1), at(1,2), at(1,3), bg(1,1), bg(1,2), bg(1,3) ) 1082 !! 1083 !! DFT section 1084 CALL qexsd_copy_dft ( output_obj%dft, nsp, atm, & 1085 dft_name, nq1, nq2, nq3, ecutfock, exx_fraction, screening_parameter, & 1086 exxdiv_treatment, x_gamma_extrapolation, ecutvcut, local_thr, & 1087 lda_plus_U, lda_plus_U_kind, U_projection, Hubbard_l, Hubbard_lmax, & 1088 Hubbard_l_back, Hubbard_l1_back, backall, Hubbard_lmax_back, Hubbard_alpha_back, & 1089 Hubbard_U, Hubbard_U_back, Hubbard_J0, Hubbard_alpha, Hubbard_beta, Hubbard_J, & 1090 vdw_corr, scal6, lon_rcut, vdw_isolated ) 1091 !! More DFT initializations 1092 CALL set_vdw_corr ( vdw_corr, llondon, ldftd3, ts_vdw, lxdm ) 1093 CALL enforce_input_dft ( dft_name, .TRUE. ) 1094 IF ( dft_is_hybrid() ) THEN 1095 ecutvcut=ecutvcut*e2 1096 ecutfock=ecutfock*e2 1097 CALL set_exx_fraction( exx_fraction ) 1098 CALL set_screening_parameter ( screening_parameter ) 1099 CALL start_exx () 1100 END IF 1101 !! Band structure section 1102 !! et and wg are allocated inside qexsd_copy_band_structure 1103 CALL qexsd_copy_band_structure( output_obj%band_structure, lsda, & 1104 nkstot, isk, natomwfc, nbnd, nupdwn(1), nupdwn(2), nelec, xk, & 1105 wk, wg, ef, ef_up, ef_dw, et ) 1106 ! convert to Ry 1107 ef = ef*e2 1108 ef_up = ef_up*e2 1109 ef_dw = ef_dw*e2 1110 two_fermi_energies = ( ef_up /= 0.0_dp ) .AND. ( ef_dw /= 0.0_dp ) 1111 et(:,:) = et(:,:)*e2 1112 ! 1113 ! ... until pools are activated, the local number of k-points nks 1114 ! ... should be equal to the global number nkstot - k-points are replicated 1115 ! 1116 nks = nkstot 1117 !! 1118 !! Magnetization section 1119 CALL qexsd_copy_magnetization ( output_obj%magnetization, lsda, noncolin,& 1120 lspinorb, domag, tot_magnetization ) 1121 ! 1122 bfield = 0.d0 1123 CALL set_spin_vars( lsda, noncolin, lspinorb, domag, & 1124 npol, nspin, nspin_lsda, nspin_mag, nspin_gga, current_spin ) 1125 !! Information for generating k-points and occupations 1126 CALL qexsd_copy_kpoints( output_obj%band_structure, & 1127 nks_start, xk_start, wk_start, nk1, nk2, nk3, k1, k2, k3, & 1128 occupations, smearing, degauss ) 1129 degauss = degauss * e2 1130 ! 1131 CALL set_occupations( occupations, smearing, degauss, & 1132 lfixed, ltetra, tetra_type, lgauss, ngauss ) 1133 IF (ltetra) ntetra = 6* nk1 * nk2 * nk3 1134 IF (lfixed) CALL errore('read_file','bad occupancies',1) 1135 IF ( lsda ) & 1136 CALL set_nelup_neldw(tot_magnetization, nelec, nelup, neldw) 1137 !! Symmetry section 1138 ALLOCATE ( irt(48,nat) ) 1139 IF ( lvalid_input ) THEN 1140 CALL qexsd_copy_symmetry ( output_obj%symmetries, & 1141 nsym, nrot, s, ft, sname, t_rev, invsym, irt, & 1142 noinv, nosym, no_t_rev, input_obj%symmetry_flags ) 1143 1144 CALL qexsd_copy_efield ( input_obj%electric_field, & 1145 tefield, dipfield, edir, emaxpos, eopreg, eamp, & 1146 gate, zgate, block, block_1, block_2, block_height, relaxz ) 1147 1148 ELSE 1149 CALL qexsd_copy_symmetry ( output_obj%symmetries, & 1150 nsym, nrot, s, ft, sname, t_rev, invsym, irt, & 1151 noinv, nosym, no_t_rev ) 1152 ENDIF 1153 !! More initialization needed for symmetry 1154 magnetic_sym = noncolin .AND. domag 1155 time_reversal = (.NOT.magnetic_sym) .AND. (.NOT.noinv) 1156 CALL inverse_s() 1157 CALL s_axis_to_cart() 1158 !! symmetry check - FIXME: must be done in a more consistent way 1159 !! IF (nat > 0) CALL checkallsym( nat, tau, ityp) 1160 !! Algorithmic info 1161 do_cutoff_2D = (output_obj%boundary_conditions%assume_isolated == "2D") 1162 CALL qexsd_copy_algorithmic_info ( output_obj%algorithmic_info, & 1163 real_space, tqr, okvan, okpaw ) 1164 ! 1165 ! ... xml data no longer needed, can be discarded 1166 ! 1167 CALL qes_reset ( output_obj ) 1168 CALL qes_reset ( geninfo_obj ) 1169 CALL qes_reset ( parinfo_obj ) 1170 IF ( TRIM(input_obj%tagname) == "input") CALL qes_reset ( input_obj) 1171 ! 1172 ! END OF READING VARIABLES FROM XML DATA FILE 1173 ! 1174 ALLOCATE( force ( 3, nat ) ) 1175 ALLOCATE( extfor( 3, nat ) ) 1176 IF ( tefield ) ALLOCATE( forcefield( 3, nat ) ) 1177 IF ( gate ) ALLOCATE( forcegate( 3, nat ) ) 1178 ! 1179 END SUBROUTINE read_xml_file 1180 ! 1181 !------------------------------------------------------------------------ 1182 SUBROUTINE read_collected_wfc ( dirname, ik, evc ) 1183 !------------------------------------------------------------------------ 1184 ! 1185 ! ... reads from directory "dirname" (new file format) for k-point "ik" 1186 ! ... wavefunctions from collected format into distributed array "evc" 1187 ! 1188 USE control_flags, ONLY : gamma_only 1189 USE lsda_mod, ONLY : nspin, isk 1190 USE noncollin_module, ONLY : noncolin, npol 1191 USE klist, ONLY : nkstot, nks, xk, ngk, igk_k 1192 USE wvfct, ONLY : npwx, g2kin, et, wg, nbnd 1193 USE gvect, ONLY : ig_l2g 1194 USE mp_bands, ONLY : root_bgrp, intra_bgrp_comm 1195 USE mp_pools, ONLY : me_pool, root_pool, & 1196 intra_pool_comm, inter_pool_comm 1197 USE mp, ONLY : mp_sum, mp_max 1198 USE io_base, ONLY : read_wfc 1199 ! 1200 IMPLICIT NONE 1201 ! 1202 CHARACTER(LEN=*), INTENT(IN) :: dirname 1203 INTEGER, INTENT(IN) :: ik 1204 COMPLEX(dp), INTENT(OUT) :: evc(:,:) 1205 ! 1206 CHARACTER(LEN=2), DIMENSION(2) :: updw = (/ 'up', 'dw' /) 1207 CHARACTER(LEN=320) :: filename, msg 1208 INTEGER :: i, ik_g, ig, ipol, ik_s 1209 INTEGER :: npol_, nbnd_ 1210 INTEGER :: nupdwn(2), ike, iks, ngk_g, npw_g, ispin 1211 INTEGER, EXTERNAL :: global_kpoint_index 1212 INTEGER, ALLOCATABLE :: mill_k(:,:) 1213 INTEGER, ALLOCATABLE :: igk_l2g(:), igk_l2g_kdip(:) 1214 LOGICAL :: opnd, ionode_k 1215 REAL(DP) :: scalef, xk_(3), b1(3), b2(3), b3(3) 1216 ! 1217 ! ... the root processor of each pool reads 1218 ! 1219 ionode_k = (me_pool == root_pool) 1220 ! 1221 iks = global_kpoint_index (nkstot, 1) 1222 ike = iks + nks - 1 1223 ! 1224 ! ik_g: index of k-point ik in the global list 1225 ! 1226 ik_g = ik + iks - 1 1227 ! 1228 ! ... the igk_l2g_kdip local-to-global map is needed to read wfcs 1229 ! 1230 ALLOCATE ( igk_l2g_kdip( npwx ) ) 1231 ! 1232 ! ... The igk_l2g array yields the correspondence between the 1233 ! ... local k+G index and the global G index - requires arrays 1234 ! ... igk_k (k+G indices) and ig_l2g (local to global G index map) 1235 ! 1236 ALLOCATE ( igk_l2g( npwx ) ) 1237 igk_l2g = 0 1238 DO ig = 1, ngk(ik) 1239 igk_l2g(ig) = ig_l2g(igk_k(ig,ik)) 1240 END DO 1241 ! 1242 ! ... npw_g: the maximum G vector index among all processors 1243 ! ... ngk_g: global number of k+G vectors for all k points 1244 ! 1245 npw_g = MAXVAL( igk_l2g(1:ngk(ik)) ) 1246 CALL mp_max( npw_g, intra_pool_comm ) 1247 ngk_g = ngk(ik) 1248 CALL mp_sum( ngk_g, intra_bgrp_comm) 1249 ! 1250 ! ... now compute the igk_l2g_kdip local-to-global map 1251 ! 1252 igk_l2g_kdip = 0 1253 CALL gk_l2gmap_kdip( npw_g, ngk_g, ngk(ik), igk_l2g, & 1254 igk_l2g_kdip ) 1255 DEALLOCATE ( igk_l2g ) 1256 ! 1257 IF ( nspin == 2 ) THEN 1258 ! 1259 ! ... LSDA: spin mapped to k-points, isk(ik) tracks up and down spin 1260 ! 1261 ik_g = MOD ( ik_g-1, nkstot/2 ) + 1 1262 ispin = isk(ik) 1263 filename = TRIM(dirname) // 'wfc' // updw(ispin) // & 1264 & TRIM(int_to_char(ik_g)) 1265 ! 1266 ELSE 1267 ! 1268 filename = TRIM(dirname) // 'wfc' // TRIM(int_to_char(ik_g)) 1269 ! 1270 ENDIF 1271 ! 1272 ! ... Miller indices are read from file (but not used) 1273 ! 1274 ALLOCATE( mill_k ( 3,npwx ) ) 1275 ! 1276 evc=(0.0_DP, 0.0_DP) 1277 ! 1278 CALL read_wfc( iunpun, filename, root_bgrp, intra_bgrp_comm, & 1279 ik_g, xk_, ispin, npol_, evc, npw_g, gamma_only, nbnd_, & 1280 igk_l2g_kdip(:), ngk(ik), b1, b2, b3, mill_k, scalef ) 1281 ! 1282 DEALLOCATE ( mill_k ) 1283 DEALLOCATE ( igk_l2g_kdip ) 1284 ! 1285 ! ... here one should check for consistency between what is read 1286 ! ... and what is expected 1287 ! 1288 IF ( nbnd_ < nbnd ) THEN 1289 WRITE (msg,'("The number of bands for this run is",I6,", but only",& 1290 & I6," bands were read from file")') nbnd, nbnd_ 1291 CALL errore ('pw_restart - read_collected_wfc', msg, 1 ) 1292 END IF 1293 ! 1294 RETURN 1295 ! 1296 END SUBROUTINE read_collected_wfc 1297 ! 1298 !------------------------------------------------------------------------ 1299 END MODULE pw_restart_new 1300