tfgg / magres-format Goto Github PK

From Paul

I'm not convinced by 'find_atom' returning NULL if it doesn't find the matching atom. But the definition doesn't make it clear what to do with an invalid atom reference. I think the best default behaviour is to skip such lines entirely.

David McKay has written such a method, sent to the CCP-NC list.

# Function to return the dihedral angle defined by four atoms

import numpy as np

import math

def my_dihedral(atoms, atom1, atom2, atom3, atom4, degrees=False):
  # Use least_mirror for positions rather than atom.position to allow for a closer image
  bond1 = np.subtract(atoms.least_mirror(atom2.position, atom1.position)[1], atom1.position)
  bond2 = np.subtract(atoms.least_mirror(atom3.position, atom2.position)[1], atom2.position)
  bond3 = np.subtract(atoms.least_mirror(atom4.position, atom3.position)[1], atom3.position)

  # Make unit vectors
  bond1 = bond1 / np.linalg.norm(bond1)
  bond2 = bond2 / np.linalg.norm(bond2)
  bond3 = bond3 / np.linalg.norm(bond3)

  norm1 = np.cross(bond1, bond2)
  norm2 = np.cross(bond2, bond3)

  # Get m, where m, norm1 and bond2 form an orthogonal frame
  m = np.cross(norm1, bond2)

  # x = cos(theta) where theta is the angle between norm1 and norm2
  x = np.dot(norm1, norm2)

  # y = cos(phi) where phi is the angle between m and norm2
  y = np.dot(m, norm2)

  # atan2(y, x) gives the theta over a range of 2pi with correct sign
  # therefore preferable over acos(x)

  if degrees:
    return np.rad2deg(math.atan2(y, x))
  else:
    return math.atan2(y, x)

From Paul

I think the "returned object" needs more thought. This is perhaps satisfactory for a "minimally conforming" parser, but a reference parser should allow the entire file to be returned or provide different read options. You should be able to write out a new file without information loss.

Documentation is needed for the Python module. Investigate how to generate it from docstrings.

From Paul

I'm wary about ignoring the units entirely. I don't think that using "non-standard" units should break major version (just not be 1.0).

From Paul

I've modified the definition of white-space; the initial version just used ' ', while the definition says anything that regex '\s' accepts. But the latter definition is ambiguous, and you can't include the line separator in the token separator. The definition doesn't actually define the line separator (I think '\n' is fine).

A problem that I've just met: shouldn't get_species(symbol) return a MagresAtoms object instead of a simple list? In this way it's possible to perform further operations like selecting only the atoms within a certain atom belonging to that specific subset.

In fact I am finding the different nature of the various return values rather confusing. I start with a MagresAtoms object; if I use a get_species, I get a list of MagresAtom. However, if I use .within, I get a list of MagresAtomImage, which in turn give me no access to the species of the original atom they are the image of, if I wanted to use an if condition to restrict the calculations only to certain species - image.atom.species() won't work because image.atom returns a string, not a MagresAtom! I find this formalism to be rather confusing. There should really be some intuitive way to do this, and I think we could do without some of the overloading of functions and members if it helped clarity.

I take it back. Just noticed all the documentation in notebooks/

I've just tried using convertoldmagres.py for the first time after installing the latest version and I got this never seen before error:

Traceback (most recent call last):
File "/usr/bin/convertoldmagres.py", line 31, in
magres_file = oldmagres_file.as_new_format()
File "/usr/lib/python2.7/site-packages/magres/oldmagres.py", line 144, in as_new_format
magres_file = format.MagresFile()
File "/usr/lib/python2.7/site-packages/magres/format.py", line 291, in init
magres.schema.validate.validate_magres(self.data_dict)
AttributeError: 'MagresFile' object has no attribute 'data_dict'

Any clues as to what's going on here? Did I do something wrong?

From Paul

I'm not sure what "Atoms can be specified out of order, but the numbers in the label should be in order." means in the definition. Does CASTEP impose restrictions on the "number in species? What should the parser do to check?

Notation like:

atoms.species('Si').set_reference(123.0)

?

Issue:

convertoldmagres.py gives no atomic positions when given old castep .magres file with ion labels.

How to replicate:

Give convertoldmagres.py a .magres file where some atoms have labels.

Fix:

Treat labels properly in parser.

From Paul

I'm uneasy about using pointers to the MagresAtom's. Using an index into the atom list would be more robust, although I'm not sure whether this is efficient with STL's std::list. The current implementation is OK if the object is being returned as const.

When a path to a file that doesn't exist is supplied to MagresAtoms.load_magres it attempts to read it as a string. This isn't a useful error if this is a mistake.

Fixing this, however, might mean removing the functionality to easily load strings as magres files. Is this important?

Following reported by Andrew Morris:

I've had magres-format generate an error.
"/home/ajm255/lib/python2.7/site-packages/magres/constants.py:17:
RuntimeWarning: invalid value encountered in double_scalars
return (eval[2], (eval[0]-eval[1])/eval[2])"

I know very little Python, but I seem to be able to generate it when
an efg tensor is exactly zero, see attached. Presumably this is
causing numpy to have problems finding its eigenvalues? In
constants.py: sorted_evals(mat)?

When looking for an atom with an element symbol, the program crashes if the element is not present into the loaded file, instead of just returning an empty selection. I suggest adding a condition to check whether a certain element is present amongst the keys of the atom dictionary.

From Paul

Note that I've written it against my own library code (the USE_LIBCMATRIX option). It wouldn't require much work to use different types, but the current code is therefore not portable.

Comparison of tests to JMol seems to get opposite sign. JMol uses the correct convention, which is (according to David):

"The convention for dihedral angles is that as you look down the B-C bond in a A-B-C-D system with A at 0 deg., a positive number is assigned to the dihedral where C is at 0-180 deg. in a clockwise direction. A negative is assigned if C is at 0-180 in the anticlockwise direction. In most software the value of the dihedral is not restricted, so if you assign an angle of +270 deg. it'll understand this as -90. Here's how IUPAC describe it: http://goldbook.iupac.org/T06406.html"

The sign should be inverted in the current code.