Unclear error associated with cell optimization

[mgoldey]

Hi everyone,

I don't quite understand why the assert in line 78 of qbox-1.60.4/src/StructureFactor.C is failing for a cell/structure optimization. Obviously, the cell has to change size, so this variable should be set to the new value at some point. Perhaps someone can enlighten me if a configuration option is setup incorrectly.

Thanks for the help,

Matthew

Input

Code: Select all

set cell 30 0 0  0 23.6 0  0 0 30
species oxygen http://fpmd.ucdavis.edu/potentials/O/O_HSCV_PBE-1.0.xml 
species carbon http://fpmd.ucdavis.edu/potentials/C/C_HSCV_PBE-1.0.xml 
species fluorine http://fpmd.ucdavis.edu/potentials/F/F_HSCV_PBE-1.0.xml 
species sulfur http://fpmd.ucdavis.edu/potentials/S/S_HSCV_PBE-1.0.xml 
species hydrogen http://fpmd.ucdavis.edu/potentials/H/H_HSCV_PBE-1.0.xml 
atom  C1 carbon     7.66922   5.26463  12.35543
atom  C2 carbon     8.24892   5.48282  15.02072
atom  S1 sulfur     8.73741   2.43299  16.52733
atom  C3 carbon     8.44649   0.98927  13.41567
atom  C4 carbon     7.62415   2.71564  11.37860
atom  H1 hydrogen   5.71819   2.20627  10.66092
atom  C5 carbon     7.27442   7.38572  10.58498
atom  C6 carbon     7.73124   9.85504  11.78565
atom  C7 carbon     8.20887  10.25310  14.46106
atom  C8 carbon     8.38861   8.04986  16.32947
atom  O1 oxygen    10.61290   8.21370  17.72765
atom  O2 oxygen     8.86714   7.11299   8.49735
atom  S2 sulfur     7.69904  12.56947   9.88816
atom  C9 carbon     8.21083  14.47816  12.61020
atom  C10 carbon    8.45086  13.01946  15.06874
atom  H2 hydrogen   8.84580  13.83159  16.96076
atom  C11 carbon    8.47087  17.26643  12.30608
atom  S3 sulfur     9.01885  19.42040  14.95935
atom  C12 carbon    8.93373  22.04954  12.75026
atom  C13 carbon    8.91570  21.21362  10.12377
atom  C14 carbon    8.41432  18.61148   9.83609
atom  H3 hydrogen   7.70276   7.14180   7.02711
atom  H4 hydrogen  10.04125   7.99976  19.50230
atom  C15 carbon    7.72399  18.05964   7.20132
atom  C16 carbon    7.82657  20.33628   5.64639
atom  F` fluorine   6.63974  15.93182   6.30401
atom  S4 sulfur     9.06398  23.04554   7.24279
atom  C17 carbon    6.70314  20.50491   3.05082
atom  O3 oxygen     6.82282  22.97567   2.12945
atom  O4 oxygen     7.95605  18.89064   1.39002
atom  H5 hydrogen   6.60215  18.23236   0.27589 
set stress ON
set ecut 60
set ecuts 60
set wf_dyn PSDA
set ecutprec 4
set xc PBE
set scf_tol 0.0001
#set debug STRESS
set cell_lock AC
set dt 3
set cell_dyn CG
randomize_wf
run 0 30
set cell_dyn CG
run 50 5

Last bit of output before crash with error message

Code: Select all

 <stress_tensor unit="GPa">
   <sigma_eks_xx>  -1.23062966 </sigma_eks_xx>
   <sigma_eks_yy>  -2.07122243 </sigma_eks_yy>
   <sigma_eks_zz>  -2.61414388 </sigma_eks_zz>
   <sigma_eks_xy>  -0.05393009 </sigma_eks_xy>
   <sigma_eks_yz>   0.43256190 </sigma_eks_yz>
   <sigma_eks_xz>  -0.20975772 </sigma_eks_xz>

   <sigma_kin_xx>   0.00000000 </sigma_kin_xx>
   <sigma_kin_yy>   0.00000000 </sigma_kin_yy>
   <sigma_kin_zz>   0.00000000 </sigma_kin_zz>
   <sigma_kin_xy>   0.00000000 </sigma_kin_xy>
   <sigma_kin_yz>   0.00000000 </sigma_kin_yz>
   <sigma_kin_xz>   0.00000000 </sigma_kin_xz>

   <sigma_ext_xx>   0.00000000 </sigma_ext_xx>
   <sigma_ext_yy>   0.00000000 </sigma_ext_yy>
   <sigma_ext_zz>   0.00000000 </sigma_ext_zz>
   <sigma_ext_xy>   0.00000000 </sigma_ext_xy>
   <sigma_ext_yz>   0.00000000 </sigma_ext_yz>
   <sigma_ext_xz>   0.00000000 </sigma_ext_xz>

   <sigma_xx>  -1.23062966 </sigma_xx>
   <sigma_yy>  -2.07122243 </sigma_yy>
   <sigma_zz>  -2.61414388 </sigma_zz>
   <sigma_xy>  -0.05393009 </sigma_xy>
   <sigma_yz>   0.43256190 </sigma_yz>
   <sigma_xz>  -0.20975772 </sigma_xz>
 </stress_tensor>
  CGCellStepper: alpha = 0.00200000
qb: StructureFactor.C:78: void StructureFactor::update(const std::vector<std::vector<double, std::allocator<double>>, std::allocator<std::vector<double, std::allocator<double>>>> &, const Basis &): Assertion `basis.localsize() == _ng' failed.

