Potential Bug: Space group not maintained during cell optimization despite KEEP_SPACE_GROUP=T

[zhangbw20]

Dear developers of cp2k, Greetings

Description

I encountered an issue where the space group was not maintained during lattice optimization despite setting "KEEP_SPACE_GROUP" to true. This issue was observed while optimizing the lattice structures of iodine (No.64,Cmce) and silicon (No.227,Fd-3m) in their primitive cell forms, where the non-zero elements of three upper-triangular elements in lattice were left unchanged during the whole process of optimization.

$ grep -A2 "Vector a" I2_N01.out
 CELL| Vector a [angstrom]:       3.763    -2.281     0.000   |a| =     4.400496
 CELL| Vector b [angstrom]:       3.763     2.281     0.000   |b| =     4.400496
 CELL| Vector c [angstrom]:       0.000     0.000     9.706   |c| =     9.706146
--
 CELL| Vector a [angstrom]:       3.802    -2.281     0.000   |a| =     4.433490
 CELL| Vector b [angstrom]:       3.802     2.325     0.000   |b| =     4.456413
 CELL| Vector c [angstrom]:       0.000     0.000     9.693   |c| =     9.692778
--
 CELL| Vector a [angstrom]:       4.028    -2.281     0.000   |a| =     4.629319
 CELL| Vector b [angstrom]:       4.028     2.575     0.000   |b| =     4.781154
 CELL| Vector c [angstrom]:       0.000    -0.000     9.604   |c| =     9.603716
--
 CELL| Vector a [angstrom]:       4.269    -2.281     0.000   |a| =     4.840494
 CELL| Vector b [angstrom]:       4.269     2.825     0.000   |b| =     5.119503
 CELL| Vector c [angstrom]:       0.000    -0.000     9.495   |c| =     9.494987
--
 CELL| Vector a [angstrom]:       4.487    -2.281     0.000   |a| =     5.033592
 CELL| Vector b [angstrom]:       4.487     3.034     0.000   |b| =     5.416581
 CELL| Vector c [angstrom]:       0.000    -0.000     9.384   |c| =     9.383704
--
...

Relevant code

Line 129~137 of gopt_f_methods.f:

            ALLOCATE (x0(3*nparticle + 6))
            CALL pack_subsys_particles(subsys=subsys, r=x0)
            idg = 3*nparticle
            DO i = 1, 3
               DO j = 1, i
                  idg = idg + 1
                  x0(idg) = cell%hmat(j, i)
               END DO
            END DO

The three upper-triangular elements were not read into x0 during the initialization of x0 which is later updated by dr for geoopt_bfgs function in bfgs_optimizer.F.
I didn't go through all the later functions to reassure whether there were other symmetry operations implemented for x0, but the output file results would show that the three upper-triangular elements were left unchanged so that the space group can't be maintained if any of the three upper-triangular elements are non-zero (which are common for primitive cell of C-centered orthorhombic lattice, FCC, etc.).

Potential Impacts

Thus, maybe cp2k in its current version (cp2k-2024.1) is malfunctioned for some primitive cell optimization. Additionally, users (like me) using primitive cell for band structure calculation might suffer from the cell relaxation, especially recently RI-HFX techniques were employed which are especially suitable for small crystal's band calculations with a great number of K points. (Certainly, this could be evaded by using a conventional cell for cell optimization and then use primitive cell for other properties calculation.)

Steps to Reproduce

Two test cases from ./QS/regtest-spgr/quartz.inp with minimal modifications are attached as postscripts.
I2_N01.inp.txt
Si_N01.inp.txt
version: cp2k-2024.1, spglib-v1.16.2 (latest as today)

Questions

Are there any known workarounds for this issue?

Notes

1. Note that the identify_space_group function (calling spgr_create) in space_groups.F uses the full hmat (all 3-by-3 elements) so that it could correctly identify the space group of a conventional unit cell (or primitive cell) by spglib (v1.16.2) in both example cases.
   Line 228 of space_groups.F:
   spgr%space_group_number = spg_get_international(spgr%international_symbol, TRANSPOSE(cell%hmat), tmp_coor, tmp_types, & spgr%n_atom_sym, eps_symmetry)

2. By the way, are there any other better ways than modifying the src code to WRITE outputs and recompiling for cp2k debug?
   I found that:
   a) the stacks of function being called are not fully written out into the output file so it would be very hard to see the layer structure and order of all the executing functions even if the TRACE tag in input file was set to T;
   b) I followed "debug build" instructions in https://manual.cp2k.org/trunk/CP2K_INPUT/DEBUG.html, but still couldn't let gdb TUI find source code;
   c) So using WRITE sentences and recompiling is my final choice to read the values of x0 to finally locate the problem, which took me quite a lot of time.

3. Earlier discussions are at https://groups.google.com/g/cp2k/c/Xwna6hLwWzk Thanks a lot for Prof. Hutter's discussions and instructions for me.