Creating internal coordinates fails in molecules with linear angles.

[S0S-90]

When adding waters with task SOLVEADD they are added as cartesian coordinates, but before they are written into the outputfile they are converted to internal coordinates and back. After this the coordinates are different than before.

I tried to track down the error and found that the number of main dihedrals which is used in the conversion from cartesian to internal is three less in SOLVEADD than when just doing the same conversion in main-function. I guess the reason for this is that the internal coordinates between the original system and the newly added waters are not taken into account because what is different are only the coordinates of the water molecules (see picture).

This leads sometimes to an atom crash (also see picture).

[S0S-90]

The different number of main dihedrals is not the reason for the error as there are less main dihedrals because the original atoms are fixed.

[JulianEr]

Maybe the following helps, though this rather makes the transformation to Cartesians a bit fishy than causing a wrong transformation: CAST does not ensure that the internal representation is even meaningful in a chemical sense. So dihedral angles are chosen somewhat arbitrary. Let's say you have a molecule A-B-C-E-D and the sequence in the XYZ-file is ABCDE (different sequence!) CAST will write a Z-Matrix describing ABC, followed by D. D hast neither A, B, or C as bonding partner (because E is the only bonding partner for D in this example), but will use one of the former atoms as first value in the Z-Matrix. This bond makes no chemical sense of course. If you look at the water molecules they all have dihedral partner in another molecule. This intermolecular internal coordinate might be hard to handle (you have at least one intermolecular bond and two angles, too, if I'm not mistaken, they won't make this easier). This should not brake the transformation as mentioned earlier but might make it less precise? Hope this is helpful in any sense. If not I'm afraid I can't help, sorry. At least it might be good to know.

[S0S-90]

I now found at least part of the error: When internal coordinates are built up three additional virtual atoms are created to define the internal coordinates of the first atom. This is necessary for correct placement of the molecule in space as the internal coordinates only give positions relative to the other atoms. If I want to solvate the molecule above my original system has 6 atoms so the internal indices of those virtual atoms are 6, 7 and 9. But after solving it there are more atoms and these indices correspond to some water atoms. That's why the Z-matrix is built up wrong. I solved this problem by clearing the bound_internals() of every atom before creating the new coordinates object (file startopt_solveadd.cc, function populate_coords()). Maybe we should keep in mind that this has to be done also in other cases when a new coordinates object is build from an old one.

However sometimes the conversion still fails (at least not always like before). I tracked the error down to the call of spherical() in function coords::Atoms::c_to_i() (file coords_atoms.cc). This is very strange as I can write out the input for this function, e.g.:

xyz[int_i] = -5.999632,1.905821,0.338626
rel_bond = -5.442202,2.144873,1.069806 rel_angle = -4.268664,2.648142,2.609134 rel_dihedral = -0.349515,-0.123188,2.609134 rel_angle - rel_bond = 1.173538,0.503269,1.539328 rel_dihedral - rel_bond = 5.092687,-2.268061,1.539328

Then the result for spherical is 0.949999,179.999969,141.839828 when just doing forth and back conversion in DEVTEST and 0.950000,180.000000,-96.671162 in task STARTOPT. The first two numbers might be rounding errors but the third is definitely something else. Very, very strange...

[S0S-90]

The problem is in file scon_vect.h, function angle_reference(): If the dot-product between a and b (line 733) is very near to 0, the atan2 is very different depending on the sign of the dot-product. So the different values can be explained by a rounding error somewhere in this dot-product. However as this results in completely different conformations (and sometimes crashing atoms) this should be considered as a bug.

[S0S-90]

The problem lies in the construction of the Z-matrix. For every atom the internal coordinates are given by a bond length, an angle and a dihedral angle. However if three of the 4 atoms that define the dihedral lie in line (angle = 180° like in the example above) it is not possible to define a dihedral angle. As linear molecules are not that rare in chemistry I guess there is a solution to this problem but at the moment I don't see one.

[S0S-90]

I have "solved" this problem now by integrating an exception in the to_internal() conversion if an angle is 180.0° (or 0.0° as this should cause the same error). Obviously this is the same way as Gaussian treats such problems (see https://www.cup.uni-muenchen.de/ch/compchem/geom/internal.html). If the cartesian to internal conversion fails, the backwards conversion is not performed but instead the original cartesian coordinates are used. This should work in most cases (and should rarely happen as mostly we don't have angles that are exactly 180°).