[fgygi]

When running variable cell simulations, a reference cell (Qbox variable name: ref_cell) must be defined in addition to the usual unit cell (Qbox variable name: cell). The reference cell must enclose the cell at all times during the simulation. Therefore the choice of the size of the reference cell requires guessing what the largest cell size will be during the simulation. It is customary to use a reference cell that is 5-10% larger than the cell. During the simulation, the cell may change size, but the reference cell remains fixed. The reference cell is used to determine the size of the plane wave basis set. Therefore, simulations are done with a fixed number of plane waves in the basis set. Since this would lead to variable resolution of the basis when the cell changes size, we compensate this by using a confinement potential in reciprocal space, which ensures that constant effective resolution is used even though the cell changes size.
The details of this approach are described in the papers cited in the Qbox User Guide under the description of the ecuts variable:

1. P. Focher, G. L. Chiarotti, M. Bernasconi, et al. Structural Phase-Transformations Via 1St-
Principles Simulation, Europhys. Lett. 26 (5): 345-351 (1994).
2. M. Bernasconi, G. L. Chiarotti, P. Focher, et al. First-Principle Constant-Pressure Molecular-
Dynamics, Journal Of Physics And Chemistry Of Solids 56 (3-4): 501-505 (1995).

When using the confinement potential, the effective plane wave cutoff used is determined by the value of the variable ecuts. It must be smaller than ecut for the confinement to be effective. In practice, it can be chosen to be 5 Ry smaller than ecut.

In the example given in your post, the reference cell was not defined. This caused the basis set to change size as soon as the cell changed, and this is caught as an error. Obviously a better error message would be needed.

If you want to perform a variable cell calculation for this problem, you can use the following script (but read the warning below):

Code: Select all

set cell 30 0 0  0 23.6 0  0 0 30
set ref_cell 32 0 0  0 25   0  0 0 32
species oxygen http://fpmd.ucdavis.edu/potentials/O/O_HSCV_PBE-1.0.xml
species carbon http://fpmd.ucdavis.edu/potentials/C/C_HSCV_PBE-1.0.xml
species fluorine http://fpmd.ucdavis.edu/potentials/F/F_HSCV_PBE-1.0.xml
species sulfur http://fpmd.ucdavis.edu/potentials/S/S_HSCV_PBE-1.0.xml
species hydrogen http://fpmd.ucdavis.edu/potentials/H/H_HSCV_PBE-1.0.xml
atom  C1 carbon     7.66922   5.26463  12.35543
atom  C2 carbon     8.24892   5.48282  15.02072
atom  S1 sulfur     8.73741   2.43299  16.52733
atom  C3 carbon     8.44649   0.98927  13.41567
atom  C4 carbon     7.62415   2.71564  11.37860
atom  H1 hydrogen   5.71819   2.20627  10.66092
atom  C5 carbon     7.27442   7.38572  10.58498
atom  C6 carbon     7.73124   9.85504  11.78565
atom  C7 carbon     8.20887  10.25310  14.46106
atom  C8 carbon     8.38861   8.04986  16.32947
atom  O1 oxygen    10.61290   8.21370  17.72765
atom  O2 oxygen     8.86714   7.11299   8.49735
atom  S2 sulfur     7.69904  12.56947   9.88816
atom  C9 carbon     8.21083  14.47816  12.61020
atom  C10 carbon    8.45086  13.01946  15.06874
atom  H2 hydrogen   8.84580  13.83159  16.96076
atom  C11 carbon    8.47087  17.26643  12.30608
atom  S3 sulfur     9.01885  19.42040  14.95935
atom  C12 carbon    8.93373  22.04954  12.75026
atom  C13 carbon    8.91570  21.21362  10.12377
atom  C14 carbon    8.41432  18.61148   9.83609
atom  H3 hydrogen   7.70276   7.14180   7.02711
atom  H4 hydrogen  10.04125   7.99976  19.50230
atom  C15 carbon    7.72399  18.05964   7.20132
atom  C16 carbon    7.82657  20.33628   5.64639
atom  F fluorine   6.63974  15.93182   6.30401
atom  S4 sulfur     9.06398  23.04554   7.24279
atom  C17 carbon    6.70314  20.50491   3.05082
atom  O3 oxygen     6.82282  22.97567   2.12945
atom  O4 oxygen     7.95605  18.89064   1.39002
atom  H5 hydrogen   6.60215  18.23236   0.27589
set stress ON
set ecut 95
set ecuts 90
set wf_dyn PSDA
set xc PBE
set scf_tol 1.e-6
randomize_wf
run 0 300
set atoms_dyn CG
run 200 20
save restart_file.xml

