diff --git a/src/Makefile.am b/src/Makefile.am index 1dcd460a1ddfc0942ccf1eaa17008b09052b2886..6f93f803025c329af577b87c285b51ad5ad96639 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -551,7 +551,8 @@ noinst_HEADERS += \ output/output_me_inc.F90 \ output/output_modelmb_inc.F90 \ output/output_mxll_inc.F90 \ - output/output_states_inc.F90 + output/output_states_inc.F90 \ + output/output_tdm_inc.F90 output_srcs = $(output_f_srcs) diff --git a/src/basic/global.F90 b/src/basic/global.F90 index a9e0ffd45a48dc2d0935fe780b93c42b57f7a4f2..3b894157b37d764eae6f0fd844ac460774ff858c 100644 --- a/src/basic/global.F90 +++ b/src/basic/global.F90 @@ -51,7 +51,7 @@ module global_oct_m integer, public, parameter :: MAX_PATH_LEN=256 - integer, public, parameter :: MAX_OUTPUT_TYPES=41 + integer, public, parameter :: MAX_OUTPUT_TYPES=42 type conf_t ! Components are public by default diff --git a/src/electrons/electrons.F90 b/src/electrons/electrons.F90 index 3dc061447ca891179a832b006c4bfc107e955ac2..574f58d1f04ddcbae5fa874000a8e5008dbfe137 100644 --- a/src/electrons/electrons.F90 +++ b/src/electrons/electrons.F90 @@ -228,7 +228,7 @@ contains call mesh_check_symmetries(this%gr%mesh, this%gr%symm, this%ions%space%periodic_dim) end if - call output_init(this%outp, this%namespace, this%space, this%st, this%st%nst, this%ks) + call output_init(this%outp, this%namespace, this%space, this%st, this%st%nst, this%ks, this%gr, this%ions, this%kpoints) call states_elec_densities_init(this%st, this%gr) call states_elec_exec_init(this%st, this%namespace, this%mc) diff --git a/src/grid/io_function.F90 b/src/grid/io_function.F90 index 6d0fcc08d861935f47a5fc821f27e33b0e1e7e4b..a3c5d7e305c4e1a4d057c3a887bf819797554ef1 100644 --- a/src/grid/io_function.F90 +++ b/src/grid/io_function.F90 @@ -340,6 +340,9 @@ contains !% Outputs the exchange-correlation torque. Only for the spinor case and in the 3D case. !%Option eigenval_kpt 41 !% Outputs the eigenvalues resolved in momentum space, with one file for each band. + !%Option tdm 42 + !% Output the transition density matrix. + !% This is used to compute the TDTDM as defined in Williams et al., JCTC 17, 1795 (2021). !%End !%Variable OutputFormat diff --git a/src/output/output.F90 b/src/output/output.F90 index 848a732efaef4904e0d314979388e7f3199c7f14..887fffe3461a4c0d7a39ca7da198cc2d0ebc7c24 100644 --- a/src/output/output.F90 +++ b/src/output/output.F90 @@ -49,6 +49,7 @@ module output_oct_m use ions_oct_m use kick_oct_m use kpoints_oct_m + use lalg_basic_oct_m use lasers_oct_m use lda_u_oct_m use lda_u_io_oct_m @@ -78,6 +79,7 @@ module output_oct_m use submesh_oct_m use symm_op_oct_m use symmetries_oct_m + use symmetrizer_oct_m use unit_oct_m use unit_system_oct_m use utils_oct_m @@ -153,24 +155,33 @@ module output_oct_m type(output_bgw_t) :: bgw !< parameters for BerkeleyGW output + integer :: supercell(3) + integer :: ip_h + integer, allocatable :: map_COM(:) + contains procedure :: what_now + final :: output_final end type output_t contains - subroutine output_init(outp, namespace, space, st, nst, ks) + subroutine output_init(outp, namespace, space, st, nst, ks, gr, ions, kpoints) type(output_t), intent(out) :: outp type(namespace_t), intent(in) :: namespace type(space_t), intent(in) :: space type(states_elec_t), intent(in) :: st integer, intent(in) :: nst type(v_ks_t), intent(inout) :: ks + type(grid_t), intent(in) :: gr + type(ions_t), intent(in) :: ions + type(kpoints_t), intent(in) :: kpoints type(block_t) :: blk - FLOAT :: norm + FLOAT :: norm, xx_h(3), xx_rel(3), xx_e(3), xx_COM(3), dmin character(len=80) :: nst_string, default + integer :: idir, ncols, rankmin, ii PUSH_SUB(output_init) outp%what = .false. @@ -381,6 +392,75 @@ contains outp%output_interval(OPTION__OUTPUT__POTENTIAL) = outp%output_interval(OPTION__OUTPUT__POTENTIAL_GRADIENT) end if + if(outp%what(OPTION__OUTPUT__TDM)) then + outp%supercell = M_ONE + if(gr%mesh%parallel_in_domains) then + call messages_not_implemented("Output = tdm with domain parallelization.") + end if + + if(kpoints%use_symmetries) then + call messages_not_implemented("Output = tdm with k-point symmetries.") + end if + + ! SupercellDimensions defined in src/utils/tdtdm.F90 + if (parse_is_defined(namespace, 'SupercellDimensions')) then + if (parse_block(namespace, 'SupercellDimensions', blk) == 0) then + ncols = parse_block_cols(blk, 0) + if (ncols /= space%dim) then + write(message(1),'(a,i3,a,i3)') 'SupercellDimensions has ', ncols, ' columns but must have ', space%dim + call messages_fatal(1, namespace=namespace) + end if + do idir = 1, space%dim + call parse_block_integer(blk, 0, idir - 1, outp%supercell(idir)) + end do + + call parse_block_end(blk) + end if + else + outp%supercell(1:space%dim) = kpoints%nik_axis(1:space%dim) + end if + + ! TDTDMHoleCoordinates defined in src/utils/tdtdm.F90 + if(parse_block(global_namespace, 'TDTDMHoleCoordinates', blk) == 0) then + if(parse_block_cols(blk,0) < space%dim) then + call messages_input_error(global_namespace, 'TDTDMHoleCoordinates') + end if + do idir = 1, space%dim + call parse_block_float(blk, 0, idir - 1, xx_h(idir), units_inp%length) + end do + call parse_block_end(blk) + else + xx_h(1:space%dim) = ions%pos(1:space%dim, 1) + end if + ! We bring back the hole into the cell + xx_h = ions%latt%fold_into_cell(xx_h) + + !At the moment, we ignore rankmin + ASSERT(.not.gr%mesh%parallel_in_domains) + outp%ip_h = mesh_nearest_point(gr%mesh, xx_h, dmin, rankmin) + write(message(1), '(a, 3(1x,f7.4,a))') "Info: Requesting the hole at (", xx_h(1), & + ",", xx_h(2), ",", xx_h(3), ")." + call mesh_r(gr%mesh, outp%ip_h, dmin, coords=xx_h) + write(message(2), '(a, 3(1x,f7.4,a))') "Info: Setting the hole at (", xx_h(1), & + ",", xx_h(2), ",", xx_h(3), ")." + + SAFE_ALLOCATE(outp%map_COM(1:gr%mesh%np)) + ! At the moment, we ignore rankmin + ASSERT(.not.gr%mesh%parallel_in_domains) + ! We use outp%ip_h to fix the center of mass for the moment + call mesh_r(gr%mesh, outp%ip_h, dmin, coords=xx_COM) + ! We find the mapping between the center of mass and the electronic coordinate + do ii = 1, gr%mesh%np + call mesh_r(gr%mesh, ii, dmin, coords=xx_e) + xx_rel(1:space%dim) = M_TWO * xx_e(1:space%dim) - xx_COM(1:space%dim) + ! We bring back the hole into the cell + xx_rel = ions%latt%fold_into_cell(xx_rel) + outp%map_COM(ii) = mesh_nearest_point(gr%mesh, xx_rel, dmin, rankmin) + end do + + end if + + !%Variable OutputDuringSCF !%Type logical @@ -559,6 +639,8 @@ contains call output_xc_torque(outp, namespace, dir, gr%mesh, hm, st, ions, ions%space) + call output_tdm(outp, namespace, ions%space, dir, st, gr, ions, hm, iter) + call profiling_out(prof) POP_SUB(output_all) end subroutine output_all @@ -1459,6 +1541,17 @@ contains end function what_now + ! --------------------------------------------------------- + subroutine output_final(this) + type(output_t), intent(inout) :: this + + PUSH_SUB(output_final) + + SAFE_DEALLOCATE_A(this%map_COM) + + POP_SUB(output_final) + end subroutine output_final + #include "output_etsf_inc.F90" #include "output_states_inc.F90" @@ -1469,6 +1562,8 @@ contains #include "output_linear_medium_inc.F90" +#include "output_tdm_inc.F90" + #include "undef.F90" #include "complex.F90" #include "output_linear_response_inc.F90" diff --git a/src/output/output_tdm_inc.F90 b/src/output/output_tdm_inc.F90 new file mode 100644 index 0000000000000000000000000000000000000000..0ba1fb84287fff7adde646d36c92afb90bbc4c5e --- /dev/null +++ b/src/output/output_tdm_inc.F90 @@ -0,0 +1,223 @@ +!! Copyright (C) 2022 N. Tancogne-Dejean +!! +!! This program is free software; you can redistribute it and/or modify +!! it under the terms of the GNU General Public License as published by +!! the Free Software Foundation; either version 2, or (at your option) +!! any later version. +!! +!! This program is distributed in the hope that it will be useful, +!! but WITHOUT ANY WARRANTY; without even the implied warranty of +!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +!! GNU General Public License for more details. +!! +!! You should have received a copy of the GNU General Public License +!! along with this program; if not, write to the Free Software +!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +!! 02110-1301, USA. +!! + +! --------------------------------------------------------- +subroutine output_tdm(outp, namespace, space, dir, st, gr, ions, hm, iter) + type(output_t), intent(in) :: outp + type(namespace_t), intent(in) :: namespace + type(space_t), intent(in) :: space + character(len=*), intent(in) :: dir + type(states_elec_t), intent(in) :: st + type(grid_t), intent(in) :: gr + type(ions_t), intent(in) :: ions + type(hamiltonian_elec_t), intent(in) :: hm + integer, intent(in) :: iter + + integer :: ii, jj, kk, ip, irep, Nreplica, ik, ikpoint, ist, iop, idim, ierr + integer :: ip_h, irep_h, iunit, rankmin + integer, allocatable :: map_COM(:) + CMPLX, allocatable, target :: psi(:,:), psi_sym_COM(:,:) + CMPLX, allocatable :: tdm(:,:), tdm_COM(:,:), tdm_1D(:,:,:,:), phase(:,:) + CMPLX, pointer :: psi_sym(:,:) + FLOAT :: weight, kpoint(3), xx_COM(3), xx_rel(3), xx_e(3), xx_h(3), dmin + FLOAT, allocatable :: centers(:,:)!, den_1D(:,:,:,:) + type(unit_t) :: fn_unit + type(symmetrizer_t) :: symmetrizer + + if (.not. outp%what_now(OPTION__OUTPUT__TDM, iter)) return + + PUSH_SUB(output_tdm) + + Nreplica = product(outp%supercell(1:space%dim)) + + ! The center of each replica of the unit cell + SAFE_ALLOCATE(centers(1:3, 1:Nreplica)) + irep = 1 + do ii = 0, outp%supercell(1)-1 + do jj = 0, outp%supercell(2)-1 + do kk = 0, outp%supercell(3)-1 + centers(1, irep) = -floor((outp%supercell(1)-1)/M_TWO)+ii + centers(2, irep) = -floor((outp%supercell(2)-1)/M_TWO)+jj + centers(3, irep) = -floor((outp%supercell(3)-1)/M_TWO)+kk + centers(1:space%dim, irep) = matmul(ions%latt%rlattice, centers(1:space%dim, irep)) + irep = irep + 1 + end do + end do + end do + + ! The phase for each center + SAFE_ALLOCATE(phase(st%d%kpt%start:st%d%kpt%end, 1:Nreplica)) + + ! Not working with k-points symmetries yet + ASSERT(.not. hm%kpoints%use_symmetries) + + do irep = 1, Nreplica + do ik = st%d%kpt%start, st%d%kpt%end + ikpoint = st%d%get_kpoint_index(ik) + kpoint(1:space%dim) = hm%kpoints%get_point(ikpoint) + phase(ik, irep) = exp(-M_zI*sum(kpoint(1:space%dim)*centers(1:space%dim, irep))) + end do + end do + + + SAFE_ALLOCATE(psi(1:gr%mesh%np, 1:st%d%dim)) + + ! Arrays needed for the center of mass frame + SAFE_ALLOCATE(psi_sym_COM(1:gr%mesh%np, 1:st%d%dim)) + + + if(hm%kpoints%use_symmetries) then + call symmetrizer_init(symmetrizer, gr%mesh, gr%symm) + SAFE_ALLOCATE(psi_sym(1:gr%mesh%np, 1:st%d%dim)) + end if + + ! TDM with a fixed hole + SAFE_ALLOCATE(tdm(1:gr%mesh%np, 1:Nreplica)) + tdm = M_z0 + + ! TDM in the center-of-mass frame (r_e+r_h) fixed + SAFE_ALLOCATE(tdm_COM(1:gr%mesh%np, 1:Nreplica)) + tdm_COM = M_z0 + + ! if(space%dim == 1) then + ! SAFE _ALLOCATE(tdm_1D(1:gr%mesh%np, 1:gr%mesh%np, 1:Nreplica, 1:Nreplica)) + ! tdm_1D = M_ZERO + ! end if + + do ik = st%d%kpt%start, st%d%kpt%end + ikpoint = st%d%get_kpoint_index(ik) + + do ist = 1, st%nst + weight = st%d%kweights(ik) * st%occ(ist, ik) / kpoints_get_num_symmetry_ops(hm%kpoints, ikpoint) + if(abs(weight) < M_EPSILON) cycle + + call states_elec_get_state(st, gr%mesh, ist, ik, psi) + if(allocated(hm%hm_base%phase)) then + call states_elec_set_phase(st%d, psi, hm%hm_base%phase(1:gr%mesh%np, ik), gr%mesh%np, .false.) + end if + + do ii = 1, kpoints_get_num_symmetry_ops(hm%kpoints, ikpoint) + iop = kpoints_get_symmetry_ops(hm%kpoints, ikpoint, ii) + + if(hm%kpoints%use_symmetries) then + do idim = 1, st%d%dim + call zsymmetrizer_apply_single(symmetrizer, gr%mesh, iop, psi(:,idim), psi_sym(:,idim)) + end do + else + psi_sym => psi + end if + + ! Change from \psi(r_e) to \psi(r_e+r_h) + do ip = 1, gr%mesh%np + psi_sym_COM(ip, 1) = psi_sym(outp%map_COM(ip),1) + end do + + ! We now compute the TDM + do irep = 1, Nreplica + call lalg_axpy(gr%mesh%np, phase(ik, irep) * weight & + * conjg(psi_sym(outp%ip_h,1)), psi_sym(:, 1), tdm(:,irep)) + ! We use outp%ip_h to fix the center of mass for the moment + call lalg_axpy(gr%mesh%np, phase(ik, irep) * weight & + * conjg(psi_sym(outp%ip_h,1)), psi_sym_COM(:, 1), tdm_COM(:,irep)) + end do + + ! if(space%dim == 1) then ! In the 1D case, we contruct the full TDTDM of r_e, r_h + ! do irep_h = 1, Nreplica + ! do irep = 1, Nreplica + ! do ip_h = 1, gr%mesh%np + ! call lalg_axpy(gr%mesh%np, phase(ik, irep) * conjg(phase(ik, irep_h)) & + ! * weight * conjg(psi_sym(ip_h,1)), psi_sym(:, 1), tdm_1D(:, ip_h, irep, irep_h)) + ! end do + ! end do + ! end do + ! end if + end do ! ii + end do + end do + + if(st%d%kpt%parallel) then + call comm_allreduce(st%d%kpt%mpi_grp, tdm) + call comm_allreduce(st%d%kpt%mpi_grp, tdm_COM) +! call comm_allreduce(st%d%kpt%mpi_grp, tdm_1D) + end if + + fn_unit = units_out%length**(-space%dim) + call io_function_output_supercell(outp%how(OPTION__OUTPUT__TDM), dir, 'tdm', gr%mesh, space, & + tdm, centers, outp%supercell, fn_unit, ierr, namespace, ions = ions, grp = st%dom_st_kpt_mpi_grp) + + call io_function_output_supercell(outp%how(OPTION__OUTPUT__TDM), dir, 'tdm_COM', gr%mesh, space, & + tdm_COM, centers, outp%supercell, fn_unit, ierr, namespace, ions = ions, grp = st%dom_st_kpt_mpi_grp) + +! if(space%dim == 1) then +! A SSERT(.not.gr%mesh%parallel_in_domains) +! +! SAFE _ALLOCATE(den_1D(1:gr%mesh%np, 1:gr%mesh%np, 1:Nreplica, 1:Nreplica)) +! do irep_h = 1, Nreplica +! do irep = 1, Nreplica +! do ip_h = 1, gr%mesh%np +! do ii = 1, gr%mesh%np +! den_1D(ii, ip_h, irep, irep_h) = TOFLOAT(tdm_1D(ii, ip_h, irep, irep_h)*conjg(tdm_1D(ii,ip_h, irep, irep_h))) +! end do +! end do +! end do +! end do + +! if (mpi_grp_is_root(mpi_world)) then +! iunit = io_open(trim(dir)//'/tdm_density', action='write') +! write(iunit, '(a)', iostat=ierr) '# r_e r_h Re(\Psi(r_e,r_h)) Im(\Psi(r_e,r_h)) |\Psi(r_e,r_h)|^2' +! +! do irep_h = 1, Nreplica +! do ip_h = 1, gr%mesh%np +! xx_h = units_from_atomic(units_out%length, mesh_x_global(gr%mesh, ip_h) + centers(1:space%dim, irep_h)) +! +! do irep = 1, Nreplica +! do ii = 1, gr%mesh%np +! xx_e = units_from_atomic(units_out%length, mesh_x_global(gr%mesh, ii) + centers(1:space%dim, irep)) +! write(iunit, '(5es23.14E3)', iostat=ierr) xx_e(1), xx_h(1), & +! TO FLOAT(units_from_atomic(fn_unit, tdm_1D(ii, ip_h, irep, irep_h))) ,& +! aimag(units_from_atomic(fn_unit, tdm_1D(ii, ip_h, irep, irep_h))), & +! units_from_atomic(fn_unit, den_1D(ii, ip_h, irep, irep_h)) +! end do +! end do +! end do +! end do +! close(iunit) +! end if +! +! SAFE _DEALLOCATE_A(den_1D) +! end if + + SAFE_DEALLOCATE_A(tdm) + SAFE_DEALLOCATE_A(tdm_COM) +! SAFE _DEALLOCATE_A(tdm_1D) + SAFE_DEALLOCATE_A(phase) + if(hm%kpoints%use_symmetries) then + SAFE_DEALLOCATE_P(psi_sym) + call symmetrizer_end(symmetrizer) + end if + SAFE_DEALLOCATE_A(psi) + SAFE_DEALLOCATE_A(psi_sym_COM) + + POP_SUB(output_tdm) + +end subroutine output_tdm + +!! Local Variables: +!! mode: f90 +!! coding: utf-8 +!! End: diff --git a/src/utils/tdtdm.F90 b/src/utils/tdtdm.F90 index b4357c02f59a70ddea7b5f212ad534d43e57f2ac..4134e12ab2205b6dd1d3558de8be2ceca76576eb 100644 --- a/src/utils/tdtdm.F90 +++ b/src/utils/tdtdm.F90 @@ -489,9 +489,9 @@ program tdtdm case(2,3) do irep = 1, Nreplica call lalg_axpy(sys%gr%mesh%np, phase(ik, irep) * weight & - * conjg(Xiak(ist,uist,ik))*conjg(psi(ip_h,1)), upsi(:, 1), tdm(:,irep)) + * conjg(Xiak(ist,uist,ik))*conjg(psi_sym(ip_h,1)), upsi_sym(:, 1), tdm(:,irep)) call lalg_axpy(sys%gr%mesh%np, phase(ik, irep) * weight & - * conjg(Yiak(ist,uist,ik))*conjg(upsi(ip_h,1)), psi(:, 1), tdm(:,irep)) + * conjg(Yiak(ist,uist,ik))*conjg(upsi_sym(ip_h,1)), psi_sym(:, 1), tdm(:,irep)) end do case(1) ! In the 1D case, we contruct the full TDTDM of r_e, r_h @@ -500,10 +500,10 @@ program tdtdm do ip_h = 1, sys%gr%mesh%np call lalg_axpy(sys%gr%mesh%np, phase(ik, irep) * conjg(phase(ik, irep_h)) & * weight * conjg(Xiak(ist,uist,ik)) * conjg(psi_sym(ip_h,1)), & - upsi(:, 1), tdm_1D(:, ip_h, irep, irep_h)) + upsi_sym(:, 1), tdm_1D(:, ip_h, irep, irep_h)) call lalg_axpy(sys%gr%mesh%np, phase(ik, irep) * conjg(phase(ik, irep_h)) & * weight * conjg(Yiak(ist,uist,ik)) * conjg(upsi_sym(ip_h,1)), & - psi(:, 1), tdm_1D(:, ip_h, irep, irep_h)) + psi_sym(:, 1), tdm_1D(:, ip_h, irep, irep_h)) end do end do end do