Symmetry analysis causing abnormal behavior in slab model relaxation calculation

[lanshuyue]

Describe the bug

When calculating the surface energy of a 7-layer bcc Mo, relaxation calculations were performed. In this process, the bottom three layers were held fixed to replicate a semi-infinite metal crystal, while the top four layers were allowed to relax.

When the symmetry setting was defaulted to 1, the calculation encountered issues as the structure optimization failed to converge. After each relaxation step, the LARGEST GRAD gradually increased, reaching up to 800 eV/Å.

In running_relax.log, it was noted that with each relaxation step, the atomic forces progressively increased. It was further observed that the forces acting on the fixed atoms and free atoms were symmetrical, which appeared to be the root cause of the abnormal behavior during the relaxation calculation process.

Upon changing the symmetry setting to 0, the structure optimization proceeded as expected.

Expected behavior

No response

To Reproduce

1. Set up a 7-layer bcc Mo structure with the specified layer relaxation conditions.

ATOMIC_SPECIES
Mo  95.95         Mo_ONCV_PBE-1.0.upf

NUMERICAL_ORBITAL
Mo_gga_7au_100Ry_4s2p2d1f.orb

LATTICE_CONSTANT
1.8897261258369282 

LATTICE_VECTORS
28.519325000      0.0000000000      0.0000000000
0.0000000000      2.2532210000      -1.593268000
0.0000000000      2.2532210000      1.5932680000

ATOMIC_POSITIONS
Direct

Mo
0.0000000000
7
0.2629800000 0.0000000000 0.0000000000 0 0 0 mag 0.0 
0.3419860000 0.5000000000 0.5000000000 0 0 0 mag 0.0 
0.4209930000 0.0000000000 0.0000000000 0 0 0 mag 0.0 
0.5000000000 0.5000000000 0.5000000000 1 1 1 mag 0.0 
0.5790070000 0.0000000000 0.0000000000 1 1 1 mag 0.0 
0.6580140000 0.5000000000 0.5000000000 1 1 1 mag 0.0 
0.7370200000 0.0000000000 0.0000000000 1 1 1 mag 0.0 

3. Use the default symmetry setting of 1 for the relaxation calculations. INPUT:

INPUT_PARAMETERS
#Parameters (1.General)
suffix                  Moslab15_fix6
calculation             relax
ntype                   1
nbands                  52
#symmetry                0
pseudo_dir              ./
cal_stress              1

#Parameters (2.Iteration)
ecutwfc                 100
scf_nmax                100
relax_nmax              100

#Parameters (3.Basis)
basis_type             lcao

#Parameters (4.Smearing)
smearing_method         mp
smearing_sigma          0.0036

#Parameters (5.Mixing)
mixing_type             broyden
mixing_beta             0.4

KPT:

K_POINTS
0
Gamma
1 7 7 1 1 1

5. Monitor the relaxation steps and observe the increasing forces on the atoms.
6. Check for symmetric forces acting on fixed and free atoms, leading to convergence failure.
7. Repeat the process after setting the symmetry to 0 for successful optimization.

Environment

abacus3.5.3

Additional Context

No response

Task list for Issue attackers (only for developers)

• [ ] Verify the issue is not a duplicate.
• [ ] Describe the bug.
• [ ] Steps to reproduce.
• [ ] Expected behavior.
• [ ] Error message.
• [ ] Environment details.
• [ ] Additional context.
• [ ] Assign a priority level (low, medium, high, urgent).
• [ ] Assign the issue to a team member.
• [ ] Label the issue with relevant tags.
• [ ] Identify possible related issues.
• [ ] Create a unit test or automated test to reproduce the bug (if applicable).
• [ ] Fix the bug.
• [ ] Test the fix.
• [ ] Update documentation (if necessary).
• [ ] Close the issue and inform the reporter (if applicable).

[kirk0830]

Noted. However I have some suggestions on your calculation setting and procedures. Does the following description match what you want to calculate?

[dyzheng]

@maki49 Hello, Do you have time to look at this issue?

[maki49]

The cause of this problem is: if in relaxation some atoms are fixed, the symmetry cannot be kept, but in ABACUS's relax procedure symmetry is analyzed only before the 1st ion step, so from the 2nd ion step the original higher symmetry is applied to the symmetrization of charge/force/stress of the current structure with lower symmetry.

@dyzheng To fix it, I will try to study and implement VASP's feature of "Analysis of constrained symmetry for selective dynamics", but it takes time. @lanshuyue You can continue your jobs with symmetry=0.

[lanshuyue]

@kirk0830 yes, i think it's quite similar.

[kirk0830]

@kirk0830 yes, i think it's quite similar.

@lanshuyue if so, I think I would suggest you should relax both side of slab, increasing the thickness of slab or say number of atomic layers, calculate surface energy till you find there is a tendency approaching a converged value, I think this way might be more appropriate. Please note if you only relax one side, then the other side surface energy will be impossible to evaluate, it is the reason why I suggest relax both side concurrently.

HTH.

[lanshuyue]

@kirk0830 yes, i think it's quite similar.

@lanshuyue if so, I think I would suggest you should relax both side of slab, increasing the thickness of slab or say number of atomic layers, calculate surface energy till you find there is a tendency approaching a converged value, I think this way might be more appropriate. Please note if you only relax one side, then the other side surface energy will be impossible to evaluate, it is the reason why I suggest relax both side concurrently.

HTH.

@kirk0830 Thank you for your suggestions!! I totally get that it's a common approach for calculating surface energies. However, my ultimate goal is to calculate adsorption energies, calculating surface energies being just one step. An asymmetric slab model aligns better with the system i'm analyzing. Plus, i'm currently replicating results from a research paper using this particular surface structure.

Thank you sooooo much for your advice and assistance!🤠