(Note that a few changes were made, in particular the number of iterations and scf_tol).

In this initial configuration, the molecule described here is not in its equilibrium geometry. The cell_dyn=CG algorithm will optimize simultaneously the atomic positions and the cell parameters. However, the calculated stress tensor will initially have a large and non-physical value due to the fact that the atoms are not in equilibrium positions. Updating the cell with this stress tensor will likely (eventually) work, but it may take a large number of iterations to converge since the cell dimension (in the y direction) is strongly coupled to the atomic positions. It would be preferable to use the following two-step approach:

1) Relax atomic positions to reach the equilibrium geometry (with the cell fixed).
2) Relax both atomic positions and cell.

Here is a script that will find equilibrium positions (with cell fixed). Note that we still use a reference cell and stress=ON to be able to follow this calculation with a cell relaxation).

Code: Select all

# set nrowmax to a value larger than the default if running on 64 or 128 tasks
set nrowmax 64 
set cell 30 0 0  0 23.6 0  0 0 30
set ref_cell 32 0 0  0 25   0  0 0 32
species oxygen http://fpmd.ucdavis.edu/potentials/O/O_HSCV_PBE-1.0.xml
species carbon http://fpmd.ucdavis.edu/potentials/C/C_HSCV_PBE-1.0.xml
species fluorine http://fpmd.ucdavis.edu/potentials/F/F_HSCV_PBE-1.0.xml
species sulfur http://fpmd.ucdavis.edu/potentials/S/S_HSCV_PBE-1.0.xml
species hydrogen http://fpmd.ucdavis.edu/potentials/H/H_HSCV_PBE-1.0.xml
atom  C1 carbon     7.66922   5.26463  12.35543
atom  C2 carbon     8.24892   5.48282  15.02072
atom  S1 sulfur     8.73741   2.43299  16.52733
atom  C3 carbon     8.44649   0.98927  13.41567
atom  C4 carbon     7.62415   2.71564  11.37860
atom  H1 hydrogen   5.71819   2.20627  10.66092
atom  C5 carbon     7.27442   7.38572  10.58498
atom  C6 carbon     7.73124   9.85504  11.78565
atom  C7 carbon     8.20887  10.25310  14.46106
atom  C8 carbon     8.38861   8.04986  16.32947
atom  O1 oxygen    10.61290   8.21370  17.72765
atom  O2 oxygen     8.86714   7.11299   8.49735
atom  S2 sulfur     7.69904  12.56947   9.88816
atom  C9 carbon     8.21083  14.47816  12.61020
atom  C10 carbon    8.45086  13.01946  15.06874
atom  H2 hydrogen   8.84580  13.83159  16.96076
atom  C11 carbon    8.47087  17.26643  12.30608
atom  S3 sulfur     9.01885  19.42040  14.95935
atom  C12 carbon    8.93373  22.04954  12.75026
atom  C13 carbon    8.91570  21.21362  10.12377
atom  C14 carbon    8.41432  18.61148   9.83609
atom  H3 hydrogen   7.70276   7.14180   7.02711
atom  H4 hydrogen  10.04125   7.99976  19.50230
atom  C15 carbon    7.72399  18.05964   7.20132
atom  C16 carbon    7.82657  20.33628   5.64639
atom  F fluorine   6.63974  15.93182   6.30401
atom  S4 sulfur     9.06398  23.04554   7.24279
atom  C17 carbon    6.70314  20.50491   3.05082
atom  O3 oxygen     6.82282  22.97567   2.12945
atom  O4 oxygen     7.95605  18.89064   1.39002
atom  H5 hydrogen   6.60215  18.23236   0.27589
set stress ON
set ecut 95
set ecuts 90
set wf_dyn PSDA
set xc PBE
set scf_tol 1.e-6
randomize_wf
run 0 300
set atoms_dyn CG
run 200 20
save sample.xml

The following script can then be used to relax the cell

Code: Select all

set nrowmax 64 
load sample.xml
set stress ON
set ecuts 90
set wf_dyn PSDA
set xc PBE
set scf_tol 1.e-7
set cell_dyn CG
set cell_lock ACS
run 200 30
save relax1.xml

Note the use of cell_lock = ACS, which keeps the length of A and C fixed and the shape of the cell rectangular.

[mgoldey]

Thank you very much for your help explaining that oddity. I'll comment further if any other problems arise. It makes sense in hindsight, but I did not intuit the use of ref_cell from the manual's description.

[fgygi]

The User Guide is clearly insufficient in terms of documentation of this feature. Additional examples should be added for variable cell calculations. This is one of the many items on the todo list