From 2f7c2b01a74b2f9f32cda40ead6c7a7d1ea7307a Mon Sep 17 00:00:00 2001 From: chpxu Date: Fri, 28 Aug 2026 14:40:02 +0100 Subject: [PATCH 1/3] Implement rotation code The previous branch f-rotate-pDOS became broken. This commit readds the necessary code. --- tools/PostProcessing/local_module.f90 | 20 +- tools/PostProcessing/process_module.f90 | 1226 ++++++++++++++++++++++- tools/PostProcessing/read_module.f90 | 105 +- 3 files changed, 1314 insertions(+), 37 deletions(-) diff --git a/tools/PostProcessing/local_module.f90 b/tools/PostProcessing/local_module.f90 index 7ec43a61..bb6e8dfa 100644 --- a/tools/PostProcessing/local_module.f90 +++ b/tools/PostProcessing/local_module.f90 @@ -3,9 +3,9 @@ module local use datatypes ! These give the number of blocks in x, y and z - integer :: nblockx, nblocky, nblockz + integer :: nblockx, nblocky, nblockz real(double) :: block_size_x, block_size_y, block_size_z, grid_x, grid_y, grid_z - + ! Processes used integer :: nprocs @@ -17,14 +17,26 @@ module local integer, dimension(:,:,:), allocatable :: band_active_kp integer, dimension(:,:), allocatable :: band_active_all real(double), dimension(:,:,:), allocatable :: current + ! pDOS rotation + logical :: flag_rotate_pdos, flag_rotate_pdos_debug + integer :: flag_rotate_pdos_mode, rotate_pdos_natoms + character(len=5) :: flag_rotate_pdos_units + real(double), dimension(:,:,:), allocatable :: U1, U2 ! Rotation matrices + real(double) :: pdos_ax(3), pdos_ay(3), pdos_az(3) ! axes to rotate pDOS into + real(double), dimension(:,:), allocatable :: euler_angles ! Euler angles to rotate pDOS by + real(double), dimension(:,:), allocatable :: axes_angles ! Axis angles to rotate pDOS by + integer, dimension(:), allocatable :: rotate_pdos_atoms ! Array of atoms to rotate (mode 1 or 3) + integer, dimension(:,:), allocatable :: find_neighbours ! mode 2 + integer, dimension(:), allocatable :: nghbr_arr ! holds nearest neighbour atoms + ! Store eigenvector coefficients complex(double_cplx), allocatable, dimension(:,:,:,:,:), save :: evec_coeff ! PAOs, atoms, bands, kpoints, spin complex(double_cplx), dimension(:,:,:,:,:), allocatable :: scaled_evec_coeff - + character(len=50) :: root_file - + real(double) :: stm_bias, fermi_offset, stm_z_min, stm_z_max, stm_x_min, stm_x_max, & stm_y_min, stm_y_max, stm_broad, gpv, E_wf_min, E_wf_max, E_procwf_min, E_procwf_max integer :: nptsx, nptsy, nptsz, nxmin, nymin, nzmin diff --git a/tools/PostProcessing/process_module.f90 b/tools/PostProcessing/process_module.f90 index 779951c3..b18473fc 100644 --- a/tools/PostProcessing/process_module.f90 +++ b/tools/PostProcessing/process_module.f90 @@ -6,7 +6,7 @@ module process ! Maximum possible number of spin components for simplicity real(double), dimension(4) :: range_offset - + contains subroutine assign_blocks @@ -14,7 +14,7 @@ subroutine assign_blocks use datatypes use local, ONLY: block_store, nprocs, block_size_x, block_size_y, block_size_z, & stm_z_min, stm_z_max, stm_x_min, stm_x_max, stm_y_min, stm_y_max - + implicit none integer :: proc, iblock, ig1, ind_group, block_x, block_y, block_z, nblock @@ -36,7 +36,7 @@ subroutine assign_blocks ! Now LHS of block and RHS of area rbx = block_size_x*real(block_store(proc)%nx(iblock)-1,double) rby = block_size_y*real(block_store(proc)%ny(iblock)-1,double) - rbz = block_size_z*real(block_store(proc)%nz(iblock)-1,double) + rbz = block_size_z*real(block_store(proc)%nz(iblock)-1,double) if(rbx<=stm_x_max.AND.rby<=stm_y_max.AND.rbz<=stm_z_max) & block_store(proc)%active(iblock) = block_store(proc)%active(iblock) + 1 !write(*,*) 'RHS: ',rbx,rby,rbz @@ -80,7 +80,7 @@ subroutine process_charge if(ispin==1) ci = TRIM(charge_stub)//"_up" if(ispin==2) ci = TRIM(charge_stub)//"_dn" end if - do proc = 1, nprocs + do proc = 1, nprocs call get_file_name(ci,nprocs,proc,filename) ! Open file open(unit=17,file=filename) @@ -303,7 +303,7 @@ subroutine process_dos use units, ONLY: HaToeV implicit none - + ! Local variables integer :: i_band, i_kp, i_spin, n_DOS_wid, n_band, n_min, n_max, i real(double) :: Ebin, dE_DOS, a, pf_DOS, spin_fac, peak_width @@ -442,10 +442,15 @@ end subroutine process_dos subroutine process_pdos use datatypes - use numbers, ONLY: zero, RD_ERR, twopi, half, one, two, four, six + use numbers, ONLY: zero, RD_ERR, pi, twopi, half, one, two, four, six use local, ONLY: eigenvalues, n_bands_total, nkp, wtk, efermi, flag_total_iDOS, & evec_coeff, scaled_evec_coeff, flag_procwf_range_Ef, flag_l_resolved, flag_lm_resolved, & band_full_to_active, n_atoms_pDOS, pDOS_atom_index + ! pDOS rotation variables + use local, ONLY: flag_rotate_pdos, flag_rotate_pdos_mode, flag_rotate_pdos_debug, & + pdos_ax, pdos_ay, pdos_az, rotate_pdos_atoms, axes_angles, euler_angles, find_neighbours,& + rotate_pdos_natoms, nghbr_arr, U1, U2 + use read, ONLY: read_eigenvalues, read_psi_coeffs, read_nprocs_from_blocks use global_module, ONLY: nspin, n_DOS, E_DOS_min, E_DOS_max, sigma_DOS, ni_in_cell, species_glob use units, ONLY: HaToeV @@ -453,9 +458,9 @@ subroutine process_pdos use pao_format, ONLY: pao implicit none - + ! Local variables - integer :: i_band, i_kp, i_spin, n_DOS_wid, n_band, n_min, n_max, i, i_atom,max_nsf, i_spec, & + integer :: i_band, i_kp, i_spin, n_DOS_wid, n_band, n_min, n_max, i, j, i_atom,max_nsf, i_spec, & i_l, nzeta, sf_offset, max_l, norbs, i_m, i_band_c, i_z real(double) :: Ebin, dE_DOS, a, pf_DOS, spin_fac, coeff, check_electrons, peak_width real(double), dimension(:,:,:), allocatable :: pDOS @@ -464,6 +469,11 @@ subroutine process_pdos real(double), dimension(:,:), allocatable :: occ real(double), dimension(:,:), allocatable :: total_electrons real(double), dimension(:,:,:), allocatable :: total_electrons_l + ! pDOS rotation variables + real(double) :: A1(3, 3), A2(5, 5), C1(3, 3), C2(5, 5), E1(3,3), E2(5,5), Qxyz(3,3) + real(double) :: rod(9), angle, axis(3), temp_matrix(3,3) + real(double), dimension(:,:), allocatable :: bond(:,:) + character(len=25) :: filename,fmt_dos complex(double_cplx),external :: zdotc @@ -500,6 +510,192 @@ subroutine process_pdos ! Read eigenvector coefficients call read_psi_coeffs("Process") allocate(occ(n_bands_total,nkp)) + + ! Call pDOS rotation subroutines if desired. Also includes debug output + if (flag_rotate_pdos) then + write(*,fmt='(2x,"Rotating wavefunction coefficients")') + call initialise_A_mat(A1, A2) + if (flag_rotate_pdos_debug .and. flag_rotate_pdos_mode /= 1) then + write(*, fmt='(/2x,"ROTATION DEBUG OUTPUT: MODE ", (I0,1x))') & + flag_rotate_pdos_mode + + write(*, fmt='(/2x,"ROTATION DEBUG OUTPUT: A^(l) MATRICES")') + write(*, fmt='(/4x, "A1: ")') + do j = 1, 3 + write(*, fmt='(/4x,3(f10.5,1X))') A1(j, :) + end do + write(*, fmt='(/4x, "A2: ")') + do j = 1, 5 + write(*, fmt='(/4x,5(f10.5,1X))') A2(j, :) + end do + write(*, fmt='(/4x, "For l = 1, orbital basis is [|y>, |z>, |x>]")') + write(*, fmt='(/4x, "For l = 2, orbital basis is [|-xy>, |yz>, |3z^2-r^2>, |xz>, |x^2-y^2>]")') + end if + if (allocated(U1)) deallocate(U1) + if (allocated(U2)) deallocate(U2) + + if (flag_rotate_pdos_mode == 0) then + write(*,fmt='(2x,"Using user input axes")') + allocate(U1(3,3,n_atoms_pDOS)) + allocate(U2(5,5,n_atoms_pDOS)) + call get_pdos_axes + + if(n_atoms_pDOS==ni_in_cell) then ! All atoms + write(*, fmt='(/2x,"New local axes for all atoms in unit cell:")') + else + write(*, fmt='(/2x,"New local axes for specified atoms: ", *(I0,1x))') & + pDOS_atom_index + end if + write(*, fmt='(/4x,"x: ",3(f10.5,1X))', advance="no") pdos_ax + write(*, fmt='(/4x,"y: ",3(f10.5,1X))', advance="no") pdos_ay + write(*, fmt='(/4x,"z: ",3(f10.5,1X))', advance="no") pdos_az + call calculate_axis_angle(pdos_ax, pdos_ay, pdos_az, axis, angle) + call construct_rodrigues(axis, angle, rod) + ! Construct all C^l matrices from rodrigues + call construct_C1(rod, C1) + call construct_C2(rod, C2) + do i = 1, n_atoms_pDOS + call construct_Ul(1, A1, C1, U1(:,:,i)) + call construct_Ul(2, A2, C2, U2(:,:,i)) + U1(:,:,i) = transpose(U1(:,:,i)) + U2(:,:,i) = transpose(U2(:,:,i)) + + if (flag_rotate_pdos_debug) then + ! call euler_from_axisangle(axis, angle) + write(*, fmt='(/4x, "C1: ")') + do j = 1, 3 + write(*, fmt='(/6x,3(f10.5,1X))', advance='no') C1(j,:) + end do + write(*, fmt='(/4x, "C2: ")') + do j = 1, 5 + write(*, fmt='(/6x,5(f10.5,1X))', advance='no') C2(j, :) + end do + call print_orbital_weights(i) + end if ! end rotation debug output + end do ! end mode 0, axes input + else if (flag_rotate_pdos_mode == 1) then + write(*,fmt='(2x,"Using extrinsic Euler angles in active zyz convention")') + allocate(U1(3,3,rotate_pdos_natoms)) + allocate(U2(5,5,rotate_pdos_natoms)) + Qxyz = 0.0 + Qxyz(1,2) = 1.0 + Qxyz(2,3) = 1.0 + Qxyz(3,1) = 1.0 + do i = 1, rotate_pdos_natoms + rod = 0.0 + write(*, fmt='(/2x,"New local axes for specified atoms: ", *(I0,1x))') & + rotate_pdos_atoms(i) + call construct_EulerMatrices(E1, E2, i) + U1(:,:,i) = (E1) + U2(:,:,i) = (E2) + + ! We only explicitly construct wavefunction rotation matrices + ! However can extrapolate new local axes using E1 since orbital basis is {y,z,x} + write(*, fmt='(/4x,"Image of x,y,z under this active rotation, in standard basis x,y,z")', advance="no") + write(*, fmt='(/4x,"x: ",3(f10.5,1X))', advance="no") & + E1(3,3), E1(1,3), E1(2,3) + write(*, fmt='(/4x,"y: ",3(f10.5,1X))', advance="no") & + E1(3,1), E1(1,1), E1(2,1) + write(*, fmt='(/4x,"z: ",3(f10.5,1X))', advance="no") & + E1(3,2), E1(1,2), E1(2,2) + temp_matrix = matmul(inv(Qxyz), matmul(E1, Qxyz)) + temp_matrix = inv(temp_matrix) + write(*, fmt='(/4x,"Equivalent basis transformation (inverse of active coordinates)")', advance="no") + write(*, fmt='(/4x,"x: ",3(f10.5,1X))', advance="no") & + temp_matrix(:,1) + write(*, fmt='(/4x,"y: ",3(f10.5,1X))', advance="no") & + temp_matrix(:,2) + write(*, fmt='(/4x,"z: ",3(f10.5,1X))', advance="no") & + temp_matrix(:,3) + if (flag_rotate_pdos_debug) then + write(*, fmt='(/4x, "alpha(z) beta(y) gamma(z): ",3(f10.5,1X))') & + euler_angles(:,i) + call print_orbital_weights(i) + end if ! end rotation debug mode + end do ! end mode 1, euler angles + else if (flag_rotate_pdos_mode == 2) then + write(*,fmt='(2x,"Using user input atom numbers and local geometry")') + allocate(U1(3,3,rotate_pdos_natoms)) + allocate(U2(5,5,rotate_pdos_natoms)) + + do i = 1, rotate_pdos_natoms + rod = 0.0 + call nearest_neighbours(find_neighbours(1, i), bond) + write(*, fmt='(/2x,"Located neighours of atom ", I0, ": ", *(I0,1x))') & + find_neighbours(1,i), nghbr_arr + call axes_from_nn(find_neighbours(1, i), bond) + write(*, fmt='(/2x,"New local axes for atom ", I0, ": ")') & + find_neighbours(1,i) + write(*, fmt='(/4x,"x: ",3(f10.5,1X))', advance="no") pdos_ax + write(*, fmt='(/4x,"y: ",3(f10.5,1X))', advance="no") pdos_ay + write(*, fmt='(/4x,"z: ",3(f10.5,1X))', advance="no") pdos_az + call get_pdos_axes + call calculate_axis_angle(pdos_ax, pdos_ay, pdos_az, axis, angle) + call construct_rodrigues(axis, angle, rod) + ! Construct all C^l matrices from rodrigues + call construct_C1(rod, C1) + call construct_C2(rod, C2) + call construct_Ul(1, A1, C1, U1(:,:,i)) + call construct_Ul(2, A2, C2, U2(:,:,i)) + U1(:,:,i) = transpose(U1(:,:,i)) + U2(:,:,i) = transpose(U2(:,:,i)) + if (flag_rotate_pdos_debug) then + ! call euler_from_axisangle(axis, angle) + write(*, fmt='(/4x, "C1: ")') + do j = 1, 3 + write(*, fmt='(/6x,3(f10.5,1X))', advance='no') C1(j,:) + end do + write(*, fmt='(/4x, "C2: ")') + do j = 1, 5 + write(*, fmt='(/6x,5(f10.5,1X))', advance='no') C2(j, :) + end do + call print_orbital_weights(i) + end if ! end rotation debug output + end do ! end mode 2, neighbour input + else if (flag_rotate_pdos_mode == 3) then + write(*,fmt='(2x,"Using user input axis and angle")') + allocate(U1(3,3,rotate_pdos_natoms)) + allocate(U2(5,5,rotate_pdos_natoms)) + + do i = 1, rotate_pdos_natoms + rod = 0.0 + axis = axes_angles(1:3,i) + angle = axes_angles(4,i) + write(*, fmt='(/2x,"User input axis: ", 3(f10.5))') axis + write(*, fmt='(/2x,"User input angle (rad, deg) ", 2(f10.5))') angle, angle * (180.0 / pi) + call construct_rodrigues(axis, angle, rod) + write(*, fmt='(/2x,"New local axes for atom ", I0, ": ")') & + rotate_pdos_atoms(i) + + ! Read off columns for local axes + write(*, fmt='(/4x,"x: ",3(f10.5,1X))', advance="no") rod(1), rod(4), rod(7) + write(*, fmt='(/4x,"y: ",3(f10.5,1X))', advance="no") rod(2), rod(5), rod(8) + write(*, fmt='(/4x,"z: ",3(f10.5,1X))', advance="no") rod(3), rod(6), rod(9) + ! Construct all C^l matrices from rodrigues + call construct_C1(rod, C1) + call construct_C2(rod, C2) + call construct_Ul(1, A1, C1, U1(:,:,i)) + call construct_Ul(2, A2, C2, U2(:,:,i)) + U1(:,:,i) = transpose(U1(:,:,i)) + U2(:,:,i) = transpose(U2(:,:,i)) + if (flag_rotate_pdos_debug) then + call euler_from_axisangle(axis, angle) + write(*, fmt='(/4x, "C1: ")') + do j = 1, 3 + write(*, fmt='(/6x,3(f10.5,1X))', advance='no') C1(j,:) + end do + write(*, fmt='(/4x, "C2: ")') + do j = 1, 5 + write(*, fmt='(/6x,5(f10.5,1X))', advance='no') C2(j, :) + end do + call print_orbital_weights(i) + end if ! end rotation debug output + end do ! end mode 3, direct axis angle + end if ! pdos rotation mode + call rotate_coefficients + deallocate(U1) + deallocate(U2) + end if ! rotate pdos ! Set up storage based on pDOS per atom, or l/lm resolved per atom if(flag_lm_resolved) then allocate(pDOS_lm(-max_l:max_l,0:max_l,n_atoms_pDOS,n_DOS,nspin)) @@ -723,6 +919,1000 @@ subroutine process_pdos return end subroutine process_pdos + ! pDOS rotation subroutines + ! misc subroutine for printing out orbital rotation matrices + subroutine print_orbital_weights(atom_index) + use local, ONLY: U1, U2 + implicit none + + integer, intent(in) :: atom_index + + integer :: i,j,k + real(double) :: pdos_weight, identity3(3,3), identity5(5,5) + character(len=30), parameter :: p_orb(3) = (/ "|y>", "|z>","|x>"/) + character(len=30), parameter :: d_orb(5) = (/ "|-xy> ", "|yz> ", & + "|3z^2-r^2>", "|xz> ", "|x^2-y^2> "/) + character(len=50) :: pdos_weight_str + character(len=100) :: line + + write(*, fmt='(/4x, "U1: ")') + do j = 1, 3 + write(*, fmt='(/6x,3(f10.5,1X))', advance='no') U1(j,:,atom_index) + end do + write(*, fmt='(/4x, "p-orbital weight decomposition")') + write(*, '(/6x, A10, 5A10)') "", (trim(p_orb(j)), j = 1, 3) + do k = 1, 3 + write(*, '(6x, A10, 5F10.5)') trim(p_orb(k)), & + (U1(k,j,atom_index)*U1(k,j,atom_index), j = 1, 3) + end do ! end printing l = 2 orbital weights + write(*, fmt='(/4x, "U2: ")') + ! Manual write because first row and column need to be negated + write(*, fmt='(/6x,5(f10.5,1X))', advance='no') U2(1,1,atom_index), -U2(1,2:,atom_index) + write(*, fmt='(/6x,5(f10.5,1X))', advance='no') -U2(2,1,atom_index), U2(2,2:,atom_index) + write(*, fmt='(/6x,5(f10.5,1X))', advance='no') -U2(3,1,atom_index), U2(3,2:,atom_index) + write(*, fmt='(/6x,5(f10.5,1X))', advance='no') -U2(4,1,atom_index), U2(4,2:,atom_index) + write(*, fmt='(/6x,5(f10.5,1X))', advance='no') -U2(5,1,atom_index), U2(5,2:,atom_index) + + write(*, fmt='(/4x, "d-orbital weight decomposition")') + + write(*, '(/6x, A10, 5A10)') "", (trim(d_orb(j)), j = 1, 5) + do k = 1, 5 + write(*, '(6x, A10, 5F10.5)') trim(d_orb(k)), & + (U2(k,j,atom_index)*U2(k,j,atom_index), j = 1, 5) + end do ! end printing l = 2 orbital weights + + ! Orthogonality check + identity3 = 0.0 + do i = 1, size(identity3(:,1)) + identity3(i,i) = 1.0 + end do + identity5 = 0.0 + do i = 1, size(identity5(:,1)) + identity5(i,i) = 1.0 + end do + if(all(abs(matmul(U1(:,:,atom_index), transpose(U1(:,:,atom_index))) - identity3) < 1e-7)) then + write(*, '(/2x, A)', advance='no') 'U1 matrix is orthogonal: PASS' + else + write(*, '(/2x, A)', advance='no') 'U1 matrix is orthogonal: FAIL' + endif + + if(all(abs(matmul(U2(:,:,atom_index), transpose(U2(:,:,atom_index))) - identity5) < 1e-7)) then + write(*, '(/2x, A)', advance='no') 'U2 matrix is orthogonal: PASS' + else + write(*, '(/2x, A)', advance='no') 'U2 matrix is orthogonal: FAIL' + endif + write(*,*) + end subroutine + + ! ----------------------------------------------------------------------------- + ! Subroutine construct_EulerMatrices + ! ----------------------------------------------------------------------------- + + !!****f* ProcModule/construct_EulerMatrices * + !! + !! NAME + !! construct_EulerMatrices + !! USAGE + !! construct_EulerMatrices(E1, E2, atom_index) + !! PURPOSE + !! Create rotation matrices using Euler angle input: right-handed, + !! active Euler angles about extrinsic (fixed) cell axes. Positive angle + !! is anticlockwise. + !! INPUTS + !! integer, intent(in) :: atom_index - atom to perform rotation for + !! USES + !! datatypes, numbers, local, GenComms + !! AUTHOR + !! C. Xu + !! CREATION DATE + !! 10/04/2026 + !! MODIFICATION HISTORY + !! 13/04/2026 - C. Xu: Add Euler matrix for d-orbitals + !! 20/07/2026 - C. Xu: Add output matrices in debugging mode + !! 30/07/2026 - C. Xu: Correct E2(3,2): alpha -> beta + !! 26/08/2026 - C. Xu: Fix orbital convention + !! SOURCE + !! Credits to R. Johnson for the mathematical expressions + !! Also see: Quantum Theory of Angular Momentum + subroutine construct_EulerMatrices(E1, E2, atom_index) + use datatypes + use numbers, ONLY: pi, two + use local, ONLY: euler_angles, flag_rotate_pdos_debug + use GenComms, ONLY: cq_abort + + implicit none + integer, intent(in) :: atom_index + real(double), intent(out) :: E1(3,3), E2(5,5) + real(double) :: ca, sa, cb, sb, cg, sg, c2a, s2a, c2b, s2b, c2g, s2g + real(double) :: euler_alpha, euler_beta, euler_gamma + real(double) :: identity3(3,3), identity5(5,5) + integer :: i + euler_alpha = euler_angles(1,atom_index) + euler_beta = euler_angles(2,atom_index) + euler_gamma = euler_angles(3,atom_index) + E1 = 0.0 + E2 = 0.0 + + ca = cos(euler_alpha) + sa = sin(euler_alpha) + cb = cos(euler_beta) + sb = sin(euler_beta) + cg = cos(euler_gamma) + sg = sin(euler_gamma) + + c2a = cos(two*euler_alpha) + s2a = sin(two*euler_alpha) + c2b = cos(two*euler_beta) + s2b = sin(two*euler_beta) + c2g = cos(two*euler_gamma) + s2g = sin(two*euler_gamma) + + ! Rotation of l = 1 coefficients + E1(1,1) = ca*cg - cb*sa*sg + E1(1,2) = sa*sb + E1(1,3) = cb*cg*sa + ca*sg + E1(2,1) = sb*sg + E1(2,2) = cb + E1(2,3) = -cg*sb + E1(3,1) = -cg*sa - ca*cb*sg + E1(3,2) = ca*sb + E1(3,3) = ca*cb*cg - sa*sg + + + ! Rotation of l = 2 coefficients + E2(1,1) = c2a*cb*c2g - sa*ca*(c2b+3.0)*sg*cg + E2(1,2) = sb*(c2a*cg - 2.0*sa*ca*cb*sg) + E2(1,3) = sqrt(3.0)*sa*ca*sb*sb + E2(1,4) = sb*(s2a*cb*cg + c2a*sg) + E2(1,5) = cb*c2a*s2g +0.25*s2a*(3.0+c2b)*c2g + + E2(2,1) = sb*(sa*cb*s2g - ca*c2g) + E2(2,2) = ca*cb*cg - sa*c2b*sg + E2(2,3) = sqrt(3.0)*sa*sb*cb + E2(2,4) = c2b*sa*cg + ca*cb*sg + E2(2,5) = -sb*(sa*cb*c2g +s2g*ca) + + E2(3,1) = -sqrt(3.0)*sb*sb*sg*cg + E2(3,2) = sqrt(3.0)*sb*cb*sg + E2(3,3) = 0.25*(3.0*c2b+1.0) + E2(3,4) = -sqrt(3.0)*sb*cb*cg + E2(3,5) = 0.5*sqrt(3.0)*sb*sb*c2g + + E2(4,1) = sb*(ca*cb*s2g + sa*c2g) + E2(4,2) = sa*-cb*cg - ca*c2b*sg + E2(4,3) = sqrt(3.0)*ca*sb*cb + E2(4,4) = ca*c2b*cg - sa*cb*sg + E2(4,5) = sb*sa*s2g - 0.5*ca*s2b*c2g + + E2(5,1) = -0.25*c2a*(c2b + 3.0)*s2g - c2g*cb*s2a + E2(5,2) = -sb*(c2a*cb*sg + s2a*cg) + E2(5,3) = 0.5*sqrt(3.0)*sb*sb*c2a + E2(5,4) = 0.5*c2a*s2b*cg - s2a*sb*sg + E2(5,5) = -s2a*cb*s2g + 0.25*c2a*(3.0 + c2b)*c2g + + if (flag_rotate_pdos_debug) then + identity3 = 0.0 + identity5 = 0.0 + do i = 1,3 + identity3(i,i) = 1.0 + end do + do i = 1,5 + identity5(i,i) = 1.0 + end do + + if(all(abs(matmul(E1, transpose(E1)) - identity3) < 1e-7)) then + write(*, '(/2x, "E1 matrix is orthogonal: PASS")', advance='no') + else + write(*, '(/2x, "E1 matrix is orthogonal: FAIL")', advance='no') + endif + if(all(abs(matmul(E2, transpose(E2)) - identity5) < 1e-7)) then + write(*, '(/2x, "E2 matrix is orthogonal: PASS")', advance='no') + else + write(*, '(/2x, "E2 matrix is orthogonal: FAIL")', advance='no') + endif + write(*,*) + end if + end subroutine construct_EulerMatrices + ! ----------------------------------------------------------------------------- + ! Subroutine get_pdos_axes + ! ----------------------------------------------------------------------------- + + !!****f* ProcModule/get_pdos_axes * + !! + !! NAME + !! get_pdos_axes - Create and normalise local rotation axes + !! USAGE + !! get_pdos_axes + !! PURPOSE + !! Create and normalise local rotation axes. Checks for mutual orthogonality. + !! INPUTS + !! NONE + !! USES + !! datatypes, local, GenComms + !! AUTHOR + !! C. Xu + !! CREATION DATE + !! 05/03/2026 + !! MODIFICATION HISTORY + !! 03/2026 - C. Xu: Orthogonality check + !! SOURCE + !! + subroutine get_pdos_axes + use datatypes + use local, ONLY: pdos_ax, pdos_ay, pdos_az + use GenComms, ONLY: cq_abort + implicit none + real(double), parameter :: tol = 1e-10 + + ! Normalise new x-axis correctly + pdos_ax = pdos_ax / norm2(pdos_ax) + ! Normalise new y-axis correctly + pdos_ay = pdos_ay / norm2(pdos_ay) + ! Normalise new z-axis correctly + pdos_az = pdos_az / norm2(pdos_az) + if (abs(dot_product(pdos_ax, pdos_ay)) > tol) & + call cq_abort("get_pdos_axes: pDOS_ax vector was not orthogonal to pDOS_ay.") + if (abs(dot_product(pdos_az, pdos_ay)) > tol) & + call cq_abort("get_pdos_axes: pDOS_az vector was not orthogonal to pDOS_ay.") + if (abs(dot_product(pdos_az, pdos_ax)) > tol) & + call cq_abort("get_pdos_axes: pDOS_az vector was not orthogonal to pDOS_ax.") + end subroutine get_pdos_axes + + subroutine initialise_A_mat(A1, A2) + use datatypes + implicit none + real(double), intent(out) :: A1(3, 3), A2(5, 5) + A1 = 0.0 + A2 = 0.0 + ! A1, FOR p-orbitals, l = 1 + A1(1, 3) = 1.0 + A1(2, 1) = 1.0 + A1(3, 2) = 1.0 + ! d-orbitals, l = 2 + A2(1, 2) = 1.0 + A2(2, 4) = 1.0 + A2(3, 1) = 1.0 + A2(4, 5) = 2.0 + A2(5, 3) = 2.0*sqrt(3.0) + end subroutine initialise_A_mat + + subroutine antisym_matrix(vector, matrix) + ! Create the antisymmetric matrix of a vector + use datatypes + implicit none + real(double), intent(in) :: vector(3) + real(double), intent(out) :: matrix(3,3) + matrix = 0.0 + matrix(1, 2) = -vector(3) + matrix(1, 3) = vector(2) + matrix(2, 1) = vector(3) + matrix(2, 3) = -vector(1) + matrix(3, 1) = -vector(2) + matrix(3, 2) = vector(1) + end subroutine antisym_matrix + ! ----------------------------------------------------------------------------- + ! Subroutine euler_from_axisangle + ! ----------------------------------------------------------------------------- + + !!****f* ProcModule/euler_from_axisangle * + !! + !! NAME + !! euler_from_axisangle + !! USAGE + !! euler_from_axisangle(axis, angle) + !! PURPOSE + !! Calculate Euler angles (alpha, beta, gamma) from a rotation axis and angle. + !! The rotation axis is assumed to be in [0, pi] + !! INPUTS + !! real(double), intent(in) :: axis(3) - axis of rotation + !! real(double), intent(in) :: angle - rotation angle in [0, pi] + !! USES + !! datatypes, numbers + !! AUTHOR + !! C. Xu + !! CREATION DATE + !! 20/08/2026 + !! MODIFICATION HISTORY + !! + !! SOURCE + subroutine euler_from_axisangle(axis, angle) + use datatypes + use numbers, ONLY: pi, one, two + implicit none + real(double), intent(in) :: axis(3), angle + + write(*, fmt= & + '(/2x,"Equivalent Euler angles. Using `Process.RotatePDOSMode 1`should give the same result.")') + write(*, fmt='(/4x,"alpha: ", (f10.5,1X))', advance="no") & + 180/pi*(datan2(axis(3)*sin(angle / two),cos(angle / two)) + datan2(axis(2), axis(1)) - (pi/two)) + write(*, fmt='(/4x,"beta: ", (f10.5,1X))', advance="no") & + 180/pi*2*asin(sqrt(axis(1)*axis(1) + axis(2)*axis(2)) * sin(angle / two)) + write(*, fmt='(/4x,"gamma: ", (f10.5,1X))', advance="no") & + 180/pi*(datan2(axis(3)*sin(angle/two),cos(angle/two)) - datan2(axis(2), axis(1)) + (pi/two)) + end subroutine + + ! ----------------------------------------------------------------------------- + ! Subroutine calculate_axis_angle + ! ----------------------------------------------------------------------------- + + !!****f* ProcModule/calculate_axis_angle * + !! + !! NAME + !! calculate_axis_angle + !! USAGE + !! call calculate_axis_angle(w1, w2, w3, axis, angle) + !! PURPOSE + !! Given a set of 3 orthonormal vectors, calculate the axis of rotation and angle needed to go from the simulation + !! cell axes to this basis (order matters). + !! This subroutine forms the basis change matrix and finds its normalised eigenvector with real eigenvalue one + !! This eigenvector is the rotation axis. The rotation angle is computed as the angle satisfying + !! cos \theta = (Tr(R) - 1) /2; sin\theta = -Tr(K_n R) / 2, \Theta = datan2(sin, cos) and shifted to the interval + !! [0, 2pi] + !! INPUTS + !! None + !! USES + !! datatypes, local, global, pao_format, GenComms + !! AUTHOR + !! C. Xu + !! CREATION DATE + !! 02/03/2026 + !! MODIFICATION HISTORY + !! 04/03/2026 - C. Xu: use LAPACK to find rotation axis + !! 05/03/2026 - C. Xu: Correct use of lapack and rotation angle + !! 09/03/2026 - C. Xu: use atan2 for rotation angle + !! 19/03/2026 - C. Xu: initial attempt to fix 90 degree rotation + !! 24/03/2026 - C. Xu: add error checking + !! 10/04/2026 - C. Xu: fixes to determinant check, negative angles and euler matrix + !! 26/05/2026, 08/06 - C. Xu: comments and clearer error messages + !! 20/08/2026 - C. Xu: Clean debug and negate axis/angle when necessary + !! SOURCE + !! + subroutine calculate_axis_angle(w1, w2, w3, axis, angle) + use datatypes + use numbers, ONLY: pi, one + use GenComms, ONLY: cq_abort + use local, ONLY: flag_rotate_pdos_debug + implicit none + external DGEEV + + ! Allow the user to define new coordinate system + ! Rotate from simulation axes (standard Cartesian) to local axes + real(double), intent(in) :: w1(3), w2(3), w3(3) + real(double), intent(out) :: axis(3), angle + real(double) :: pos_diff(3), det, sin_angle, cos_angle, angle2 + real(double) :: basis_matrix(3,3), temp(3,3), antisym_axis(3,3) + ! Local variables + integer :: i, j, N, LDA, LDVL, LDVR, INFO, LDWORK + real(double), allocatable:: A(:,:), WR(:), WI(:), VL(:,:), VR(:,:), WORK(:), tra + real(double), parameter :: tol = 1e-10 + real(double) :: real_eigenvalue,imag_eigenvalue + N = 3 + LDA = N + LDVL = LDA + LDVR = LDA + allocate(A(N,N)) + allocate(WR(N)) + allocate(WI(N)) + allocate(VL(LDVL, N)) + allocate(VR(LDVR, N)) + allocate(WORK(5*N)) + LDWORK = 5*N + ! Define change of basis matrix; + ! Columns of new basis in terms of old basis + basis_matrix(:,1) = w1 / norm2(w1) + basis_matrix(:,2) = w2 / norm2(w2) + basis_matrix(:,3) = w3 / norm2(w3) + ! Make copy because DGEEV destroys matrix + basis_matrix = (basis_matrix) + A = basis_matrix + ! If determinant flips sign, orientation of the coordinate system has changed + det = A(1,1)*(A(2,2)*A(3,3) - A(2,3)*A(3,2)) & + - A(1,2)*(A(2,1)*A(3,3) - A(2,3)*A(3,1)) & + + A(1,3)*(A(2,1)*A(3,2) - A(2,2)*A(3,1)) + + if (det < 1.0 - tol .or. det > 1.0 + tol) & + call cq_abort("calculate_axis_angle: determinant is not +1. Make sure input axes form a right-handed basis.") + call DGEEV("N", "V", N, A, LDA, WR, WI, VL, LDVL, VR, LDVR, & + WORK, LDWORK, INFO) + if (INFO /= 0) then + write(*, fmt='(A,I0)') 'DGEEV INFO = ', INFO + call cq_abort('calculate_axis_angle: DGEEV, Eigenvalues of basis change matrix could not be found!') + end if + ! Require eigenvector with eigenvalue 1, no complex, to find axis + do i = 1, N + if (abs(WR(i) - one) < tol .and. abs(WI(i)) < tol) then + axis = VR(:, i) ! Column i of VR is the i-th right eigenvector + exit + end if + end do + if (flag_rotate_pdos_debug) then + real_eigenvalue = 0.0 + imag_eigenvalue = 0.0 + write(*, '(/2x, "Eigenvalues/Eigenvectors of change of base matrix: ")') + do j = 1, N + ! Print eigenvalue + if (abs(WI(j)) < 1e-7) then + write(*,'(/4x, A,I3,A,F12.6)') 'Eigenvalue ', j, ': ', WR(j) + else + write(*,'(/4x, A,I3,A,F12.6,SP,F12.6,A)') 'Eigenvalue ', j, ': ', WR(j), WI(j), 'i' + end if + write(*,'(/6x, 2X,A)', advance="no") 'Eigenvector:' + do i = 1, N + if (WI(j) == 0.0_double) then + write(*,'(/7x, 4X,A,I3,A,F10.5)', advance="no") 'component ', i, ': ', VR(i,j) + else if (WI(j) > 0.0d0) then + ! First complex conjugate pair: v = VR(:,j) + i*VR(:,j+1) + write(*,'(/7x, 4X,A,I3,A,F10.5,SP,F10.5,A)', advance="no") & + 'component ', i, ': ', VR(i,j), VR(i,j+1), 'i' + else + ! Second complex conjugate pair: v = VR(:,j-1) - i*VR(:,j) + write(*,'(/7x, 4XA,I3,A,F10.5,SP,F10.5,A)', advance="no") & + 'component ', i, ': ', VR(i,j-1), -VR(i,j), 'i' + end if + end do + write(*,*) + end do + end if ! End debug: print all eigenvalues/vectors + if (any(axis /= axis)) & + call cq_abort("calculate_axis_angle: NaN in rotation axis.") + if (all(abs(axis) < tol)) & + call cq_abort("calculate_axis_angle: Rotation axis was (0,0,0). Cannot perform rotation.") + tra = basis_matrix(1,1) + basis_matrix(2,2) + basis_matrix(3,3) + antisym_axis = 0.0 + axis = axis / norm2(axis) + cos_angle = (tra - 1.0) + call antisym_matrix(axis, antisym_axis) + temp = matmul(antisym_axis, basis_matrix) + sin_angle = -(temp(1,1) + temp(2,2) + temp(3,3)) + if (abs(sin_angle) < tol .and. abs(cos_angle) < tol) & + call cq_abort("calculate_axis_angle: Both arguments to datan2 are zero.") + angle = datan2(sin_angle, cos_angle) + if (angle /= angle) & + call cq_abort("calculate_axis_angle: NaN rotation angle.") + if (flag_rotate_pdos_debug) & + write(*, fmt='(/4x,"Rotation angle (rad/deg) [-pi, pi]: ",2(f10.5,1X))') & + angle, angle * 180/pi + write(*, fmt='(/4x,"Rotation axis from LAPACK : ",3(f10.5,1X))') & + axis + if (angle < 0.0) then + ! Keep angle between [0, pi] and negate both axis and angle + angle = -angle + axis = -axis !R(n, theta) = R(-n, -theta) + end if + write(*, fmt='(/4x,"Rotation angle (rad/deg) [0, pi]: ",2(f10.5,1X))') & + angle, angle * 180/pi + write(*, fmt='(/4x,"Rotation axis: ", 3(f10.5,1X))') axis + end subroutine calculate_axis_angle + subroutine construct_rodrigues(axis, angle, matrix) + use datatypes + implicit none +! We will compute the rodrigues matrix, specifically R^T(axis, -angle) + real(double), intent(in) :: axis(3), angle + real(double), intent(out) :: matrix(9) + real(double) :: K(3, 3), KT(3, 3), identity(3, 3), temp(3, 3) + real(double) :: norm_axis(3) + K = 0.0 + identity = 0.0 + identity(1, 1) = 1.0 + identity(2, 2) = 1.0 + identity(3, 3) = 1.0 + norm_axis = axis / norm2(axis) + ! Define K + call antisym_matrix(norm_axis, K) + + KT = transpose(K) + ! Rodrigues rotation matrix but with -angle + temp = identity - (sin(angle)*KT) + (1 - cos(angle))*matmul(KT, KT) + ! Use transpose to restore matrix order in maths + matrix = reshape(transpose(temp), shape=(/9/)) + end subroutine construct_rodrigues + + subroutine construct_C1(rod, C1) + use datatypes + use local, ONLY: flag_rotate_pdos_debug + implicit none + real(double), intent(in) :: rod(9) + real(double), intent(out) :: C1(3, 3) + real(double) :: identity(3,3) + integer :: i + C1(1, 1) = rod(1) + C1(1, 2) = rod(2) + C1(1, 3) = rod(3) + C1(2, 1) = rod(4) + C1(2, 2) = rod(5) + C1(2, 3) = rod(6) + C1(3, 1) = rod(7) + C1(3, 2) = rod(8) + C1(3, 3) = rod(9) + if (flag_rotate_pdos_debug) then + identity = 0.0 + identity(1, 1) = 1.0 + identity(2, 2) = 1.0 + identity(3, 3) = 1.0 + if(all(abs(matmul(C1, transpose(C1)) - identity) < 1e-7)) then + write(*, '(/2x, "C1 matrix is orthogonal: PASS")', advance='no') + else + write(*, '(/2x, "C1 matrix is orthogonal: FAIL")', advance='no') + endif + end if + end subroutine construct_C1 + subroutine construct_C2(r, C2) + use datatypes + use local, ONLY: flag_rotate_pdos_debug + implicit none + real(double), intent(in) :: r(9) + real(double), intent(out) :: C2(5, 5) + real(double) :: identity(5,5) + + C2(1, 1) = r(6)*r(8) + r(5)*r(9) + C2(1, 2) = r(6)*r(7) + r(4)*r(9) + C2(1, 3) = r(5)*r(7) + r(4)*r(8) + C2(1, 5) = r(6)*r(9)*0.5 + C2(1, 4) = r(4)*r(7) + r(6)*r(9)*0.5 +! Row 2 + C2(2, 1) = r(3)*r(8) + r(2)*r(9) + C2(2, 2) = r(3)*r(7) + r(1)*r(9) + C2(2, 3) = r(2)*r(7) + r(1)*r(8) + C2(2, 5) = r(3)*r(9)*0.5 + C2(2, 4) = r(1)*r(7) + r(3)*r(9)*0.5 +! Row 3 + C2(3, 1) = r(3)*r(5) + r(2)*r(6) + C2(3, 2) = r(3)*r(4) + r(1)*r(6) + C2(3, 3) = r(2)*r(4) + r(1)*r(5) + C2(3, 5) = r(3)*r(6)*0.5 + C2(3, 4) = r(1)*r(4) + r(3)*r(6)*0.5 + +! Row 4 + C2(4, 1) = 2*(r(2)*r(3) - r(5)*r(6)) + C2(4, 2) = 2*(r(1)*r(3) - r(4)*r(6)) + C2(4, 3) = 2*(r(1)*r(2) - r(4)*r(5)) + C2(4, 5) = 0.5*(r(3)*r(3) - r(6)*r(6)) + C2(4, 4) = r(1)*r(1) - r(4)*r(4) + 0.5*(r(3)*r(3) - r(6)*r(6)) + +! Row 5 + C2(5, 1) = (4*r(8)*r(9)) - 2*(r(2)*r(3) + r(5)*r(6)) + C2(5, 2) = (4*r(7)*r(9)) - 2*(r(1)*r(3) + r(4)*r(6)) + C2(5, 3) = (4*r(7)*r(8)) - 2*(r(1)*r(2) + r(4)*r(5)) + C2(5, 5) = r(9)*r(9) - 0.5*(r(3)*r(3) + r(6)*r(6)) + C2(5, 4) = C2(5, 5) - r(1)*r(1) - r(4)*r(4) + 2*r(7)*r(7) + end subroutine construct_C2 + function inv(A) result(Ainv) + use datatypes + use GenComms, ONLY: cq_abort + implicit none + + real(double), dimension(:, :), intent(in) :: A + real(double), dimension(size(A, 1), size(A, 2)) :: Ainv + + real(double), dimension(size(A, 1)) :: work ! work array for LAPACK + integer, dimension(size(A, 1)) :: ipiv ! pivot indices + integer :: n, info + + external DGETRF + external DGETRI + + Ainv = A + n = size(A, 1) + + ! DGETRF computes an LU factorization of a general M-by-N matrix A + call DGETRF(n, n, Ainv, n, ipiv, info) + + if (info /= 0) then + call cq_abort("inv: Matrix is numerically singular!") + end if + + ! DGETRI computes the inverse using the LU factorization by DGETRF. + call DGETRI(n, Ainv, n, ipiv, work, n, info) + + if (info /= 0) then + call cq_abort("inv: Matrix inversion failed!") + end if + end function inv + + subroutine construct_Ul(angmom, Al, Cl, Ul) + use datatypes + implicit none + + integer, intent(in) :: angmom + real(double), intent(in) :: Al(2*angmom + 1, 2*angmom + 1), Cl(2*angmom + 1, 2*angmom + 1) + real(double), intent(out) :: Ul(:,:) + real(double) :: Al_inv(2*angmom + 1, 2*angmom + 1) + Al_inv = inv(Al) + Ul = matmul(Al_inv, matmul(Cl, Al)) + + end subroutine construct_Ul + ! ----------------------------------------------------------------------------- + ! Subroutine rotate_coefficients + ! ----------------------------------------------------------------------------- + + !!****f* ProcModule/rotate_coefficients * + !! + !! NAME + !! rotate_coefficients + !! USAGE + !! call rotate_coefficients + !! PURPOSE + !! Apply rotation matrices to wavefunction coefficients. Formally, for each l, each 2l + 1 coefficients + !! corresponding to a (m, zeta) are multiplied by the corresponding U^{(l)} rotation matrix. + !! As different rotation modes use different arrays and lookups, these are taken care of at the start of the + !! subroutine + !! The loop over all values is pasted and modified from process_pdos. + !! INPUTS + !! None + !! USES + !! datatypes, local, global, pao_format, GenComms + !! AUTHOR + !! C. Xu + !! CREATION DATE + !! 02/03/2026 + !! MODIFICATION HISTORY + !! 14/04/2026 - C. Xu: loop over active bands only + !! 03/05/2026 - C. Xu: correct g_atom_lookup + !! 12/05/2026 - C. Xu: correct de/allocation behaviour. Edit desc + !! 27/08/2026 - C. Xu: Similarity transform on U2 to account for orbital basis + !! SOURCE + !! + subroutine rotate_coefficients + use datatypes + use local, ONLY: n_bands_total, nkp, n_atoms_pDOS, evec_coeff, scaled_evec_coeff, & + pDOS_atom_index, band_full_to_active, rotate_pdos_natoms, find_neighbours, & + flag_rotate_pdos_mode, rotate_pdos_atoms, U1, U2 + use global_module, ONLY: nspin, species_glob + use pao_format, ONLY: pao + use GenComms, ONLY: cq_abort + + implicit none + ! Local variables + integer :: i_atom, i_spec, i_band, i_band_c, i_kp, i_spin, g_atom, i + integer :: i_l, i_z, nzeta, norbs, sf_offset, rotate_counter + integer, dimension(:), allocatable :: g_atom_lookup + ! Currently input axes depends on n_atoms_pDOS + if (flag_rotate_pdos_mode == 0) then + rotate_counter = n_atoms_pDOS + else + rotate_counter = rotate_pdos_natoms + end if + if (.not. allocated(g_atom_lookup)) allocate(g_atom_lookup(rotate_counter)) + ! Define the lookups correctly for each mode + if (flag_rotate_pdos_mode == 2) then + ! Atom counter should be the order the user input pDOSNeighbours block + g_atom_lookup = find_neighbours(1, :) + else if (flag_rotate_pdos_mode == 1 .or. flag_rotate_pdos_mode == 3) then + ! Atom counter should be the order the user input other blocks + g_atom_lookup = rotate_pdos_atoms + else + g_atom_lookup = pDOS_atom_index + end if + do i = 1, size(g_atom_lookup) + ! As in CONQUEST, m = -2 = -xy NOT xy, need to negate first row and column, except U2(1,1) + U2(1,2:,i) = -U2(1,2:,i) + U2(2:,1,i) = -U2(2:,1,i) + end do + + do i_spin = 1, nspin + do i_kp = 1, nkp + do i_band = 1, n_bands_total + i_band_c = band_full_to_active(i_band) + if(i_band_c > 0) then + do i_atom = 1, rotate_counter + ! Get global atom number from input + g_atom = g_atom_lookup(i_atom) + i_spec = species_glob(g_atom) + sf_offset = 0 + ! Include l = 0 to correctly calculate offset + do i_l = 0, pao(i_spec)%greatest_angmom + nzeta = pao(i_spec)%angmom(i_l)%n_zeta_in_angmom + norbs = 2*i_l + 1 + do i_z = 1, nzeta + select case(i_l) + case(1) + evec_coeff(sf_offset+1:sf_offset+norbs, g_atom, i_band_c, i_kp, i_spin) = & + matmul(U1(:,:,i_atom), evec_coeff(sf_offset+1:sf_offset+norbs, g_atom, i_band_c, i_kp, i_spin)) + scaled_evec_coeff(sf_offset+1:sf_offset+norbs, g_atom, i_band_c, i_kp, i_spin) = & + matmul(U1(:,:,i_atom), scaled_evec_coeff(sf_offset+1:sf_offset+norbs, g_atom, i_band_c, i_kp, i_spin)) + case(2) + evec_coeff(sf_offset+1:sf_offset+norbs, g_atom, i_band_c, i_kp, i_spin) = & + matmul(U2(:,:,i_atom), evec_coeff(sf_offset+1:sf_offset+norbs, g_atom, i_band_c, i_kp, i_spin)) + scaled_evec_coeff(sf_offset+1:sf_offset+norbs, g_atom, i_band_c, i_kp, i_spin) = & + matmul(U2(:,:,i_atom), scaled_evec_coeff(sf_offset+1:sf_offset+norbs, g_atom, i_band_c, i_kp, i_spin)) + end select + sf_offset = sf_offset + norbs + end do ! i_z + end do ! i_l + end do ! i_atom + end if ! if band is active + end do ! i_band + end do ! i_kp + end do ! i_spin + deallocate(g_atom_lookup) + end subroutine rotate_coefficients + + ! ----------------------------------------------------------------------------- + ! Subroutine nearest_neighbours + ! ----------------------------------------------------------------------------- + + !!****f* ProcModule/nearest_neighbours * + !! + !! NAME + !! nearest_neighbours - Find nearest neighbours and their bond vectors + !! USAGE + !! nearest_neighbours(atomno, bond) + !! PURPOSE + !! Evaluates the nearest neighbours using periodic boundary conditions + !! and returns their bond vectors + !! + !! INPUTS + !! integer, intent(in) :: atomno !atom to find neighbours of + !! real(double), intent(out), allocatable :: bond(:,:) ! array to hold bond + !! information regarding an atom + !! USES + !! datatypes, dimens, local, global, GenComms + !! AUTHOR + !! C. Xu + !! CREATION DATE + !! 15/04/2026 + !! MODIFICATION HISTORY + !! 12/05/2026 - C. Xu: correct de/allocation behaviour. Edit desc + !! SOURCE + !! + subroutine nearest_neighbours(atomno, bond) + use datatypes + use dimens, ONLY: r_super_x, r_super_y, r_super_z + use local, ONLY: find_neighbours, nghbr_arr, flag_rotate_pdos_debug + use global_module, ONLY: ni_in_cell, atom_coord, species_glob + use GenComms, ONLY: cq_abort + use species_module, ONLY: species_label + + implicit none + + integer, intent(in) :: atomno + real(double), intent(out), allocatable :: bond(:,:) + + + + + ! Local variables + real(double), parameter :: err = 1e-5 + real(double) :: distances(1:ni_in_cell), temp(1:ni_in_cell) + real(double) ::cx, cy, cz, dx, dy, dz, atomno_pos(3) + integer :: i, j, min_idx, ibond_max_len, idx_direction + ! Reset allocatables as this subroutine is called in a loop + ! and populated with new data every time + ! and populated with new data every time + if (allocated(nghbr_arr)) deallocate(nghbr_arr) + if (allocated(bond)) deallocate(bond) + + cx = atom_coord(1, atomno) + cy = atom_coord(2, atomno) + cz = atom_coord(3, atomno) + ! ACCOUNT FOR PERIODICITY + do i = 1, ni_in_cell + + if (i == atomno) then + distances(i) = huge(1.0d0) + cycle + end if + + dx = atom_coord(1,i) - cx + dy = atom_coord(2,i) - cy + dz = atom_coord(3,i) - cz + + dx = dx - r_super_x * nint(dx / r_super_x) + dy = dy - r_super_y * nint(dy / r_super_y) + dz = dz - r_super_z * nint(dz / r_super_z) + + distances(i) = dx*dx + dy*dy + dz*dz + end do + + temp = distances + ! Eliminate zeros - self-distance + where (abs(temp) < err) temp = huge(1.0) + select case (find_neighbours(2, findloc(find_neighbours(1,:), atomno, dim=1))) + case (0) + allocate(nghbr_arr(4)) + case (1) + allocate(nghbr_arr(6)) + case default + call cq_abort("nearest_neighbours: Did not correctly allocate nghbr_arr.") + end select + do j = 1, size(nghbr_arr) + min_idx = minloc(temp, 1) + nghbr_arr(j) = min_idx + temp(min_idx) = huge(1.0) + end do + + allocate(bond(3, size(nghbr_arr))) + + atomno_pos = (/cx, cy, cz/) + do i = 1, size(nghbr_arr) + bond(1, i) = atom_coord(1, nghbr_arr(i)) - atomno_pos(1) + bond(2, i) = atom_coord(2, nghbr_arr(i)) - atomno_pos(2) + bond(3, i) = atom_coord(3, nghbr_arr(i)) - atomno_pos(3) + + bond(1, i) = bond(1, i) - r_super_x * nint(bond(1, i) / r_super_x) + bond(2, i) = bond(2, i) - r_super_y * nint(bond(2, i) / r_super_y) + bond(3, i) = bond(3, i) - r_super_z * nint(bond(3, i) / r_super_z) + + end do + + if (flag_rotate_pdos_debug) then + write(*, fmt='(/2x,"Unit cell dimensions in Bohr (x,y,z)", 3(f10.5))', advance="no") & + r_super_x, r_super_y, r_super_z + write(*, fmt='(/2x,"Outputting cell-periodic distances (Bohr)", & + "of all atoms relative to atom", 1X,I0)', advance="no") atomno + do i = 1, ni_in_cell + if (i /= atomno) then + write(*, fmt='(/4x, I0, A, f10.5)', advance="no") & + i, " " // species_label(species_glob(i)), sqrt(distances(i)) + end if + end do + write(*, fmt='(/2x,"Outputting cell-periodic unnormalised bond lengths")') + + do i = 1, size(nghbr_arr) + write(*, fmt='(/4x, I0, A, 3f10.5)', advance="no") & + nghbr_arr(i), " " // species_label(species_glob(nghbr_arr(i))), bond(:,i) + end do + + end if + + end subroutine + ! ----------------------------------------------------------------------------- + ! Subroutine axes_from_nn + ! ----------------------------------------------------------------------------- + + !!****f* ProcModule/axes_from_nn * + !! + !! NAME + !! axes_from_nn - Gets rotation axes from bond vectors + !! USAGE + !! axes_from_nn(atomno, bond) + !! PURPOSE + !! This subroutine will use computed nearest_neighbours + !! This subroutine will use computed nearest_neighbours + !! and specified local geometry to construct a new set of local + !! axes to rotate the pDOS into + !! + !! For square-planar geometry: + !! + !! For square-planar geometry: + !! Longest bond is chosen as x_hat. + !! Bonds which are closest to orthogonality are + !! projected onto a plane defined by x_hat + !! y_hat is chosen from the projected direction + !! which has the minimal difference to its unprojected bond + !! z_hat is computed as the cross-product from these 2 directions, and thus + !! z_hat is computed as the cross-product from these 2 directions, and thus + !! point perpendicular to the planar geometry + !! For octahedra: choose z_hat as longest or shortest bond + !! Use this direction to define a plane - project the 4 atoms onto it + !! Call y_hat the bond which is the closest to the plane + !! Find x_hat = z cross y + !! INPUTS + !! integer, intent(in) :: atomno - atom to find neighbours of + !! real(double), intent(in) :: bond(:,:) - array to hold bond vectors + !! USES + !! datatypes, local, global, GenComms + !! AUTHOR + !! C. Xu + !! CREATION DATE + !! 15/04/2026 + !! MODIFICATION HISTORY + !! 20/04/2026 C. Xu - allow user to specify neighbour as secondary direction + !! SOURCE + !! + subroutine axes_from_nn(atomno, bond) + use datatypes + use local, ONLY: find_neighbours, nghbr_arr, pdos_ax, pdos_ay, pdos_az + use global_module, ONLY: ni_in_cell, atom_coord + use GenComms, ONLY: cq_abort + implicit none + ! This subroutine will use computed nearest_neighbours + ! and specified local geometry to construct a new set of local + ! axes to rotate the pDOS into + ! + ! For square-planar geometry: + ! Longest bond is chosen as x_hat. + ! Bonds which are closest to orthogonality are + ! projected onto a plane defined by x_hat + ! y_hat is chosen from the projected direction + ! which has the minimal difference to its unprojected bond + ! z_hat is computed as the cross-product from these 2 directions, and thus + ! point perpendicular to the planar geometry + ! For octahedra: choose z_hat as longest or shortest bond + ! Use this direction to define a plane - project the 4 atoms onto it + ! Call y_hat the bond which is the closest to the plane + ! Find x_hat by x_hat = z cross y + integer, intent(in) :: atomno + real(double), intent(in) :: bond(:,:) + ! Local variables + real(double) ::atomno_pos(3), plane_normal(3), proj_vector(3), cos_angle + real(double), allocatable :: dots(:), bond_lengths(:), norm_bond(:,:) + integer :: i, ibond_principal, ibond_sec, idx_direction, minormax, find_atomno + + allocate(dots(size(bond(1,:)))) + allocate(norm_bond(3, size(bond(1,:)))) + allocate(bond_lengths(size(bond(1,:)))) + do i = 1, size(bond_lengths) + bond_lengths(i) = sqrt(dot_product(bond(:, i), bond(:, i))) + if (bond_lengths(i) > 1e-10) then + norm_bond(:, i) = bond(:, i) / bond_lengths(i) + else + call cq_abort("axes_from_nn: Zero-length bond") + end if + end do + ! find_atomno: index of user supplied central atom, atomno + find_atomno = findloc(find_neighbours(1,:), atomno, dim=1) + minormax = find_neighbours(3, find_atomno) + if (minormax < 0) then + ibond_principal = minloc(bond_lengths,1) ! gets neighbour with min bond length + else if (minormax .eq. 0) then + ibond_principal = maxloc(bond_lengths,1) ! gets neighbour with maximum bond length + else + ibond_principal = findloc(nghbr_arr, minormax, dim=1) + end if + if (ibond_principal == 0) & + call cq_abort('axes_from_nn: principal neighbour not found for atom ', find_atomno) + + write(*,fmt='(/2x,"Principal neighbour found: atom ", (I0,1x))') nghbr_arr(ibond_principal) + dots = matmul(transpose(norm_bond), norm_bond(:, ibond_principal)) ! dot prod with chosen normal + ! select index of closest to perpendicular bond - there may be two + idx_direction = minloc(abs(dots), 1) + if (find_neighbours(4, find_atomno) == 0) then + ! If 0, we choose second direction by closest projection + call project_onto_plane(norm_bond(:, ibond_principal),norm_bond(:, idx_direction), proj_vector) + else + ! This entry > 0 -> corresponds to neighbour + ibond_sec = findloc(nghbr_arr, find_neighbours(4, find_atomno), dim=1) + ! Handle when input neighbour is not found + if (ibond_sec == 0) & + call cq_abort('axes_from_nn: second neighbour not found for atom ', find_atomno) + + if (abs(dots(ibond_sec)) > 0.5) & + print *, "WARNING: Chosen secondary neighbour appears to not be very perpendicular to principal." + write(*,fmt='(/2x,"Secondary neighbour found: atom ", (I0,1x))') nghbr_arr(ibond_sec) + + call project_onto_plane(norm_bond(:, ibond_principal),norm_bond(:, ibond_sec), proj_vector) + end if ! choice of 2nd direction + if ( find_neighbours(1, find_atomno) == find_neighbours(4, find_atomno)) & + call cq_abort("axes_from_nn: cannot have the principal neighbour also as the second direction") + select case (find_neighbours(2, find_atomno)) + case (0) + pdos_ax = norm_bond(:, ibond_principal) + pdos_ay = proj_vector + pdos_az = cross_product(pdos_ax, pdos_ay) + case (1) + pdos_ay = proj_vector + pdos_az = norm_bond(:, ibond_principal) + pdos_ax = cross_product(pdos_ay, pdos_az) + end select + pdos_ax = pdos_ax / norm2(pdos_ax) + pdos_ay = pdos_ay / norm2(pdos_ay) + pdos_az = pdos_az / norm2(pdos_az) + end subroutine axes_from_nn + function cross_product(a, b) result(cross) + use datatypes + implicit none + real(double), dimension(3), intent(in) :: a, b + real(double), dimension(3) :: cross + + cross(1) = a(2) * b(3) - a(3) * b(2) + cross(2) = a(3) * b(1) - a(1) * b(3) + cross(3) = a(1) * b(2) - a(2) * b(1) + end function cross_product + + subroutine project_onto_plane(plane_normal, vector, proj_vector) + use datatypes + implicit none + real(double), intent(in) :: plane_normal(3), vector(3) + real(double), intent(out) :: proj_vector(3) + real(double) :: norm_plane_normal(3) + ! Subroutine to project vector onto plane defined by plane_normal vector + norm_plane_normal = plane_normal / norm2(plane_normal) + + proj_vector = vector - ((dot_product(vector, norm_plane_normal))*norm_plane_normal) + end subroutine project_onto_plane + ! End pDOS rotation subroutines + subroutine process_band_structure use datatypes @@ -734,7 +1924,7 @@ subroutine process_band_structure use units, ONLY: HaToeV implicit none - + ! Local variables integer :: i_band, i_kp, i_spin, n_DOS_wid, n_band, n_min, n_max, i real(double) :: Ebin, dE_DOS, a, pf_DOS, spin_fac, dE @@ -868,7 +2058,7 @@ subroutine pao_dpao_to_grid(i_band, i_kp, i_spin, psi, dpsi) real(double) :: a, b, c, d, r1, r2, r3, r4, rr, kr, krx, kry, krz real(double), dimension(3) :: dsph_rl, dg complex(double_cplx) :: phase, phase_shift - + psi = zero dpsi = zero ! Grid spacing @@ -881,11 +2071,11 @@ subroutine pao_dpao_to_grid(i_band, i_kp, i_spin, psi, dpsi) if(atom_coord(3, i_atom) + RadiusAtomf(i_spec) >= stm_z_min) then ! Is the atom in STM region? kr = kx(i_kp)*atom_coord(1, i_atom) + ky(i_kp)*atom_coord(2, i_atom) + kz(i_kp)*atom_coord(3, i_atom) ! Find grid limits - minx = floor( (atom_coord(1, i_atom) - RadiusAtomf(i_spec))/dg(1) ) + minx = floor( (atom_coord(1, i_atom) - RadiusAtomf(i_spec))/dg(1) ) maxx = floor( (atom_coord(1, i_atom) + RadiusAtomf(i_spec))/dg(1) ) + 1 - miny = floor( (atom_coord(2, i_atom) - RadiusAtomf(i_spec))/dg(2) ) + miny = floor( (atom_coord(2, i_atom) - RadiusAtomf(i_spec))/dg(2) ) maxy = floor( (atom_coord(2, i_atom) + RadiusAtomf(i_spec))/dg(2) ) + 1 - minz = floor( (atom_coord(3, i_atom) - RadiusAtomf(i_spec))/dg(3) ) + minz = floor( (atom_coord(3, i_atom) - RadiusAtomf(i_spec))/dg(3) ) maxz = floor( (atom_coord(3, i_atom) + RadiusAtomf(i_spec))/dg(3) ) + 1 if(i_job==4.or.i_job==5) then ! STM not band density, so no z periodicity if(stm_z_min>zero) then @@ -998,7 +2188,7 @@ subroutine pao_to_grid(i_band, i_kp, i_spin, psi) real(double) :: a, b, c, d, r1, r2, r3, r4, rr, kr, krx, kry, krz real(double), dimension(3) :: dsph_rl, dg complex(double_cplx) :: phase, phase_shift - + psi = zero ! Grid spacing dg(1) = grid_x!/BohrToAng @@ -1011,11 +2201,11 @@ subroutine pao_to_grid(i_band, i_kp, i_spin, psi) kr = kx(i_kp)*atom_coord(1, i_atom) + ky(i_kp)*atom_coord(2, i_atom) + kz(i_kp)*atom_coord(3, i_atom) !phase = cmplx(cos(kr),sin(kr)) ! Find grid limits - minx = floor( (atom_coord(1, i_atom) - RadiusAtomf(i_spec))/dg(1) ) + minx = floor( (atom_coord(1, i_atom) - RadiusAtomf(i_spec))/dg(1) ) maxx = floor( (atom_coord(1, i_atom) + RadiusAtomf(i_spec))/dg(1) ) + 1 - miny = floor( (atom_coord(2, i_atom) - RadiusAtomf(i_spec))/dg(2) ) + miny = floor( (atom_coord(2, i_atom) - RadiusAtomf(i_spec))/dg(2) ) maxy = floor( (atom_coord(2, i_atom) + RadiusAtomf(i_spec))/dg(2) ) + 1 - minz = floor( (atom_coord(3, i_atom) - RadiusAtomf(i_spec))/dg(3) ) + minz = floor( (atom_coord(3, i_atom) - RadiusAtomf(i_spec))/dg(3) ) maxz = floor( (atom_coord(3, i_atom) + RadiusAtomf(i_spec))/dg(3) ) + 1 ! Account for STM limits if(i_job==4.or.i_job==5) then ! STM not band density @@ -1095,7 +2285,7 @@ subroutine read_domain(lun,proc,data) use numbers use local, ONLY: block_store, nxmin, nymin, nzmin, current, nptsx, nptsy, nptsz use block_module, only: n_pts_in_block, in_block_x,in_block_y,in_block_z - + implicit none ! Passed diff --git a/tools/PostProcessing/read_module.f90 b/tools/PostProcessing/read_module.f90 index e60d2bae..41545c62 100644 --- a/tools/PostProcessing/read_module.f90 +++ b/tools/PostProcessing/read_module.f90 @@ -4,7 +4,7 @@ module read implicit none character(len=80) :: block_file - + contains ! Read Conquest_input file for parameters from simulation, output parameters and coordinates @@ -24,7 +24,7 @@ subroutine read_input use pseudo_tm_info, only: setup_pseudo_info use GenComms, only: cq_abort use pseudopotential_common, only: pseudo_type, ABINIT, OLDPS, SIESTA - + implicit none character(len=80) :: input_string, proc_coords, tmp @@ -56,7 +56,7 @@ subroutine read_input ! Grid spacing n_grid_x = fdf_integer('Grid.PointsAlongX',0) n_grid_y = fdf_integer('Grid.PointsAlongY',0) - n_grid_z = fdf_integer('Grid.PointsAlongZ',0) + n_grid_z = fdf_integer('Grid.PointsAlongZ',0) if(n_grid_x>0.AND.n_grid_y>0.AND.n_grid_z>0) then dk = pi/min(n_grid_x, n_grid_y, n_grid_z) GridCutoff = half*dk*dk @@ -150,7 +150,7 @@ subroutine read_input flag_proc_band_str = 4 end if end if - ! + ! charge_stub = fdf_string(80,'Process.ChargeStub','chden') ! STM parameters ! NB Bias will be in volts @@ -315,7 +315,7 @@ subroutine read_input else allocate(pDOS_atom_index(n_atoms_pDOS)) if(fdf_block('pDOS_atoms')) then - if(1+block_end-block_start 3 reads + read (unit=input_array(block_start),fmt=*) pdos_ax + read (unit=input_array(block_start+1),fmt=*) pdos_ay + read (unit=input_array(block_start+2),fmt=*) pdos_az + call fdf_endblock + else if (fdf_block('pDOSEuler') .and. flag_rotate_pdos_mode == 1) then + if(1+block_end-block_start 2) & + call cq_abort("Local geometry flag in block pDOSNeighbours was not 0 or 1: ",& + 1+block_end-block_start,2) + if (find_neighbours(3,i) < -1 .or. find_neighbours(3,i) > ni_in_cell) & + call cq_abort("Input for principal axis must be -1 (shortest bond), 0 (longest bond) or neighbour in cell",& + 1+block_end-block_start,3) + if (find_neighbours(4,i) < 0 .or. find_neighbours(4,i) > ni_in_cell) & + call cq_abort("Input for second axis must be 0 or neighbour in cell",& + 1+block_end-block_start,4) + end do + call fdf_endblock + else if (fdf_block('pDOSAxisAngle') .and. flag_rotate_pdos_mode == 3) then + if(1+block_end-block_start Date: Fri, 28 Aug 2026 14:42:52 +0100 Subject: [PATCH 2/3] Add rotation code documentation --- docs/post-proc.rst | 95 +++++++++++++++++++++++++++++++++++++++++++-- docs/references.bib | 28 ++++++++++++- 2 files changed, 117 insertions(+), 6 deletions(-) diff --git a/docs/post-proc.rst b/docs/post-proc.rst index 0770e000..e3c6b0f3 100644 --- a/docs/post-proc.rst +++ b/docs/post-proc.rst @@ -75,7 +75,7 @@ and ``hilbert_make_blk.dat`` or ``raster_make_blk.dat``. Parameters that can be set include: :: - + Process.ChargeStub string (default: chden) The ChargeStub simply defines the filename which will be read, and @@ -267,7 +267,7 @@ Once the files have been generated by CONQUEST, the output can be processed by setting the output tag: :: - + Process.Job pdos This is all that is needed for the simplest output. The number of @@ -277,7 +277,7 @@ bins and smearing of the peaks can be set using: Process.sigma_DOS 0.002 Process.n_DOS 10001 - + To resolve the DOS by angular momentum as well as by atom, then the following flags can be set: @@ -296,7 +296,7 @@ The energy range for the projected DOS can also be specified: :: - + Process.min_DOS_E -0.35 Process.max_DOS_E 0.35 Process.WFRangeRelative T @@ -318,6 +318,93 @@ in the block ``pDOS_atoms``: Go to :ref:`top `. + +Rotated pDOS +------------ + +CONQUEST, by default, projects the orbitals along the simulation cell axes. Sometimes it is convenient to be able to project onto a coordinate system defined by an atom's local environment, e.g. bonds, to facilitate chemical analysis. CONQUEST supports rotating the wavefunction coefficients either by inputting the basis of the final coordinate system, using the algorithm implemented by Maintz *et al* and Romanowski *et al* [:cite:`pp-maintz2016`, :cite:`pp-romanowski2008`] or active Euler angles in the extrinsic :math:`zyz` convention. + +There are 4 operating modes that are possible: + +#. Supply a set of axes to apply to a set of (or all) atoms (mode 0) +#. Supply 3 Euler angles to apply to a set of (or all) atoms (mode 1) +#. Supply specific atoms and information about their neighbours to construct local axes for each atom (mode 2) +#. Supply specific atoms each with a rotation axis and angle (mode 3) + + +To enable any form of pDOS rotation, the flag ``Process.RotatePDOS T`` must be set in ``Conquest_input``. It defaults to ``F``, which will use the standard projection onto the simulation cell axes. + +To use mode 0, have the following snippet in ``Conquest_input``: +:: + + Process.RotatePDOS T (default: F) + Process.RotatePDOSMode 0 (default) + %block pDOSAxes + x1 y1 z1 ! new x + x2 y2 z2 ! new y + x3 y3 z3 ! new z + %endblock + Process.n_atoms_pDOS 2 + %block pDOS_atoms + 1 + 12 + %endblock + +The above snippet says to use the user-input axes defined in the ``pDOSAxes`` block. The first line defines the new :math:`x` direction along vector :math:`(x_1, y_1, z_1)`, then :math:`y` along :math:`(x_2, y_2, z_2)` and :math:`z` along :math:`(x_3, y_3, z_3)`. In this mode, all atoms will have this rotation applied to them by default. The ``pDOS_atoms`` block states which of these are written out to disk. The code will abort if the vectors provided do not form a right-handed coordinate system, or if the vectors are not orthogonal. It is recommended to supply the directions as vectors of integers - the code will take care of normalisation. + +If Euler angles are desired, set ``Process.RotatePDOSMode 1`` and in ``Conquest_input`` use the following lines: + +:: + + Process.RotatePDOS T + Process.RotatePDOSMode 1 + Process.RotatePDOS.NumAtoms int (required) + Process.RotatePDOSAngle deg (default: deg, optional: rad) + %block pDOSEuler + A alpha beta gamma !n entries, equal to Process.RotatePDOS.NumAtoms + %endblock pDOSEuler + + +Here, ``A``: the atom number to assign the Euler angles to. The default angle units is in degrees: ``Process.RotatePDOSAngle deg`` and units cannot be mixed. In the extrinsic :math:`zyz` convention, 3 consecutive rotations happen about a fixed set of coordinates, i.e. the CONQUEST cell simulation axes (as CONQUEST only supports orthorhombic cells), first about the :math:`z`-axis by :math:`\gamma` , then :math:`y`-axis by :math:`\beta` and then about :math:`z`-axis by :math:`\alpha` in a right-handed sense (i.e., looking down any of the axes towards the origin will mean the positive direction of rotation is anticlockwise). + +Determining local axes is possible depending on local geometries. + +:: + + Process.RotatePDOS T + Process.RotatePDOSMode 2 + Process.RotatePDOS.NumAtoms int (required) + %block pDOSNeighbours + A B C D !n entries, same number as in Process.RotatePDOS.NumAtoms + %endblock pDOSNeighbours + +where ``A``, ``B``, ``C``, ``D`` are integers with different conditions: + +* ``A``: the atom number to perform the rotation for +* ``B``: if it is ``0``, assume square planar geometry, if ``1`` assume octahedral. This controls nearest-neighbour searching. This wil search for the nearest 4 or 6 neighbours respectively. +* ``C``: if ``-1``, set the principal direction along the shortest bond. If ``0``, set along the longest bond. If it is a positive integer, it is interpreted as an atom number which must be a neighbour of the atom specified in ``A``. All runs will output the atom neighbours, so setting this to ``0`` or ``-1`` as a first run is recommended unless the neighbours are known beforehand +* ``D``: if ``0``, the second direction is chosen by the bond which changes the least under projection to the plane defined by the bond vector calculated from the neighbour specified in ``C``. If ``D`` is a positive integer, then it is interpreted as an atom number which must be a neighbour of the atom specified in ``A``. The code will error out if it is the same as ``C``. A warning will be displayed if the direction towards neighbour ``D`` is determined to deviate significantly away from orthogonal. + + +If the rotation axis and angle are known beforehand, using mode 3 can be convenient. + +:: + + Process.RotatePDOS T + Process.RotatePDOSMode 3 + Process.RotatePDOS.NumAtoms int (required) + %block pDOSAxisAngle + atom_number x y z angle !n entries, same number as in Process.RotatePDOS.NumAtoms + %endblock + +The atom number that this defines a local axis for must be specified. Then it is followed by three floats specifying the rotation axis in the basis of the simulation cell coordinate system. Finally, this is followed by a rotation angle which is positive when looking down the rotation axis and rotating anticlockwise. + +For modes 1, 2 and 3, setting the block ``pDOS_atoms`` is optional since only the specified atoms will be rotated whilst the rest will be in the simulation cell axes. + +For users interested in verbose output, set ``Process.RotatePDOSDebug T`` (default ``F``) in ``Conquest_input`` which will output detailed information about rotation matrices and orbital weights. + +Go to :ref:`top `. + .. _pp_band_str Band structure diff --git a/docs/references.bib b/docs/references.bib index e39ca8fa..d8705a21 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -1,10 +1,10 @@ %% This BibTeX bibliography file was created using BibDesk. %% https://bibdesk.sourceforge.io/ -%% Created for David Bowler at 2026-02-11 11:22:57 +0000 +%% Created for David Bowler at 2026-02-11 11:22:57 +0000 -%% Saved with string encoding Unicode (UTF-8) +%% Saved with string encoding Unicode (UTF-8) @@ -722,4 +722,28 @@ @article{Cococcioni:2005aa publisher = {American Physical Society}, doi = {10.1103/PhysRevB.71.035105}, url = {https://link.aps.org/doi/10.1103/PhysRevB.71.035105} +} + +@article{maintz2016, + author = {Maintz, S. and Esser, M. and Dronskowski, R.}, + journal = {Acta Physica Polonica B}, + month = {1}, + number = {4}, + pages = {1165}, + title = {{Efficient rotation of local basis functions using real spherical harmonics}}, + volume = {47}, + year = {2016}, + doi = {10.5506/aphyspolb.47.1165}, + url = {https://doi.org/10.5506/aphyspolb.47.1165}, +} +@article{romanowski2008, + author = {Romanowski, Zbigniew and Krukowski, Stanisław and Jalbout, Abraham F.}, + journal = {Acta Physica Polonica B}, + month = {8}, + number = {8}, + pages = {1985--1991}, + title = {{Transformation of Real Spherical Harmonics under Rotations}}, + volume = {39}, + year = {2008}, + url = {https://www.actaphys.uj.edu.pl/fulltext?series=Reg&vol=39&page=1985}, } \ No newline at end of file From f6921e521ccc2420d847b7faf2fece94649080a9 Mon Sep 17 00:00:00 2001 From: David Bowler Date: Fri, 28 Aug 2026 15:52:05 +0100 Subject: [PATCH 3/3] Mainly formatting updates to pDOS rotation code --- tools/PostProcessing/process_module.f90 | 2361 ++++++++++++----------- tools/PostProcessing/read_module.f90 | 105 +- 2 files changed, 1258 insertions(+), 1208 deletions(-) diff --git a/tools/PostProcessing/process_module.f90 b/tools/PostProcessing/process_module.f90 index b18473fc..253566fc 100644 --- a/tools/PostProcessing/process_module.f90 +++ b/tools/PostProcessing/process_module.f90 @@ -513,188 +513,187 @@ subroutine process_pdos ! Call pDOS rotation subroutines if desired. Also includes debug output if (flag_rotate_pdos) then - write(*,fmt='(2x,"Rotating wavefunction coefficients")') - call initialise_A_mat(A1, A2) - if (flag_rotate_pdos_debug .and. flag_rotate_pdos_mode /= 1) then - write(*, fmt='(/2x,"ROTATION DEBUG OUTPUT: MODE ", (I0,1x))') & - flag_rotate_pdos_mode - - write(*, fmt='(/2x,"ROTATION DEBUG OUTPUT: A^(l) MATRICES")') - write(*, fmt='(/4x, "A1: ")') - do j = 1, 3 - write(*, fmt='(/4x,3(f10.5,1X))') A1(j, :) - end do - write(*, fmt='(/4x, "A2: ")') - do j = 1, 5 - write(*, fmt='(/4x,5(f10.5,1X))') A2(j, :) - end do - write(*, fmt='(/4x, "For l = 1, orbital basis is [|y>, |z>, |x>]")') - write(*, fmt='(/4x, "For l = 2, orbital basis is [|-xy>, |yz>, |3z^2-r^2>, |xz>, |x^2-y^2>]")') - end if - if (allocated(U1)) deallocate(U1) - if (allocated(U2)) deallocate(U2) - - if (flag_rotate_pdos_mode == 0) then - write(*,fmt='(2x,"Using user input axes")') - allocate(U1(3,3,n_atoms_pDOS)) - allocate(U2(5,5,n_atoms_pDOS)) - call get_pdos_axes - - if(n_atoms_pDOS==ni_in_cell) then ! All atoms - write(*, fmt='(/2x,"New local axes for all atoms in unit cell:")') - else - write(*, fmt='(/2x,"New local axes for specified atoms: ", *(I0,1x))') & - pDOS_atom_index - end if - write(*, fmt='(/4x,"x: ",3(f10.5,1X))', advance="no") pdos_ax - write(*, fmt='(/4x,"y: ",3(f10.5,1X))', advance="no") pdos_ay - write(*, fmt='(/4x,"z: ",3(f10.5,1X))', advance="no") pdos_az - call calculate_axis_angle(pdos_ax, pdos_ay, pdos_az, axis, angle) - call construct_rodrigues(axis, angle, rod) - ! Construct all C^l matrices from rodrigues - call construct_C1(rod, C1) - call construct_C2(rod, C2) - do i = 1, n_atoms_pDOS - call construct_Ul(1, A1, C1, U1(:,:,i)) - call construct_Ul(2, A2, C2, U2(:,:,i)) - U1(:,:,i) = transpose(U1(:,:,i)) - U2(:,:,i) = transpose(U2(:,:,i)) - - if (flag_rotate_pdos_debug) then - ! call euler_from_axisangle(axis, angle) - write(*, fmt='(/4x, "C1: ")') - do j = 1, 3 - write(*, fmt='(/6x,3(f10.5,1X))', advance='no') C1(j,:) - end do - write(*, fmt='(/4x, "C2: ")') - do j = 1, 5 - write(*, fmt='(/6x,5(f10.5,1X))', advance='no') C2(j, :) - end do - call print_orbital_weights(i) - end if ! end rotation debug output + write(*,fmt='(2x,"Rotating wavefunction coefficients")') + call initialise_A_mat(A1, A2) + if (flag_rotate_pdos_debug .and. flag_rotate_pdos_mode /= 1) then + write(*, fmt='(/2x,"ROTATION DEBUG OUTPUT: MODE ", (I0,1x))') & + flag_rotate_pdos_mode + + write(*, fmt='(/2x,"ROTATION DEBUG OUTPUT: A^(l) MATRICES")') + write(*, fmt='(/4x, "A1: ")') + do j = 1, 3 + write(*, fmt='(/4x,3(f10.5,1X))', advance="no") A1(j, :) + end do + write(*, fmt='(/4x, "A2: ")') + do j = 1, 5 + write(*, fmt='(/4x,5(f10.5,1X))', advance="no") A2(j, :) + end do + write(*, fmt='(/4x, "For l = 1, orbital basis is [|y>, |z>, |x>]")') + write(*, fmt='(/4x, "For l = 2, orbital basis is [|-xy>, |yz>, |3z^2-r^2>, |xz>, |x^2-y^2>]")') + end if + if (allocated(U1)) deallocate(U1) + if (allocated(U2)) deallocate(U2) + + if (flag_rotate_pdos_mode == 0) then + write(*,fmt='(2x,"Using user input axes")') + allocate(U1(3,3,n_atoms_pDOS)) + allocate(U2(5,5,n_atoms_pDOS)) + call get_pdos_axes + + if(n_atoms_pDOS==ni_in_cell) then ! All atoms + write(*, fmt='(/2x,"New local axes for all atoms in unit cell:")') + else + write(*, fmt='(/2x,"New local axes for specified atoms: ", *(I0,1x))') & + pDOS_atom_index + end if + write(*, fmt='(/4x,"x: ",3(f10.5,1X))', advance="no") pdos_ax + write(*, fmt='(/4x,"y: ",3(f10.5,1X))', advance="no") pdos_ay + write(*, fmt='(/4x,"z: ",3(f10.5,1X))', advance="no") pdos_az + call calculate_axis_angle(pdos_ax, pdos_ay, pdos_az, axis, angle) + call construct_rodrigues(axis, angle, rod) + ! Construct all C^l matrices from rodrigues + call construct_C1(rod, C1) + call construct_C2(rod, C2) + do i = 1, n_atoms_pDOS + call construct_Ul(1, A1, C1, U1(:,:,i)) + call construct_Ul(2, A2, C2, U2(:,:,i)) + U1(:,:,i) = transpose(U1(:,:,i)) + U2(:,:,i) = transpose(U2(:,:,i)) + + if (flag_rotate_pdos_debug) then + ! call euler_from_axisangle(axis, angle) + write(*, fmt='(/4x, "C1: ")') + do j = 1, 3 + write(*, fmt='(/6x,3(f10.5,1X))', advance='no') C1(j,:) + end do + write(*, fmt='(/4x, "C2: ")') + do j = 1, 5 + write(*, fmt='(/6x,5(f10.5,1X))', advance='no') C2(j, :) + end do + call print_orbital_weights(i) + end if ! end rotation debug output end do ! end mode 0, axes input - else if (flag_rotate_pdos_mode == 1) then - write(*,fmt='(2x,"Using extrinsic Euler angles in active zyz convention")') - allocate(U1(3,3,rotate_pdos_natoms)) - allocate(U2(5,5,rotate_pdos_natoms)) - Qxyz = 0.0 - Qxyz(1,2) = 1.0 - Qxyz(2,3) = 1.0 - Qxyz(3,1) = 1.0 - do i = 1, rotate_pdos_natoms - rod = 0.0 - write(*, fmt='(/2x,"New local axes for specified atoms: ", *(I0,1x))') & - rotate_pdos_atoms(i) - call construct_EulerMatrices(E1, E2, i) - U1(:,:,i) = (E1) - U2(:,:,i) = (E2) - - ! We only explicitly construct wavefunction rotation matrices - ! However can extrapolate new local axes using E1 since orbital basis is {y,z,x} - write(*, fmt='(/4x,"Image of x,y,z under this active rotation, in standard basis x,y,z")', advance="no") - write(*, fmt='(/4x,"x: ",3(f10.5,1X))', advance="no") & - E1(3,3), E1(1,3), E1(2,3) - write(*, fmt='(/4x,"y: ",3(f10.5,1X))', advance="no") & - E1(3,1), E1(1,1), E1(2,1) - write(*, fmt='(/4x,"z: ",3(f10.5,1X))', advance="no") & - E1(3,2), E1(1,2), E1(2,2) - temp_matrix = matmul(inv(Qxyz), matmul(E1, Qxyz)) - temp_matrix = inv(temp_matrix) + else if (flag_rotate_pdos_mode == 1) then + write(*,fmt='(2x,"Using extrinsic Euler angles in active zyz convention")') + allocate(U1(3,3,rotate_pdos_natoms)) + allocate(U2(5,5,rotate_pdos_natoms)) + Qxyz = zero + Qxyz(1,2) = one + Qxyz(2,3) = one + Qxyz(3,1) = one + do i = 1, rotate_pdos_natoms + write(*, fmt='(/2x,"New local axes for specified atoms: ", *(I0,1x))') & + rotate_pdos_atoms(i) + call construct_EulerMatrices(E1, E2, i) + U1(:,:,i) = (E1) + U2(:,:,i) = (E2) + + ! We only explicitly construct wavefunction rotation matrices + ! However can extrapolate new local axes using E1 since orbital basis is {y,z,x} + write(*, fmt='(/4x,"Image of x,y,z under this active rotation, in standard basis x,y,z")', advance="no") + write(*, fmt='(/4x,"x: ",3(f10.5,1X))', advance="no") & + E1(3,3), E1(1,3), E1(2,3) + write(*, fmt='(/4x,"y: ",3(f10.5,1X))', advance="no") & + E1(3,1), E1(1,1), E1(2,1) + write(*, fmt='(/4x,"z: ",3(f10.5,1X))', advance="no") & + E1(3,2), E1(1,2), E1(2,2) + temp_matrix = matmul(inv(Qxyz), matmul(E1, Qxyz)) + temp_matrix = inv(temp_matrix) write(*, fmt='(/4x,"Equivalent basis transformation (inverse of active coordinates)")', advance="no") write(*, fmt='(/4x,"x: ",3(f10.5,1X))', advance="no") & - temp_matrix(:,1) - write(*, fmt='(/4x,"y: ",3(f10.5,1X))', advance="no") & - temp_matrix(:,2) - write(*, fmt='(/4x,"z: ",3(f10.5,1X))', advance="no") & - temp_matrix(:,3) - if (flag_rotate_pdos_debug) then - write(*, fmt='(/4x, "alpha(z) beta(y) gamma(z): ",3(f10.5,1X))') & - euler_angles(:,i) - call print_orbital_weights(i) - end if ! end rotation debug mode - end do ! end mode 1, euler angles - else if (flag_rotate_pdos_mode == 2) then - write(*,fmt='(2x,"Using user input atom numbers and local geometry")') - allocate(U1(3,3,rotate_pdos_natoms)) - allocate(U2(5,5,rotate_pdos_natoms)) - - do i = 1, rotate_pdos_natoms - rod = 0.0 - call nearest_neighbours(find_neighbours(1, i), bond) - write(*, fmt='(/2x,"Located neighours of atom ", I0, ": ", *(I0,1x))') & - find_neighbours(1,i), nghbr_arr - call axes_from_nn(find_neighbours(1, i), bond) - write(*, fmt='(/2x,"New local axes for atom ", I0, ": ")') & - find_neighbours(1,i) - write(*, fmt='(/4x,"x: ",3(f10.5,1X))', advance="no") pdos_ax - write(*, fmt='(/4x,"y: ",3(f10.5,1X))', advance="no") pdos_ay - write(*, fmt='(/4x,"z: ",3(f10.5,1X))', advance="no") pdos_az - call get_pdos_axes - call calculate_axis_angle(pdos_ax, pdos_ay, pdos_az, axis, angle) - call construct_rodrigues(axis, angle, rod) - ! Construct all C^l matrices from rodrigues - call construct_C1(rod, C1) - call construct_C2(rod, C2) - call construct_Ul(1, A1, C1, U1(:,:,i)) - call construct_Ul(2, A2, C2, U2(:,:,i)) - U1(:,:,i) = transpose(U1(:,:,i)) - U2(:,:,i) = transpose(U2(:,:,i)) - if (flag_rotate_pdos_debug) then - ! call euler_from_axisangle(axis, angle) - write(*, fmt='(/4x, "C1: ")') - do j = 1, 3 - write(*, fmt='(/6x,3(f10.5,1X))', advance='no') C1(j,:) - end do - write(*, fmt='(/4x, "C2: ")') - do j = 1, 5 - write(*, fmt='(/6x,5(f10.5,1X))', advance='no') C2(j, :) - end do - call print_orbital_weights(i) - end if ! end rotation debug output - end do ! end mode 2, neighbour input - else if (flag_rotate_pdos_mode == 3) then - write(*,fmt='(2x,"Using user input axis and angle")') - allocate(U1(3,3,rotate_pdos_natoms)) - allocate(U2(5,5,rotate_pdos_natoms)) - - do i = 1, rotate_pdos_natoms - rod = 0.0 - axis = axes_angles(1:3,i) - angle = axes_angles(4,i) - write(*, fmt='(/2x,"User input axis: ", 3(f10.5))') axis - write(*, fmt='(/2x,"User input angle (rad, deg) ", 2(f10.5))') angle, angle * (180.0 / pi) - call construct_rodrigues(axis, angle, rod) - write(*, fmt='(/2x,"New local axes for atom ", I0, ": ")') & - rotate_pdos_atoms(i) - - ! Read off columns for local axes - write(*, fmt='(/4x,"x: ",3(f10.5,1X))', advance="no") rod(1), rod(4), rod(7) - write(*, fmt='(/4x,"y: ",3(f10.5,1X))', advance="no") rod(2), rod(5), rod(8) - write(*, fmt='(/4x,"z: ",3(f10.5,1X))', advance="no") rod(3), rod(6), rod(9) - ! Construct all C^l matrices from rodrigues - call construct_C1(rod, C1) - call construct_C2(rod, C2) - call construct_Ul(1, A1, C1, U1(:,:,i)) - call construct_Ul(2, A2, C2, U2(:,:,i)) - U1(:,:,i) = transpose(U1(:,:,i)) - U2(:,:,i) = transpose(U2(:,:,i)) - if (flag_rotate_pdos_debug) then - call euler_from_axisangle(axis, angle) - write(*, fmt='(/4x, "C1: ")') - do j = 1, 3 - write(*, fmt='(/6x,3(f10.5,1X))', advance='no') C1(j,:) - end do - write(*, fmt='(/4x, "C2: ")') - do j = 1, 5 - write(*, fmt='(/6x,5(f10.5,1X))', advance='no') C2(j, :) - end do - call print_orbital_weights(i) - end if ! end rotation debug output - end do ! end mode 3, direct axis angle - end if ! pdos rotation mode - call rotate_coefficients - deallocate(U1) - deallocate(U2) + temp_matrix(:,1) + write(*, fmt='(/4x,"y: ",3(f10.5,1X))', advance="no") & + temp_matrix(:,2) + write(*, fmt='(/4x,"z: ",3(f10.5,1X))', advance="no") & + temp_matrix(:,3) + if (flag_rotate_pdos_debug) then + write(*, fmt='(/4x, "alpha(z) beta(y) gamma(z): ",3(f10.5,1X))') & + euler_angles(:,i) + call print_orbital_weights(i) + end if ! end rotation debug mode + end do ! end mode 1, euler angles + else if (flag_rotate_pdos_mode == 2) then + write(*,fmt='(2x,"Using user input atom numbers and local geometry")') + allocate(U1(3,3,rotate_pdos_natoms)) + allocate(U2(5,5,rotate_pdos_natoms)) + + do i = 1, rotate_pdos_natoms + rod = zero + call nearest_neighbours(find_neighbours(1, i), bond) + write(*, fmt='(/2x,"Located neighours of atom ", I0, ": ", *(I0,1x))') & + find_neighbours(1,i), nghbr_arr + call axes_from_nn(find_neighbours(1, i), bond) + write(*, fmt='(/2x,"New local axes for atom ", I0, ": ")') & + find_neighbours(1,i) + write(*, fmt='(/4x,"x: ",3(f10.5,1X))', advance="no") pdos_ax + write(*, fmt='(/4x,"y: ",3(f10.5,1X))', advance="no") pdos_ay + write(*, fmt='(/4x,"z: ",3(f10.5,1X))', advance="no") pdos_az + call get_pdos_axes + call calculate_axis_angle(pdos_ax, pdos_ay, pdos_az, axis, angle) + call construct_rodrigues(axis, angle, rod) + ! Construct all C^l matrices from rodrigues + call construct_C1(rod, C1) + call construct_C2(rod, C2) + call construct_Ul(1, A1, C1, U1(:,:,i)) + call construct_Ul(2, A2, C2, U2(:,:,i)) + U1(:,:,i) = transpose(U1(:,:,i)) + U2(:,:,i) = transpose(U2(:,:,i)) + if (flag_rotate_pdos_debug) then + ! call euler_from_axisangle(axis, angle) + write(*, fmt='(/4x, "C1: ")') + do j = 1, 3 + write(*, fmt='(/6x,3(f10.5,1X))', advance='no') C1(j,:) + end do + write(*, fmt='(/4x, "C2: ")') + do j = 1, 5 + write(*, fmt='(/6x,5(f10.5,1X))', advance='no') C2(j, :) + end do + call print_orbital_weights(i) + end if ! end rotation debug output + end do ! end mode 2, neighbour input + else if (flag_rotate_pdos_mode == 3) then + write(*,fmt='(2x,"Using user input axis and angle")') + allocate(U1(3,3,rotate_pdos_natoms)) + allocate(U2(5,5,rotate_pdos_natoms)) + + do i = 1, rotate_pdos_natoms + rod = zero + axis = axes_angles(1:3,i) + angle = axes_angles(4,i) + write(*, fmt='(/2x,"User input axis: ", 3(f10.5))') axis + write(*, fmt='(/2x,"User input angle (rad, deg) ", 2(f10.5))') angle, angle * (180.0_double / pi) + call construct_rodrigues(axis, angle, rod) + write(*, fmt='(/2x,"New local axes for atom ", I0, ": ")') & + rotate_pdos_atoms(i) + + ! Read off columns for local axes + write(*, fmt='(/4x,"x: ",3(f10.5,1X))', advance="no") rod(1), rod(4), rod(7) + write(*, fmt='(/4x,"y: ",3(f10.5,1X))', advance="no") rod(2), rod(5), rod(8) + write(*, fmt='(/4x,"z: ",3(f10.5,1X))', advance="no") rod(3), rod(6), rod(9) + ! Construct all C^l matrices from rodrigues + call construct_C1(rod, C1) + call construct_C2(rod, C2) + call construct_Ul(1, A1, C1, U1(:,:,i)) + call construct_Ul(2, A2, C2, U2(:,:,i)) + U1(:,:,i) = transpose(U1(:,:,i)) + U2(:,:,i) = transpose(U2(:,:,i)) + if (flag_rotate_pdos_debug) then + call euler_from_axisangle(axis, angle) + write(*, fmt='(/4x, "C1: ")') + do j = 1, 3 + write(*, fmt='(/6x,3(f10.5,1X))', advance='no') C1(j,:) + end do + write(*, fmt='(/4x, "C2: ")') + do j = 1, 5 + write(*, fmt='(/6x,5(f10.5,1X))', advance='no') C2(j, :) + end do + call print_orbital_weights(i) + end if ! end rotation debug output + end do ! end mode 3, direct axis angle + end if ! pdos rotation mode + call rotate_coefficients + deallocate(U1) + deallocate(U2) end if ! rotate pdos ! Set up storage based on pDOS per atom, or l/lm resolved per atom if(flag_lm_resolved) then @@ -920,998 +919,1056 @@ subroutine process_pdos end subroutine process_pdos ! pDOS rotation subroutines - ! misc subroutine for printing out orbital rotation matrices + ! misc subroutine for printing out orbital rotation matrices subroutine print_orbital_weights(atom_index) - use local, ONLY: U1, U2 - implicit none - - integer, intent(in) :: atom_index - - integer :: i,j,k - real(double) :: pdos_weight, identity3(3,3), identity5(5,5) - character(len=30), parameter :: p_orb(3) = (/ "|y>", "|z>","|x>"/) - character(len=30), parameter :: d_orb(5) = (/ "|-xy> ", "|yz> ", & - "|3z^2-r^2>", "|xz> ", "|x^2-y^2> "/) - character(len=50) :: pdos_weight_str - character(len=100) :: line - - write(*, fmt='(/4x, "U1: ")') - do j = 1, 3 - write(*, fmt='(/6x,3(f10.5,1X))', advance='no') U1(j,:,atom_index) - end do - write(*, fmt='(/4x, "p-orbital weight decomposition")') - write(*, '(/6x, A10, 5A10)') "", (trim(p_orb(j)), j = 1, 3) - do k = 1, 3 - write(*, '(6x, A10, 5F10.5)') trim(p_orb(k)), & + + use datatypes + use numbers, only : zero, one + use local, only: U1, U2 + + implicit none + + integer, intent(in) :: atom_index + + integer :: i,j,k + real(double) :: pdos_weight, identity3(3,3), identity5(5,5) + character(len=30), parameter :: p_orb(3) = (/ "|y>", "|z>","|x>"/) + character(len=30), parameter :: d_orb(5) = (/ "|-xy> ", "|yz> ", & + "|3z^2-r^2>", "|xz> ", "|x^2-y^2> "/) + character(len=50) :: pdos_weight_str + character(len=100) :: line + + write(*, fmt='(/4x, "U1: ")') + do j = 1, 3 + write(*, fmt='(/6x,3(f10.5,1X))', advance='no') U1(j,:,atom_index) + end do + write(*, fmt='(/4x, "p-orbital weight decomposition")') + write(*, '(/6x, A10, 5A10)') "", (trim(p_orb(j)), j = 1, 3) + do k = 1, 3 + write(*, '(6x, A10, 5F10.5)') trim(p_orb(k)), & (U1(k,j,atom_index)*U1(k,j,atom_index), j = 1, 3) - end do ! end printing l = 2 orbital weights - write(*, fmt='(/4x, "U2: ")') - ! Manual write because first row and column need to be negated - write(*, fmt='(/6x,5(f10.5,1X))', advance='no') U2(1,1,atom_index), -U2(1,2:,atom_index) - write(*, fmt='(/6x,5(f10.5,1X))', advance='no') -U2(2,1,atom_index), U2(2,2:,atom_index) - write(*, fmt='(/6x,5(f10.5,1X))', advance='no') -U2(3,1,atom_index), U2(3,2:,atom_index) - write(*, fmt='(/6x,5(f10.5,1X))', advance='no') -U2(4,1,atom_index), U2(4,2:,atom_index) - write(*, fmt='(/6x,5(f10.5,1X))', advance='no') -U2(5,1,atom_index), U2(5,2:,atom_index) - - write(*, fmt='(/4x, "d-orbital weight decomposition")') - - write(*, '(/6x, A10, 5A10)') "", (trim(d_orb(j)), j = 1, 5) - do k = 1, 5 - write(*, '(6x, A10, 5F10.5)') trim(d_orb(k)), & + end do ! end printing l = 2 orbital weights + write(*, fmt='(/4x, "U2: ")') + ! Manual write because first row and column need to be negated + write(*, fmt='(/6x,5(f10.5,1X))', advance='no') U2(1,1,atom_index), -U2(1,2:,atom_index) + write(*, fmt='(/6x,5(f10.5,1X))', advance='no') -U2(2,1,atom_index), U2(2,2:,atom_index) + write(*, fmt='(/6x,5(f10.5,1X))', advance='no') -U2(3,1,atom_index), U2(3,2:,atom_index) + write(*, fmt='(/6x,5(f10.5,1X))', advance='no') -U2(4,1,atom_index), U2(4,2:,atom_index) + write(*, fmt='(/6x,5(f10.5,1X))', advance='no') -U2(5,1,atom_index), U2(5,2:,atom_index) + + write(*, fmt='(/4x, "d-orbital weight decomposition")') + + write(*, '(/6x, A10, 5A10)') "", (trim(d_orb(j)), j = 1, 5) + do k = 1, 5 + write(*, '(6x, A10, 5F10.5)') trim(d_orb(k)), & (U2(k,j,atom_index)*U2(k,j,atom_index), j = 1, 5) - end do ! end printing l = 2 orbital weights - - ! Orthogonality check - identity3 = 0.0 - do i = 1, size(identity3(:,1)) - identity3(i,i) = 1.0 - end do - identity5 = 0.0 - do i = 1, size(identity5(:,1)) - identity5(i,i) = 1.0 - end do - if(all(abs(matmul(U1(:,:,atom_index), transpose(U1(:,:,atom_index))) - identity3) < 1e-7)) then - write(*, '(/2x, A)', advance='no') 'U1 matrix is orthogonal: PASS' - else - write(*, '(/2x, A)', advance='no') 'U1 matrix is orthogonal: FAIL' - endif - - if(all(abs(matmul(U2(:,:,atom_index), transpose(U2(:,:,atom_index))) - identity5) < 1e-7)) then - write(*, '(/2x, A)', advance='no') 'U2 matrix is orthogonal: PASS' - else - write(*, '(/2x, A)', advance='no') 'U2 matrix is orthogonal: FAIL' - endif - write(*,*) - end subroutine - - ! ----------------------------------------------------------------------------- - ! Subroutine construct_EulerMatrices - ! ----------------------------------------------------------------------------- - - !!****f* ProcModule/construct_EulerMatrices * - !! - !! NAME - !! construct_EulerMatrices - !! USAGE - !! construct_EulerMatrices(E1, E2, atom_index) - !! PURPOSE - !! Create rotation matrices using Euler angle input: right-handed, - !! active Euler angles about extrinsic (fixed) cell axes. Positive angle - !! is anticlockwise. - !! INPUTS - !! integer, intent(in) :: atom_index - atom to perform rotation for - !! USES - !! datatypes, numbers, local, GenComms - !! AUTHOR - !! C. Xu - !! CREATION DATE - !! 10/04/2026 - !! MODIFICATION HISTORY - !! 13/04/2026 - C. Xu: Add Euler matrix for d-orbitals - !! 20/07/2026 - C. Xu: Add output matrices in debugging mode - !! 30/07/2026 - C. Xu: Correct E2(3,2): alpha -> beta - !! 26/08/2026 - C. Xu: Fix orbital convention - !! SOURCE - !! Credits to R. Johnson for the mathematical expressions - !! Also see: Quantum Theory of Angular Momentum - subroutine construct_EulerMatrices(E1, E2, atom_index) - use datatypes - use numbers, ONLY: pi, two - use local, ONLY: euler_angles, flag_rotate_pdos_debug - use GenComms, ONLY: cq_abort - - implicit none - integer, intent(in) :: atom_index - real(double), intent(out) :: E1(3,3), E2(5,5) - real(double) :: ca, sa, cb, sb, cg, sg, c2a, s2a, c2b, s2b, c2g, s2g - real(double) :: euler_alpha, euler_beta, euler_gamma - real(double) :: identity3(3,3), identity5(5,5) - integer :: i - euler_alpha = euler_angles(1,atom_index) - euler_beta = euler_angles(2,atom_index) - euler_gamma = euler_angles(3,atom_index) - E1 = 0.0 - E2 = 0.0 - - ca = cos(euler_alpha) - sa = sin(euler_alpha) - cb = cos(euler_beta) - sb = sin(euler_beta) - cg = cos(euler_gamma) - sg = sin(euler_gamma) - - c2a = cos(two*euler_alpha) - s2a = sin(two*euler_alpha) - c2b = cos(two*euler_beta) - s2b = sin(two*euler_beta) - c2g = cos(two*euler_gamma) - s2g = sin(two*euler_gamma) - - ! Rotation of l = 1 coefficients - E1(1,1) = ca*cg - cb*sa*sg - E1(1,2) = sa*sb - E1(1,3) = cb*cg*sa + ca*sg - E1(2,1) = sb*sg - E1(2,2) = cb - E1(2,3) = -cg*sb - E1(3,1) = -cg*sa - ca*cb*sg - E1(3,2) = ca*sb - E1(3,3) = ca*cb*cg - sa*sg - - - ! Rotation of l = 2 coefficients - E2(1,1) = c2a*cb*c2g - sa*ca*(c2b+3.0)*sg*cg - E2(1,2) = sb*(c2a*cg - 2.0*sa*ca*cb*sg) - E2(1,3) = sqrt(3.0)*sa*ca*sb*sb - E2(1,4) = sb*(s2a*cb*cg + c2a*sg) - E2(1,5) = cb*c2a*s2g +0.25*s2a*(3.0+c2b)*c2g - - E2(2,1) = sb*(sa*cb*s2g - ca*c2g) - E2(2,2) = ca*cb*cg - sa*c2b*sg - E2(2,3) = sqrt(3.0)*sa*sb*cb - E2(2,4) = c2b*sa*cg + ca*cb*sg - E2(2,5) = -sb*(sa*cb*c2g +s2g*ca) - - E2(3,1) = -sqrt(3.0)*sb*sb*sg*cg - E2(3,2) = sqrt(3.0)*sb*cb*sg - E2(3,3) = 0.25*(3.0*c2b+1.0) - E2(3,4) = -sqrt(3.0)*sb*cb*cg - E2(3,5) = 0.5*sqrt(3.0)*sb*sb*c2g - - E2(4,1) = sb*(ca*cb*s2g + sa*c2g) - E2(4,2) = sa*-cb*cg - ca*c2b*sg - E2(4,3) = sqrt(3.0)*ca*sb*cb - E2(4,4) = ca*c2b*cg - sa*cb*sg - E2(4,5) = sb*sa*s2g - 0.5*ca*s2b*c2g - - E2(5,1) = -0.25*c2a*(c2b + 3.0)*s2g - c2g*cb*s2a - E2(5,2) = -sb*(c2a*cb*sg + s2a*cg) - E2(5,3) = 0.5*sqrt(3.0)*sb*sb*c2a - E2(5,4) = 0.5*c2a*s2b*cg - s2a*sb*sg - E2(5,5) = -s2a*cb*s2g + 0.25*c2a*(3.0 + c2b)*c2g - - if (flag_rotate_pdos_debug) then - identity3 = 0.0 - identity5 = 0.0 - do i = 1,3 - identity3(i,i) = 1.0 - end do - do i = 1,5 - identity5(i,i) = 1.0 - end do - - if(all(abs(matmul(E1, transpose(E1)) - identity3) < 1e-7)) then - write(*, '(/2x, "E1 matrix is orthogonal: PASS")', advance='no') - else - write(*, '(/2x, "E1 matrix is orthogonal: FAIL")', advance='no') - endif - if(all(abs(matmul(E2, transpose(E2)) - identity5) < 1e-7)) then - write(*, '(/2x, "E2 matrix is orthogonal: PASS")', advance='no') - else - write(*, '(/2x, "E2 matrix is orthogonal: FAIL")', advance='no') - endif - write(*,*) - end if - end subroutine construct_EulerMatrices - ! ----------------------------------------------------------------------------- - ! Subroutine get_pdos_axes - ! ----------------------------------------------------------------------------- - - !!****f* ProcModule/get_pdos_axes * - !! - !! NAME - !! get_pdos_axes - Create and normalise local rotation axes - !! USAGE - !! get_pdos_axes - !! PURPOSE - !! Create and normalise local rotation axes. Checks for mutual orthogonality. - !! INPUTS - !! NONE - !! USES - !! datatypes, local, GenComms - !! AUTHOR - !! C. Xu - !! CREATION DATE - !! 05/03/2026 - !! MODIFICATION HISTORY - !! 03/2026 - C. Xu: Orthogonality check - !! SOURCE - !! + end do ! end printing l = 2 orbital weights + + ! Orthogonality check + identity3 = zero + do i = 1, size(identity3(:,1)) + identity3(i,i) = one + end do + identity5 = zero + do i = 1, size(identity5(:,1)) + identity5(i,i) = one + end do + if(all(abs(matmul(U1(:,:,atom_index), transpose(U1(:,:,atom_index))) - identity3) < 1e-7)) then + write(*, '(/2x, A)', advance='no') 'U1 matrix is orthogonal: PASS' + else + write(*, '(/2x, A)', advance='no') 'U1 matrix is orthogonal: FAIL' + endif + + if(all(abs(matmul(U2(:,:,atom_index), transpose(U2(:,:,atom_index))) - identity5) < 1e-7)) then + write(*, '(/2x, A)', advance='no') 'U2 matrix is orthogonal: PASS' + else + write(*, '(/2x, A)', advance='no') 'U2 matrix is orthogonal: FAIL' + endif + write(*,*) + end subroutine print_orbital_weights + + ! ----------------------------------------------------------------------------- + ! Subroutine construct_EulerMatrices + ! ----------------------------------------------------------------------------- + + !!****f* ProcModule/construct_EulerMatrices * + !! + !! NAME + !! construct_EulerMatrices + !! USAGE + !! construct_EulerMatrices(E1, E2, atom_index) + !! PURPOSE + !! Create rotation matrices using Euler angle input: right-handed, + !! active Euler angles about extrinsic (fixed) cell axes. Positive angle + !! is anticlockwise. + !! INPUTS + !! integer, intent(in) :: atom_index - atom to perform rotation for + !! USES + !! datatypes, numbers, local, GenComms + !! AUTHOR + !! C. Xu + !! CREATION DATE + !! 10/04/2026 + !! MODIFICATION HISTORY + !! 13/04/2026 - C. Xu: Add Euler matrix for d-orbitals + !! 20/07/2026 - C. Xu: Add output matrices in debugging mode + !! 30/07/2026 - C. Xu: Correct E2(3,2): alpha -> beta + !! 26/08/2026 - C. Xu: Fix orbital convention + !! SOURCE + !! Credits to R. Johnson for the mathematical expressions + !! Also see: Quantum Theory of Angular Momentum + subroutine construct_EulerMatrices(E1, E2, atom_index) + use datatypes + use numbers, ONLY: pi, zero, one, two, three, quarter, half + use local, ONLY: euler_angles, flag_rotate_pdos_debug + use GenComms, ONLY: cq_abort + + implicit none + integer, intent(in) :: atom_index + real(double), intent(out) :: E1(3,3), E2(5,5) + real(double) :: ca, sa, cb, sb, cg, sg, c2a, s2a, c2b, s2b, c2g, s2g + real(double) :: euler_alpha, euler_beta, euler_gamma + real(double) :: identity3(3,3), identity5(5,5) + integer :: i + euler_alpha = euler_angles(1,atom_index) + euler_beta = euler_angles(2,atom_index) + euler_gamma = euler_angles(3,atom_index) + E1 = zero + E2 = zero + + ca = cos(euler_alpha) + sa = sin(euler_alpha) + cb = cos(euler_beta) + sb = sin(euler_beta) + cg = cos(euler_gamma) + sg = sin(euler_gamma) + + c2a = cos(two*euler_alpha) + s2a = sin(two*euler_alpha) + c2b = cos(two*euler_beta) + s2b = sin(two*euler_beta) + c2g = cos(two*euler_gamma) + s2g = sin(two*euler_gamma) + + ! Rotation of l = 1 coefficients + E1(1,1) = ca*cg - cb*sa*sg + E1(1,2) = sa*sb + E1(1,3) = cb*cg*sa + ca*sg + E1(2,1) = sb*sg + E1(2,2) = cb + E1(2,3) = -cg*sb + E1(3,1) = -cg*sa - ca*cb*sg + E1(3,2) = ca*sb + E1(3,3) = ca*cb*cg - sa*sg + + + ! Rotation of l = 2 coefficients + E2(1,1) = c2a*cb*c2g - sa*ca*(c2b+three)*sg*cg + E2(1,2) = sb*(c2a*cg - two*sa*ca*cb*sg) + E2(1,3) = sqrt(three)*sa*ca*sb*sb + E2(1,4) = sb*(s2a*cb*cg + c2a*sg) + E2(1,5) = cb*c2a*s2g +quarter*s2a*(three+c2b)*c2g + + E2(2,1) = sb*(sa*cb*s2g - ca*c2g) + E2(2,2) = ca*cb*cg - sa*c2b*sg + E2(2,3) = sqrt(three)*sa*sb*cb + E2(2,4) = c2b*sa*cg + ca*cb*sg + E2(2,5) = -sb*(sa*cb*c2g +s2g*ca) + + E2(3,1) = -sqrt(three)*sb*sb*sg*cg + E2(3,2) = sqrt(three)*sb*cb*sg + E2(3,3) = quarter*(three*c2b+one) + E2(3,4) = -sqrt(three)*sb*cb*cg + E2(3,5) = half*sqrt(three)*sb*sb*c2g + + E2(4,1) = sb*(ca*cb*s2g + sa*c2g) + E2(4,2) = -sa*cb*cg - ca*c2b*sg + E2(4,3) = sqrt(three)*ca*sb*cb + E2(4,4) = ca*c2b*cg - sa*cb*sg + E2(4,5) = sb*sa*s2g - half*ca*s2b*c2g + + E2(5,1) = -quarter*c2a*(c2b + three)*s2g - c2g*cb*s2a + E2(5,2) = -sb*(c2a*cb*sg + s2a*cg) + E2(5,3) = half*sqrt(three)*sb*sb*c2a + E2(5,4) = half*c2a*s2b*cg - s2a*sb*sg + E2(5,5) = -s2a*cb*s2g + quarter*c2a*(three + c2b)*c2g + + if (flag_rotate_pdos_debug) then + identity3 = zero + identity5 = zero + do i = 1,3 + identity3(i,i) = one + end do + do i = 1,5 + identity5(i,i) = one + end do + + if(all(abs(matmul(E1, transpose(E1)) - identity3) < 1e-7)) then + write(*, '(/2x, "E1 matrix is orthogonal: PASS")', advance='no') + else + write(*, '(/2x, "E1 matrix is orthogonal: FAIL")', advance='no') + endif + if(all(abs(matmul(E2, transpose(E2)) - identity5) < 1e-7)) then + write(*, '(/2x, "E2 matrix is orthogonal: PASS")', advance='no') + else + write(*, '(/2x, "E2 matrix is orthogonal: FAIL")', advance='no') + endif + write(*,*) + end if + end subroutine construct_EulerMatrices + ! ----------------------------------------------------------------------------- + ! Subroutine get_pdos_axes + ! ----------------------------------------------------------------------------- + + !!****f* ProcModule/get_pdos_axes * + !! + !! NAME + !! get_pdos_axes - Create and normalise local rotation axes + !! USAGE + !! get_pdos_axes + !! PURPOSE + !! Create and normalise local rotation axes. Checks for mutual orthogonality. + !! INPUTS + !! NONE + !! USES + !! datatypes, local, GenComms + !! AUTHOR + !! C. Xu + !! CREATION DATE + !! 05/03/2026 + !! MODIFICATION HISTORY + !! 03/2026 - C. Xu: Orthogonality check + !! SOURCE + !! subroutine get_pdos_axes - use datatypes - use local, ONLY: pdos_ax, pdos_ay, pdos_az - use GenComms, ONLY: cq_abort - implicit none - real(double), parameter :: tol = 1e-10 - - ! Normalise new x-axis correctly - pdos_ax = pdos_ax / norm2(pdos_ax) - ! Normalise new y-axis correctly - pdos_ay = pdos_ay / norm2(pdos_ay) - ! Normalise new z-axis correctly - pdos_az = pdos_az / norm2(pdos_az) - if (abs(dot_product(pdos_ax, pdos_ay)) > tol) & - call cq_abort("get_pdos_axes: pDOS_ax vector was not orthogonal to pDOS_ay.") - if (abs(dot_product(pdos_az, pdos_ay)) > tol) & - call cq_abort("get_pdos_axes: pDOS_az vector was not orthogonal to pDOS_ay.") - if (abs(dot_product(pdos_az, pdos_ax)) > tol) & - call cq_abort("get_pdos_axes: pDOS_az vector was not orthogonal to pDOS_ax.") + + use datatypes + use numbers, only: RD_ERR + use local, only: pdos_ax, pdos_ay, pdos_az + use GenComms, only: cq_abort + + implicit none + + ! Normalise new x-axis correctly + pdos_ax = pdos_ax / norm2(pdos_ax) + ! Normalise new y-axis correctly + pdos_ay = pdos_ay / norm2(pdos_ay) + ! Normalise new z-axis correctly + pdos_az = pdos_az / norm2(pdos_az) + if (abs(dot_product(pdos_ax, pdos_ay)) > RD_ERR) & + call cq_abort("get_pdos_axes: pDOS_ax vector was not orthogonal to pDOS_ay.") + if (abs(dot_product(pdos_az, pdos_ay)) > RD_ERR) & + call cq_abort("get_pdos_axes: pDOS_az vector was not orthogonal to pDOS_ay.") + if (abs(dot_product(pdos_az, pdos_ax)) > RD_ERR) & + call cq_abort("get_pdos_axes: pDOS_az vector was not orthogonal to pDOS_ax.") end subroutine get_pdos_axes - subroutine initialise_A_mat(A1, A2) - use datatypes - implicit none - real(double), intent(out) :: A1(3, 3), A2(5, 5) - A1 = 0.0 - A2 = 0.0 - ! A1, FOR p-orbitals, l = 1 - A1(1, 3) = 1.0 - A1(2, 1) = 1.0 - A1(3, 2) = 1.0 - ! d-orbitals, l = 2 - A2(1, 2) = 1.0 - A2(2, 4) = 1.0 - A2(3, 1) = 1.0 - A2(4, 5) = 2.0 - A2(5, 3) = 2.0*sqrt(3.0) - end subroutine initialise_A_mat - - subroutine antisym_matrix(vector, matrix) - ! Create the antisymmetric matrix of a vector - use datatypes - implicit none - real(double), intent(in) :: vector(3) - real(double), intent(out) :: matrix(3,3) - matrix = 0.0 - matrix(1, 2) = -vector(3) - matrix(1, 3) = vector(2) - matrix(2, 1) = vector(3) - matrix(2, 3) = -vector(1) - matrix(3, 1) = -vector(2) - matrix(3, 2) = vector(1) - end subroutine antisym_matrix - ! ----------------------------------------------------------------------------- - ! Subroutine euler_from_axisangle - ! ----------------------------------------------------------------------------- - - !!****f* ProcModule/euler_from_axisangle * - !! - !! NAME - !! euler_from_axisangle - !! USAGE - !! euler_from_axisangle(axis, angle) - !! PURPOSE - !! Calculate Euler angles (alpha, beta, gamma) from a rotation axis and angle. - !! The rotation axis is assumed to be in [0, pi] - !! INPUTS - !! real(double), intent(in) :: axis(3) - axis of rotation - !! real(double), intent(in) :: angle - rotation angle in [0, pi] - !! USES - !! datatypes, numbers - !! AUTHOR - !! C. Xu - !! CREATION DATE - !! 20/08/2026 - !! MODIFICATION HISTORY - !! - !! SOURCE - subroutine euler_from_axisangle(axis, angle) - use datatypes - use numbers, ONLY: pi, one, two - implicit none - real(double), intent(in) :: axis(3), angle - - write(*, fmt= & + subroutine initialise_A_mat(A1, A2) + + use datatypes + use numbers, only: zero, one, two, three + + implicit none + + real(double), intent(out) :: A1(3, 3), A2(5, 5) + + A1 = zero + A2 = zero + ! A1, FOR p-orbitals, l = 1 + A1(1, 3) = one + A1(2, 1) = one + A1(3, 2) = one + ! d-orbitals, l = 2 + A2(1, 2) = one + A2(2, 4) = one + A2(3, 1) = one + A2(4, 5) = two + A2(5, 3) = two*sqrt(three) + end subroutine initialise_A_mat + + subroutine antisym_matrix(vector, matrix) + ! Create the antisymmetric matrix of a vector + + use datatypes + use numbers, only: zero + + implicit none + + real(double), intent(in) :: vector(3) + real(double), intent(out) :: matrix(3,3) + + matrix = zero + matrix(1, 2) = -vector(3) + matrix(1, 3) = vector(2) + matrix(2, 1) = vector(3) + matrix(2, 3) = -vector(1) + matrix(3, 1) = -vector(2) + matrix(3, 2) = vector(1) + end subroutine antisym_matrix + + ! ----------------------------------------------------------------------------- + ! Subroutine euler_from_axisangle + ! ----------------------------------------------------------------------------- + + !!****f* ProcModule/euler_from_axisangle * + !! + !! NAME + !! euler_from_axisangle + !! USAGE + !! euler_from_axisangle(axis, angle) + !! PURPOSE + !! Calculate Euler angles (alpha, beta, gamma) from a rotation axis and angle. + !! The rotation axis is assumed to be in [0, pi] + !! INPUTS + !! real(double), intent(in) :: axis(3) - axis of rotation + !! real(double), intent(in) :: angle - rotation angle in [0, pi] + !! USES + !! datatypes, numbers + !! AUTHOR + !! C. Xu + !! CREATION DATE + !! 20/08/2026 + !! MODIFICATION HISTORY + !! + !! SOURCE + subroutine euler_from_axisangle(axis, angle) + + use datatypes + use numbers, ONLY: pi, one, two + + implicit none + + real(double), intent(in) :: axis(3), angle + + write(*, fmt= & '(/2x,"Equivalent Euler angles. Using `Process.RotatePDOSMode 1`should give the same result.")') - write(*, fmt='(/4x,"alpha: ", (f10.5,1X))', advance="no") & - 180/pi*(datan2(axis(3)*sin(angle / two),cos(angle / two)) + datan2(axis(2), axis(1)) - (pi/two)) - write(*, fmt='(/4x,"beta: ", (f10.5,1X))', advance="no") & - 180/pi*2*asin(sqrt(axis(1)*axis(1) + axis(2)*axis(2)) * sin(angle / two)) - write(*, fmt='(/4x,"gamma: ", (f10.5,1X))', advance="no") & - 180/pi*(datan2(axis(3)*sin(angle/two),cos(angle/two)) - datan2(axis(2), axis(1)) + (pi/two)) - end subroutine - - ! ----------------------------------------------------------------------------- - ! Subroutine calculate_axis_angle - ! ----------------------------------------------------------------------------- - - !!****f* ProcModule/calculate_axis_angle * - !! - !! NAME - !! calculate_axis_angle - !! USAGE - !! call calculate_axis_angle(w1, w2, w3, axis, angle) - !! PURPOSE - !! Given a set of 3 orthonormal vectors, calculate the axis of rotation and angle needed to go from the simulation - !! cell axes to this basis (order matters). - !! This subroutine forms the basis change matrix and finds its normalised eigenvector with real eigenvalue one - !! This eigenvector is the rotation axis. The rotation angle is computed as the angle satisfying - !! cos \theta = (Tr(R) - 1) /2; sin\theta = -Tr(K_n R) / 2, \Theta = datan2(sin, cos) and shifted to the interval - !! [0, 2pi] - !! INPUTS - !! None - !! USES - !! datatypes, local, global, pao_format, GenComms - !! AUTHOR - !! C. Xu - !! CREATION DATE - !! 02/03/2026 - !! MODIFICATION HISTORY - !! 04/03/2026 - C. Xu: use LAPACK to find rotation axis - !! 05/03/2026 - C. Xu: Correct use of lapack and rotation angle - !! 09/03/2026 - C. Xu: use atan2 for rotation angle - !! 19/03/2026 - C. Xu: initial attempt to fix 90 degree rotation - !! 24/03/2026 - C. Xu: add error checking - !! 10/04/2026 - C. Xu: fixes to determinant check, negative angles and euler matrix - !! 26/05/2026, 08/06 - C. Xu: comments and clearer error messages - !! 20/08/2026 - C. Xu: Clean debug and negate axis/angle when necessary - !! SOURCE - !! - subroutine calculate_axis_angle(w1, w2, w3, axis, angle) - use datatypes - use numbers, ONLY: pi, one - use GenComms, ONLY: cq_abort - use local, ONLY: flag_rotate_pdos_debug - implicit none - external DGEEV - - ! Allow the user to define new coordinate system - ! Rotate from simulation axes (standard Cartesian) to local axes - real(double), intent(in) :: w1(3), w2(3), w3(3) - real(double), intent(out) :: axis(3), angle - real(double) :: pos_diff(3), det, sin_angle, cos_angle, angle2 - real(double) :: basis_matrix(3,3), temp(3,3), antisym_axis(3,3) - ! Local variables - integer :: i, j, N, LDA, LDVL, LDVR, INFO, LDWORK - real(double), allocatable:: A(:,:), WR(:), WI(:), VL(:,:), VR(:,:), WORK(:), tra - real(double), parameter :: tol = 1e-10 - real(double) :: real_eigenvalue,imag_eigenvalue - N = 3 - LDA = N - LDVL = LDA - LDVR = LDA - allocate(A(N,N)) - allocate(WR(N)) - allocate(WI(N)) - allocate(VL(LDVL, N)) - allocate(VR(LDVR, N)) - allocate(WORK(5*N)) - LDWORK = 5*N - ! Define change of basis matrix; - ! Columns of new basis in terms of old basis - basis_matrix(:,1) = w1 / norm2(w1) - basis_matrix(:,2) = w2 / norm2(w2) - basis_matrix(:,3) = w3 / norm2(w3) - ! Make copy because DGEEV destroys matrix - basis_matrix = (basis_matrix) - A = basis_matrix - ! If determinant flips sign, orientation of the coordinate system has changed - det = A(1,1)*(A(2,2)*A(3,3) - A(2,3)*A(3,2)) & - - A(1,2)*(A(2,1)*A(3,3) - A(2,3)*A(3,1)) & - + A(1,3)*(A(2,1)*A(3,2) - A(2,2)*A(3,1)) - - if (det < 1.0 - tol .or. det > 1.0 + tol) & + write(*, fmt='(/4x,"alpha: ", (f10.5,1X))', advance="no") & + 180.0_double/pi*(datan2(axis(3)*sin(angle / two),cos(angle / two)) + datan2(axis(2), axis(1)) - (pi/two)) + write(*, fmt='(/4x,"beta: ", (f10.5,1X))', advance="no") & + 180.0_double/pi*2*asin(sqrt(axis(1)*axis(1) + axis(2)*axis(2)) * sin(angle / two)) + write(*, fmt='(/4x,"gamma: ", (f10.5,1X))', advance="no") & + 180.0_double/pi*(datan2(axis(3)*sin(angle/two),cos(angle/two)) - datan2(axis(2), axis(1)) + (pi/two)) + end subroutine euler_from_axisangle + + ! ----------------------------------------------------------------------------- + ! Subroutine calculate_axis_angle + ! ----------------------------------------------------------------------------- + + !!****f* ProcModule/calculate_axis_angle * + !! + !! NAME + !! calculate_axis_angle + !! USAGE + !! call calculate_axis_angle(w1, w2, w3, axis, angle) + !! PURPOSE + !! Given a set of 3 orthonormal vectors, calculate the axis of rotation and angle needed to go from the simulation + !! cell axes to this basis (order matters). + !! This subroutine forms the basis change matrix and finds its normalised eigenvector with real eigenvalue one + !! This eigenvector is the rotation axis. The rotation angle is computed as the angle satisfying + !! cos \theta = (Tr(R) - 1) /2; sin\theta = -Tr(K_n R) / 2, \Theta = datan2(sin, cos) and shifted to the interval + !! [0, 2pi] + !! INPUTS + !! None + !! USES + !! datatypes, local, global, pao_format, GenComms + !! AUTHOR + !! C. Xu + !! CREATION DATE + !! 02/03/2026 + !! MODIFICATION HISTORY + !! 04/03/2026 - C. Xu: use LAPACK to find rotation axis + !! 05/03/2026 - C. Xu: Correct use of lapack and rotation angle + !! 09/03/2026 - C. Xu: use atan2 for rotation angle + !! 19/03/2026 - C. Xu: initial attempt to fix 90 degree rotation + !! 24/03/2026 - C. Xu: add error checking + !! 10/04/2026 - C. Xu: fixes to determinant check, negative angles and euler matrix + !! 26/05/2026, 08/06 - C. Xu: comments and clearer error messages + !! 20/08/2026 - C. Xu: Clean debug and negate axis/angle when necessary + !! SOURCE + !! + subroutine calculate_axis_angle(w1, w2, w3, axis, angle) + + use datatypes + use numbers, ONLY: zero, pi, one, RD_ERR + use GenComms, ONLY: cq_abort + use local, ONLY: flag_rotate_pdos_debug + + implicit none + + external DGEEV + + ! Allow the user to define new coordinate system + ! Rotate from simulation axes (standard Cartesian) to local axes + real(double), intent(in) :: w1(3), w2(3), w3(3) + real(double), intent(out) :: axis(3), angle + real(double) :: pos_diff(3), det, sin_angle, cos_angle, angle2 + real(double) :: basis_matrix(3,3), temp(3,3), antisym_axis(3,3) + ! Local variables + integer :: i, j, N, LDA, LDVL, LDVR, INFO, LDWORK + real(double), allocatable:: A(:,:), WR(:), WI(:), VL(:,:), VR(:,:), WORK(:), tra + real(double) :: real_eigenvalue,imag_eigenvalue + + N = 3 + LDA = N + LDVL = LDA + LDVR = LDA + allocate(A(N,N)) + allocate(WR(N)) + allocate(WI(N)) + allocate(VL(LDVL, N)) + allocate(VR(LDVR, N)) + allocate(WORK(5*N)) + LDWORK = 5*N + ! Define change of basis matrix; + ! Columns of new basis in terms of old basis + basis_matrix(:,1) = w1 / norm2(w1) + basis_matrix(:,2) = w2 / norm2(w2) + basis_matrix(:,3) = w3 / norm2(w3) + ! Make copy because DGEEV destroys matrix + A = basis_matrix + ! If determinant flips sign, orientation of the coordinate system has changed + det = A(1,1)*(A(2,2)*A(3,3) - A(2,3)*A(3,2)) & + - A(1,2)*(A(2,1)*A(3,3) - A(2,3)*A(3,1)) & + + A(1,3)*(A(2,1)*A(3,2) - A(2,2)*A(3,1)) + + if (det < one - RD_ERR .or. det > one + RD_ERR) & call cq_abort("calculate_axis_angle: determinant is not +1. Make sure input axes form a right-handed basis.") - call DGEEV("N", "V", N, A, LDA, WR, WI, VL, LDVL, VR, LDVR, & - WORK, LDWORK, INFO) - if (INFO /= 0) then - write(*, fmt='(A,I0)') 'DGEEV INFO = ', INFO - call cq_abort('calculate_axis_angle: DGEEV, Eigenvalues of basis change matrix could not be found!') - end if - ! Require eigenvector with eigenvalue 1, no complex, to find axis - do i = 1, N - if (abs(WR(i) - one) < tol .and. abs(WI(i)) < tol) then - axis = VR(:, i) ! Column i of VR is the i-th right eigenvector - exit - end if - end do - if (flag_rotate_pdos_debug) then - real_eigenvalue = 0.0 - imag_eigenvalue = 0.0 - write(*, '(/2x, "Eigenvalues/Eigenvectors of change of base matrix: ")') - do j = 1, N - ! Print eigenvalue - if (abs(WI(j)) < 1e-7) then - write(*,'(/4x, A,I3,A,F12.6)') 'Eigenvalue ', j, ': ', WR(j) - else - write(*,'(/4x, A,I3,A,F12.6,SP,F12.6,A)') 'Eigenvalue ', j, ': ', WR(j), WI(j), 'i' - end if - write(*,'(/6x, 2X,A)', advance="no") 'Eigenvector:' - do i = 1, N - if (WI(j) == 0.0_double) then - write(*,'(/7x, 4X,A,I3,A,F10.5)', advance="no") 'component ', i, ': ', VR(i,j) - else if (WI(j) > 0.0d0) then - ! First complex conjugate pair: v = VR(:,j) + i*VR(:,j+1) - write(*,'(/7x, 4X,A,I3,A,F10.5,SP,F10.5,A)', advance="no") & - 'component ', i, ': ', VR(i,j), VR(i,j+1), 'i' - else - ! Second complex conjugate pair: v = VR(:,j-1) - i*VR(:,j) - write(*,'(/7x, 4XA,I3,A,F10.5,SP,F10.5,A)', advance="no") & - 'component ', i, ': ', VR(i,j-1), -VR(i,j), 'i' - end if - end do - write(*,*) - end do - end if ! End debug: print all eigenvalues/vectors - if (any(axis /= axis)) & + call DGEEV("N", "V", N, A, LDA, WR, WI, VL, LDVL, VR, LDVR, & + WORK, LDWORK, INFO) + if (INFO /= 0) then + write(*, fmt='(A,I0)') 'DGEEV INFO = ', INFO + call cq_abort('calculate_axis_angle: DGEEV, Eigenvalues of basis change matrix could not be found!') + end if + ! Require eigenvector with eigenvalue 1, no complex, to find axis + do i = 1, N + if (abs(WR(i) - one) < RD_ERR .and. abs(WI(i)) < RD_ERR) then + axis = VR(:, i) ! Column i of VR is the i-th right eigenvector + exit + end if + end do + if (flag_rotate_pdos_debug) then + real_eigenvalue = zero + imag_eigenvalue = zero + write(*, '(/2x, "Eigenvalues/Eigenvectors of change of base matrix: ")') + do j = 1, N + ! Print eigenvalue + if (abs(WI(j)) < RD_ERR) then + write(*,'(/4x, A,I3,A,F12.6)') 'Eigenvalue ', j, ': ', WR(j) + else + write(*,'(/4x, A,I3,A,F12.6,SP,F12.6,A)') 'Eigenvalue ', j, ': ', WR(j), WI(j), 'i' + end if + write(*,'(/6x, 2X,A)', advance="no") 'Eigenvector:' + do i = 1, N + if (abs(WI(j)) < RD_ERR) then + write(*,'(/7x, 4X,A,I3,A,F10.5)', advance="no") 'component ', i, ': ', VR(i,j) + else if (WI(j) > zero) then + ! First complex conjugate pair: v = VR(:,j) + i*VR(:,j+1) + write(*,'(/7x, 4X,A,I3,A,F10.5,SP,F10.5,A)', advance="no") & + 'component ', i, ': ', VR(i,j), VR(i,j+1), 'i' + else + ! Second complex conjugate pair: v = VR(:,j-1) - i*VR(:,j) + write(*,'(/7x,4X,A,I3,A,F10.5,SP,F10.5,A)', advance="no") & + 'component ', i, ': ', VR(i,j-1), -VR(i,j), 'i' + end if + end do + write(*,*) + end do + end if ! End debug: print all eigenvalues/vectors + ! Check this ! + if (any(axis /= axis)) & call cq_abort("calculate_axis_angle: NaN in rotation axis.") - if (all(abs(axis) < tol)) & + if (all(abs(axis) < RD_ERR)) & call cq_abort("calculate_axis_angle: Rotation axis was (0,0,0). Cannot perform rotation.") - tra = basis_matrix(1,1) + basis_matrix(2,2) + basis_matrix(3,3) - antisym_axis = 0.0 - axis = axis / norm2(axis) - cos_angle = (tra - 1.0) - call antisym_matrix(axis, antisym_axis) - temp = matmul(antisym_axis, basis_matrix) - sin_angle = -(temp(1,1) + temp(2,2) + temp(3,3)) - if (abs(sin_angle) < tol .and. abs(cos_angle) < tol) & - call cq_abort("calculate_axis_angle: Both arguments to datan2 are zero.") - angle = datan2(sin_angle, cos_angle) - if (angle /= angle) & - call cq_abort("calculate_axis_angle: NaN rotation angle.") - if (flag_rotate_pdos_debug) & - write(*, fmt='(/4x,"Rotation angle (rad/deg) [-pi, pi]: ",2(f10.5,1X))') & - angle, angle * 180/pi - write(*, fmt='(/4x,"Rotation axis from LAPACK : ",3(f10.5,1X))') & - axis - if (angle < 0.0) then - ! Keep angle between [0, pi] and negate both axis and angle - angle = -angle - axis = -axis !R(n, theta) = R(-n, -theta) - end if - write(*, fmt='(/4x,"Rotation angle (rad/deg) [0, pi]: ",2(f10.5,1X))') & - angle, angle * 180/pi - write(*, fmt='(/4x,"Rotation axis: ", 3(f10.5,1X))') axis - end subroutine calculate_axis_angle - subroutine construct_rodrigues(axis, angle, matrix) - use datatypes - implicit none -! We will compute the rodrigues matrix, specifically R^T(axis, -angle) - real(double), intent(in) :: axis(3), angle - real(double), intent(out) :: matrix(9) - real(double) :: K(3, 3), KT(3, 3), identity(3, 3), temp(3, 3) - real(double) :: norm_axis(3) - K = 0.0 - identity = 0.0 - identity(1, 1) = 1.0 - identity(2, 2) = 1.0 - identity(3, 3) = 1.0 - norm_axis = axis / norm2(axis) - ! Define K - call antisym_matrix(norm_axis, K) - - KT = transpose(K) - ! Rodrigues rotation matrix but with -angle - temp = identity - (sin(angle)*KT) + (1 - cos(angle))*matmul(KT, KT) - ! Use transpose to restore matrix order in maths - matrix = reshape(transpose(temp), shape=(/9/)) - end subroutine construct_rodrigues - - subroutine construct_C1(rod, C1) - use datatypes - use local, ONLY: flag_rotate_pdos_debug - implicit none - real(double), intent(in) :: rod(9) - real(double), intent(out) :: C1(3, 3) - real(double) :: identity(3,3) - integer :: i - C1(1, 1) = rod(1) - C1(1, 2) = rod(2) - C1(1, 3) = rod(3) - C1(2, 1) = rod(4) - C1(2, 2) = rod(5) - C1(2, 3) = rod(6) - C1(3, 1) = rod(7) - C1(3, 2) = rod(8) - C1(3, 3) = rod(9) - if (flag_rotate_pdos_debug) then - identity = 0.0 - identity(1, 1) = 1.0 - identity(2, 2) = 1.0 - identity(3, 3) = 1.0 - if(all(abs(matmul(C1, transpose(C1)) - identity) < 1e-7)) then - write(*, '(/2x, "C1 matrix is orthogonal: PASS")', advance='no') - else - write(*, '(/2x, "C1 matrix is orthogonal: FAIL")', advance='no') - endif - end if - end subroutine construct_C1 - subroutine construct_C2(r, C2) - use datatypes - use local, ONLY: flag_rotate_pdos_debug - implicit none - real(double), intent(in) :: r(9) - real(double), intent(out) :: C2(5, 5) - real(double) :: identity(5,5) - - C2(1, 1) = r(6)*r(8) + r(5)*r(9) - C2(1, 2) = r(6)*r(7) + r(4)*r(9) - C2(1, 3) = r(5)*r(7) + r(4)*r(8) - C2(1, 5) = r(6)*r(9)*0.5 - C2(1, 4) = r(4)*r(7) + r(6)*r(9)*0.5 -! Row 2 - C2(2, 1) = r(3)*r(8) + r(2)*r(9) - C2(2, 2) = r(3)*r(7) + r(1)*r(9) - C2(2, 3) = r(2)*r(7) + r(1)*r(8) - C2(2, 5) = r(3)*r(9)*0.5 - C2(2, 4) = r(1)*r(7) + r(3)*r(9)*0.5 -! Row 3 - C2(3, 1) = r(3)*r(5) + r(2)*r(6) - C2(3, 2) = r(3)*r(4) + r(1)*r(6) - C2(3, 3) = r(2)*r(4) + r(1)*r(5) - C2(3, 5) = r(3)*r(6)*0.5 - C2(3, 4) = r(1)*r(4) + r(3)*r(6)*0.5 - -! Row 4 - C2(4, 1) = 2*(r(2)*r(3) - r(5)*r(6)) - C2(4, 2) = 2*(r(1)*r(3) - r(4)*r(6)) - C2(4, 3) = 2*(r(1)*r(2) - r(4)*r(5)) - C2(4, 5) = 0.5*(r(3)*r(3) - r(6)*r(6)) - C2(4, 4) = r(1)*r(1) - r(4)*r(4) + 0.5*(r(3)*r(3) - r(6)*r(6)) - -! Row 5 - C2(5, 1) = (4*r(8)*r(9)) - 2*(r(2)*r(3) + r(5)*r(6)) - C2(5, 2) = (4*r(7)*r(9)) - 2*(r(1)*r(3) + r(4)*r(6)) - C2(5, 3) = (4*r(7)*r(8)) - 2*(r(1)*r(2) + r(4)*r(5)) - C2(5, 5) = r(9)*r(9) - 0.5*(r(3)*r(3) + r(6)*r(6)) - C2(5, 4) = C2(5, 5) - r(1)*r(1) - r(4)*r(4) + 2*r(7)*r(7) - end subroutine construct_C2 - function inv(A) result(Ainv) - use datatypes - use GenComms, ONLY: cq_abort - implicit none - - real(double), dimension(:, :), intent(in) :: A - real(double), dimension(size(A, 1), size(A, 2)) :: Ainv - - real(double), dimension(size(A, 1)) :: work ! work array for LAPACK - integer, dimension(size(A, 1)) :: ipiv ! pivot indices - integer :: n, info - - external DGETRF - external DGETRI - - Ainv = A - n = size(A, 1) - - ! DGETRF computes an LU factorization of a general M-by-N matrix A - call DGETRF(n, n, Ainv, n, ipiv, info) - - if (info /= 0) then - call cq_abort("inv: Matrix is numerically singular!") - end if - - ! DGETRI computes the inverse using the LU factorization by DGETRF. - call DGETRI(n, Ainv, n, ipiv, work, n, info) - - if (info /= 0) then - call cq_abort("inv: Matrix inversion failed!") - end if - end function inv - - subroutine construct_Ul(angmom, Al, Cl, Ul) - use datatypes - implicit none - - integer, intent(in) :: angmom - real(double), intent(in) :: Al(2*angmom + 1, 2*angmom + 1), Cl(2*angmom + 1, 2*angmom + 1) - real(double), intent(out) :: Ul(:,:) - real(double) :: Al_inv(2*angmom + 1, 2*angmom + 1) - Al_inv = inv(Al) - Ul = matmul(Al_inv, matmul(Cl, Al)) - - end subroutine construct_Ul - ! ----------------------------------------------------------------------------- - ! Subroutine rotate_coefficients - ! ----------------------------------------------------------------------------- - - !!****f* ProcModule/rotate_coefficients * - !! - !! NAME - !! rotate_coefficients - !! USAGE - !! call rotate_coefficients - !! PURPOSE - !! Apply rotation matrices to wavefunction coefficients. Formally, for each l, each 2l + 1 coefficients - !! corresponding to a (m, zeta) are multiplied by the corresponding U^{(l)} rotation matrix. - !! As different rotation modes use different arrays and lookups, these are taken care of at the start of the - !! subroutine - !! The loop over all values is pasted and modified from process_pdos. - !! INPUTS - !! None - !! USES - !! datatypes, local, global, pao_format, GenComms - !! AUTHOR - !! C. Xu - !! CREATION DATE - !! 02/03/2026 - !! MODIFICATION HISTORY - !! 14/04/2026 - C. Xu: loop over active bands only - !! 03/05/2026 - C. Xu: correct g_atom_lookup - !! 12/05/2026 - C. Xu: correct de/allocation behaviour. Edit desc - !! 27/08/2026 - C. Xu: Similarity transform on U2 to account for orbital basis - !! SOURCE - !! - subroutine rotate_coefficients - use datatypes - use local, ONLY: n_bands_total, nkp, n_atoms_pDOS, evec_coeff, scaled_evec_coeff, & - pDOS_atom_index, band_full_to_active, rotate_pdos_natoms, find_neighbours, & - flag_rotate_pdos_mode, rotate_pdos_atoms, U1, U2 - use global_module, ONLY: nspin, species_glob - use pao_format, ONLY: pao - use GenComms, ONLY: cq_abort - - implicit none - ! Local variables - integer :: i_atom, i_spec, i_band, i_band_c, i_kp, i_spin, g_atom, i - integer :: i_l, i_z, nzeta, norbs, sf_offset, rotate_counter - integer, dimension(:), allocatable :: g_atom_lookup - ! Currently input axes depends on n_atoms_pDOS - if (flag_rotate_pdos_mode == 0) then - rotate_counter = n_atoms_pDOS - else - rotate_counter = rotate_pdos_natoms - end if - if (.not. allocated(g_atom_lookup)) allocate(g_atom_lookup(rotate_counter)) - ! Define the lookups correctly for each mode - if (flag_rotate_pdos_mode == 2) then - ! Atom counter should be the order the user input pDOSNeighbours block - g_atom_lookup = find_neighbours(1, :) - else if (flag_rotate_pdos_mode == 1 .or. flag_rotate_pdos_mode == 3) then - ! Atom counter should be the order the user input other blocks - g_atom_lookup = rotate_pdos_atoms - else - g_atom_lookup = pDOS_atom_index - end if - do i = 1, size(g_atom_lookup) + tra = basis_matrix(1,1) + basis_matrix(2,2) + basis_matrix(3,3) + antisym_axis = zero + axis = axis / norm2(axis) + cos_angle = (tra - one) + call antisym_matrix(axis, antisym_axis) + temp = matmul(antisym_axis, basis_matrix) + sin_angle = -(temp(1,1) + temp(2,2) + temp(3,3)) + if (abs(sin_angle) < RD_ERR .and. abs(cos_angle) < RD_ERR) & + call cq_abort("calculate_axis_angle: Both arguments to datan2 are zero.") + angle = datan2(sin_angle, cos_angle) + if (angle /= angle) & + call cq_abort("calculate_axis_angle: NaN rotation angle.") + if (flag_rotate_pdos_debug) & + write(*, fmt='(/4x,"Rotation angle (rad/deg) [-pi, pi]: ",2(f10.5,1X))') & + angle, angle * 180/pi + write(*, fmt='(/4x,"Rotation axis from LAPACK : ",3(f10.5,1X))') & + axis + if (angle < zero) then + ! Keep angle between [0, pi] and negate both axis and angle + angle = -angle + axis = -axis !R(n, theta) = R(-n, -theta) + end if + write(*, fmt='(/4x,"Rotation angle (rad/deg) [0, pi]: ",2(f10.5,1X))') & + angle, angle * 180/pi + write(*, fmt='(/4x,"Rotation axis: ", 3(f10.5,1X))') axis + end subroutine calculate_axis_angle + + subroutine construct_rodrigues(axis, angle, matrix) + + use datatypes + use numbers, ONLY: zero, one + + implicit none + + ! We will compute the rodrigues matrix, specifically R^T(axis, -angle) + real(double), intent(in) :: axis(3), angle + real(double), intent(out) :: matrix(9) + real(double) :: K(3, 3), KT(3, 3), identity(3, 3), temp(3, 3) + real(double) :: norm_axis(3) + + K = zero + identity = zero + identity(1, 1) = one + identity(2, 2) = one + identity(3, 3) = one + norm_axis = axis / norm2(axis) + ! Define K + call antisym_matrix(norm_axis, K) + KT = transpose(K) + ! Rodrigues rotation matrix but with -angle + temp = identity - (sin(angle)*KT) + (one - cos(angle))*matmul(KT, KT) + ! Use transpose to restore matrix order in maths + matrix = reshape(transpose(temp), shape=(/9/)) + end subroutine construct_rodrigues + + subroutine construct_C1(rod, C1) + + use datatypes + use numbers, ONLY: zero, one, RD_ERR + use local, ONLY: flag_rotate_pdos_debug + + implicit none + + real(double), intent(in) :: rod(9) + real(double), intent(out) :: C1(3, 3) + real(double) :: identity(3,3) + integer :: i + C1(1, 1) = rod(1) + C1(1, 2) = rod(2) + C1(1, 3) = rod(3) + C1(2, 1) = rod(4) + C1(2, 2) = rod(5) + C1(2, 3) = rod(6) + C1(3, 1) = rod(7) + C1(3, 2) = rod(8) + C1(3, 3) = rod(9) + if (flag_rotate_pdos_debug) then + identity = zero + identity(1, 1) = one + identity(2, 2) = one + identity(3, 3) = one + if(all(abs(matmul(C1, transpose(C1)) - identity) < RD_ERR)) then + write(*, '(/2x, "C1 matrix is orthogonal: PASS")', advance='no') + else + write(*, '(/2x, "C1 matrix is orthogonal: FAIL")', advance='no') + endif + end if + end subroutine construct_C1 + + subroutine construct_C2(r, C2) + + use datatypes + use numbers, ONLY: half, two, four + use local, ONLY: flag_rotate_pdos_debug + + implicit none + + real(double), intent(in) :: r(9) + real(double), intent(out) :: C2(5, 5) + real(double) :: identity(5,5) + + C2(1, 1) = r(6)*r(8) + r(5)*r(9) + C2(1, 2) = r(6)*r(7) + r(4)*r(9) + C2(1, 3) = r(5)*r(7) + r(4)*r(8) + C2(1, 5) = r(6)*r(9)*half + C2(1, 4) = r(4)*r(7) + r(6)*r(9)*half + ! Row 2 + C2(2, 1) = r(3)*r(8) + r(2)*r(9) + C2(2, 2) = r(3)*r(7) + r(1)*r(9) + C2(2, 3) = r(2)*r(7) + r(1)*r(8) + C2(2, 5) = r(3)*r(9)*half + C2(2, 4) = r(1)*r(7) + r(3)*r(9)*half + ! Row 3 + C2(3, 1) = r(3)*r(5) + r(2)*r(6) + C2(3, 2) = r(3)*r(4) + r(1)*r(6) + C2(3, 3) = r(2)*r(4) + r(1)*r(5) + C2(3, 5) = r(3)*r(6)*half + C2(3, 4) = r(1)*r(4) + r(3)*r(6)*half + + ! Row 4 + C2(4, 1) = two*(r(2)*r(3) - r(5)*r(6)) + C2(4, 2) = two*(r(1)*r(3) - r(4)*r(6)) + C2(4, 3) = two*(r(1)*r(2) - r(4)*r(5)) + C2(4, 5) = half*(r(3)*r(3) - r(6)*r(6)) + C2(4, 4) = r(1)*r(1) - r(4)*r(4) + half*(r(3)*r(3) - r(6)*r(6)) + + ! Row 5 + C2(5, 1) = (four*r(8)*r(9)) - two*(r(2)*r(3) + r(5)*r(6)) + C2(5, 2) = (four*r(7)*r(9)) - two*(r(1)*r(3) + r(4)*r(6)) + C2(5, 3) = (four*r(7)*r(8)) - two*(r(1)*r(2) + r(4)*r(5)) + C2(5, 5) = r(9)*r(9) - half*(r(3)*r(3) + r(6)*r(6)) + C2(5, 4) = C2(5, 5) - r(1)*r(1) - r(4)*r(4) + two*r(7)*r(7) + end subroutine construct_C2 + + function inv(A) result(Ainv) + + use datatypes + use GenComms, ONLY: cq_abort + + implicit none + + real(double), dimension(:, :), intent(in) :: A + real(double), dimension(size(A, 1), size(A, 2)) :: Ainv + + real(double), dimension(size(A, 1)) :: work ! work array for LAPACK + integer, dimension(size(A, 1)) :: ipiv ! pivot indices + integer :: n, info + + external DGETRF + external DGETRI + + Ainv = A + n = size(A, 1) + + ! DGETRF computes an LU factorization of a general M-by-N matrix A + call DGETRF(n, n, Ainv, n, ipiv, info) + + if (info /= 0) then + call cq_abort("inv: Matrix is numerically singular!") + end if + + ! DGETRI computes the inverse using the LU factorization by DGETRF. + call DGETRI(n, Ainv, n, ipiv, work, n, info) + + if (info /= 0) then + call cq_abort("inv: Matrix inversion failed!") + end if + end function inv + + subroutine construct_Ul(angmom, Al, Cl, Ul) + + use datatypes + + implicit none + + integer, intent(in) :: angmom + real(double), intent(in) :: Al(2*angmom + 1, 2*angmom + 1), Cl(2*angmom + 1, 2*angmom + 1) + real(double), intent(out) :: Ul(:,:) + real(double) :: Al_inv(2*angmom + 1, 2*angmom + 1) + + Al_inv = inv(Al) + Ul = matmul(Al_inv, matmul(Cl, Al)) + + end subroutine construct_Ul + + ! ----------------------------------------------------------------------------- + ! Subroutine rotate_coefficients + ! ----------------------------------------------------------------------------- + + !!****f* ProcModule/rotate_coefficients * + !! + !! NAME + !! rotate_coefficients + !! USAGE + !! call rotate_coefficients + !! PURPOSE + !! Apply rotation matrices to wavefunction coefficients. Formally, for each l, each 2l + 1 coefficients + !! corresponding to a (m, zeta) are multiplied by the corresponding U^{(l)} rotation matrix. + !! As different rotation modes use different arrays and lookups, these are taken care of at the start of the + !! subroutine + !! The loop over all values is pasted and modified from process_pdos. + !! INPUTS + !! None + !! USES + !! datatypes, local, global, pao_format, GenComms + !! AUTHOR + !! C. Xu + !! CREATION DATE + !! 02/03/2026 + !! MODIFICATION HISTORY + !! 14/04/2026 - C. Xu: loop over active bands only + !! 03/05/2026 - C. Xu: correct g_atom_lookup + !! 12/05/2026 - C. Xu: correct de/allocation behaviour. Edit desc + !! 27/08/2026 - C. Xu: Similarity transform on U2 to account for orbital basis + !! SOURCE + !! + subroutine rotate_coefficients + + use datatypes + use local, ONLY: n_bands_total, nkp, n_atoms_pDOS, evec_coeff, scaled_evec_coeff, & + pDOS_atom_index, band_full_to_active, rotate_pdos_natoms, find_neighbours, & + flag_rotate_pdos_mode, rotate_pdos_atoms, U1, U2 + use global_module, ONLY: nspin, species_glob + use pao_format, ONLY: pao + use GenComms, ONLY: cq_abort + + implicit none + + ! Local variables + integer :: i_atom, i_spec, i_band, i_band_c, i_kp, i_spin, g_atom, i + integer :: i_l, i_z, nzeta, norbs, sf_offset, rotate_counter + integer, dimension(:), allocatable :: g_atom_lookup + + ! Currently input axes depends on n_atoms_pDOS + if (flag_rotate_pdos_mode == 0) then + rotate_counter = n_atoms_pDOS + else + rotate_counter = rotate_pdos_natoms + end if + if (.not. allocated(g_atom_lookup)) allocate(g_atom_lookup(rotate_counter)) + ! Define the lookups correctly for each mode + if (flag_rotate_pdos_mode == 2) then + ! Atom counter should be the order the user input pDOSNeighbours block + g_atom_lookup = find_neighbours(1, :) + else if (flag_rotate_pdos_mode == 1 .or. flag_rotate_pdos_mode == 3) then + ! Atom counter should be the order the user input other blocks + g_atom_lookup = rotate_pdos_atoms + else + g_atom_lookup = pDOS_atom_index + end if + do i = 1, size(g_atom_lookup) ! As in CONQUEST, m = -2 = -xy NOT xy, need to negate first row and column, except U2(1,1) - U2(1,2:,i) = -U2(1,2:,i) - U2(2:,1,i) = -U2(2:,1,i) - end do - - do i_spin = 1, nspin - do i_kp = 1, nkp - do i_band = 1, n_bands_total - i_band_c = band_full_to_active(i_band) - if(i_band_c > 0) then - do i_atom = 1, rotate_counter - ! Get global atom number from input - g_atom = g_atom_lookup(i_atom) - i_spec = species_glob(g_atom) - sf_offset = 0 - ! Include l = 0 to correctly calculate offset - do i_l = 0, pao(i_spec)%greatest_angmom - nzeta = pao(i_spec)%angmom(i_l)%n_zeta_in_angmom - norbs = 2*i_l + 1 - do i_z = 1, nzeta - select case(i_l) - case(1) - evec_coeff(sf_offset+1:sf_offset+norbs, g_atom, i_band_c, i_kp, i_spin) = & - matmul(U1(:,:,i_atom), evec_coeff(sf_offset+1:sf_offset+norbs, g_atom, i_band_c, i_kp, i_spin)) - scaled_evec_coeff(sf_offset+1:sf_offset+norbs, g_atom, i_band_c, i_kp, i_spin) = & - matmul(U1(:,:,i_atom), scaled_evec_coeff(sf_offset+1:sf_offset+norbs, g_atom, i_band_c, i_kp, i_spin)) - case(2) - evec_coeff(sf_offset+1:sf_offset+norbs, g_atom, i_band_c, i_kp, i_spin) = & - matmul(U2(:,:,i_atom), evec_coeff(sf_offset+1:sf_offset+norbs, g_atom, i_band_c, i_kp, i_spin)) - scaled_evec_coeff(sf_offset+1:sf_offset+norbs, g_atom, i_band_c, i_kp, i_spin) = & - matmul(U2(:,:,i_atom), scaled_evec_coeff(sf_offset+1:sf_offset+norbs, g_atom, i_band_c, i_kp, i_spin)) - end select - sf_offset = sf_offset + norbs - end do ! i_z - end do ! i_l - end do ! i_atom - end if ! if band is active - end do ! i_band - end do ! i_kp - end do ! i_spin - deallocate(g_atom_lookup) - end subroutine rotate_coefficients - - ! ----------------------------------------------------------------------------- - ! Subroutine nearest_neighbours - ! ----------------------------------------------------------------------------- - - !!****f* ProcModule/nearest_neighbours * - !! - !! NAME - !! nearest_neighbours - Find nearest neighbours and their bond vectors - !! USAGE - !! nearest_neighbours(atomno, bond) - !! PURPOSE - !! Evaluates the nearest neighbours using periodic boundary conditions - !! and returns their bond vectors - !! - !! INPUTS - !! integer, intent(in) :: atomno !atom to find neighbours of - !! real(double), intent(out), allocatable :: bond(:,:) ! array to hold bond - !! information regarding an atom - !! USES - !! datatypes, dimens, local, global, GenComms - !! AUTHOR - !! C. Xu - !! CREATION DATE - !! 15/04/2026 - !! MODIFICATION HISTORY - !! 12/05/2026 - C. Xu: correct de/allocation behaviour. Edit desc - !! SOURCE - !! - subroutine nearest_neighbours(atomno, bond) - use datatypes - use dimens, ONLY: r_super_x, r_super_y, r_super_z - use local, ONLY: find_neighbours, nghbr_arr, flag_rotate_pdos_debug - use global_module, ONLY: ni_in_cell, atom_coord, species_glob - use GenComms, ONLY: cq_abort - use species_module, ONLY: species_label - - implicit none - - integer, intent(in) :: atomno - real(double), intent(out), allocatable :: bond(:,:) - - - - - ! Local variables - real(double), parameter :: err = 1e-5 - real(double) :: distances(1:ni_in_cell), temp(1:ni_in_cell) - real(double) ::cx, cy, cz, dx, dy, dz, atomno_pos(3) - integer :: i, j, min_idx, ibond_max_len, idx_direction - ! Reset allocatables as this subroutine is called in a loop - ! and populated with new data every time - ! and populated with new data every time - if (allocated(nghbr_arr)) deallocate(nghbr_arr) - if (allocated(bond)) deallocate(bond) - - cx = atom_coord(1, atomno) - cy = atom_coord(2, atomno) - cz = atom_coord(3, atomno) - ! ACCOUNT FOR PERIODICITY - do i = 1, ni_in_cell - - if (i == atomno) then - distances(i) = huge(1.0d0) - cycle - end if - - dx = atom_coord(1,i) - cx - dy = atom_coord(2,i) - cy - dz = atom_coord(3,i) - cz - - dx = dx - r_super_x * nint(dx / r_super_x) - dy = dy - r_super_y * nint(dy / r_super_y) - dz = dz - r_super_z * nint(dz / r_super_z) - - distances(i) = dx*dx + dy*dy + dz*dz - end do - - temp = distances - ! Eliminate zeros - self-distance - where (abs(temp) < err) temp = huge(1.0) - select case (find_neighbours(2, findloc(find_neighbours(1,:), atomno, dim=1))) - case (0) - allocate(nghbr_arr(4)) - case (1) - allocate(nghbr_arr(6)) - case default - call cq_abort("nearest_neighbours: Did not correctly allocate nghbr_arr.") - end select - do j = 1, size(nghbr_arr) - min_idx = minloc(temp, 1) - nghbr_arr(j) = min_idx - temp(min_idx) = huge(1.0) - end do - - allocate(bond(3, size(nghbr_arr))) - - atomno_pos = (/cx, cy, cz/) - do i = 1, size(nghbr_arr) - bond(1, i) = atom_coord(1, nghbr_arr(i)) - atomno_pos(1) - bond(2, i) = atom_coord(2, nghbr_arr(i)) - atomno_pos(2) - bond(3, i) = atom_coord(3, nghbr_arr(i)) - atomno_pos(3) - - bond(1, i) = bond(1, i) - r_super_x * nint(bond(1, i) / r_super_x) - bond(2, i) = bond(2, i) - r_super_y * nint(bond(2, i) / r_super_y) - bond(3, i) = bond(3, i) - r_super_z * nint(bond(3, i) / r_super_z) - - end do - - if (flag_rotate_pdos_debug) then - write(*, fmt='(/2x,"Unit cell dimensions in Bohr (x,y,z)", 3(f10.5))', advance="no") & + U2(1,2:,i) = -U2(1,2:,i) + U2(2:,1,i) = -U2(2:,1,i) + end do + + do i_spin = 1, nspin + do i_kp = 1, nkp + do i_band = 1, n_bands_total + i_band_c = band_full_to_active(i_band) + if(i_band_c > 0) then + do i_atom = 1, rotate_counter + ! Get global atom number from input + g_atom = g_atom_lookup(i_atom) + i_spec = species_glob(g_atom) + sf_offset = 0 + ! Include l = 0 to correctly calculate offset + do i_l = 0, pao(i_spec)%greatest_angmom + nzeta = pao(i_spec)%angmom(i_l)%n_zeta_in_angmom + norbs = 2*i_l + 1 + do i_z = 1, nzeta + select case(i_l) + case(1) + evec_coeff(sf_offset+1:sf_offset+norbs, g_atom, i_band_c, i_kp, i_spin) = & + matmul(U1(:,:,i_atom), evec_coeff(sf_offset+1:sf_offset+norbs, g_atom, & + i_band_c, i_kp, i_spin)) + scaled_evec_coeff(sf_offset+1:sf_offset+norbs, g_atom, i_band_c, i_kp, i_spin) = & + matmul(U1(:,:,i_atom), scaled_evec_coeff(sf_offset+1:sf_offset+norbs, g_atom, & + i_band_c, i_kp, i_spin)) + case(2) + evec_coeff(sf_offset+1:sf_offset+norbs, g_atom, i_band_c, i_kp, i_spin) = & + matmul(U2(:,:,i_atom), evec_coeff(sf_offset+1:sf_offset+norbs, g_atom, & + i_band_c, i_kp, i_spin)) + scaled_evec_coeff(sf_offset+1:sf_offset+norbs, g_atom, i_band_c, i_kp, i_spin) = & + matmul(U2(:,:,i_atom), scaled_evec_coeff(sf_offset+1:sf_offset+norbs, g_atom, & + i_band_c, i_kp, i_spin)) + end select + sf_offset = sf_offset + norbs + end do ! i_z + end do ! i_l + end do ! i_atom + end if ! if band is active + end do ! i_band + end do ! i_kp + end do ! i_spin + deallocate(g_atom_lookup) + end subroutine rotate_coefficients + + ! ----------------------------------------------------------------------------- + ! Subroutine nearest_neighbours + ! ----------------------------------------------------------------------------- + + !!****f* ProcModule/nearest_neighbours * + !! + !! NAME + !! nearest_neighbours - Find nearest neighbours and their bond vectors + !! USAGE + !! nearest_neighbours(atomno, bond) + !! PURPOSE + !! Evaluates the nearest neighbours using periodic boundary conditions + !! and returns their bond vectors + !! + !! INPUTS + !! integer, intent(in) :: atomno !atom to find neighbours of + !! real(double), intent(out), allocatable :: bond(:,:) ! array to hold bond + !! information regarding an atom + !! USES + !! datatypes, dimens, local, global, GenComms + !! AUTHOR + !! C. Xu + !! CREATION DATE + !! 15/04/2026 + !! MODIFICATION HISTORY + !! 12/05/2026 - C. Xu: correct de/allocation behaviour. Edit desc + !! SOURCE + !! + subroutine nearest_neighbours(atomno, bond) + + use datatypes + use numbers, only: BIG, RD_ERR + use dimens, ONLY: r_super_x, r_super_y, r_super_z + use local, ONLY: find_neighbours, nghbr_arr, flag_rotate_pdos_debug + use global_module, ONLY: ni_in_cell, atom_coord, species_glob + use GenComms, ONLY: cq_abort + use species_module, ONLY: species_label + + implicit none + + integer, intent(in) :: atomno + real(double), intent(out), allocatable :: bond(:,:) + + ! Local variables + real(double) :: distances(1:ni_in_cell), temp(1:ni_in_cell) + real(double) ::cx, cy, cz, dx, dy, dz, atomno_pos(3) + integer :: i, j, min_idx, ibond_max_len, idx_direction + + ! Reset allocatables as this subroutine is called in a loop + ! and populated with new data every time + ! and populated with new data every time + if (allocated(nghbr_arr)) deallocate(nghbr_arr) + if (allocated(bond)) deallocate(bond) + + cx = atom_coord(1, atomno) + cy = atom_coord(2, atomno) + cz = atom_coord(3, atomno) + ! ACCOUNT FOR PERIODICITY + do i = 1, ni_in_cell + + if (i == atomno) then + distances(i) = BIG + cycle + end if + + dx = atom_coord(1,i) - cx + dy = atom_coord(2,i) - cy + dz = atom_coord(3,i) - cz + + dx = dx - r_super_x * nint(dx / r_super_x) + dy = dy - r_super_y * nint(dy / r_super_y) + dz = dz - r_super_z * nint(dz / r_super_z) + + distances(i) = dx*dx + dy*dy + dz*dz + end do + temp = distances + ! Eliminate zeros - self-distance + where (abs(temp) < RD_ERR) temp = BIG + select case (find_neighbours(2, findloc(find_neighbours(1,:), atomno, dim=1))) + case (0) + allocate(nghbr_arr(4)) + case (1) + allocate(nghbr_arr(6)) + case default + call cq_abort("nearest_neighbours: Did not correctly allocate nghbr_arr.") + end select + do j = 1, size(nghbr_arr) + min_idx = minloc(temp, 1) + nghbr_arr(j) = min_idx + temp(min_idx) = BIG + end do + + allocate(bond(3, size(nghbr_arr))) + + atomno_pos = (/cx, cy, cz/) + do i = 1, size(nghbr_arr) + bond(1, i) = atom_coord(1, nghbr_arr(i)) - atomno_pos(1) + bond(2, i) = atom_coord(2, nghbr_arr(i)) - atomno_pos(2) + bond(3, i) = atom_coord(3, nghbr_arr(i)) - atomno_pos(3) + + bond(1, i) = bond(1, i) - r_super_x * nint(bond(1, i) / r_super_x) + bond(2, i) = bond(2, i) - r_super_y * nint(bond(2, i) / r_super_y) + bond(3, i) = bond(3, i) - r_super_z * nint(bond(3, i) / r_super_z) + + end do + + if (flag_rotate_pdos_debug) then + write(*, fmt='(/2x,"Unit cell dimensions in Bohr (x,y,z)", 3(f10.5))', advance="no") & r_super_x, r_super_y, r_super_z - write(*, fmt='(/2x,"Outputting cell-periodic distances (Bohr)", & + write(*, fmt='(/2x,"Outputting cell-periodic distances (Bohr)", & "of all atoms relative to atom", 1X,I0)', advance="no") atomno - do i = 1, ni_in_cell - if (i /= atomno) then - write(*, fmt='(/4x, I0, A, f10.5)', advance="no") & + do i = 1, ni_in_cell + if (i /= atomno) then + write(*, fmt='(/4x, I0, A, f10.5)', advance="no") & i, " " // species_label(species_glob(i)), sqrt(distances(i)) - end if - end do - write(*, fmt='(/2x,"Outputting cell-periodic unnormalised bond lengths")') + end if + end do + write(*, fmt='(/2x,"Outputting cell-periodic unnormalised bond lengths")') - do i = 1, size(nghbr_arr) - write(*, fmt='(/4x, I0, A, 3f10.5)', advance="no") & + do i = 1, size(nghbr_arr) + write(*, fmt='(/4x, I0, A, 3f10.5)', advance="no") & nghbr_arr(i), " " // species_label(species_glob(nghbr_arr(i))), bond(:,i) - end do - - end if - - end subroutine - ! ----------------------------------------------------------------------------- - ! Subroutine axes_from_nn - ! ----------------------------------------------------------------------------- - - !!****f* ProcModule/axes_from_nn * - !! - !! NAME - !! axes_from_nn - Gets rotation axes from bond vectors - !! USAGE - !! axes_from_nn(atomno, bond) - !! PURPOSE - !! This subroutine will use computed nearest_neighbours - !! This subroutine will use computed nearest_neighbours - !! and specified local geometry to construct a new set of local - !! axes to rotate the pDOS into - !! - !! For square-planar geometry: - !! - !! For square-planar geometry: - !! Longest bond is chosen as x_hat. - !! Bonds which are closest to orthogonality are - !! projected onto a plane defined by x_hat - !! y_hat is chosen from the projected direction - !! which has the minimal difference to its unprojected bond - !! z_hat is computed as the cross-product from these 2 directions, and thus - !! z_hat is computed as the cross-product from these 2 directions, and thus - !! point perpendicular to the planar geometry - !! For octahedra: choose z_hat as longest or shortest bond - !! Use this direction to define a plane - project the 4 atoms onto it - !! Call y_hat the bond which is the closest to the plane - !! Find x_hat = z cross y - !! INPUTS - !! integer, intent(in) :: atomno - atom to find neighbours of - !! real(double), intent(in) :: bond(:,:) - array to hold bond vectors - !! USES - !! datatypes, local, global, GenComms - !! AUTHOR - !! C. Xu - !! CREATION DATE - !! 15/04/2026 - !! MODIFICATION HISTORY - !! 20/04/2026 C. Xu - allow user to specify neighbour as secondary direction - !! SOURCE - !! - subroutine axes_from_nn(atomno, bond) - use datatypes - use local, ONLY: find_neighbours, nghbr_arr, pdos_ax, pdos_ay, pdos_az - use global_module, ONLY: ni_in_cell, atom_coord - use GenComms, ONLY: cq_abort - implicit none - ! This subroutine will use computed nearest_neighbours - ! and specified local geometry to construct a new set of local - ! axes to rotate the pDOS into - ! - ! For square-planar geometry: - ! Longest bond is chosen as x_hat. - ! Bonds which are closest to orthogonality are - ! projected onto a plane defined by x_hat - ! y_hat is chosen from the projected direction - ! which has the minimal difference to its unprojected bond - ! z_hat is computed as the cross-product from these 2 directions, and thus - ! point perpendicular to the planar geometry - ! For octahedra: choose z_hat as longest or shortest bond - ! Use this direction to define a plane - project the 4 atoms onto it - ! Call y_hat the bond which is the closest to the plane - ! Find x_hat by x_hat = z cross y - integer, intent(in) :: atomno - real(double), intent(in) :: bond(:,:) - ! Local variables - real(double) ::atomno_pos(3), plane_normal(3), proj_vector(3), cos_angle - real(double), allocatable :: dots(:), bond_lengths(:), norm_bond(:,:) - integer :: i, ibond_principal, ibond_sec, idx_direction, minormax, find_atomno - - allocate(dots(size(bond(1,:)))) - allocate(norm_bond(3, size(bond(1,:)))) - allocate(bond_lengths(size(bond(1,:)))) - do i = 1, size(bond_lengths) - bond_lengths(i) = sqrt(dot_product(bond(:, i), bond(:, i))) - if (bond_lengths(i) > 1e-10) then - norm_bond(:, i) = bond(:, i) / bond_lengths(i) - else - call cq_abort("axes_from_nn: Zero-length bond") - end if - end do - ! find_atomno: index of user supplied central atom, atomno - find_atomno = findloc(find_neighbours(1,:), atomno, dim=1) - minormax = find_neighbours(3, find_atomno) - if (minormax < 0) then - ibond_principal = minloc(bond_lengths,1) ! gets neighbour with min bond length - else if (minormax .eq. 0) then - ibond_principal = maxloc(bond_lengths,1) ! gets neighbour with maximum bond length - else - ibond_principal = findloc(nghbr_arr, minormax, dim=1) - end if - if (ibond_principal == 0) & - call cq_abort('axes_from_nn: principal neighbour not found for atom ', find_atomno) - - write(*,fmt='(/2x,"Principal neighbour found: atom ", (I0,1x))') nghbr_arr(ibond_principal) - dots = matmul(transpose(norm_bond), norm_bond(:, ibond_principal)) ! dot prod with chosen normal - ! select index of closest to perpendicular bond - there may be two - idx_direction = minloc(abs(dots), 1) - if (find_neighbours(4, find_atomno) == 0) then - ! If 0, we choose second direction by closest projection - call project_onto_plane(norm_bond(:, ibond_principal),norm_bond(:, idx_direction), proj_vector) - else - ! This entry > 0 -> corresponds to neighbour - ibond_sec = findloc(nghbr_arr, find_neighbours(4, find_atomno), dim=1) - ! Handle when input neighbour is not found - if (ibond_sec == 0) & + end do + + end if + + end subroutine nearest_neighbours + + ! ----------------------------------------------------------------------------- + ! Subroutine axes_from_nn + ! ----------------------------------------------------------------------------- + + !!****f* ProcModule/axes_from_nn * + !! + !! NAME + !! axes_from_nn - Gets rotation axes from bond vectors + !! USAGE + !! axes_from_nn(atomno, bond) + !! PURPOSE + !! This subroutine will use computed nearest_neighbours + !! This subroutine will use computed nearest_neighbours + !! and specified local geometry to construct a new set of local + !! axes to rotate the pDOS into + !! + !! For square-planar geometry: + !! + !! For square-planar geometry: + !! Longest bond is chosen as x_hat. + !! Bonds which are closest to orthogonality are + !! projected onto a plane defined by x_hat + !! y_hat is chosen from the projected direction + !! which has the minimal difference to its unprojected bond + !! z_hat is computed as the cross-product from these 2 directions, and thus + !! z_hat is computed as the cross-product from these 2 directions, and thus + !! point perpendicular to the planar geometry + !! For octahedra: choose z_hat as longest or shortest bond + !! Use this direction to define a plane - project the 4 atoms onto it + !! Call y_hat the bond which is the closest to the plane + !! Find x_hat = z cross y + !! INPUTS + !! integer, intent(in) :: atomno - atom to find neighbours of + !! real(double), intent(in) :: bond(:,:) - array to hold bond vectors + !! USES + !! datatypes, local, global, GenComms + !! AUTHOR + !! C. Xu + !! CREATION DATE + !! 15/04/2026 + !! MODIFICATION HISTORY + !! 20/04/2026 C. Xu - allow user to specify neighbour as secondary direction + !! SOURCE + !! + subroutine axes_from_nn(atomno, bond) + + use datatypes + use local, ONLY: find_neighbours, nghbr_arr, pdos_ax, pdos_ay, pdos_az + use global_module, ONLY: ni_in_cell, atom_coord + use GenComms, ONLY: cq_abort + + implicit none + + ! This subroutine will use computed nearest_neighbours + ! and specified local geometry to construct a new set of local + ! axes to rotate the pDOS into + ! + ! For square-planar geometry: + ! Longest bond is chosen as x_hat. + ! Bonds which are closest to orthogonality are + ! projected onto a plane defined by x_hat + ! y_hat is chosen from the projected direction + ! which has the minimal difference to its unprojected bond + ! z_hat is computed as the cross-product from these 2 directions, and thus + ! point perpendicular to the planar geometry + ! For octahedra: choose z_hat as longest or shortest bond + ! Use this direction to define a plane - project the 4 atoms onto it + ! Call y_hat the bond which is the closest to the plane + ! Find x_hat by x_hat = z cross y + + integer, intent(in) :: atomno + real(double), intent(in) :: bond(:,:) + + ! Local variables + real(double) ::atomno_pos(3), plane_normal(3), proj_vector(3), cos_angle + real(double), allocatable :: dots(:), bond_lengths(:), norm_bond(:,:) + integer :: i, ibond_principal, ibond_sec, idx_direction, minormax, find_atomno + + allocate(dots(size(bond(1,:)))) + allocate(norm_bond(3, size(bond(1,:)))) + allocate(bond_lengths(size(bond(1,:)))) + do i = 1, size(bond_lengths) + bond_lengths(i) = sqrt(dot_product(bond(:, i), bond(:, i))) + if (bond_lengths(i) > 1e-10) then + norm_bond(:, i) = bond(:, i) / bond_lengths(i) + else + call cq_abort("axes_from_nn: Zero-length bond") + end if + end do + ! find_atomno: index of user supplied central atom, atomno + find_atomno = findloc(find_neighbours(1,:), atomno, dim=1) + minormax = find_neighbours(3, find_atomno) + if (minormax < 0) then + ibond_principal = minloc(bond_lengths,1) ! gets neighbour with min bond length + else if (minormax .eq. 0) then + ibond_principal = maxloc(bond_lengths,1) ! gets neighbour with maximum bond length + else + ibond_principal = findloc(nghbr_arr, minormax, dim=1) + end if + if (ibond_principal == 0) & + call cq_abort('axes_from_nn: principal neighbour not found for atom ', find_atomno) + + write(*,fmt='(/2x,"Principal neighbour found: atom ", (I0,1x))') nghbr_arr(ibond_principal) + dots = matmul(transpose(norm_bond), norm_bond(:, ibond_principal)) ! dot prod with chosen normal + ! select index of closest to perpendicular bond - there may be two + idx_direction = minloc(abs(dots), 1) + if (find_neighbours(4, find_atomno) == 0) then + ! If 0, we choose second direction by closest projection + call project_onto_plane(norm_bond(:, ibond_principal),norm_bond(:, idx_direction), proj_vector) + else + ! This entry > 0 -> corresponds to neighbour + ibond_sec = findloc(nghbr_arr, find_neighbours(4, find_atomno), dim=1) + ! Handle when input neighbour is not found + if (ibond_sec == 0) & call cq_abort('axes_from_nn: second neighbour not found for atom ', find_atomno) - if (abs(dots(ibond_sec)) > 0.5) & + if (abs(dots(ibond_sec)) > 0.5) & print *, "WARNING: Chosen secondary neighbour appears to not be very perpendicular to principal." - write(*,fmt='(/2x,"Secondary neighbour found: atom ", (I0,1x))') nghbr_arr(ibond_sec) - - call project_onto_plane(norm_bond(:, ibond_principal),norm_bond(:, ibond_sec), proj_vector) - end if ! choice of 2nd direction - if ( find_neighbours(1, find_atomno) == find_neighbours(4, find_atomno)) & - call cq_abort("axes_from_nn: cannot have the principal neighbour also as the second direction") - select case (find_neighbours(2, find_atomno)) - case (0) - pdos_ax = norm_bond(:, ibond_principal) - pdos_ay = proj_vector - pdos_az = cross_product(pdos_ax, pdos_ay) - case (1) - pdos_ay = proj_vector - pdos_az = norm_bond(:, ibond_principal) - pdos_ax = cross_product(pdos_ay, pdos_az) - end select - pdos_ax = pdos_ax / norm2(pdos_ax) - pdos_ay = pdos_ay / norm2(pdos_ay) - pdos_az = pdos_az / norm2(pdos_az) - end subroutine axes_from_nn - function cross_product(a, b) result(cross) - use datatypes - implicit none - real(double), dimension(3), intent(in) :: a, b - real(double), dimension(3) :: cross - - cross(1) = a(2) * b(3) - a(3) * b(2) - cross(2) = a(3) * b(1) - a(1) * b(3) - cross(3) = a(1) * b(2) - a(2) * b(1) - end function cross_product - - subroutine project_onto_plane(plane_normal, vector, proj_vector) - use datatypes - implicit none - real(double), intent(in) :: plane_normal(3), vector(3) - real(double), intent(out) :: proj_vector(3) - real(double) :: norm_plane_normal(3) - ! Subroutine to project vector onto plane defined by plane_normal vector - norm_plane_normal = plane_normal / norm2(plane_normal) - - proj_vector = vector - ((dot_product(vector, norm_plane_normal))*norm_plane_normal) - end subroutine project_onto_plane - ! End pDOS rotation subroutines + write(*,fmt='(/2x,"Secondary neighbour found: atom ", (I0,1x))') nghbr_arr(ibond_sec) + + call project_onto_plane(norm_bond(:, ibond_principal),norm_bond(:, ibond_sec), proj_vector) + end if ! choice of 2nd direction + if ( find_neighbours(1, find_atomno) == find_neighbours(4, find_atomno)) & + call cq_abort("axes_from_nn: cannot have the principal neighbour also as the second direction") + select case (find_neighbours(2, find_atomno)) + case (0) + pdos_ax = norm_bond(:, ibond_principal) + pdos_ay = proj_vector + pdos_az = cross_product(pdos_ax, pdos_ay) + case (1) + pdos_ay = proj_vector + pdos_az = norm_bond(:, ibond_principal) + pdos_ax = cross_product(pdos_ay, pdos_az) + end select + pdos_ax = pdos_ax / norm2(pdos_ax) + pdos_ay = pdos_ay / norm2(pdos_ay) + pdos_az = pdos_az / norm2(pdos_az) + end subroutine axes_from_nn + function cross_product(a, b) result(cross) + use datatypes + implicit none + real(double), dimension(3), intent(in) :: a, b + real(double), dimension(3) :: cross + + cross(1) = a(2) * b(3) - a(3) * b(2) + cross(2) = a(3) * b(1) - a(1) * b(3) + cross(3) = a(1) * b(2) - a(2) * b(1) + end function cross_product + + subroutine project_onto_plane(plane_normal, vector, proj_vector) + + use datatypes + + implicit none + + real(double), intent(in) :: plane_normal(3), vector(3) + real(double), intent(out) :: proj_vector(3) + real(double) :: norm_plane_normal(3) + ! Subroutine to project vector onto plane defined by plane_normal vector + norm_plane_normal = plane_normal / norm2(plane_normal) + + proj_vector = vector - ((dot_product(vector, norm_plane_normal))*norm_plane_normal) + end subroutine project_onto_plane + ! End pDOS rotation subroutines subroutine process_band_structure diff --git a/tools/PostProcessing/read_module.f90 b/tools/PostProcessing/read_module.f90 index 41545c62..ce5d33d7 100644 --- a/tools/PostProcessing/read_module.f90 +++ b/tools/PostProcessing/read_module.f90 @@ -332,75 +332,68 @@ subroutine read_input flag_rotate_pdos_mode = fdf_integer('Process.RotatePDOSMode',0) rotate_pdos_natoms = fdf_integer('Process.RotatePDOS.NumAtoms',1) flag_rotate_pdos_units = fdf_string(7, 'Process.RotatePDOSAngle',"deg") ! deg or rad - if (flag_rotate_pDOS) then - if(fdf_block('pDOSAxes') .and. flag_rotate_pdos_mode == 0) then - if(1+block_end-block_start<3) & - call cq_abort("Too few vectors in pDOS_axes: ",& - 1+block_end-block_start,3) - ! Expect exactly 3 vectors -> 3 reads - read (unit=input_array(block_start),fmt=*) pdos_ax - read (unit=input_array(block_start+1),fmt=*) pdos_ay - read (unit=input_array(block_start+2),fmt=*) pdos_az - call fdf_endblock - else if (fdf_block('pDOSEuler') .and. flag_rotate_pdos_mode == 1) then - if(1+block_end-block_start 3 reads + read (unit=input_array(block_start),fmt=*) pdos_ax + read (unit=input_array(block_start+1),fmt=*) pdos_ay + read (unit=input_array(block_start+2),fmt=*) pdos_az + call fdf_endblock + else if (fdf_block('pDOSEuler') .and. flag_rotate_pdos_mode == 1) then + if(1+block_end-block_start 2) & - call cq_abort("Local geometry flag in block pDOSNeighbours was not 0 or 1: ",& + if (find_neighbours(2,i) < 0 .or. find_neighbours(2,i) > 2) & + call cq_abort("Local geometry flag in block pDOSNeighbours was not 0 or 1: ",& 1+block_end-block_start,2) - if (find_neighbours(3,i) < -1 .or. find_neighbours(3,i) > ni_in_cell) & - call cq_abort("Input for principal axis must be -1 (shortest bond), 0 (longest bond) or neighbour in cell",& + if (find_neighbours(3,i) < -1 .or. find_neighbours(3,i) > ni_in_cell) & + call cq_abort("Input for principal axis must be -1 (shortest bond), 0 (longest bond) or neighbour in cell",& 1+block_end-block_start,3) - if (find_neighbours(4,i) < 0 .or. find_neighbours(4,i) > ni_in_cell) & - call cq_abort("Input for second axis must be 0 or neighbour in cell",& + if (find_neighbours(4,i) < 0 .or. find_neighbours(4,i) > ni_in_cell) & + call cq_abort("Input for second axis must be 0 or neighbour in cell",& 1+block_end-block_start,4) end do - call fdf_endblock + call fdf_endblock else if (fdf_block('pDOSAxisAngle') .and. flag_rotate_pdos_mode == 3) then - if(1+block_end-block_start