1! Copyright (C) 2019 Quantum ESPRESSO foundation 2! This file is distributed under the terms of the 3! GNU General Public License. See the file `License' 4! in the root directory of the present distribution, 5! or http://www.gnu.org/copyleft/gpl.txt . 6! 7!---------------------------------------------------------------------------- 8MODULE qexsd_init 9 !---------------------------------------------------------------------------- 10 ! 11 ! This module contains some common subroutines used to copy data used by 12 ! the Quantum ESPRESSO package into XML format 13 ! 14 ! Written by Paolo Giannozzi, building upon pre-existing code qexsd.f90 15 ! 16 ! 17 USE kinds, ONLY : DP 18 ! 19 USE qes_types_module 20 USE qes_reset_module, ONLY: qes_reset 21 USE qes_init_module, ONLY: qes_init 22 ! FIXME: none of the following modules should be used here 23 USE constants, ONLY : e2 24 USE mp_world, ONLY : nproc 25 USE mp_images, ONLY : nimage,nproc_image 26 USE mp_pools, ONLY : npool 27 USE mp_bands, ONLY : ntask_groups, nproc_bgrp, nbgrp 28 ! 29 IMPLICIT NONE 30 ! 31 PRIVATE 32 SAVE 33 ! 34 ! type of objects 35 ! 36 TYPE (berryPhaseOutput_type), TARGET :: qexsd_bp_obj 37 TYPE (k_points_IBZ_type) :: qexsd_start_k_obj 38 TYPE (occupations_type) :: qexsd_occ_obj 39 ! 40 PUBLIC :: qexsd_bp_obj, qexsd_start_k_obj, qexsd_occ_obj 41 ! 42 ! public subroutines. They all work in the same way: 43 ! call qexsd_init_*( xml object, list of QE variables) 44 ! copies QE variables into the xml object 45 ! 46 PUBLIC :: qexsd_init_convergence_info, qexsd_init_algorithmic_info, & 47 qexsd_init_atomic_species, qexsd_init_atomic_structure, & 48 qexsd_init_symmetries, qexsd_init_basis_set, qexsd_init_dft, & 49 qexsd_init_magnetization, qexsd_init_band_structure, & 50 qexsd_init_total_energy, qexsd_init_forces, qexsd_init_stress, & 51 qexsd_init_dipole_info, qexsd_init_outputElectricField, & 52 qexsd_init_outputPBC, qexsd_init_gate_info, qexsd_init_hybrid, & 53 qexsd_init_dftU, qexsd_init_vdw, qexsd_init_berryPhaseOutput 54 ! 55CONTAINS 56 ! 57 ! 58 !------------------------------------------------------------------------ 59 SUBROUTINE qexsd_init_convergence_info(obj, n_scf_steps, scf_has_converged, scf_error, & 60 optimization_has_converged, n_opt_steps, grad_norm ) 61 !------------------------------------------------------------------------ 62 IMPLICIT NONE 63 ! 64 TYPE(convergence_info_type) :: obj 65 INTEGER, INTENT(IN) :: n_scf_steps 66 LOGICAL, INTENT(IN) :: scf_has_converged 67 REAL(DP), INTENT(IN) :: scf_error 68 LOGICAL, OPTIONAL, INTENT(IN) :: optimization_has_converged 69 INTEGER, OPTIONAL, INTENT(in) :: n_opt_steps 70 REAL(DP),OPTIONAL, INTENT(IN) :: grad_norm 71 ! 72 CHARACTER(27) :: subname="qexsd_init_convergence_info" 73 TYPE(scf_conv_type) :: scf_conv 74 TYPE(opt_conv_type),POINTER :: opt_conv 75 ! 76 NULLIFY(opt_conv) 77 call qes_init (scf_conv, "scf_conv", scf_has_converged, n_scf_steps, scf_error) 78 ! 79 IF ( PRESENT(optimization_has_converged )) THEN 80 ! 81 IF ( .NOT. PRESENT(n_opt_steps) ) CALL errore(subname,"n_opt_steps not present",10) 82 IF ( .NOT. PRESENT(grad_norm) ) CALL errore(subname,"grad_norm not present",10) 83 ALLOCATE ( opt_conv) 84 ! 85 call qes_init (opt_conv, "opt_conv", optimization_has_converged, n_opt_steps, grad_norm) 86 ENDIF 87 ! 88 call qes_init (obj, "convergence_info", scf_conv, opt_conv) 89 ! 90 call qes_reset (scf_conv) 91 IF (ASSOCIATED(opt_conv)) THEN 92 CALL qes_reset (opt_conv) 93 NULLIFY ( opt_conv) 94 END IF 95 ! 96 END SUBROUTINE qexsd_init_convergence_info 97 ! 98 ! 99 !------------------------------------------------------------------------ 100 SUBROUTINE qexsd_init_algorithmic_info(obj, real_space_beta, real_space_q, uspp, paw ) 101 !------------------------------------------------------------------------ 102 IMPLICIT NONE 103 ! 104 TYPE(algorithmic_info_type) :: obj 105 LOGICAL, INTENT(IN) :: real_space_beta, real_space_q, uspp, paw 106 ! 107 CALL qes_init (obj, "algorithmic_info", REAL_SPACE_Q = real_space_q, & 108 REAL_SPACE_BETA = real_space_beta, USPP = uspp, PAW = paw) 109 ! 110 END SUBROUTINE qexsd_init_algorithmic_info 111 ! 112 ! 113 !------------------------------------------------------------------------ 114 SUBROUTINE qexsd_init_atomic_species(obj, nsp, atm, psfile, amass, starting_magnetization,& 115 angle1,angle2) 116 !------------------------------------------------------------------------ 117 IMPLICIT NONE 118 ! 119 TYPE(atomic_species_type) :: obj 120 INTEGER, INTENT(IN) :: nsp 121 CHARACTER(len=*), INTENT(IN) :: atm(:) 122 CHARACTER(len=*), INTENT(IN) :: psfile(:) 123 REAL(DP), OPTIONAL,TARGET, INTENT(IN) :: amass(:) 124 REAL(DP), OPTIONAL,TARGET, INTENT(IN) :: starting_magnetization(:) 125 REAL(DP), OPTIONAL,TARGET, INTENT(IN) :: angle1(:),angle2(:) 126 ! 127 TYPE(species_type), ALLOCATABLE :: species(:) 128 REAL(DP),POINTER :: amass_ 129 REAL(DP),POINTER :: start_mag_ 130 REAL(DP),POINTER :: spin_teta 131 REAL(DP),POINTER :: spin_phi 132 INTEGER :: i 133 134 ALLOCATE(species(nsp)) 135 NULLIFY ( amass_, start_mag_, spin_teta, spin_phi) 136 ! 137 DO i = 1, nsp 138 ! 139 IF ( PRESENT(amass) ) THEN 140 IF (amass(i) .GT. 0._DP) amass_=>amass(i) 141 END IF 142 IF ( PRESENT(starting_magnetization) ) THEN 143 IF (ANY( starting_magnetization(1:nsp) /= 0.0_DP)) start_mag_ => starting_magnetization(i) 144 END IF 145 IF ( PRESENT( angle1 ) ) THEN 146 IF (ANY ( angle1(1:nsp) /= 0.0_DP)) spin_teta => angle1(i) 147 END IF 148 IF ( PRESENT( angle2 ) ) THEN 149 IF (ANY(angle2(1:nsp) /= 0.0_DP)) spin_phi => angle2(i) 150 END IF 151 ! 152 CALL qes_init ( species(i), "species", NAME = TRIM(atm(i)), PSEUDO_FILE = TRIM(psfile(i)), MASS = amass_, & 153 STARTING_MAGNETIZATION = start_mag_, SPIN_TETA = spin_teta, SPIN_PHI = spin_phi ) 154 ENDDO 155 ! 156 CALL qes_init (obj, "atomic_species", nsp, species) 157 ! 158 DO i = 1, nsp 159 CALL qes_reset (species(i)) 160 ENDDO 161 DEALLOCATE(species) 162 ! 163 END SUBROUTINE qexsd_init_atomic_species 164 ! 165 ! 166 !------------------------------------------------------------------------ 167 SUBROUTINE qexsd_init_atomic_structure(obj, nsp, atm, ityp, nat, tau, & 168 alat, a1, a2, a3, ibrav) 169 !------------------------------------------------------------------------ 170 IMPLICIT NONE 171 ! 172 TYPE(atomic_structure_type) :: obj 173 INTEGER, INTENT(IN) :: nsp, nat 174 INTEGER, INTENT(in) :: ityp(:) 175 CHARACTER(LEN=*), INTENT(in) :: atm(:) 176 REAL(DP), INTENT(IN) :: tau(3,*)! cartesian atomic positions, a.u. 177 REAL(DP), INTENT(IN) :: alat 178 REAL(DP), INTENT(IN) :: a1(:), a2(:), a3(:) 179 INTEGER, INTENT(IN) :: ibrav 180 ! 181 INTEGER :: ia 182 TYPE(atom_type), ALLOCATABLE :: atom(:) 183 TYPE(cell_type) :: cell 184 TYPE(atomic_positions_type) :: atomic_pos 185 TYPE(wyckoff_positions_type) :: wyckoff_pos 186 REAL(DP) :: new_alat 187 INTEGER,TARGET :: ibrav_tgt 188 INTEGER,POINTER :: ibrav_ptr 189 CHARACTER(LEN=256),POINTER :: use_alt_axes_ 190 CHARACTER(LEN=256),TARGET :: use_alt_axes 191 ! 192 ! atomic positions 193 ! 194 NULLIFY(use_alt_axes_, ibrav_ptr) 195 IF ( ibrav .ne. 0 ) THEN 196 ibrav_tgt = abs(ibrav) 197 ibrav_ptr => ibrav_tgt 198 use_alt_axes_ => use_alt_axes 199 SELECT CASE(ibrav) 200 CASE(-3) 201 use_alt_axes="b:a-b+c:-c" 202 CASE(-5) 203 use_alt_axes="3fold-111" 204 CASE(-9) 205 use_alt_axes="-b:a:c" 206 CASE (91) 207 ibrav_tgt = 9 208 use_alt_axes ="bcoA-type" 209 CASE(-12,-13) 210 use_alt_axes="unique-axis-b" 211 CASE default 212 NULLIFY (use_alt_axes_) 213 END SELECT 214 END IF 215 216 ! 217 ALLOCATE(atom(nat)) 218 DO ia = 1, nat 219 CALL qes_init ( atom(ia), "atom", name=trim(atm(ityp(ia))), atom=tau(1:3,ia), index = ia ) 220 ENDDO 221 ! 222 CALL qes_init (atomic_pos, "atomic_positions", atom) 223 ! 224 DO ia = 1, nat 225 CALL qes_reset ( atom(ia) ) 226 ENDDO 227 DEALLOCATE(atom) 228 ! 229 ! cell 230 ! 231 CALL qes_init (cell, "cell", a1, a2, a3) 232 ! 233 ! global init 234 ! 235 CALL qes_init (obj, "atomic_structure", NAT=nat, ALAT=alat, & 236 ATOMIC_POSITIONS=atomic_pos, CELL=cell , & 237 BRAVAIS_INDEX=ibrav_ptr, ALTERNATIVE_AXES = use_alt_axes_ ) 238 ! 239 ! cleanup 240 ! 241 CALL qes_reset (atomic_pos) 242 CALL qes_reset (cell) 243 ! 244 END SUBROUTINE qexsd_init_atomic_structure 245 ! 246 ! 247 !------------------------------------------------------------------------ 248 SUBROUTINE qexsd_init_symmetries(obj, nsym, nrot, space_group, s, ft, sname, t_rev, nat, irt, & 249 class_names, verbosity, noncolin) 250 !------------------------------------------------------------------------ 251 IMPLICIT NONE 252 ! 253 TYPE(symmetries_type) :: obj 254 INTEGER, INTENT(IN) :: nsym, nrot, nat 255 INTEGER, INTENT(IN) :: space_group 256 INTEGER, INTENT(IN) :: s(:,:,:), irt(:,:) 257 REAL(DP), INTENT(IN) :: ft(:,:) 258 INTEGER, INTENT(IN) :: t_rev(:) 259 CHARACTER(LEN=*), INTENT(IN) :: sname(:), verbosity 260 CHARACTER(LEN=15),INTENT(IN) :: class_names(:) 261 LOGICAL,INTENT(IN) :: noncolin 262 ! 263 TYPE(symmetry_type), ALLOCATABLE :: symm(:) 264 TYPE(equivalent_atoms_type) :: equiv_atm 265 TYPE(info_type) :: info 266 TYPE(matrix_type) :: matrix 267 CHARACTER(LEN=15),POINTER :: classname 268 CHARACTER(LEN=256) :: la_info 269 LOGICAL :: class_ispresent = .FALSE., time_reversal_ispresent = .FALSE. 270 INTEGER :: i 271 REAL(DP) :: mat_(3,3) 272 LOGICAL :: true_=.TRUE., false_ = .FALSE. 273 LOGICAL,POINTER :: trev 274 TARGET :: class_names, true_, false_ 275 ALLOCATE(symm(nrot)) 276 NULLIFY( classname, trev) 277 ! 278 IF ( TRIM(verbosity) .EQ. 'high' .OR. TRIM(verbosity) .EQ. 'medium') class_ispresent= .TRUE. 279 IF ( noncolin ) time_reversal_ispresent = .TRUE. 280 DO i = 1, nrot 281 ! 282 IF (class_ispresent ) classname => class_names(i) 283 IF (time_reversal_ispresent) THEN 284 SELECT CASE (t_rev(i)) 285 CASE (1) 286 trev => true_ 287 CASE default 288 trev => false_ 289 END SELECT 290 END IF 291 IF ( i .LE. nsym ) THEN 292 la_info = "crystal_symmetry" 293 ELSE 294 la_info = "lattice_symmetry" 295 END IF 296 CALL qes_init (info, "info", name=sname(i), class=classname, time_reversal= trev, INFO= TRIM(la_info) ) 297 ! 298 mat_ = real(s(:,:,i),DP) 299 CALL qes_init (matrix, "rotation", DIMS=[3,3], mat=mat_ ) 300 ! 301 IF ( i .LE. nsym ) THEN 302 CALL qes_init (equiv_atm, "equivalent_atoms", nat=nat, equivalent_atoms = irt(i,1:nat) ) 303 ! 304 CALL qes_init (symm(i),"symmetry", info=info, rotation=matrix, fractional_translation=ft(:,i), & 305 equivalent_atoms=equiv_atm) 306 ELSE 307 CALL qes_init ( symm(i), "symmetry", INFO = info, ROTATION = matrix ) 308 END IF 309 ! 310 CALL qes_reset (info) 311 CALL qes_reset (matrix) 312 IF ( i .LT. nsym ) THEN 313 CALL qes_reset ( equiv_atm ) 314 ELSE IF ( i .EQ. nrot ) THEN 315 CALL qes_reset ( equiv_atm ) 316 END IF 317 ! 318 ENDDO 319 ! 320 CALL qes_init (obj,"symmetries",NSYM = nsym, NROT=nrot, SPACE_GROUP = space_group, SYMMETRY=symm ) 321 ! 322 DO i = 1, nsym 323 CALL qes_reset (symm(i)) 324 ENDDO 325 DEALLOCATE(symm) 326 ! 327 END SUBROUTINE qexsd_init_symmetries 328 ! 329 ! 330 !------------------------------------------------------------------------ 331 SUBROUTINE qexsd_init_basis_set(obj, gamma_only, ecutwfc, ecutrho, & 332 nr1, nr2, nr3, nr1s, nr2s, nr3s, & 333 fft_box_ispresent, nr1b, nr2b, nr3b, & 334 ngm, ngms, npwx, b1, b2, b3 ) 335 !------------------------------------------------------------------------ 336 IMPLICIT NONE 337 ! 338 TYPE(basis_set_type) :: obj 339 LOGICAL, INTENT(IN) :: gamma_only 340 INTEGER, INTENT(IN) :: nr1, nr2, nr3 341 INTEGER, INTENT(IN) :: nr1s, nr2s, nr3s 342 LOGICAL, INTENT(IN) :: fft_box_ispresent 343 INTEGER, INTENT(IN) :: nr1b, nr2b, nr3b 344 INTEGER, INTENT(IN) :: ngm, ngms, npwx 345 REAL(DP), INTENT(IN) :: ecutwfc, ecutrho 346 REAL(DP), INTENT(IN) :: b1(3), b2(3), b3(3) 347 ! 348 TYPE(basisSetItem_type) :: fft_grid 349 TYPE(basisSetItem_type) :: fft_smooth 350 TYPE(basisSetItem_type) :: fft_box 351 TYPE(reciprocal_lattice_type) :: recipr_latt 352 353 CALL qes_init (fft_grid, "fft_grid", nr1, nr2, nr3, "") 354 CALL qes_init (fft_smooth, "fft_smooth", nr1s, nr2s, nr3s, "") 355 CALL qes_init (fft_box, "fft_box", nr1b, nr2b, nr3b, "" ) 356 CALL qes_init (recipr_latt, "reciprocal_lattice", b1, b2, b3) 357 358 CALL qes_init (obj, "basis_set", GAMMA_ONLY=gamma_only, ECUTWFC=ecutwfc, ECUTRHO=ecutrho, FFT_GRID=fft_grid, & 359 FFT_SMOOTH=fft_smooth, FFT_BOX=fft_box, NGM=ngm, NGMS=ngms, NPWX=npwx, & 360 RECIPROCAL_LATTICE=recipr_latt ) 361 ! 362 CALL qes_reset(fft_grid) 363 CALL qes_reset(fft_smooth) 364 CALL qes_reset(fft_box) 365 CALL qes_reset(recipr_latt) 366 ! 367 END SUBROUTINE qexsd_init_basis_set 368 ! 369 ! 370 SUBROUTINE qexsd_init_dft (obj, functional, hybrid_, vdW_, dftU_) 371 IMPLICIT NONE 372 TYPE (dft_type),INTENT(INOUT) :: obj 373 CHARACTER(LEN=*),INTENT(IN) :: functional 374 TYPE(hybrid_type),OPTIONAL,INTENT(IN) :: hybrid_ 375 TYPE(vdW_type),OPTIONAL,INTENT(IN) :: vdW_ 376 TYPE(dftU_type),OPTIONAL,INTENT(IN) :: dftU_ 377 ! 378 CALL qes_init(obj, 'dft', functional, hybrid_, dftU_, vdW_) 379 END SUBROUTINE qexsd_init_dft 380 381 !------------------------------------------------------------------------ 382 SUBROUTINE qexsd_init_hybrid ( obj, dft_is_hybrid, nq1, nq2, nq3, ecutfock, exx_fraction, screening_parameter,& 383 exxdiv_treatment, x_gamma_extrapolation, ecutvcut, local_thr ) 384 IMPLICIT NONE 385 TYPE (hybrid_type),INTENT(INOUT) :: obj 386 LOGICAL,INTENT(IN) :: dft_is_hybrid 387 INTEGER,OPTIONAL, INTENT(IN) :: nq1, nq2, nq3 388 REAL(DP),OPTIONAL,INTENT(IN) :: ecutfock, exx_fraction, screening_parameter, ecutvcut,& 389 local_thr 390 CHARACTER(LEN=*), INTENT(IN) :: exxdiv_treatment 391 LOGICAL,OPTIONAL,INTENT(IN) :: x_gamma_extrapolation 392 ! 393 TYPE (qpoint_grid_type),TARGET :: qpoint_grid 394 TYPE (qpoint_grid_type),POINTER :: qpoint_grid_opt 395 ! 396 NULLIFY ( qpoint_grid_opt) 397 IF (.NOT. dft_is_hybrid) RETURN 398 IF (PRESENT(nq1) .AND. PRESENT(nq2) .AND. PRESENT(nq3) ) THEN 399 qpoint_grid_opt => qpoint_grid 400 CALL qes_init (qpoint_grid, "qpoint_grid", nq1, nq2, nq3, "") 401 END IF 402 ! 403 CALL qes_init ( obj, "hybrid", qpoint_grid_opt, ecutfock, exx_fraction, & 404 screening_parameter, exxdiv_treatment, x_gamma_extrapolation, ecutvcut,& 405 local_thr ) 406 ! 407 IF (ASSOCIATED (qpoint_grid_opt)) CALL qes_reset (qpoint_grid_opt) 408 ! 409 END SUBROUTINE qexsd_init_hybrid 410 ! 411 SUBROUTINE qexsd_init_dftU (obj, nsp, psd, species, ityp, is_hubbard, & 412 is_hubbard_back, backall, hubb_l_back, hubb_l1_back, & 413 noncolin, lda_plus_u_kind, U_projection_type, U, U_back, J0, J, & 414 alpha, beta, alpha_back, starting_ns, Hub_ns, Hub_ns_nc ) 415 IMPLICIT NONE 416 TYPE(dftU_type),INTENT(INOUT) :: obj 417 INTEGER,INTENT(IN) :: nsp 418 CHARACTER(LEN=*),INTENT(IN) :: psd(nsp) 419 CHARACTER(LEN=*),INTENT(IN) :: species(nsp) 420 INTEGER,INTENT(IN) :: ityp(:) 421 LOGICAL,INTENT(IN) :: is_hubbard(nsp) 422 LOGICAL,OPTIONAL,INTENT(IN) :: is_hubbard_back(nsp) 423 LOGICAL,OPTIONAL,INTENT(IN) :: backall(nsp) 424 INTEGER,OPTIONAL,INTENT(IN) :: hubb_l_back(nsp) 425 INTEGER,OPTIONAL,INTENT(IN) :: hubb_l1_back(nsp) 426 INTEGER,INTENT(IN) :: lda_plus_u_kind 427 CHARACTER(LEN=*),INTENT(IN) :: U_projection_type 428 LOGICAL,OPTIONAL,INTENT(IN) :: noncolin 429 REAL(DP),OPTIONAL,INTENT(IN) :: U(:), U_back(:), J0(:), alpha(:), alpha_back(:), & 430 beta(:), J(:,:) 431 REAL(DP),OPTIONAL,INTENT(IN) :: starting_ns(:,:,:), Hub_ns(:,:,:,:) 432 COMPLEX(DP),OPTIONAL,INTENT(IN) :: Hub_ns_nc(:,:,:,:) 433 ! 434 CHARACTER(10), ALLOCATABLE :: label(:) 435 TYPE(HubbardCommon_type),ALLOCATABLE :: U_(:), U_back_(:), J0_(:), alpha_(:), & 436 alpha_back_(:), beta_(:) 437 TYPE(HubbardJ_type),ALLOCATABLE :: J_(:) 438 TYPE(starting_ns_type),ALLOCATABLE :: starting_ns_(:) 439 TYPE(Hubbard_ns_type),ALLOCATABLE :: Hubbard_ns_(:), Hubbard_ns_nc_(:) 440 TYPE(HubbardBack_type),ALLOCATABLE :: Hub_back_(:) 441 LOGICAL :: noncolin_ =.FALSE. 442 ! 443 CALL set_labels () 444 IF ( PRESENT(noncolin)) noncolin_ = noncolin 445 ! 446 IF (PRESENT(U)) CALL init_hubbard_commons(U, U_, label, "Hubbard_U") 447 IF (PRESENT(U_back)) CALL init_hubbard_commons(U_back, U_back_, label, "Hubbard_U_back") 448 IF (PRESENT(J0)) CALL init_hubbard_commons(J0, J0_, label, "Hubbard_J0" ) 449 IF (PRESENT(alpha)) CALL init_hubbard_commons(alpha, alpha_,label, "Hubbard_alpha") 450 IF (PRESENT(alpha_back)) CALL init_hubbard_commons(alpha_back, alpha_back_,label, "Hubbard_alpha_back") 451 IF (PRESENT(beta)) CALL init_hubbard_commons(beta, beta_, label, "Hubbard_beta") 452 IF (PRESENT(J)) CALL init_hubbard_J (J, J_, label, "Hubbard_J" ) 453 IF (PRESENT(starting_ns)) CALL init_starting_ns(starting_ns_ , label) 454 IF (PRESENT(Hub_ns)) CALL init_Hubbard_ns(Hubbard_ns_ , label) 455 IF (PRESENT(Hub_ns_nc)) CALL init_Hubbard_ns(Hubbard_ns_nc_ , label) 456 IF (ANY(is_hubbard_back) .AND. PRESENT(hubb_l_back)) & 457 CALL init_Hubbard_back(is_hubbard_back, Hub_back_, hubb_l_back, backall, hubb_l1_back) 458 IF (ANY(is_hubbard_back) .AND. .NOT. PRESENT (hubb_l_back)) & 459 CALL errore('qexsd_init_dft:',& 460 'internal error background is set to true but hubb_l_back is not present',1) 461 ! 462 CALL qes_init (obj, "dftU", lda_plus_u_kind, U_, J0_, alpha_, beta_, J_, starting_ns_, Hubbard_ns_, & 463 U_projection_type, Hub_back_, U_back_, alpha_back_, Hubbard_ns_nc_) 464 ! 465 CALL reset_hubbard_commons(U_) 466 CALL reset_hubbard_commons(U_back_) 467 CALL reset_hubbard_commons(beta_) 468 CALL reset_hubbard_commons(J0_) 469 CALL reset_hubbard_commons(alpha_) 470 CALL reset_hubbard_commons(alpha_back_) 471 CALL reset_hubbard_J(J_) 472 CALL reset_starting_ns(starting_ns_) 473 CALL reset_Hubbard_ns(Hubbard_ns_) 474 ! 475 CONTAINS 476 SUBROUTINE set_labels() 477 IMPLICIT NONE 478 CHARACTER :: hubbard_shell(4)=['s','p','d','f'] 479 INTEGER,EXTERNAL :: set_hubbard_l,set_hubbard_n 480 INTEGER,EXTERNAL :: set_hubbard_l_back,set_hubbard_n_back 481 INTEGER :: i, hubb_l, hubb_n 482 ! 483 ALLOCATE(label(nsp)) 484 DO i = 1, nsp 485 IF (is_hubbard(i)) THEN 486 hubb_l=set_hubbard_l(psd(i)) 487 hubb_n=set_hubbard_n(psd(i)) 488 WRITE (label(i),'(I0,A)') hubb_n,hubbard_shell(hubb_l+1) 489 ELSE 490 label(i)="no Hubbard" 491 ENDIF 492 ENDDO 493 END SUBROUTINE set_labels 494 495 SUBROUTINE init_hubbard_commons(dati, objs, labs, tag) 496 IMPLICIT NONE 497 REAL(DP) :: dati(:) 498 TYPE(HubbardCommon_type),ALLOCATABLE :: objs(:) 499 CHARACTER(LEN=*) :: labs(:), tag 500 INTEGER :: i 501 ! 502 ALLOCATE (objs(nsp)) 503 DO i = 1, nsp 504 CALL qes_init( objs(i), TRIM(tag), TRIM(species(i)), dati(i), TRIM(labs(i))) 505 IF (TRIM(labs(i)) =='no Hubbard') objs(i)%lwrite = .FALSE. 506 END DO 507 END SUBROUTINE init_hubbard_commons 508 ! 509 SUBROUTINE init_hubbard_J(dati, objs, labs, tag) 510 IMPLICIT NONE 511 REAL(DP) :: dati(:,:) 512 TYPE(HubbardJ_type),ALLOCATABLE :: objs(:) 513 CHARACTER(LEN=*) :: labs(:), tag 514 INTEGER :: i 515 ! 516 IF ( SIZE(dati,2) .LE. 0 ) RETURN 517 ALLOCATE (objs(nsp)) 518 DO i = 1, nsp 519 CALL qes_init( objs(i), TRIM(tag), TRIM(species(i)), HubbardJ = dati(1:3,i), LABEL = TRIM(labs(i))) 520 IF (TRIM(labs(i)) =='no Hubbard') objs(i)%lwrite = .FALSE. 521 END DO 522 END SUBROUTINE init_hubbard_J 523 ! 524 SUBROUTINE reset_hubbard_commons(objs) 525 IMPLICIT NONE 526 TYPE(HubbardCommon_type),ALLOCATABLE :: objs(:) 527 INTEGER :: i 528 IF (.NOT. ALLOCATED(objs)) RETURN 529 DO i =1, SIZE(objs) 530 CALL qes_reset(objs(i)) 531 END DO 532 DEALLOCATE(objs) 533 END SUBROUTINE reset_hubbard_commons 534 ! 535 SUBROUTINE reset_hubbard_J(objs) 536 IMPLICIT NONE 537 TYPE(HubbardJ_type),ALLOCATABLE :: objs(:) 538 INTEGER :: i 539 IF (.NOT. ALLOCATED(objs)) RETURN 540 DO i =1, SIZE(objs) 541 CALL qes_reset(objs(i)) 542 END DO 543 DEALLOCATE(objs) 544 END SUBROUTINE reset_hubbard_J 545 ! 546 SUBROUTINE init_starting_ns(objs, labs ) 547 IMPLICIT NONE 548 TYPE(starting_ns_type), ALLOCATABLE :: objs(:) 549 CHARACTER(len=*) :: labs(nsp) 550 INTEGER :: i, is, ind, llmax, nspin 551 ! 552 IF ( .NOT. PRESENT(starting_ns)) RETURN 553 554 IF (noncolin_) THEN 555 llmax = SIZE(starting_ns,1) 556 nspin = 1 557 ALLOCATE(objs(nsp)) 558 DO i = 1, nsp 559 IF (.NOT. ANY(starting_ns(1:2*llmax,1,i)>0.d0)) CYCLE 560 ind = ind + 1 561 CALL qes_init(objs(ind),"starting_ns", TRIM(species(i)), TRIM(labs(i)), 1, & 562 MAX(starting_ns(1:2*llmax,1,i),0._DP)) 563 END DO 564 RETURN 565 ELSE 566 llmax = SIZE (starting_ns, 1) 567 nspin = SIZE(starting_ns, 2) 568 ALLOCATE(objs(nspin*nsp)) 569 ind = 0 570 DO is = 1, nspin 571 DO i = 1, nsp 572 IF (.NOT. ANY (starting_ns(1:llmax,is,i) > 0.d0)) CYCLE 573 ind = ind + 1 574 CALL qes_init(objs(ind), "starting_ns", TRIM(species(i)), TRIM (labs(i)), & 575 is,MAX(starting_ns(1:llmax,is,i),0._DP)) 576 END DO 577 END DO 578 RETURN 579 END IF 580 END SUBROUTINE init_starting_ns 581 ! 582 SUBROUTINE init_Hubbard_ns(objs, labs ) 583 IMPLICIT NONE 584 TYPE (Hubbard_ns_type),ALLOCATABLE :: objs(:) 585 CHARACTER(LEN=*) :: labs(nsp) 586 ! 587 REAL(DP), ALLOCATABLE :: Hubb_occ_aux(:,:) 588 INTEGER :: i, is,ind, ldim, m1, m2, llmax, nat, nspin 589 ! 590 IF (PRESENT(Hub_ns_nc )) THEN 591 llmax = SIZE ( Hub_ns_nc, 1) 592 nat = size(Hub_ns_nc,4) 593 ALLOCATE (objs(nat)) 594 ldim = SIZE(Hub_ns_nc,1) 595 ALLOCATE (Hubb_occ_aux(2*ldim, 2*ldim)) 596 DO i =1, nat 597 Hubb_occ_aux = 0._DP 598 DO m2 = 1, ldim 599 DO m1 =1, ldim 600 Hubb_occ_aux(m1,m2)=SQRT(DCONJG(Hub_ns_nc(m1,m2,1,i))*Hub_ns_nc(m1,m2,1,i)) 601 Hubb_occ_aux(m1,ldim+m2)=SQRT(DCONJG(Hub_ns_nc(m1,m2,2,i))*Hub_ns_nc(m1,m2,2,i)) 602 Hubb_occ_aux(ldim+m1,m2)=SQRT(DCONJG(Hub_ns_nc(m1,m2,3,i))*Hub_ns_nc(m1,m2,3,i)) 603 Hubb_occ_aux(ldim+m1,ldim+m2)=SQRT(DCONJG(Hub_ns_nc(m1,m2,4,i))*Hub_ns_nc(m1,m2,4,i)) 604 END DO 605 END DO 606 CALL qes_init (objs(i), TAGNAME = "Hubbard_ns_mod", SPECIE = TRIM(species(ityp(i))), & 607 LABEL = TRIM(labs(ityp(i))), SPIN =1, INDEX = i,ORDER ='F',Hubbard_NS = Hubb_occ_aux) 608 IF (TRIM(labs(ityp(i))) == 'no Hubbard') objs(i)%lwrite = .FALSE. 609 END DO 610 RETURN 611 ELSE IF (PRESENT (Hub_ns)) THEN 612 llmax = SIZE ( Hub_ns,1) 613 nat = size(Hub_ns,4) 614 nspin = size(Hub_ns,3) 615 ALLOCATE( objs(nspin*nat) ) 616 ind = 0 617 DO i = 1, nat 618 DO is = 1, nspin 619 ind = ind+1 620 CALL qes_init(objs(ind),"Hubbard_ns", SPECIE = TRIM(species(ityp(i))), SPIN = is, & 621 ORDER = 'F', INDEX = ind, LABEL = TRIM(labs(ityp(i))), Hubbard_NS = Hub_ns(:,:,is,i)) 622 IF (TRIM(labs(ityp(i))) =='no Hubbard' ) objs(ind)%lwrite=.FALSE. 623 END DO 624 END DO 625 END IF 626 RETURN 627 ! 628 END SUBROUTINE init_Hubbard_ns 629 630 SUBROUTINE init_Hubbard_back(is_back, objs, l_back, backall_, l1_back) 631 IMPLICIT NONE 632 LOGICAL, INTENT(IN) :: is_back(nsp) 633 INTEGER, INTENT(IN) :: l_back(nsp) 634 TYPE(HubbardBack_type),ALLOCATABLE,INTENT(INOUT) :: objs(:) 635 LOGICAL,OPTIONAL,INTENT(IN) :: backall_(nsp) 636 INTEGER,OPTIONAL,INTENT(IN) :: l1_back(nsp) 637 ! 638 INTEGER :: isp, il, ndimbackL 639 LOGICAL,ALLOCATABLE :: temp(:) 640 TYPE(backL_type) :: backL_objs(2) 641 CHARACTER(LEN=16) :: backchar 642 ! 643 ALLOCATE(objs(nsp), temp(nsp)) 644 IF (PRESENT(backall_)) THEN 645 temp(1:nsp) = backall_(1:nsp) 646 ELSE 647 temp(1:nsp) = .FALSE. 648 END IF 649 DO isp =1, nsp 650 CALL qes_init(backL_objs(1), "l_number", l_index=0, backL = l_back(isp)) 651 ndimbackL = 1 652 IF (temp(isp) .AND. PRESENT(l1_back) ) THEN 653 IF (l1_back(isp) >=0) THEN 654 ndimbackL=2 655 CALL qes_init(backL_objs(2), "l_number", l_index=1, backL = l1_back(isp)) 656 ELSE 657 CALL errore ('qexsd_init_dftU:', 'internal error: l1_back < 0',1) 658 END IF 659 ELSEIF (temp(isp) .AND. .NOT.PRESENT(l1_back) ) THEN 660 CALL errore ('qexsd_init_dftU:', 'internal error: backall is true but l1_back is not present',1) 661 END IF 662 IF (temp(isp)) THEN 663 backchar = 'two_orbitals' 664 ELSE 665 backchar = 'one_orbital' 666 END IF 667 CALL qes_init(objs(isp), "Hubbard_back", SPECIES = TRIM(species(ityp(isp))), & 668 background=TRIM(backchar), l_number = backL_objs(1:ndimbackL)) 669 IF (.NOT. is_back(isp)) objs(isp)%lwrite = .FALSE. 670 DO il = 1, ndimbackL 671 CALL qes_reset(backL_objs(il)) 672 END DO 673 END DO 674 END SUBROUTINE init_Hubbard_back 675 676 677 SUBROUTINE reset_Hubbard_ns(objs) 678 IMPLICIT NONE 679 ! 680 TYPE(hubbard_ns_type),OPTIONAL :: objs(:) 681 INTEGER :: i_ 682 683 IF ( .NOT. PRESENT (objs)) RETURN 684 DO i_ = 1, SIZE(objs) 685 CALL qes_reset(objs(i_)) 686 END DO 687 END SUBROUTINE reset_Hubbard_ns 688 689 SUBROUTINE reset_starting_ns(obj) 690 IMPLICIT NONE 691 TYPE (starting_ns_type), OPTIONAL :: obj(:) 692 INTEGER :: i 693 IF ( .NOT. PRESENT(obj) ) RETURN 694 DO i = 1, SIZE(obj) 695 CALL qes_reset(obj(i)) 696 END DO 697 END SUBROUTINE reset_starting_ns 698 ! 699 700 END SUBROUTINE qexsd_init_dftU 701 ! 702 ! 703 SUBROUTINE qexsd_init_vdw(obj, non_local_term, vdw_corr, vdw_term, ts_thr, ts_isol,& 704 london_s6, london_c6, london_rcut, species, xdm_a1, xdm_a2,& 705 dftd3_version, dftd3_threebody ) 706 IMPLICIT NONE 707 TYPE(vdW_type) :: obj 708 CHARACTER(LEN=*),OPTIONAL,INTENT(IN) :: non_local_term, vdw_corr 709 REAL(DP),OPTIONAL,INTENT(IN) :: vdw_term, london_c6(:), london_rcut, xdm_a1, xdm_a2, ts_thr,& 710 london_s6 711 INTEGER,OPTIONAL,INTENT(IN) :: dftd3_version 712 CHARACTER(LEN=*),OPTIONAL :: species(:) 713 LOGICAL,OPTIONAL,INTENT(IN) :: ts_isol, dftd3_threebody 714 ! 715 LOGICAL :: empirical_vdw = .FALSE. , dft_is_vdw = .FALSE. 716 TYPE(HubbardCommon_type),ALLOCATABLE :: london_c6_obj(:) 717 INTEGER :: isp 718 ! 719 empirical_vdw = PRESENT(vdw_corr) 720 dft_is_vdw = PRESENT(non_local_term) 721 IF ( .NOT. (dft_is_vdW .OR. empirical_vdw)) RETURN 722 IF ( PRESENT (london_c6)) CALL init_londonc6(london_c6, london_c6_obj) 723 CALL qes_init (obj, "vdW", VDW_CORR = vdw_corr, NON_LOCAL_TERM = non_local_term,& 724 TOTAL_ENERGY_TERM = vdw_term, LONDON_S6 = london_s6,& 725 TS_VDW_ECONV_THR = ts_thr, TS_VDW_ISOLATED = ts_isol, LONDON_RCUT = london_rcut, & 726 XDM_A1 = xdm_a1, XDM_A2 = xdm_a2, LONDON_C6 = london_c6_obj, & 727 DFTD3_VERSION = dftd3_version, DFTD3_THREEBODY = dftd3_threebody) 728 ! 729 IF (ALLOCATED(london_c6_obj)) THEN 730 DO isp=1, SIZE(london_c6_obj,1) 731 CALL qes_reset(london_c6_obj(isp)) 732 END DO 733 END IF 734 CONTAINS 735 ! 736 SUBROUTINE init_londonc6(c6data, c6objs ) 737 USE constants, ONLY: eps16 738 IMPLICIT NONE 739 REAL(DP),INTENT(IN) :: c6data(:) 740 TYPE(HubbardCommon_type),ALLOCATABLE,INTENT(INOUT) :: c6objs(:) 741 ! 742 INTEGER :: ndim_london_c6, isp, ind, nsp 743 ! 744 IF (.NOT. PRESENT ( species)) RETURN 745 nsp = SIZE(c6data) 746 ndim_london_c6 = COUNT ( c6data .GT. -eps16) 747 IF ( ndim_london_c6 .GT. 0 ) THEN 748 ALLOCATE (c6objs(ndim_london_c6)) 749 DO isp = 1, nsp 750 IF ( c6data(isp) .GT. -eps16 ) THEN 751 ind = ind + 1 752 CALL qes_init(c6objs(ind ), "london_c6", SPECIE = TRIM(species(isp)), HUBBARDCOMMON = c6data(isp)) 753 END IF 754 END DO 755 END IF 756 END SUBROUTINE init_londonc6 757 ! 758 END SUBROUTINE qexsd_init_vdw 759 !-------------------------------------------------------------------------------------------- 760 SUBROUTINE qexsd_init_outputPBC(obj,assume_isolated) 761 !-------------------------------------------------------------------------------------------- 762 ! 763 IMPLICIT NONE 764 ! 765 TYPE (outputPBC_type) :: obj 766 CHARACTER(LEN=*),INTENT(IN) :: assume_isolated 767 CHARACTER(LEN=*),PARAMETER :: TAGNAME="boundary_conditions" 768 ! 769 CALL qes_init (obj,TAGNAME,ASSUME_ISOLATED =assume_isolated) 770 END SUBROUTINE qexsd_init_outputPBC 771 ! 772 ! 773 !--------------------------------------------------------------------------------------- 774 SUBROUTINE qexsd_init_magnetization(obj, lsda, noncolin, spinorbit, total_mag, total_mag_nc, & 775 absolute_mag, do_magnetization) 776 !------------------------------------------------------------------------------------ 777 IMPLICIT NONE 778 ! 779 TYPE(magnetization_type) :: obj 780 LOGICAL, INTENT(IN) :: lsda, noncolin, spinorbit 781 REAL(DP), INTENT(IN) :: total_mag, absolute_mag 782 REAL(DP), INTENT(IN) :: total_mag_nc(3) 783 LOGICAL, INTENT(IN) :: do_magnetization 784 ! 785 CALL qes_init(obj, "magnetization", lsda, noncolin, spinorbit, total_mag, absolute_mag, & 786 do_magnetization) 787 ! 788 END SUBROUTINE qexsd_init_magnetization 789 ! 790 ! 791 !--------------------------------------------------------------------------------------- 792 SUBROUTINE qexsd_init_band_structure(obj, lsda, noncolin, lspinorb, nelec, n_wfc_at, et, wg, nks, xk, ngk, wk, & 793 starting_kpoints, occupations_kind, wf_collected, & 794 smearing, nbnd, nbnd_up, nbnd_dw, fermi_energy, ef_updw, homo, lumo) 795 !---------------------------------------------------------------------------------------- 796 IMPLICIT NONE 797 ! 798 TYPE(band_structure_type) :: obj 799 CHARACTER(LEN=*), PARAMETER :: TAGNAME="band_structure" 800 LOGICAL,INTENT(IN) :: lsda, noncolin, lspinorb 801 INTEGER,INTENT(IN) :: nks, n_wfc_at 802 INTEGER,OPTIONAL,INTENT(IN) :: nbnd, nbnd_up, nbnd_dw 803 REAL(DP),INTENT(IN) :: nelec 804 REAL(DP),OPTIONAL,INTENT(IN) :: fermi_energy, ef_updw(2), homo, lumo 805 REAL(DP),DIMENSION(:,:),INTENT(IN) :: et, wg, xk 806 REAL(DP),DIMENSION(:),INTENT(IN) :: wk 807 INTEGER,DIMENSION(:),INTENT(IN) :: ngk 808 TYPE(k_points_IBZ_type),INTENT(IN) :: starting_kpoints 809 TYPE(occupations_type), INTENT(IN) :: occupations_kind 810 TYPE(smearing_type),OPTIONAL,INTENT(IN) :: smearing 811 LOGICAL,INTENT(IN) :: wf_collected 812 ! 813 LOGICAL :: n_wfc_at_ispresent = .TRUE. 814 INTEGER :: ndim_ks_energies, ik 815 INTEGER,TARGET :: nbnd_, nbnd_up_, nbnd_dw_ 816 INTEGER,POINTER :: nbnd_opt, nbnd_up_opt, nbnd_dw_opt 817 TYPE(k_point_type) :: kp_obj 818 TYPE(ks_energies_type),ALLOCATABLE :: ks_objs(:) 819 TYPE (k_points_IBZ_type) :: starting_k_points_ 820 REAL(DP),DIMENSION(:),ALLOCATABLE :: eigenvalues, occupations 821 TYPE (smearing_type) :: smearing_ 822 ! 823 ! 824 ndim_ks_energies=nks 825 ! 826 NULLIFY( nbnd_opt, nbnd_up_opt, nbnd_dw_opt) 827 IF ( lsda ) THEN 828 ndim_ks_energies=ndim_ks_energies/2 829 nbnd_up_opt => nbnd_up_ 830 nbnd_dw_opt => nbnd_dw_ 831 IF ( PRESENT(nbnd_up) .AND. PRESENT(nbnd_dw) ) THEN 832 nbnd_ = nbnd_up+nbnd_dw 833 nbnd_up_ = nbnd_up 834 nbnd_dw_ = nbnd_dw 835 ELSE IF ( PRESENT (nbnd) ) THEN 836 nbnd_ = 2*nbnd 837 nbnd_up_ = nbnd 838 nbnd_dw_ = nbnd 839 ELSE 840 CALL errore ( "qexsd:qexsd_init_band_structure: ", & 841 "in case of lsda nbnd_up+nbnd_dw or nbnd must be givens as arguments", 10) 842 END IF 843 ELSE 844 IF (.NOT. PRESENT(nbnd) ) & 845 CALL errore ("qexsd:qexsd_init_band_structure:", "lsda is false but needed nbnd argument is missing", 10) 846 nbnd_=nbnd 847 nbnd_opt => nbnd_ 848 END IF 849 ! 850 ! 851 ALLOCATE(eigenvalues(nbnd_),occupations(nbnd_)) 852 ALLOCATE(ks_objs(ndim_ks_energies)) 853 ! 854 ks_objs%tagname="ks_energies" 855 DO ik=1,ndim_ks_energies 856 CALL qes_init(kp_obj,"k_point",WEIGHT = wk(ik), K_POINT = xk(:,ik)) 857 IF ( lsda ) THEN 858 eigenvalues(1:nbnd_up_)=et(1:nbnd_up_,ik)/e2 859 eigenvalues(nbnd_up_+1:nbnd_)=et(1:nbnd_dw_,ndim_ks_energies+ik)/e2 860 ELSE 861 eigenvalues(1:nbnd_)= et(1:nbnd_,ik)/e2 862 END IF 863 ! 864 ! 865 IF (lsda) THEN 866 IF ( ABS(wk(ik)).GT.1.d-10) THEN 867 occupations(1:nbnd_up_)=wg(1:nbnd_up_,ik)/wk(ik) 868 occupations(nbnd_up_+1:nbnd_)=wg(1:nbnd_dw_,ndim_ks_energies+ik)/wk(ndim_ks_energies+ik) 869 ELSE 870 occupations(1:nbnd_up_)=wg(1:nbnd_up_,ik) 871 occupations(nbnd_up_+1:nbnd_)=wg(1:nbnd_dw_,ik) 872 END IF 873 ELSE 874 IF (ABS(wk(ik)).GT.1.d-10) THEN 875 occupations(1:nbnd_)=wg(1:nbnd_,ik)/wk(ik) 876 ELSE 877 occupations(1:nbnd_)=wg(1:nbnd_,ik) 878 END IF 879 END IF 880 ! 881 ! 882 ks_objs(ik)%k_point = kp_obj 883 ks_objs(ik)%npw = ngk(ik) 884 CALL qes_init(ks_objs(ik)%eigenvalues, "eigenvalues",eigenvalues) 885 CALL qes_init(ks_objs(ik)%occupations, "occupations",occupations) 886 ! 887 eigenvalues=0.d0 888 occupations=0.d0 889 CALL qes_reset(kp_obj) 890 END DO 891 ks_objs%lwrite = .TRUE. 892 ks_objs%lread = .TRUE. 893 ! 894 IF ( PRESENT(smearing) ) smearing_ = smearing 895! 896 starting_k_points_ = starting_kpoints 897 starting_k_points_%tagname = "starting_k_points" 898! 899! 900 CALL qes_init (obj, TAGNAME, LSDA = lsda, NONCOLIN = noncolin, SPINORBIT = lspinorb, NBND = nbnd_opt, & 901 NELEC = nelec, WF_COLLECTED = wf_collected, STARTING_K_POINTS = starting_k_points_, & 902 NKS = ndim_ks_energies, OCCUPATIONS_KIND = occupations_kind, KS_ENERGIES = ks_objs, & 903 NBND_UP = nbnd_up_opt, NBND_DW = nbnd_dw_opt, NUM_OF_ATOMIC_WFC = n_wfc_at, & 904 FERMI_ENERGY = fermi_energy, HIGHESTOCCUPIEDLEVEL = homo, TWO_FERMI_ENERGIES = ef_updw, & 905 SMEARING = smearing, LOWESTUNOCCUPIEDLEVEL = lumo) 906 DO ik=1,ndim_ks_energies 907 CALL qes_reset(ks_objs(ik)) 908 END DO 909 CALL qes_reset( starting_k_points_ ) 910 DEALLOCATE (ks_objs,eigenvalues, occupations) 911 END SUBROUTINE qexsd_init_band_structure 912 ! 913 ! 914 !--------------------------------------------------------------------------------------- 915 SUBROUTINE qexsd_init_total_energy(obj, etot, eband, ehart, vtxc, etxc, ewald, degauss, demet, & 916 electric_field_corr, potentiostat_contr, gate_contribution, dispersion_contribution ) 917 !---------------------------------------------------------------------------------------- 918 ! 919 ! 920 IMPLICIT NONE 921 ! 922 TYPE (total_energy_type) :: obj 923 REAL(DP),INTENT(IN) :: etot, ehart,vtxc,etxc 924 REAL(DP),OPTIONAL,INTENT(IN) :: ewald,demet, eband, degauss 925 REAL(DP),OPTIONAL :: electric_field_corr 926 REAL(DP),OPTIONAL :: potentiostat_contr 927 REAL(DP),OPTIONAL :: gate_contribution 928 REAL(DP),OPTIONAL :: dispersion_contribution 929 ! 930 CHARACTER(LEN=*),PARAMETER :: TAGNAME="total_energy" 931 ! 932 CALL qes_init (obj, TAGNAME, ETOT = etot, EBAND = eband, EHART = ehart, VTXC = vtxc, ETXC = etxc, & 933 EWALD = ewald, DEMET = demet, EFIELDCORR = electric_field_corr, POTENTIOSTAT_CONTR = potentiostat_contr, & 934 GATEFIELD_CONTR = gate_contribution, vdW_term = dispersion_contribution ) 935 936 END SUBROUTINE qexsd_init_total_energy 937 ! 938 ! 939 !-------------------------------------------------------------------------------------------------------- 940 SUBROUTINE qexsd_init_forces(obj,nat,forces,tprnfor) 941 !-------------------------------------------------------------------------------------------------------- 942 ! 943 IMPLICIT NONE 944 ! 945 TYPE(matrix_type) :: obj 946 INTEGER,INTENT(IN) :: nat 947 REAL(DP),DIMENSION(:,:),INTENT(IN) :: forces 948 LOGICAL,INTENT(IN) :: tprnfor 949 ! 950 CHARACTER(LEN=*),PARAMETER :: TAGNAME="forces" 951 REAL(DP),DIMENSION(:,:),ALLOCATABLE :: forces_aux 952 ! 953 IF (.NOT. tprnfor) THEN 954 obj%lwrite=.FALSE. 955 obj%lread =.FALSE. 956 RETURN 957 END IF 958 ! 959 ALLOCATE (forces_aux(3,nat)) 960 forces_aux(1:3,1:nat)=forces(1:3,1:nat)/e2 961 ! 962 CALL qes_init(obj,TAGNAME,[3,nat],forces_aux ) 963 ! 964 DEALLOCATE (forces_aux) 965 ! 966 END SUBROUTINE qexsd_init_forces 967 ! 968 ! 969 !--------------------------------------------------------------------------------------------- 970 SUBROUTINE qexsd_init_stress(obj,stress,tstress) 971 !--------------------------------------------------------------------------------------------- 972 ! 973 IMPLICIT NONE 974 TYPE( matrix_type) :: obj 975 REAL(DP),DIMENSION(3,3),INTENT(IN) :: stress 976 LOGICAL,INTENT(IN) :: tstress 977 ! 978 CHARACTER(LEN=*),PARAMETER :: TAGNAME="stress" 979 REAL(DP),DIMENSION(3,3) :: stress_aux 980 981 IF ( .NOT. tstress ) THEN 982 obj%lwrite = .FALSE. 983 obj%lread = .FALSE. 984 stress_aux = 0.d0 985 RETURN 986 END IF 987 ! 988 stress_aux=stress/e2 989 CALL qes_init(obj,TAGNAME,[3,3],stress_aux ) 990 ! 991 END SUBROUTINE qexsd_init_stress 992 ! 993 ! 994 !------------------------------------------------------------------------------------------------ 995 SUBROUTINE qexsd_init_dipole_info (dipole_info, el_dipole, ion_dipole, edir, eamp, emaxpos, eopreg) 996 !------------------------------------------------------------------------------------------------ 997 ! 998 USE kinds, ONLY : DP 999 USE constants, ONLY : e2, fpi 1000 USE qes_types_module,ONLY : dipoleOutput_type, scalarQuantity_type 1001 USE qes_libs_module, ONLY : qes_init 1002 USE cell_base, ONLY : alat, at, omega 1003 ! 1004 IMPLICIT NONE 1005 ! 1006 TYPE ( dipoleOutput_type ), INTENT(OUT) :: dipole_info 1007 REAL(DP),INTENT(IN) :: el_dipole, ion_dipole, eamp, emaxpos, eopreg 1008 INTEGER , INTENT(IN) :: edir 1009 ! 1010 REAL(DP) :: tot_dipole, length, vamp, fac 1011 TYPE ( scalarQuantity_type) :: temp_qobj 1012 ! 1013 tot_dipole = -el_dipole+ion_dipole 1014 ! 1015 dipole_info%idir = edir 1016 fac=omega/fpi 1017 dipole_info%tagname = "dipoleInfo" 1018 dipole_info%lwrite = .TRUE. 1019 dipole_info%lread = .TRUE. 1020 CALL qes_init (dipole_info%ion_dipole, "ion_dipole" , units="Atomic Units", & 1021 scalarQuantity= ion_dipole*fac) 1022 CALL qes_init (dipole_info%elec_dipole,"elec_dipole" , units="Atomic Units",& 1023 scalarQuantity= el_dipole*fac) 1024 CALL qes_init (dipole_info%dipole,"dipole" , units="Atomic Units", & 1025 scalarQuantity= tot_dipole*fac) 1026 CALL qes_init (dipole_info%dipoleField,"dipoleField" , units="Atomic Units", & 1027 scalarQuantity= tot_dipole) 1028 ! 1029 length=(1._DP-eopreg)*(alat*SQRT(at(1,edir)**2+at(2,edir)**2+at(3,edir)**2)) 1030 vamp=e2*(eamp-tot_dipole)*length 1031 ! 1032 CALL qes_init (dipole_info%potentialAmp,"potentialAmp" , units="Atomic Units",& 1033 scalarQuantity= vamp) 1034 CALL qes_init (dipole_info%totalLength, "totalLength", units = "Bohr",& 1035 scalarQuantity = length ) 1036 1037 END SUBROUTINE qexsd_init_dipole_info 1038 !--------------------------------------------------------------------------------------------- 1039 SUBROUTINE qexsd_init_outputElectricField(obj, lelfield, tefield, ldipole, lberry, bp_obj, el_pol, & 1040 ion_pol, dipole_obj , gateInfo) 1041 !--------------------------------------------------------------------------------------------- 1042 ! 1043 IMPLICIT NONE 1044 ! 1045 TYPE(outputElectricField_type) :: obj 1046 ! 1047 LOGICAL,INTENT(IN) :: lberry, lelfield, tefield, ldipole 1048 REAL(DP),OPTIONAL,INTENT(IN) :: el_pol(:), ion_pol(:) 1049 TYPE(berryPhaseOutput_type),OPTIONAL,INTENT(IN) :: bp_obj 1050 TYPE ( dipoleOutput_type ),OPTIONAL, INTENT(IN) :: dipole_obj 1051 TYPE ( gateInfo_type),OPTIONAL,INTENT(IN) :: gateInfo 1052 ! 1053 CHARACTER(LEN=*),PARAMETER :: TAGNAME="electric_field" 1054 TYPE ( berryPhaseOutput_type ) :: bp_loc_obj 1055 TYPE ( dipoleOutput_type ) :: dip_loc_obj 1056 TYPE ( finiteFieldOut_type ) :: finiteField_obj 1057 LOGICAL :: bp_is = .FALSE. , finfield_is = .FALSE. , & 1058 dipo_is = .FALSE. 1059 ! 1060 1061 IF (lberry .AND. PRESENT ( bp_obj)) THEN 1062 bp_is = .TRUE. 1063 bp_loc_obj = bp_obj 1064 END IF 1065 IF ( lelfield .AND. PRESENT(el_pol) .AND. PRESENT (ion_pol ) ) THEN 1066 finfield_is=.TRUE. 1067 CALL qes_init (finiteField_obj, "finiteElectricFieldInfo", el_pol, ion_pol) 1068 END IF 1069 IF ( ldipole .AND. PRESENT( dipole_obj ) ) THEN 1070 dipo_is = .TRUE. 1071 dip_loc_obj=dipole_obj 1072 END IF 1073 CALL qes_init (obj, TAGNAME, BerryPhase = bp_obj, & 1074 finiteElectricFieldInfo = finiteField_obj, & 1075 dipoleInfo = dipole_obj, & 1076 GATEINFO = gateInfo ) 1077 IF ( finfield_is) CALL qes_reset ( finiteField_obj) 1078 ! 1079 END SUBROUTINE qexsd_init_outputElectricField 1080 ! 1081 !------------------------------------------------------------------------------------------------- 1082 SUBROUTINE qexsd_init_berryPhaseOutput( obj, gpar, gvec, nppstr, nkort, xk, pdl_ion, & 1083 mod_ion, pdl_ion_tot, mod_ion_tot, nstring, pdl_elec, & 1084 mod_elec, wstring, pdl_elec_up, mod_elec_up, pdl_elec_dw,& 1085 mod_elec_dw, pdl_elec_tot,mod_elec_tot, pdl_tot, mod_tot,& 1086 upol, rmod) 1087 !--------------------------------------------------------------------------------------------------- 1088 ! 1089 USE ions_base, ONLY: nat, tau, atm, zv, ityp 1090 USE cell_base, ONLY: omega 1091 USE noncollin_module, ONLY: nspin_lsda 1092 IMPLICIT NONE 1093 ! 1094 TYPE (berryPhaseOutput_type) :: obj 1095 REAL(DP),INTENT(IN) :: gpar(3), gvec, pdl_ion(nat), pdl_ion_tot, xk(3,*) 1096 1097 REAL(DP),INTENT(IN) :: pdl_elec(:), pdl_elec_up, pdl_elec_dw, pdl_elec_tot, & 1098 pdl_tot, upol(3), rmod 1099 ! 1100 INTEGER,INTENT(IN) :: mod_ion(nat), mod_ion_tot, mod_elec(:), mod_elec_up, & 1101 mod_elec_dw, mod_elec_tot, mod_tot, nppstr, nkort, nstring 1102 ! 1103 REAL(DP),INTENT(IN) :: wstring(nstring) 1104 ! 1105 CHARACTER(LEN=*),PARAMETER :: TAGNAME = "BerryPhase" 1106 TYPE ( polarization_type) :: tot_pol_obj 1107 ! 1108 TYPE ( electronicPolarization_type),ALLOCATABLE :: str_pol_obj(:) 1109 TYPE ( ionicPolarization_type ), ALLOCATABLE :: ion_pol_obj(:) 1110 TYPE ( k_point_type ) :: kp_obj 1111 TYPE ( phase_type) :: el_phase, ion_phase, tot_phase 1112 TYPE ( atom_type ) :: atom_obj 1113 TYPE ( scalarQuantity_type ) :: pol_val 1114 INTEGER :: iat, istring, indstring 1115 INTEGER,POINTER :: ispin 1116 INTEGER, TARGET :: spin_val 1117 CHARACTER(10) :: mod_string 1118 LOGICAL :: spin_is = .FALSE. 1119 ! 1120 ALLOCATE (ion_pol_obj(nat)) 1121 ALLOCATE (str_pol_obj(nstring)) 1122 NULLIFY(ispin) 1123 DO iat =1, nat 1124 WRITE(mod_string,'("(mod" ,I1,")")') mod_ion(iat) 1125 CALL qes_init (ion_phase,"phase", modulus = TRIM(mod_string), phase = pdl_ion(iat) ) 1126 CALL qes_init (atom_obj,"ion",name=TRIM(atm(ityp(iat))),atom = tau(:,iat)) 1127 CALL qes_init (ion_pol_obj(iat), "ionicPolarization", atom_obj, zv(ityp(iat)), ion_phase ) 1128 CALL qes_reset (ion_phase) 1129 CALL qes_reset (atom_obj) 1130 END DO 1131 ! 1132 IF ( nspin_lsda .EQ. 2 ) ispin => spin_val 1133 DO istring= 1, nstring 1134 indstring = 1+(istring-1)*nppstr 1135 WRITE(mod_string,'("(mod ",I1,")")') mod_elec(istring) 1136 CALL qes_init(el_phase, "phase", modulus = TRIM (mod_string), phase = pdl_elec(istring) ) 1137 IF ( istring .LE. nstring/nspin_lsda ) THEN 1138 spin_val = 1 1139 ELSE 1140 spin_val = 2 1141 END IF 1142 CALL qes_init(kp_obj, "firstKeyPoint", WEIGHT = wstring(istring), K_POINT = xk(:,indstring)) 1143 CALL qes_init(str_pol_obj(istring),"electronicPolarization", kp_obj, el_phase, ispin ) 1144 CALL qes_reset( el_phase ) 1145 CALL qes_reset(kp_obj) 1146 END DO 1147 ! 1148 WRITE(mod_string,'("(mod ",I1,")")') mod_tot 1149 CALL qes_init (tot_phase, "totalPhase", IONIC = pdl_ion_tot, ELECTRONIC = pdl_elec_tot, & 1150 MODULUS = TRIM(mod_string), PHASE = pdl_tot) 1151 ! 1152 CALL qes_init ( pol_val, "polarization", Units="e/bohr^2", scalarQuantity=(rmod/omega)*pdl_tot ) 1153 ! 1154 CALL qes_init (tot_pol_obj, "totalPolarization", pol_val, modulus = (rmod/omega)*dble(mod_tot), & 1155 direction = upol ) 1156 ! 1157 CALL qes_init ( obj, TAGNAME, totalPolarization = tot_pol_obj, totalPhase = tot_phase, & 1158 ionicPolarization = ion_pol_obj, electronicPolarization = str_pol_obj ) 1159 ! 1160 DO istring=1,nstring 1161 CALL qes_reset (str_pol_obj(istring)) 1162 END DO 1163 DEALLOCATE (str_pol_obj) 1164 DO iat=1, nat 1165 CALL qes_reset (ion_pol_obj(iat)) 1166 END DO 1167 DEALLOCATE (ion_pol_obj) 1168 CALL qes_reset (tot_pol_obj) 1169 CALL qes_reset (pol_val) 1170 CALL qes_reset (tot_phase) 1171 ! 1172 END SUBROUTINE qexsd_init_berryPhaseOutput 1173 ! 1174!----------------------------------------------------------------------------------- 1175SUBROUTINE qexsd_init_gate_info(obj, tagname, gatefield_en, zgate_, nelec_, alat_, at_, bg_, zv_, ityp_) 1176 !-------------------------------------------------------------------------------- 1177 USE kinds, ONLY : DP 1178 USE constants, ONLY : tpi 1179 ! 1180 IMPLICIT NONE 1181 TYPE (gateInfo_type),INTENT(INOUT) :: obj; 1182 CHARACTER(LEN=*) :: tagname 1183 REAL(DP), INTENT(IN) :: gatefield_en, zgate_, alat_, at_(3,3), bg_(3,3), zv_(:), nelec_ 1184 INTEGER,INTENT(IN) :: ityp_(:) 1185 ! 1186 REAL(DP) :: bmod, area, ionic_charge, gateamp, gate_gate_term 1187 ! 1188 bmod=SQRT(bg_(1,3)**2+bg_(2,3)**2+bg_(3,3)**2) 1189 ionic_charge = SUM( zv_(ityp_(:)) ) 1190 area = ABS((at_(1,1)*at_(2,2)-at_(2,1)*at_(1,2))*alat_**2) 1191 gateamp = (-(nelec_-ionic_charge)/area*tpi) 1192 gate_gate_term = (- (nelec_-ionic_charge) * gateamp * (alat_/bmod) / 6.0) 1193 obj = gateInfo_type( TAGNAME = TRIM(tagname), lwrite = .TRUE., lread = .FALSE., POT_PREFACTOR = gateamp, & 1194 GATE_ZPOS = zgate_, GATE_GATE_TERM = gate_gate_term, GATEFIELDENERGY = gatefield_en) 1195 ! 1196END SUBROUTINE qexsd_init_gate_info 1197 1198 1199 1200 END MODULE qexsd_init 1201