Error during Symmetrization with custom gradient.

[JonathanSchmidt1]

I was running some relaxations of a ternary structure with one species fixed using the custom gradient function from the FixAtoms example (The relaxations runs fine without the custom gradient). The calculations seem to run fine for a while until this error pops up:

  self.dyn.structure.coords[:,:] = new_struct
Traceback (most recent call last):
  File "SSCHA_relax_fixed_atom.py", line 213, in <module>
    relax.relax(ensemble_loc = "ensemble/")
  File "/opt/conda/lib/python3.7/site-packages/sscha/Relax.py", line 365, in relax
    custom_function_gradient = self.__cfg__)
  File "/opt/conda/lib/python3.7/site-packages/sscha/SchaMinimizer.py", line 1163, in run
    self.minimization_step(custom_function_gradient)
  File "/opt/conda/lib/python3.7/site-packages/sscha/SchaMinimizer.py", line 541, in minimization_step
    self.dyn.Symmetrize(use_spglib = self.use_spglib)
  File "/opt/conda/lib/python3.7/site-packages/cellconstructor/Phonons.py", line 3181, in Symmetrize
    qe_sym.SymmetrizeFCQ(fcq, self.q_stars, asr = asr, verbose = verbose)
  File "/opt/conda/lib/python3.7/site-packages/cellconstructor/symmetries.py", line 900, in SymmetrizeFCQ
    self.ApplyQStar(fcq[q0_index : q0_index + q_len, :,:], np.array(q_stars[i]))
  File "/opt/conda/lib/python3.7/site-packages/cellconstructor/symmetries.py", line 759, in ApplyQStar
    count))
ValueError: Error, the vector (-0.061, -0.061, 0.043) has 0 identification in the star
double free or corruption (!prev)

I am running a rather old version of SSCHA so I will see if an update helps. Besides that anyone got any idea what's happening?
thank you very much for you help.
Jonathan

[mesonepigreco]

Hi Jonathan,
This error means there is a change of the symmetries during the minimization, and the star of q points changed (that explains why the code is complaining that it cannot find a q vector in the symmetry-generated star of q).

Using a custom gradient function like FixAtoms could introduce a symmetry breaking.
My suggestion is to replace the fixed atoms with the constraints on the modes if you can obtain a similar effect (for example, if you have heavy atoms that you want to constrain in the simulation, you can achieve the same by constraining all small frequency modes).

Can you share the snippet of the code you employ to fix the atoms?
All the best,
Lorenzo

[JonathanSchmidt1]

Thank you for the quick reply.
I was just using the snippet from the example:

fix_struct = np.array([x_at == FIX_ATOM_TYPE for x_at in dyn.structure.atoms])

# Obtain an array that can be applied also on the dyn gradient
# By repeating each value 3 times (the xyz cartesian coordinates)
fix_dyn = np.tile(fix_struct, (3, 1)).T.ravel()
all_grads = []

def fix_atoms(gradient_dyn, gradient_struct):
    # gradient_struct is a (n_at, 3) shaped array,
    # that contains the forces on atoms at each step.
    # We put to zero this force on the atoms we want
    #gradient_struct[fix_struct, :] = 0
    
    # Now the gradient violates translations
    # Compute and subtract the average force
    av_force = np.mean(gradient_struct, axis = 0)
    gradient_struct[:,:] -= np.tile(av_force, (dyn.structure.N_atoms, 1))

    all_grads.append(gradient_struct)

    
    # gradient_dyn is a (nq, 3*n_at, 3*n_at) shaped array.
    # nq is the q points, then there is the dynamical matrix.
    # Here cartesian and atomic indices are contracted.
    # For this reason we created a fix_dyn mask.
    gradient_dyn[:, :, fix_dyn] = 0
    gradient_dyn[:, fix_dyn, :] = 0

    # We now must impose the acoustic sum rule
    # We violated it because we fixed some atoms.
    # It can be imposed in the gamma point
    CC.symmetries.CustomASR(gradient_dyn[0, :, :])

and than adding the the function to the relaxation setup:
relax.setup_custom_functions(custom_function_gradient = fix_atoms, custom_function_post = custom_functions)
I have a SrTiO3 structure and I was fixing the Ti positions but I guess there would be little harm in my case to fix the Sr positions as well. However, I would be surprised if the mass difference between oxygen and titanium is sufficiently large. What heuristic do you use to select the modes you want to fix?