Gradients on the cell not available for most calculations
paulsaxe opened this issue · 2 comments
Is your feature request related to a problem? Please describe.
The gradients of the translation vectors are not available in the AUX file except for 1SCF calculations, which means it is impossible to calculate the stress.
Describe the solution you'd like
Add the gradients of the translation vectors to the AUX file win a consistent fashion. Currently they are appended to the atomic gradients for 1SCF calculations -- 9 extra values for a 3-D periodic system. They could be appended in other calculations -- optimization, FORCE, etc. -- or it might be cleaner to add them as a separate field, as the translation vectors are currently:
INITIAL_TRANS_VECTS:ANGSTROMS[9]=
5.4309 0.0000 0.0000
0.0000 5.4309 0.0000
0.0000 0.0000 5.4309
They are not appended to the coordinates, so it would be more consistent to not append the gradients (and also might break fewer codes!)
It would also be easy and valuable to calculate the complete stress tensor and make it available, probably in Voigt notation. The documentation provides the formula in the section on the pressure on the faces of cells:
The scalar, f1, of the component of the gradient vector, g1, in the direction of the translation vector, Tv1, is given by:
f1 = (g1(x)·Tv1(x) + g1(y)·Tv1(y) + g1(z)·Tv1(z)).
Using the volume of the unit cell, V = |Tv1·(Tv2⊗Tv3)| or:
V = |Tv3(x)*(Tv1(y)·Tv2(z) - Tv1(z)·Tv2(y)) + Tv3(y)*(Tv1(z)·Tv2(x) - Tv1(x)·Tv2(z)) + Tv3(z)*(Tv1(x)·Tv2(y) - Tv1(y)·Tv2(x))|,
and conversion factors to go from kcals and moles to ergs and unit cells, the pressure, P, in Pascals, can be calculated using:
P = -4184·1030/(6.022·1023)·f·V-1
For a 3-D system there are 3 translation vectors and their corresponding gradients. If we label them 1, 2, 3 or x, y, z then the stress tensor following the equation for P (which is really stress) is:
S11 S12 S13
S21 S22 S23
S31 S32 S33
The Voigt notation is a 6-vector: Sxx Syy Szz Syz Sxz Sxy (note the unusual order!). The terms are:
Sxx = S11
Syy = S22
Szz = S33
Syz = (S23 + S32) / 2
Sxz = (S13 + S32) / 2
Sxy = (S12 + S21) / 2
Which accounts for the small symmetry breaking in MOPAC.
These would be useful both in the output and in the AUX file.
Sorry about the delay in addressing this issue, which on the surface was a very straightforward request. I started working on it multiple times this year, but ended up repeatedly getting confused and stopping. After finally committing to work through the confusion, I think the major points of this request are now resolved. There are some lingering issues, but they will require a more thoughtful approach when I focus on interface revisions in upcoming versions.
About the stress tensor, it is now implemented and tested against a finite-difference reference. There were a lot of minor confusions and problems along the way. First, the sign of the pressure is internally flipped in MOPAC. Second, there isn't universal agreement on the sign of the stress tensor itself in electronic structure codes, although I'm using the more standard engineering convention here (i.e. the diagonals of the stress tensor are the negative of pressure). Third, your assertion that Jimmy's note on "Pressure on the faces of a unit cell" was directly related to stress tensor calculations was not correct, and I needed to more broadly refresh my knowledge of the topic. Jimmy is calculating a pressure on a crystal face in a very literal way, which is in the direction of the translation vectors, whereas the stress tensor indices are the Cartesian coordinates (i.e. they are distinct derivatives if the translation vectors are not aligned with Cartesian coordinates). Also, proper calculations of stress tensors are calculating a uniform deformation of space, so the gradients of atomic coordinates also contribute to the stress tensor (see DOI:10.1016/j.cpc.2015.01.003 for a thorough discussion). Finally, the finite-pressure gradients in MOPAC were never implemented correctly for non-orthorhombic crystals, and I ended up having to debug and re-implement them. This bug might have something to do with "the small symmetry breaking in MOPAC" you are referring to.
About printing gradients and stresses more consistently in the AUX file, I've added them in whenever updated coordinates are printed (i.e. geometry optimization and reaction coordinate methods). While I've added in more print statements, they will only trigger if the gradients are calculated. I am not printing gradient information with FORCE, because it requires the gradients to be very close to zero to function correctly, so gradient information should be useless when that feature is used as intended.
About separating atomic-coordinate and translation-vector derivatives, there are some subtle problems with changing the existing approach. MOPAC has 3 overlapping systems of bookkeeping - the "atoms" (which include things like sparkles and translation vectors), the "real atoms" (which are specifically atoms with atomic orbitals attached to them), and the optimized degrees of freedom (which include the x, y, and z coordinates of each atom and translation vector). If a particular coordinate of a particular atom is not set to be optimized, then it won't be on the list of optimized degrees of freedom and its derivative won't be calculated. MOPAC only prints the gradient components that it calculates. I can't guarantee that a list of gradients on all of the atoms or translation vectors doesn't have gaps because some components weren't calculated. I could retain the current approach and also print the set of atomic and translation vector gradients separately when they are all available, but that would be redundant and conditional. This is a case where MOPAC's flexibility complicates its interface, no matter which approach you take. A simpler interface would involve restricting some of that flexibility (e.g. forcing all coordinates to be optimized).
Weeeelllll, I didn't realize it would be that complicated! Sorry! But thanks for the tremendous effort. This will be very useful and important for periodic systems. As for the limitations you mention at the end, I would think that in essentially all cases, no parameters will be fixed, with the possible exception of lattice parameters that are zero by symmetry, which should be simple to handle externally. For example, in an orthorhombic cell you might choose to fix the angles at 90º so that noise doesn't slip into them, optimizing only the three cell lengths.