quantum-espresso/Modules/qexsd.f90

1253 lines
50 KiB
Fortran

!
! Copyright (C) 2003-2015 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
#ifdef __XSD
!----------------------------------------------------------------------------
MODULE qexsd_module
!----------------------------------------------------------------------------
!
! This module contains some common subroutines used to read and write
! in XML format the data produced by Quantum ESPRESSO package.
!
! Written by Giovanni Borghi, A. Ferretti, ... (2015).
!
! Based on the qexml.f90 routine:
! Written by Andrea Ferretti (2006).
! Initial work by Carlo Sbraccia (xml_io_base.f90)
! Modified by Simone Ziraldo (2013).
!
USE kinds, ONLY : DP
USE input_parameters, ONLY : input_xml_schema_file, title
USE mp_world, ONLY : nproc
USE mp_images, ONLY : nimage,nproc_image
USE mp_pools, ONLY : npool
USE mp_bands, ONLY : ntask_groups, nproc_bgrp, nbgrp
USE global_version, ONLY: version_number, svn_revision
!
USE iotk_base, ONLY : iotk_indent, iotk_maxindent
USE constants, ONLY : e2
USE iotk_module
USE qes_module
!
IMPLICIT NONE
!
PRIVATE
SAVE
!
! definitions for the fmt
!
CHARACTER(5), PARAMETER :: fmt_name = "QEXSD"
CHARACTER(5), PARAMETER :: fmt_version = "0.1.0"
!
! some default for kinds
!
!INTEGER, PARAMETER :: DP = selected_real_kind( 14, 200 )
!
! internal data to be set
!
CHARACTER(256) :: datadir_in, datadir_out
INTEGER :: iunit, ounit
!
! vars to manage back compatibility
!
CHARACTER(10) :: qexsd_current_version = " "
CHARACTER(10) :: qexsd_default_version = trim( fmt_version )
LOGICAL :: qexsd_current_version_init = .FALSE.
!
LOGICAL :: qexsd_use_large_indent = .FALSE.
!
CHARACTER(iotk_attlenx) :: attr
!
TYPE (input_type) :: qexsd_input_obj
TYPE (general_info_type) :: general_info
TYPE (parallel_info_type) :: parallel_info
TYPE (berryPhaseOutput_type) :: qexsd_bp_obj
TYPE (dipoleOutput_type ) :: qexsd_dipol_obj
TYPE ( step_type),ALLOCATABLE :: steps(:)
TYPE ( status_type ) :: exit_status
TYPE ( closed_type ) :: qexsd_closed_element
INTEGER :: step_counter
!
! end of declarations
!
PUBLIC :: qexsd_current_version, qexsd_default_version
PUBLIC :: qexsd_current_version_init
!
PUBLIC :: qexsd_input_obj
!
PUBLIC :: qexsd_init_schema, qexsd_openschema, qexsd_closeschema
!
PUBLIC :: qexsd_init_convergence_info, qexsd_init_algorithmic_info, &
qexsd_init_atomic_species, qexsd_init_atomic_structure, &
qexsd_init_symmetries, qexsd_init_basis_set, qexsd_init_dft, &
qexsd_init_magnetization, qexsd_init_band_structure, &
qexsd_init_total_energy, qexsd_init_forces, qexsd_init_stress, &
qexsd_init_outputElectricField
!
PUBLIC :: qexsd_step_addstep, qexsd_set_status
!
PUBLIC :: qexsd_init_berryPhaseOutput, qexsd_bp_obj, qexsd_dipol_obj
CONTAINS
!
!-------------------------------------------
! ... basic (public) subroutines
!-------------------------------------------
!
!------------------------------------------------------------------------
SUBROUTINE qexsd_init_schema( unit_in, unit_out )
!------------------------------------------------------------------------
!
! just init module data
!
IMPLICIT NONE
INTEGER, INTENT(in) :: unit_in
INTEGER, OPTIONAL, INTENT(in) :: unit_out
!
iunit = unit_in
ounit = unit_in
IF ( present( unit_out ) ) ounit = unit_out
!
! increase the xml indentation for better reading
IF ( qexsd_use_large_indent ) THEN
iotk_indent = 8
iotk_maxindent = 64
ENDIF
!
END SUBROUTINE qexsd_init_schema
!
!
!------------------------------------------------------------------------
FUNCTION check_file_exst( filename )
!------------------------------------------------------------------------
!
IMPLICIT NONE
!
LOGICAL :: check_file_exst
CHARACTER(len=*) :: filename
!
LOGICAL :: lexists
!
INQUIRE( FILE = trim( filename ), EXIST = lexists )
!
check_file_exst = lexists
RETURN
!
END FUNCTION check_file_exst
!
!
!------------------------------------------------------------------------
SUBROUTINE qexsd_openschema( filename )
!------------------------------------------------------------------------
!
IMPLICIT NONE
!
CHARACTER(iotk_attlenx) :: attr
CHARACTER(len=*), INTENT(IN) :: filename
CHARACTER(len=16) :: subname = 'qexsd_openschema'
INTEGER :: ierr, len_steps, i_step
!
! we need a qes-version number here
CALL iotk_write_attr(attr,"xmlns:qes","http://www.quantum-espresso.org/ns/qes/qes-1.0",FIRST=.true.)
!
! we need a schema location number (put a string in fortran-input)?
CALL iotk_write_attr(attr,"xmlns:xsi","http://www.w3.org/2001/XMLSchema-instance",newline=.true.)
!
CALL iotk_write_attr(attr,"xsi:schemaLocation","http://www.quantum-espresso.org/ns/qes/qes-1.0 espresso.xsd",&
newline=.true.)
!
CALL iotk_open_write(ounit,FILE=filename , root="qes:espresso",attr=attr,binary=.false., &
skip_head=.true.,IERR=ierr)
!
IF (ierr /= 0) call errore(subname, 'opening xml input file', ierr)
!
IF ( input_xml_schema_file == ' ' ) CALL get_environment_VARIABLE ( 'QEXSD', input_xml_schema_file )
!
! the input file is mandatory to have a validating schema
! here an error should be issued, instead
!
CALL qexsd_init_general_info(general_info)
CALL qes_write_general_info(ounit,general_info)
CALL qes_reset_general_info(general_info)
!
CALL qexsd_init_parallel_info(parallel_info)
CALL qes_write_parallel_info(ounit,parallel_info)
CALL qes_reset_parallel_info(parallel_info)
IF ( check_file_exst(input_xml_schema_file) ) THEN
CALL qexsd_cp_line_by_line(ounit,input_xml_schema_file, spec_tag="input")
ELSE IF ( TRIM(qexsd_input_obj%tagname) == "input") THEN
CALL qes_write_input(ounit,qexsd_input_obj)
END IF
IF (ALLOCATED(steps) ) THEN
len_steps= step_counter
DO i_step = 1, len_steps
CALL qes_write_step(ounit, steps(i_step) )
END DO
END IF
!
END SUBROUTINE qexsd_openschema
!
!
!---------------------------------------------------------------------------------------
SUBROUTINE qexsd_init_general_info(obj)
!---------------------------------------------------------------------------------------
IMPLICIT NONE
!
TYPE( general_info_type ) :: obj
CHARACTER(LEN=*),PARAMETER :: TAGNAME="general_info"
TYPE( creator_type ) :: creator_obj
TYPE( created_type ) :: created_obj
TYPE( xml_format_type) :: xml_fmt_obj
CHARACTER(LEN=256) :: version
CHARACTER(9) :: cdate, ctime
CHARACTER(60) :: timestamp
!
version=TRIM(version_number)
IF (svn_revision .NE. "unknown") version=TRIM(version)//&
" (svn rev. " // TRIM (svn_revision) // ")"
CALL qes_init_creator(creator_obj, "creator", "PWSCF", version, "XML file generated by PWSCF")
!
CALL date_and_tim(cdate, ctime)
timestamp = 'This run was terminated on: ' // ctime // ' ' // cdate(1:2) // &
' '//cdate(3:5) // ' '// cdate (6:9)
CALL qes_init_created (created_obj, "created", cdate, ctime, timestamp )
!
CALL qes_init_xml_format(xml_fmt_obj, "xml_format", fmt_name, fmt_version, fmt_name//"_"//fmt_version)
!
CALL qes_init_general_info( obj, TAGNAME, xml_fmt_obj, creator = creator_obj, created = created_obj,&
job=title)
!
CALL qes_reset_creator(creator_obj)
CALL qes_reset_created(created_obj)
CALL qes_reset_xml_format(xml_fmt_obj)
END SUBROUTINE qexsd_init_general_info
!
!---------------------------------------------------------------------------------------------
SUBROUTINE qexsd_init_parallel_info(obj)
!---------------------------------------------------------------------------------------------
IMPLICIT NONE
!
TYPE ( parallel_info_type ) :: obj
!
INTEGER :: nthreads=1
#ifdef __OMP
INTEGER,EXTERNAL :: omp_get_max
!
nthreads = omp_get_max()
#endif
CALL qes_init_parallel_info(obj, "parallel_info", nproc, nthreads, ntask_groups, &
nbgrp, npool, nproc_bgrp)
END SUBROUTINE qexsd_init_parallel_info
!
!------------------------------------------------------------------------
SUBROUTINE qexsd_closeschema()
!------------------------------------------------------------------------
USE mytime, ONLY: nclock, clock_label
IMPLICIT NONE
REAL(DP),EXTERNAL :: get_clock
!
CHARACTER(len=17) :: subname = 'qexsd_closeschema'
INTEGER :: ierr
!
IF (TRIM(exit_status%tagname) == "status") THEN
CALL qes_write_status(ounit, exit_status)
CALL qexsd_set_closed()
CALL iotk_write_begin(ounit, "cputime", attr="",new_line=.FALSE.)
!
WRITE(ounit, '(I0)', advance = 'no' ) nint(get_clock('PWSCF'))
CALL iotk_write_end(ounit, "cputime",indentation=.FALSE.)
CALL qes_write_closed(ounit, qexsd_closed_element)
CALL iotk_close_write(ounit, IERR=ierr)
!
CALL errore(subname, 'closing xml input file', ierr)
!
END IF
END SUBROUTINE qexsd_closeschema
!
!
!------------------------------------------------------------------------
SUBROUTINE qexsd_cp_line_by_line(iun_out,filename,spec_tag)
!------------------------------------------------------------------------
implicit none
!
integer, intent(in) :: iun_out
character(*), intent(in) :: filename
character(*), optional, intent(in) :: spec_tag
!
integer :: iun, ierr
character(256) :: str
logical :: icopy, exists
call iotk_free_unit(iun)
!
INQUIRE(FILE=trim(filename), EXIST=exists)
!
IF(.not.exists) THEN
CALL errore('qexsd_cp_line_by_line', 'input xml file "' // &
& TRIM(filename) // '" not found', 1)
ENDIF
!
open(iun,FILE=trim(filename),status="old", IOSTAT=ierr)
!
icopy=.false.
copy_loop: do
!
read(iun,"(a256)",iostat=ierr) str
if (ierr<0) exit copy_loop
if (present(spec_tag)) then
!
if (index(str,"<"//trim(adjustl(spec_tag))//">")/=0) then
!
icopy=.true.
!
endif
!
else
!
icopy=.true.
!
endif
!
! filtering
!
if ( index(str,"<Root>")/=0 .or. index(str,"<Root>")/=0 .or. &
index(str,"<?")/=0 .or. .not.icopy) cycle copy_loop
!
write(iun_out,"(a)") trim(str)
!
if (present(spec_tag)) then
if (index(str,"</input>")/=0) icopy=.false.
endif
!
enddo copy_loop
!
close(iun)
!
END SUBROUTINE qexsd_cp_line_by_line
!
!
!-------------------------------------------
! ... write subroutines
!-------------------------------------------
!
!
!------------------------------------------------------------------------
SUBROUTINE qexsd_init_convergence_info(obj, n_scf_steps, scf_error, &
opt_conv_ispresent, n_opt_steps, grad_norm )
!------------------------------------------------------------------------
IMPLICIT NONE
!
TYPE(convergence_info_type) :: obj
INTEGER, INTENT(IN) :: n_scf_steps
REAL(DP), INTENT(IN) :: scf_error
LOGICAL, INTENT(IN) :: opt_conv_ispresent
INTEGER, OPTIONAL, INTENT(in) :: n_opt_steps
REAL(DP),OPTIONAL, INTENT(IN) :: grad_norm
!
CHARACTER(27) :: subname="qexsd_init_convergence_info"
TYPE(scf_conv_type) :: scf_conv
TYPE(opt_conv_type) :: opt_conv
!
call qes_init_scf_conv(scf_conv, "scf_conv", n_scf_steps, scf_error)
!
IF ( opt_conv_ispresent ) THEN
!
IF ( .NOT. PRESENT(n_opt_steps) ) CALL errore(subname,"n_opt_steps not present",10)
IF ( .NOT. PRESENT(grad_norm) ) CALL errore(subname,"grad_norm not present",10)
!
call qes_init_opt_conv(opt_conv, "opt_conv", n_opt_steps, grad_norm)
ENDIF
!
call qes_init_convergence_info(obj, "convergence_info", scf_conv, opt_conv_ispresent, opt_conv)
!
call qes_reset_scf_conv(scf_conv)
call qes_reset_opt_conv(opt_conv)
!
END SUBROUTINE qexsd_init_convergence_info
!
!
!------------------------------------------------------------------------
SUBROUTINE qexsd_init_algorithmic_info(obj, real_space_q, uspp, paw )
!------------------------------------------------------------------------
IMPLICIT NONE
!
TYPE(algorithmic_info_type) :: obj
LOGICAL, INTENT(IN) :: real_space_q, uspp, paw
!
CALL qes_init_algorithmic_info(obj, "algorithmic_info", real_space_q, uspp, paw)
!
END SUBROUTINE qexsd_init_algorithmic_info
!
!
!------------------------------------------------------------------------
SUBROUTINE qexsd_init_atomic_species(obj, nsp, atm, psfile, amass, starting_magnetization,&
angle1,angle2)
!------------------------------------------------------------------------
IMPLICIT NONE
!
TYPE(atomic_species_type) :: obj
INTEGER, INTENT(IN) :: nsp
CHARACTER(len=*), INTENT(IN) :: atm(:)
CHARACTER(len=*), INTENT(IN) :: psfile(:)
REAL(DP), OPTIONAL, INTENT(IN) :: amass(:)
REAL(DP), OPTIONAL, INTENT(IN) :: starting_magnetization(:)
REAL(DP), OPTIONAL, INTENT(IN) :: angle1(:),angle2(:)
!
TYPE(species_type), ALLOCATABLE :: species(:)
REAL(DP) :: amass_ = 0.0d0
REAL(DP) :: start_mag_ = 0.0d0
REAL(DP) :: spin_teta = 0.0d0
REAL(DP) :: spin_phi = 0.0d0
INTEGER :: i
ALLOCATE(species(nsp))
!
DO i = 1, nsp
!
IF ( PRESENT(amass) ) amass_=amass(i)
IF ( PRESENT(starting_magnetization) ) start_mag_=starting_magnetization(i)
IF ( PRESENT( angle1 ) ) spin_teta =angle1(i)
IF ( PRESENT( angle2 ) ) spin_phi = angle2(i)
!
CALL qes_init_species( species(i), "species", TRIM(atm(i)),PRESENT(amass),amass_, &
TRIM(psfile(i)), PRESENT(starting_magnetization), start_mag_,&
PRESENT(angle1),spin_teta,PRESENT(angle2),spin_phi)
ENDDO
!
CALL qes_init_atomic_species(obj, "atomic_species", nsp, SIZE(species), species)
!
DO i = 1, nsp
CALL qes_reset_species(species(i))
ENDDO
DEALLOCATE(species)
!
END SUBROUTINE qexsd_init_atomic_species
!
!
!------------------------------------------------------------------------
SUBROUTINE qexsd_init_atomic_structure(obj, nsp, atm, ityp, nat, tau, tau_units, &
alat, a1, a2, a3)
!------------------------------------------------------------------------
IMPLICIT NONE
!
TYPE(atomic_structure_type) :: obj
INTEGER, INTENT(IN) :: nsp, nat
INTEGER, INTENT(in) :: ityp(:)
CHARACTER(LEN=*), INTENT(in) :: atm(:)
REAL(DP), INTENT(IN) :: tau(3,*)
CHARACTER(LEN=*), INTENT(IN) :: tau_units
REAL(DP), INTENT(IN) :: alat
REAL(DP), INTENT(IN) :: a1(:), a2(:), a3(:)
!
INTEGER :: ia
TYPE(atom_type), ALLOCATABLE :: atom(:)
TYPE(cell_type) :: cell
TYPE(atomic_positions_type) :: atomic_pos
TYPE(wyckoff_positions_type) :: wyckoff_pos
REAL(DP) :: new_alat
!
! atomic positions
!
ALLOCATE(atom(nat))
DO ia = 1, nat
CALL qes_init_atom( atom(ia), "atom", name=trim(atm(ityp(ia))), &
position="", position_ispresent=.FALSE., atom=tau(1:3,ia))
ENDDO
!
CALL qes_init_atomic_positions(atomic_pos, "atomic_positions", SIZE(atom), atom)
!
DO ia = 1, nat
CALL qes_reset_atom( atom(ia) )
ENDDO
DEALLOCATE(atom)
!
! cell
!
CALL qes_init_cell(cell, "cell", a1, a2, a3)
!
! global init
!
CALL qes_init_atomic_structure(obj, "atomic_structure", nat=nat, nat_ispresent=.TRUE., &
alat=alat, alat_ispresent=.TRUE., atomic_positions_ispresent=.TRUE., &
atomic_positions=atomic_pos, wyckoff_positions_ispresent=.FALSE., &
wyckoff_positions=wyckoff_pos, cell=cell )
!
! cleanup
!
CALL qes_reset_atomic_positions(atomic_pos)
CALL qes_reset_wyckoff_positions(wyckoff_pos)
CALL qes_reset_cell(cell)
!
END SUBROUTINE qexsd_init_atomic_structure
!
!
!------------------------------------------------------------------------
SUBROUTINE qexsd_init_symmetries(obj, nsym, nrot, space_group, s, ft, sname, t_rev, nat, irt, class_names,&
verbosity, noncolin)
!------------------------------------------------------------------------
IMPLICIT NONE
!
TYPE(symmetries_type) :: obj
INTEGER, INTENT(IN) :: nsym, nrot, nat
INTEGER, INTENT(IN) :: space_group
INTEGER, INTENT(IN) :: s(:,:,:), irt(:,:)
REAL(DP), INTENT(IN) :: ft(:,:)
INTEGER, INTENT(IN) :: t_rev(:)
CHARACTER(LEN=*), INTENT(IN) :: sname(:), verbosity
CHARACTER(LEN=15),INTENT(IN) :: class_names(:)
LOGICAL,INTENT(IN) :: noncolin
!
TYPE(symmetry_type), ALLOCATABLE :: symm(:)
TYPE(equivalent_atoms_type) :: equiv_atm
TYPE(info_type) :: info
TYPE(matrix_type) :: matrix
CHARACTER(LEN=15) :: classname
LOGICAL :: class_ispresent = .FALSE., time_reversal_ispresent = .FALSE.
INTEGER :: i
ALLOCATE(symm(nsym))
!
IF ( TRIM(verbosity) .EQ. 'high' .OR. TRIM(verbosity) .EQ. 'medium') class_ispresent= .TRUE.
IF ( noncolin ) time_reversal_ispresent = .TRUE.
DO i = 1, nsym
!
classname = class_names(i)
CALL qes_init_info(info, "info", name=sname(i), name_ispresent=.TRUE., &
class=classname, class_ispresent = class_ispresent, &
time_reversal=(t_rev(i)==1), time_reversal_ispresent = time_reversal_ispresent, &
info= '')
!
CALL qes_init_matrix(matrix, "rotation", ndim1_mat=3, ndim2_mat=3, mat=real(s(:,:,i),DP))
!
CALL qes_init_equivalent_atoms(equiv_atm, "equivalent_atoms", nat=nat, ndim_index_list=nat, &
index_list=irt(i,1:nat) )
!
CALL qes_init_symmetry(symm(i),"symmetry", info=info, rotation=matrix, &
fractional_translation_ispresent=.TRUE., fractional_translation=ft(:,i), &
equivalent_atoms_ispresent=.TRUE., equivalent_atoms=equiv_atm)
!
CALL qes_reset_info(info)
CALL qes_reset_matrix(matrix)
CALL qes_reset_equivalent_atoms(equiv_atm)
!
ENDDO
!
CALL qes_init_symmetries(obj,"symmetries",nsym=nsym, nrot=nrot, space_group=space_group, &
ndim_symmetry=SIZE(symm), symmetry=symm )
!
DO i = 1, nsym
CALL qes_reset_symmetry(symm(i))
ENDDO
DEALLOCATE(symm)
!
END SUBROUTINE qexsd_init_symmetries
!
!
!------------------------------------------------------------------------
SUBROUTINE qexsd_init_basis_set(obj, gamma_only, ecutwfc, ecutrho, &
nr1, nr2, nr3, nr1s, nr2s, nr3s, &
fft_box_ispresent, nr1b, nr2b, nr3b, &
ngm, ngms, npwx, b1, b2, b3 )
!------------------------------------------------------------------------
IMPLICIT NONE
!
TYPE(basis_set_type) :: obj
LOGICAL, INTENT(IN) :: gamma_only
INTEGER, INTENT(IN) :: nr1, nr2, nr3
INTEGER, INTENT(IN) :: nr1s, nr2s, nr3s
LOGICAL, INTENT(IN) :: fft_box_ispresent
INTEGER, INTENT(IN) :: nr1b, nr2b, nr3b
INTEGER, INTENT(IN) :: ngm, ngms, npwx
REAL(DP), INTENT(IN) :: ecutwfc, ecutrho
REAL(DP), INTENT(IN) :: b1(3), b2(3), b3(3)
!
TYPE(basisSetItem_type) :: fft_grid
TYPE(basisSetItem_type) :: fft_smooth
TYPE(basisSetItem_type) :: fft_box
TYPE(reciprocal_lattice_type) :: recipr_latt
CALL qes_init_basisSetItem(fft_grid, "fft_grid", nr1, nr2, nr3, "")
CALL qes_init_basisSetItem(fft_smooth, "fft_smooth", nr1s, nr2s, nr3s, "")
CALL qes_init_basisSetItem(fft_box, "fft_box", nr1b, nr2b, nr3b, "" )
CALL qes_init_reciprocal_lattice(recipr_latt, "reciprocal_lattice", b1, b2, b3)
CALL qes_init_basis_set(obj, "basis_set", gamma_only_ispresent=.TRUE., gamma_only=gamma_only, &
ecutwfc=ecutwfc, ecutrho_ispresent=.TRUE., ecutrho=ecutrho, fft_grid=fft_grid, &
fft_smoooth_ispresent=.TRUE., fft_smoooth=fft_smooth, &
fft_box_ispresent=fft_box_ispresent, fft_box=fft_box, ngm=ngm, &
ngms_ispresent=.TRUE., ngms=ngms, npwx=npwx, reciprocal_lattice=recipr_latt)
!
CALL qes_reset_basisSetItem(fft_grid)
CALL qes_reset_basisSetItem(fft_smooth)
CALL qes_reset_basisSetItem(fft_box)
CALL qes_reset_reciprocal_lattice(recipr_latt)
!
END SUBROUTINE qexsd_init_basis_set
!
!
!------------------------------------------------------------------------
SUBROUTINE qexsd_init_dft(obj, functional, &
dft_is_hybrid, nqx1, nqx2, nqx3, ecutfock, exx_fraction, screening_parameter, &
exxdiv_treatment, x_gamma_extrapolation, ecutvcut, &
dft_is_lda_plus_U, lda_plus_U_kind, llmax, nspin, nsp, ldim, nat, species, ityp, Hubbard_U, Hubbard_J0, &
Hubbard_alpha, Hubbard_beta, Hubbard_J, starting_ns, Hubbard_ns, Hubbard_ns_nc, &
U_projection_type, &
dft_is_vdW, vdw_corr, london_s6, london_rcut, xdm_a1, xdm_a2 ,is_hubbard,psd)
!------------------------------------------------------------------------
USE parameters, ONLY: lqmax
USE input_parameters, ONLY: nspinx
IMPLICIT NONE
!
TYPE(dft_type) :: obj
CHARACTER(len=*), INTENT(IN) :: functional
LOGICAL, INTENT(IN) :: dft_is_hybrid
INTEGER, INTENT(IN) :: nqx1, nqx2, nqx3
REAL(DP), INTENT(IN) :: ecutfock
REAL(DP), INTENT(IN) :: exx_fraction
REAL(DP), INTENT(IN) :: screening_parameter
CHARACTER(len=*), INTENT(IN) :: exxdiv_treatment
LOGICAL, INTENT(IN) :: x_gamma_extrapolation
REAL(DP), INTENT(IN) :: ecutvcut
!
LOGICAL, INTENT(IN) :: dft_is_lda_plus_U
INTEGER, INTENT(IN) :: lda_plus_U_kind
INTEGER, INTENT(IN) :: llmax, nspin, nsp, ldim, nat
CHARACTER(len=*), INTENT(IN) :: species(nsp)
INTEGER, INTENT(IN) :: ityp(nat)
REAL(DP), INTENT(IN) :: Hubbard_U(nsp)
REAL(DP), INTENT(IN) :: Hubbard_J0(nsp)
REAL(DP), INTENT(IN) :: Hubbard_alpha(nsp)
REAL(DP), INTENT(IN) :: Hubbard_beta(nsp)
REAL(DP), INTENT(IN) :: Hubbard_J(3,nsp)
REAL(DP), INTENT(IN) :: starting_ns(lqmax,nspinx,nsp)
REAL(DP), INTENT(IN) :: Hubbard_ns(ldim,ldim,nspin,nat)
COMPLEX(DP), INTENT(IN) :: Hubbard_ns_nc(ldim,ldim,nspin,nat)
CHARACTER(len=*), INTENT(IN) :: U_projection_type
LOGICAL,INTENT(IN) :: is_hubbard(nsp)
CHARACTER(LEN=2),INTENT(IN) :: psd(nsp)
!
LOGICAL, INTENT(IN) :: dft_is_vdW
CHARACTER(len=*), INTENT(IN) :: vdw_corr
REAL(DP), INTENT(IN) :: london_s6
REAL(DP), INTENT(IN) :: london_rcut
REAL(DP), INTENT(IN) :: xdm_a1
REAL(DP), INTENT(IN) :: xdm_a2
!
INTEGER :: i, is, ind,hubb_l,hubb_n
TYPE(hybrid_type) :: hybrid
TYPE(qpoint_grid_type) :: qpoint_grid
TYPE(dftU_type) :: dftU
TYPE(vdW_type) :: vdW
TYPE(HubbardCommon_type), ALLOCATABLE :: Hubbard_U_(:)
TYPE(HubbardCommon_type), ALLOCATABLE :: Hubbard_J0_(:)
TYPE(HubbardCommon_type), ALLOCATABLE :: Hubbard_alpha_(:)
TYPE(HubbardCommon_type), ALLOCATABLE :: Hubbard_beta_(:)
TYPE(HubbardJ_type), ALLOCATABLE :: Hubbard_J_(:)
TYPE(starting_ns_type), ALLOCATABLE :: starting_ns_(:)
TYPE(Hubbard_ns_type), ALLOCATABLE :: Hubbard_ns_(:)
LOGICAL :: Hubbard_U_ispresent
LOGICAL :: Hubbard_J0_ispresent
LOGICAL :: Hubbard_alpha_ispresent
LOGICAL :: Hubbard_beta_ispresent
LOGICAL :: Hubbard_J_ispresent
LOGICAL :: starting_ns_ispresent
LOGICAL :: Hubbard_ns_ispresent
CHARACTER(10), ALLOCATABLE :: label(:)
CHARACTER :: hubbard_shell
INTEGER,EXTERNAL :: set_hubbard_l,set_hubbard_n
!
!
IF ( dft_is_hybrid ) THEN
!
CALL qes_init_qpoint_grid(qpoint_grid, "qpoint_grid", nqx1, nqx2, nqx3, "")
!
CALL qes_init_hybrid(hybrid, "hybrid", qpoint_grid, ecutfock, exx_fraction, &
screening_parameter, exxdiv_treatment, x_gamma_extrapolation, ecutvcut)
!
CALL qes_reset_qpoint_grid(qpoint_grid)
!
ENDIF
!
IF ( dft_is_lda_plus_U ) THEN
!
ALLOCATE(label(nsp))
DO i = 1, nsp
IF (is_hubbard(i)) THEN
hubb_l=set_hubbard_l(psd(i))
hubb_n=set_hubbard_n(psd(i))
SELECT CASE ( hubb_l )
CASE ( 0)
hubbard_shell='s'
CASE ( 1 )
hubbard_shell='p'
CASE( 2 )
hubbard_shell='d'
CASE( 3 )
hubbard_shell='f'
END SELECT
WRITE (label(i),'(I0,A)') hubb_n,hubbard_shell
ELSE
label(i)="no Hubbard"
END IF
END DO
!
Hubbard_U_ispresent = (SIZE(Hubbard_U)>0)
Hubbard_J0_ispresent = (SIZE(Hubbard_J0)>0)
Hubbard_alpha_ispresent = (SIZE(Hubbard_alpha)>0)
Hubbard_beta_ispresent = (SIZE(Hubbard_beta)>0)
Hubbard_J_ispresent = (SIZE(Hubbard_J)>0)
Hubbard_ns_ispresent = (SIZE(Hubbard_ns)>0)
starting_ns_ispresent = (SIZE(starting_ns)>0)
!
ALLOCATE( Hubbard_U_(nsp) )
ALLOCATE( Hubbard_J0_(nsp) )
ALLOCATE( Hubbard_alpha_(nsp) )
ALLOCATE( Hubbard_beta_(nsp) )
ALLOCATE( Hubbard_J_(nsp) )
!
ALLOCATE( starting_ns_(min(nspin,nspinx)*nsp) )
ALLOCATE( Hubbard_ns_(nspin*nat) )
!
DO i = 1, nsp
CALL qes_init_HubbardCommon(Hubbard_U_(i),"Hubbard_U",TRIM(species(i)),TRIM(label(i)),Hubbard_U(i))
CALL qes_init_HubbardCommon(Hubbard_J0_(i),"Hubbard_J0",TRIM(species(i)),TRIM(label(i)),Hubbard_J0(i))
CALL qes_init_HubbardCommon(Hubbard_alpha_(i),"Hubbard_alpha",TRIM(species(i)),TRIM(label(i)),&
Hubbard_alpha(i))
CALL qes_init_HubbardCommon(Hubbard_beta_(i),"Hubbard_beta",TRIM(species(i)),TRIM(label(i)),&
Hubbard_beta(i))
CALL qes_init_HubbardJ(Hubbard_J_(i),"Hubbard_J",TRIM(species(i)),TRIM(label(i)),Hubbard_J(1:3,i))
ENDDO
!
ind = 0
DO is = 1, MIN(nspin,nspinx)
DO i = 1, nsp
ind = ind+1
CALL qes_init_starting_ns(starting_ns_(ind),"starting_ns",TRIM(species(i)),TRIM(label(i)), &
is, llmax, starting_ns(1:llmax,is,i) )
ENDDO
ENDDO
!
ind = 0
DO i = 1, nat
DO is = 1, nspin
ind = ind+1
CALL qes_init_Hubbard_ns(Hubbard_ns_(ind),"Hubbard_ns", TRIM(species(ityp(i))),TRIM(label(ityp(i))), &
is, ldim, ldim, Hubbard_ns(:,:,is,i) )
ENDDO
ENDDO
!
! main init
CALL qes_init_dftU(dftU, "dftU", .TRUE., lda_plus_u_kind, &
Hubbard_U_ispresent, SIZE(Hubbard_U_), Hubbard_U_, &
Hubbard_J0_ispresent, SIZE(Hubbard_J0_), Hubbard_J0_, &
Hubbard_alpha_ispresent, SIZE(Hubbard_alpha_), Hubbard_alpha_, &
Hubbard_beta_ispresent, SIZE(Hubbard_beta_), Hubbard_beta_, &
Hubbard_J_ispresent, SIZE(Hubbard_J_), Hubbard_J_, &
starting_ns_ispresent, SIZE(starting_ns_), starting_ns_, &
Hubbard_ns_ispresent, SIZE(Hubbard_ns_), Hubbard_ns_, &
.TRUE., U_projection_type)
!
DO i = 1, nsp
CALL qes_reset_HubbardCommon(Hubbard_U_(i))
CALL qes_reset_HubbardCommon(Hubbard_J0_(i))
CALL qes_reset_HubbardCommon(Hubbard_alpha_(i))
CALL qes_reset_HubbardCommon(Hubbard_beta_(i))
CALL qes_reset_HubbardJ(Hubbard_J_(i))
ENDDO
!
DEALLOCATE(Hubbard_U_)
DEALLOCATE(Hubbard_J0_)
DEALLOCATE(Hubbard_alpha_)
DEALLOCATE(Hubbard_beta_)
DEALLOCATE(Hubbard_J_)
!
DO i = 1, SIZE(starting_ns_)
CALL qes_reset_starting_ns(starting_ns_(i))
ENDDO
DEALLOCATE(starting_ns_)
!
DO i = 1, SIZE(Hubbard_ns_)
CALL qes_reset_Hubbard_ns(Hubbard_ns_(i))
ENDDO
DEALLOCATE(Hubbard_ns_)
!
DEALLOCATE(label)
!
ENDIF
!
IF ( dft_is_vdW ) THEN
!
CALL qes_init_vdW(vdW, "vdW", vdw_corr, london_s6, london_rcut, xdm_a1, xdm_a2 )
!
ENDIF
CALL qes_init_dft(obj, "dft", functional, dft_is_hybrid, hybrid, &
dft_is_lda_plus_U, dftU, dft_is_vdW, vdW)
!
IF (dft_is_hybrid) CALL qes_reset_hybrid(hybrid)
IF (dft_is_lda_plus_U) CALL qes_reset_dftU(dftU)
IF (dft_is_vdW) CALL qes_reset_vdW(vdW)
!
END SUBROUTINE qexsd_init_dft
!
!
!---------------------------------------------------------------------------------------
SUBROUTINE qexsd_init_magnetization(obj, lsda, noncolin, spinorbit, total_mag, total_mag_nc, &
absolute_mag, do_magnetization)
!------------------------------------------------------------------------------------
IMPLICIT NONE
!
TYPE(magnetization_type) :: obj
LOGICAL, INTENT(IN) :: lsda, noncolin, spinorbit
REAL(DP), INTENT(IN) :: total_mag, absolute_mag
REAL(DP), INTENT(IN) :: total_mag_nc(3)
LOGICAL, INTENT(IN) :: do_magnetization
!
CALL qes_init_magnetization(obj, "magnetization", lsda, noncolin, spinorbit, total_mag, absolute_mag, &
do_magnetization)
!
END SUBROUTINE qexsd_init_magnetization
!
!
!---------------------------------------------------------------------------------------
SUBROUTINE qexsd_init_band_structure(obj, lsda,noncolin, lspinorb, nbnd, nelec, &
fermi_energy, et, wg, nks, xk, ngk, wk)
!----------------------------------------------------------------------------------------
IMPLICIT NONE
!
TYPE(band_structure_type) :: obj
CHARACTER(LEN=*), PARAMETER :: TAGNAME="band_structure"
LOGICAL,INTENT(IN) :: lsda, noncolin, lspinorb
INTEGER,INTENT(IN) :: nbnd, nks
REAL(DP),INTENT(IN) :: nelec, fermi_energy
REAL(DP),DIMENSION(:,:),INTENT(IN) :: et, wg, xk
REAL(DP),DIMENSION(:),INTENT(IN) :: wk
INTEGER,DIMENSION(:),INTENT(IN) :: ngk
!
LOGICAL :: nbnd_up_ispresent, nbnd_dw_ispresent, &
fermi_energy_ispresent
INTEGER :: nbnd_up,nbnd_dw
INTEGER :: ndim_ks_energies,nbnd_tot,ik
TYPE(k_point_type) :: kp_obj
TYPE(ks_energies_type),ALLOCATABLE :: ks_objs(:)
REAL(DP),DIMENSION(:),ALLOCATABLE :: eigenvalues,occupations
!
!
ndim_ks_energies=nks
nbnd_tot=nbnd
!
IF ( lsda ) THEN
ndim_ks_energies=ndim_ks_energies/2
nbnd_up=nbnd
nbnd_dw=nbnd
nbnd_tot=nbnd_up+nbnd_dw
nbnd_up_ispresent=.true.
nbnd_dw_ispresent=.true.
ELSE
nbnd_up=0
nbnd_dw=0
nbnd_up_ispresent=.false.
nbnd_dw_ispresent=.false.
END IF
IF (fermi_energy.GT.-1.D6) THEN
fermi_energy_ispresent=.TRUE.
ELSE
fermi_energy_ispresent=.FALSE.
END IF
!
!
ALLOCATE(eigenvalues(nbnd_tot),occupations(nbnd_tot))
ALLOCATE(ks_objs(ndim_ks_energies))
!
DO ik=1,ndim_ks_energies
CALL qes_init_k_point(kp_obj,"k_point",wk(ik),.true.,xk(:,ik))
IF ( lsda ) THEN
eigenvalues(1:nbnd_up)=et(1:nbnd_up,ik)
eigenvalues(nbnd_up+1:nbnd_tot)=et(1:nbnd_dw,ndim_ks_energies+ik)/e2
ELSE
eigenvalues(1:nbnd_tot)= et(1:nbnd_tot,ik)/e2
END IF
!
!
IF (lsda) THEN
IF ( ABS(wk(ik)).GT.1.d-10) THEN
occupations(1:nbnd_up)=wg(1:nbnd_up,ik)/wk(ik)
occupations(nbnd_up+1:nbnd_tot)=wg(1:nbnd_dw,ik)/wk(ik)
ELSE
occupations(1:nbnd_up)=wg(1:nbnd_up,ik)
occupations(nbnd_up+1:nbnd_tot)=wg(1:nbnd_dw,ik)
END IF
ELSE
IF (ABS(wk(ik)).GT.1.d-10) THEN
occupations(1:nbnd_tot)=wg(1:nbnd_tot,ik)/wk(ik)
ELSE
occupations(1:nbnd_tot)=wg(1:nbnd_tot,ik)
END IF
END IF
!
!
CALL qes_init_ks_energies(ks_objs(ik),"ks_energies",kp_obj,ngk(ik),nbnd_tot,eigenvalues,&
nbnd_tot,occupations)
!
eigenvalues=0.d0
occupations=0.d0
CALL qes_reset_k_point(kp_obj)
END DO
!
CALL qes_init_band_structure(obj,TAGNAME,lsda,noncolin,lspinorb,nbnd_tot,nbnd_up_ispresent,&
nbnd_up,nbnd_dw_ispresent,nbnd_dw,nelec,fermi_energy_ispresent,&
fermi_energy/e2, ndim_ks_energies,ndim_ks_energies,ks_objs)
DO ik=1,ndim_ks_energies
CALL qes_reset_ks_energies(ks_objs(ik))
END DO
DEALLOCATE (ks_objs,eigenvalues,occupations)
END SUBROUTINE qexsd_init_band_structure
!
!
!---------------------------------------------------------------------------------------
SUBROUTINE qexsd_init_total_energy(obj,etot,eband,ehart,vtxc,etxc,ewald,degauss,demet,electric_field_corr)
!----------------------------------------------------------------------------------------
!
!
IMPLICIT NONE
!
TYPE (total_energy_type) :: obj
REAL(DP),INTENT(IN) :: etot,eband,ehart,vtxc,etxc,ewald,demet
REAL(DP),INTENT(IN) :: degauss
REAL(DP),OPTIONAL,INTENT(IN) :: electric_field_corr
!
LOGICAL :: demet_ispresent
CHARACTER(LEN=*),PARAMETER :: TAGNAME="total_energy"
REAL(DP) :: etot_har,eband_har,vtxc_har,etxc_har,ewald_har,&
demet_har,ehart_har,efield_corr
etot_har = etot/e2
eband_har = etot/e2
vtxc_har = vtxc/e2
etxc_har = etxc/e2
ewald_har = ewald/e2
ehart_har = ehart/e2
IF (PRESENT(electric_field_corr)) THEN
efield_corr=electric_field_corr/e2
ELSE
efield_corr=0.d0
END IF
IF (degauss .GT. 0.D0) THEN
demet_ispresent=.TRUE.
demet_har=demet/e2
ELSE
demet_ispresent=.FALSE.
demet_har=0.d0
ENDIF
CALL qes_init_total_energy(obj,TAGNAME,etot_har,eband_ispresent=.TRUE.,eband=eband_har,&
ehart_ispresent=.TRUE., ehart=ehart_har,vtxc_ispresent=.TRUE.,&
vtxc=vtxc_har,etxc_ispresent=.TRUE., etxc=etxc_har,ewald_ispresent=.TRUE.,&
ewald=ewald_har, demet_ispresent=demet_ispresent,demet=demet_har, &
efieldcorr_ispresent=PRESENT(electric_field_corr), efieldcorr=efield_corr)
END SUBROUTINE qexsd_init_total_energy
!
!
!--------------------------------------------------------------------------------------------------------
SUBROUTINE qexsd_init_forces(obj,nat,forces,tprnfor)
!--------------------------------------------------------------------------------------------------------
!
IMPLICIT NONE
!
TYPE(matrix_type) :: obj
INTEGER,INTENT(IN) :: nat
REAL(DP),DIMENSION(:,:),INTENT(IN) :: forces
LOGICAL,INTENT(IN) :: tprnfor
!
CHARACTER(LEN=*),PARAMETER :: TAGNAME="forces"
REAL(DP),DIMENSION(:,:),ALLOCATABLE :: forces_aux
!
IF (.NOT. tprnfor) THEN
obj%lwrite=.FALSE.
obj%lread =.FALSE.
RETURN
END IF
!
ALLOCATE (forces_aux(3,nat))
forces_aux(1:3,1:nat)=forces(1:3,1:nat)/e2
!
CALL qes_init_matrix(obj,TAGNAME,3,nat,forces_aux)
!
DEALLOCATE (forces_aux)
!
END SUBROUTINE qexsd_init_forces
!
!
!---------------------------------------------------------------------------------------------
SUBROUTINE qexsd_init_stress(obj,stress,tstress)
!---------------------------------------------------------------------------------------------
!
IMPLICIT NONE
TYPE( matrix_type) :: obj
REAL(DP),DIMENSION(3,3),INTENT(IN) :: stress
LOGICAL,INTENT(IN) :: tstress
!
CHARACTER(LEN=*),PARAMETER :: TAGNAME="stress"
REAL(DP),DIMENSION(3,3) :: stress_aux
IF ( .NOT. tstress ) THEN
obj%lwrite = .FALSE.
obj%lread = .FALSE.
stress_aux = 0.d0
RETURN
END IF
!
stress_aux=stress/e2
CALL qes_init_matrix(obj,TAGNAME,3,3,stress_aux)
!
END SUBROUTINE qexsd_init_stress
!
!
!---------------------------------------------------------------------------------------------
SUBROUTINE qexsd_init_outputElectricField(obj, lelfield, tefield, ldipole, lberry, bp_obj, el_pol, &
ion_pol, dipole_obj )
!---------------------------------------------------------------------------------------------
!
IMPLICIT NONE
!
TYPE(outputElectricField_type) :: obj
!
LOGICAL,INTENT(IN) :: lberry, lelfield, tefield, ldipole
REAL(DP),OPTIONAL,DIMENSION(3),INTENT(IN) :: el_pol, ion_pol
TYPE(berryPhaseOutput_type),OPTIONAL,INTENT(IN) :: bp_obj
TYPE ( dipoleOutput_type ),OPTIONAL, INTENT(IN) :: dipole_obj
!
CHARACTER(LEN=*),PARAMETER :: TAGNAME="electric_field"
TYPE ( berryPhaseOutput_type ) :: bp_loc_obj
TYPE ( dipoleOutput_type ) :: dip_loc_obj
TYPE ( finiteFieldOut_type ) :: finiteField_obj
LOGICAL :: bp_is = .FALSE. , finfield_is = .FALSE. , &
dipo_is = .FALSE.
!
IF (lberry .AND. PRESENT ( bp_obj)) THEN
bp_is = .TRUE.
bp_loc_obj = bp_obj
END IF
IF ( lelfield .AND. PRESENT(el_pol) .AND. PRESENT (ion_pol ) ) THEN
finfield_is=.TRUE.
CALL qes_init_finiteFieldOut (finiteField_obj, "finiteElectricFieldInfo", el_pol, ion_pol)
END IF
IF ( ldipole .AND. PRESENT( dipole_obj ) ) THEN
dipo_is = .TRUE.
dip_loc_obj=dipole_obj
END IF
CALL qes_init_outputElectricField(obj, TAGNAME, BerryPhase_ispresent = bp_is, &
BerryPhase = bp_loc_obj, &
finiteElectricFieldInfo_ispresent = finfield_is, &
finiteElectricFieldInfo = finiteField_obj, &
dipoleInfo_ispresent = dipo_is, dipoleInfo = dip_loc_obj)
IF (dipo_is) CALL qes_reset_dipoleOutput( dip_loc_obj )
IF ( bp_is ) CALL qes_reset_berryPhaseOutput( bp_loc_obj )
!
END SUBROUTINE qexsd_init_outputElectricField
!
!----------------------------------------------------------------------------------------
SUBROUTINE qexsd_step_addstep( i_step, max_steps, ntyp, atm, ityp, nat,&
tau, alat, a1, a2, a3, etot, eband, ehart, vtxc, etxc,&
ewald, degauss, demet, forces, stress, n_scf_steps, scf_error)
!-----------------------------------------------------------------------------------------
IMPLICIT NONE
!
INTEGER ,INTENT(IN) :: i_step, max_steps, ntyp, nat, n_scf_steps, ityp(:)
REAL(DP),INTENT(IN) :: tau(3,nat), alat, a1(3), a2(3), a3(3), etot, eband, ehart, vtxc, &
etxc, ewald, degauss, demet, scf_error, forces(3,nat), stress(3,3)
CHARACTER(LEN=*),INTENT(IN) :: atm(:)
TYPE (step_type) :: step_obj
TYPE ( scf_conv_type ) :: scf_conv_obj
TYPE ( atomic_structure_type ) :: atomic_struct_obj
TYPE ( total_energy_type ) :: tot_en_obj
TYPE ( matrix_type ) :: mat_forces, mat_stress
!
IF ( i_step .EQ. 1 ) THEN
ALLOCATE (steps(max_steps))
step_counter = 0
END IF
step_counter = step_counter+1
PRINT *,"step counter value is ",step_counter
!
step_obj%tagname="step"
step_obj%n_step = step_counter
!
CALL qes_init_scf_conv( scf_conv_obj,"scf_conv", n_scf_steps, scf_error )
!
step_obj%scf_conv = scf_conv_obj
CALL qes_reset_scf_conv(scf_conv_obj)
!
CALL qexsd_init_atomic_structure(atomic_struct_obj, ntyp, atm, ityp, nat, tau, "bohr", &
alat, a1, a2, a3)
step_obj%atomic_structure=atomic_struct_obj
CALL qes_reset_atomic_structure( atomic_struct_obj )
!
CALL qexsd_init_total_energy ( tot_en_obj, etot, eband, ehart, vtxc, etxc, ewald, degauss, demet )
step_obj%total_energy=tot_en_obj
CALL qes_reset_total_energy( tot_en_obj )
!
CALL qes_init_matrix( mat_forces, "forces", 3, nat, forces )
step_obj%forces=mat_forces
CALL qes_reset_matrix ( mat_forces )
!
CALL qes_init_matrix( mat_stress, "stress", 3, 3, stress )
step_obj%stress = mat_stress
CALL qes_reset_matrix ( mat_stress )
!
!
steps(step_counter) = step_obj
call qes_reset_step(step_obj)
END SUBROUTINE qexsd_step_addstep
!-------------------------------------------------------------------------------------------------
SUBROUTINE qexsd_init_berryPhaseOutput( obj, gpar, gvec, nppstr, nkort, xk, pdl_ion, &
mod_ion, pdl_ion_tot, mod_ion_tot, nstring, pdl_elec, &
mod_elec, wstring, pdl_elec_up, mod_elec_up, pdl_elec_dw,&
mod_elec_dw, pdl_elec_tot,mod_elec_tot, pdl_tot, mod_tot,&
upol, rmod)
!---------------------------------------------------------------------------------------------------
!
USE ions_base, ONLY: nat, tau, atm, zv, ityp
USE cell_base, ONLY: omega
USE noncollin_module, ONLY : noncolin, nspin_lsda
IMPLICIT NONE
!
TYPE (berryPhaseOutput_type) :: obj
REAL(DP),INTENT(IN) :: gpar(3), gvec, pdl_ion(nat), pdl_ion_tot, xk(3,*)
REAL(DP),INTENT(IN) :: pdl_elec(:), pdl_elec_up, pdl_elec_dw, pdl_elec_tot, &
pdl_tot, upol(3), rmod
!
INTEGER,INTENT(IN) :: mod_ion(nat), mod_ion_tot, mod_elec(:), mod_elec_up, &
mod_elec_dw, mod_elec_tot, mod_tot, nppstr, nkort, nstring
!
REAL(DP),INTENT(IN) :: wstring(nstring)
!
CHARACTER(LEN=*),PARAMETER :: TAGNAME = "BerryPhase"
TYPE ( polarization_type) :: tot_pol_obj
!
TYPE ( electronicPolarization_type),ALLOCATABLE :: str_pol_obj(:)
TYPE ( ionicPolarization_type ), ALLOCATABLE :: ion_pol_obj(:)
TYPE ( k_point_type ) :: kp_obj
TYPE ( phase_type) :: el_phase, ion_phase, tot_phase
TYPE ( atom_type ) :: atom_obj
TYPE ( scalarQuantity_type ) :: pol_val
INTEGER :: iat, istring, indstring, ispin
CHARACTER(10) :: mod_string
LOGICAL :: spin_is = .FALSE.
!
ALLOCATE (ion_pol_obj(nat), str_pol_obj(nat))
DO iat =1, nat
WRITE(mod_string,'("(mod" ,I1,")")') mod_ion(iat)
CALL qes_init_phase(ion_phase,"phase", 0.d0,.FALSE.,0.d0,.FALSE.,TRIM(mod_string),.TRUE., pdl_ion(iat) )
CALL qes_init_atom(atom_obj,"ion",name=TRIM(atm(ityp(iat))),position_ispresent=.FALSE.,atom = tau(:,iat))
CALL qes_init_ionicPolarization(ion_pol_obj(iat), "ionicPolarization", atom_obj, zv(ityp(iat)), ion_phase )
CALL qes_reset_phase(ion_phase)
CALL qes_reset_atom(atom_obj)
END DO
!
IF ( nspin_lsda .EQ. 2 ) spin_is = .TRUE.
DO istring= 1, nstring
indstring = 1+(istring-1)*nppstr
WRITE(mod_string,'("(mod ",I1,")")') mod_elec(istring)
CALL qes_init_phase(el_phase, "phase", 0.d0, .FALSE., 0.d0, .FALSE., TRIM (mod_string), .TRUE., &
pdl_elec(iat) )
IF (istring .LE. nstring/nspin_lsda) THEN
ispin = 1
ELSE
ispin = 2
END IF
CALL qes_init_k_point(kp_obj, "firstKeyPoint", wstring(istring), .TRUE., xk(:,indstring))
CALL qes_init_electronicPolarization(str_pol_obj(istring),"electronicPolarization", kp_obj, spin_is, ispin, &
el_phase )
CALL qes_reset_phase ( el_phase )
CALL qes_reset_k_point(kp_obj)
END DO
!
WRITE(mod_string,'("(mod ",I1,")")') mod_tot
CALL qes_init_phase(tot_phase, "totalPhase", pdl_ion_tot, .TRUE. , pdl_elec_tot, .TRUE., TRIM(mod_string), &
.TRUE., pdl_tot)
!
CALL qes_init_scalarQuantity ( pol_val, "polarization", Units="e/bohr^2", scalarQuantity=(rmod/omega)*pdl_tot )
!
CALL qes_init_polarization(tot_pol_obj, "polarization", pol_val, modulus = (rmod/omega)*dble(mod_tot), &
direction = upol )
!
CALL qes_init_berryPhaseOutput( obj, TAGNAME, tot_pol_obj, tot_phase, nat, ion_pol_obj, nstring, str_pol_obj )
!
DO istring=1,nstring
CALL qes_reset_electronicPolarization(str_pol_obj(istring))
END DO
DEALLOCATE (str_pol_obj)
DO iat=1, nat
CALL qes_reset_ionicPolarization(ion_pol_obj(iat))
END DO
DEALLOCATE (ion_pol_obj)
CALL qes_reset_polarization(tot_pol_obj)
CALL qes_reset_scalarQuantity(pol_val)
CALL qes_reset_phase(tot_phase)
!
END SUBROUTINE qexsd_init_berryPhaseOutput
!
!-------------------------------------------------------------------------------------------------
SUBROUTINE qexsd_set_status(status_int)
!-------------------------------------------------------------------------------------------------
IMPLICIT NONE
!
INTEGER :: status_int
CALL qes_init_status( exit_status, "status", status_int)
END SUBROUTINE qexsd_set_status
!
!--------------------------------------------------------------------------------------------------
SUBROUTINE qexsd_set_closed()
!
IMPLICIT NONE
CHARACTER(LEN=9) :: cdate, time_string
CHARACTER(LEN=12) :: date_string
!
CALL date_and_tim( cdate, time_string )
date_string = cdate(1:2) // ' ' // cdate(3:5) // ' ' // cdate (6:9)
CALL qes_init_closed (qexsd_closed_element, "closed", date_string, time_string,&
"")
END SUBROUTINE qexsd_set_closed
!-------------------------------------------------------------------------
!
!-------------------------------------------
! ... read subroutines
!-------------------------------------------
!
!
END MODULE qexsd_module
!
!----------------
! ... dummy defs
!----------------
!
#else
!
MODULE qexsd_module
IMPLICIT NONE
INTEGER :: dummy__
END MODULE qexsd_module
!
#endif