Add local orbital rotation in CONQUEST for projecting pDOS along any axes

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 <post-proc>`.
 
+
+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 <post-proc>`.
+
 .. _pp_band_str
 
 Band structure

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},
 }

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