RESCU forces in unitcell vs in supercell in respect to grid resolution

shwhyshwhy

I have a question regarding the forces in unitcell and its own supercell under different grid resolutions.

I need minimal forces on all the atoms in order to get a correct phonon band structure. And I also need to use grids as large as they can be to reduce the computational cost since I need the phonon band structure on a large supercell with many atoms.

So I did several tests with domain.lowres = 0.2, 0.25, 0.3, 0.35

I found that when domain.lowres = 0.2, the unitcell and its counterpart supercell have very small forces.
However, when domain.lowres = 0.35, the unitcell still has very small forces, but its counterpart supercell has large forces.

I don't understand this discrepancy. The supercell is created by option.replicateCell = 3*[1,1,1] and the unicell has kpoint.gridn = [3 3 3], while the supercell has kpoint.gridn = [1 1 1]. So I think the two systems should be equivalent? Why the force is small in the unitcell, but becomes large on the supercell for grid size 0.35?

These are the full test results I have:

where E is obtained by: a=load('nacl_scf.mat'), a.energy.Etot, F is obtained by: a=load('nacl_scf.mat'),max(max(a.force.total))

The input file is:

info.calculationType = 'self-consistent'
info.savepath = './results/nacl_scf'
element(1).species = 'Na'
element(1).path = './Na_LDA_TZP.mat'
element(2).species = 'Cl'
element(2).path = './Cl_LDA_TZP.mat'
atom.element = [1 1 1 1 2 2 2 2]
atom.xyz = 10.348949044610158*...
[[0.000000 0.000000 0.000000]
[0.000000 0.500000 0.500000]
[0.500000 0.000000 0.500000]
[0.500000 0.500000 0.000000]
[0.500000 0.000000 0.000000]
[0.500000 0.500000 0.500000]
[0.000000 0.000000 0.500000]
[0.000000 0.500000 0.000000]]
option.replicateCell = 3*[1,1,1] % supercell
domain.latvec = 10.348949044610158*[1 0 0; 0 1 0; 0 0 1]
domain.lowres = 0.35
functional.libxc = true
functional.list = {'XC_LDA_X','XC_LDA_C_PW'}
%kpoint.gridn = [1 1 1] % kpoint.gridn = [3 3 3] % for unitcell 
kpoint.sigma = 0            % advised for dfpt non-metallic systems
mixing.tol = 1e-8*[1 1]     % advised for dfpt
option.saveWavefunction = 1 % required for dfpt
diffop.method = 'fft'       % required for dfpt
domain.fourierInit = false  % required for dfpt

The resculog.out comparison:

Could you please give me some guidance on this? Thank you so much!

vincentm

Dear shwhyshwhy ,

This is a good question. My guess would be that for certain resolutions, the grid of the supercell is not a strict superset of the original cell. For instance, your output shows that the original cell has cgridn = [30, 30, 30] while the supercell has cgridn = [89, 89, 89]. Not only is it not a multiple, but it is a somewhat large prime number, possibly breaking several symmetries. I would try forcing cgridn = [90, 90, 90] in that case and see whether the discrepancy resolves itself.

shwhyshwhy

Dear vincentm

I just did a test. For the case of grids=0.35, adding domain.cgridn = [90, 90, 90], the force indeed reduce from 0.018 to 8.5663e-13. That's amazing! Such a huge difference. You mentioned this is because of symmetries. If you could share more or point to some references I would love to learn more!

And for other materials, should I also manually set this domain.cgridn or let the program automatically choose? I found out that is also the reason that I get some abrupt changes when running EOS calculations, as shown below:
The Low resolution N just changes by itself when I loop on different lattices. So I should fix domain.cgridn as well for these cases, right? But how can I determine the values of cgridn for the correct symmetry for other materials? e.g. whether I should go with [30,30,30] or [31,31,31] (automatically chosen for larger lattices by the program but with less symmetry)?

Thank you so much!

vincentm

So I should fix domain.cgridn as well for these cases, right? But how can I determine the values of cgridn for the correct symmetry for other materials? e.g. whether I should go with [30,30,30] or [31,31,31] (automatically chosen for larger lattices by the program but with less symmetry)?

It is difficult for users to think of grid symmetries in general. At high enough resolution, these effects should go away indeed. The effects are more pronounced when using real space grids because the method is less translation invariant than planewaves (option.precision = 'pw-high') and there is some discrete jumps in resolution when you change the unit cell size, as you pointed out. You could fix the grid to avoid the jumps, but then make sure to try various grids since the large unit cell limit will have lower resolution than the small unit cell limit, which can modify the shape of the EOS curve.

shwhyshwhy

Dear vincentm

Thank you for your detailed reply! I try to use the grid to give relatively less energy （seems the number of more prime factors）and fix them during the loop. I got a much better EOS curve.

However, although the total energy for fixed configuration seems already converged for grids=0.35, the equilibrium lattice parameter is quite different when applying different grids resolution, is this normal? or is it converged? And should I use the lattice structure “self-consistently”? e.g. if I choose to use grids=0.35, I also use the lattice parameter optimized for grids=0.35 for all calculations, instead of using the lattice parameter relaxed with finer grids?

Thank you so much!

vincentm

However, although the total energy for fixed configuration seems already converged for grids=0.35, the equilibrium lattice parameter is quite different when applying different grids resolution, is this normal?

Yes, basically you should increase the resolution until convergence.

And should I use the lattice structure “self-consistently”? e.g. if I choose to use grids=0.35, I also use the lattice parameter optimized for grids=0.35 for all calculations, instead of using the lattice parameter relaxed with finer grids?

Yes, it is best to have all parameters as consistent as possible in a workflow.