symmetry.c

/*
 * symmetry.c
 *
 * Symmetry
 *
 * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
 *                  a research centre of the Helmholtz Association.
 *
 * Authors:
 *   2010-2012 Thomas White <taw@physics.org>
 *
 * This file is part of CrystFEL.
 *
 * CrystFEL is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * CrystFEL is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <assert.h>

#include "symmetry.h"
#include "utils.h"
#include "integer_matrix.h"


/**
 * SECTION:symmetry
 * @short_description: Point symmetry handling
 * @title: Symmetry
 * @section_id:
 * @see_also:
 * @include: "symmetry.h"
 * @Image:
 *
 * Routines to handle point symmetry.
 */


struct _symoplist
{
	IntegerMatrix **ops;
	int n_ops;
	int max_ops;
	char *name;
	int num_equivs;
};


struct _symopmask
{
	const SymOpList *list;
	int *mask;
};



static void alloc_ops(SymOpList *ops)
{
	ops->ops = realloc(ops->ops, ops->max_ops*sizeof(IntegerMatrix *));
}


/**
 * new_symopmask:
 * @list: A %SymOpList
 *
 * Returns: a new %SymOpMask, which you can use when filtering out special
 * reflections.
 **/
SymOpMask *new_symopmask(const SymOpList *list)
{
	SymOpMask *m;
	int i;

	m = malloc(sizeof(struct _symopmask));
	if ( m == NULL ) return NULL;

	m->list = list;
	m->mask = malloc(sizeof(int)*list->n_ops);
	if ( m->mask == NULL ) {
		free(m);
		return NULL;
	}

	for ( i=0; i<list->n_ops; i++ ) {
		m->mask[i] = 1;
	}

	return m;
}


/* Creates a new SymOpList */
static SymOpList *new_symoplist()
{
	SymOpList *new;
	new = malloc(sizeof(SymOpList));
	if ( new == NULL ) return NULL;
	new->max_ops = 16;
	new->n_ops = 0;
	new->ops = NULL;
	new->name = NULL;
	new->num_equivs = 1;
	alloc_ops(new);
	return new;
}


/**
 * free_symoplist:
 * @ops: A %SymOpList to free
 *
 * Frees a %SymOpList and all associated resources.
 **/
void free_symoplist(SymOpList *ops)
{
	int i;

	if ( ops == NULL ) return;
	for ( i=0; i<ops->n_ops; i++ ) {
		intmat_free(ops->ops[i]);
	}
	if ( ops->ops != NULL ) free(ops->ops);
	if ( ops->name != NULL ) free(ops->name);
	free(ops);
}

/**
 * free_symopmask:
 * @m: A %SymOpMask to free
 *
 * Frees a %SymOpMask and all associated resources.
 **/
void free_symopmask(SymOpMask *m)
{
	if ( m == NULL ) return;
	free(m->mask);
	free(m);
}


/* This returns the number of operations in "ops".  This might be different
 * to num_equivs() if the point group is being constructed. */
static int num_ops(const SymOpList *ops)
{
	return ops->n_ops;
}


/**
 * add_symop:
 * @ops: A %SymOpList
 * @m: An %IntegerMatrix
 *
 * Adds @m to @ops.
 **/
void add_symop(SymOpList *ops, IntegerMatrix *m)
{
	if ( ops->n_ops == ops->max_ops ) {
		ops->max_ops += 16;
		alloc_ops(ops);
	}

	ops->ops[ops->n_ops++] = m;
}


/* Add a operation to a SymOpList, starting from v(..) */
static void add_symop_v(SymOpList *ops,
                        signed int *h, signed int *k, signed int *l)
{
	IntegerMatrix *m;
	int i;

	m = intmat_new(3, 3);
	assert(m != NULL);

	for ( i=0; i<3; i++ ) intmat_set(m, 0, i, h[i]);
	for ( i=0; i<3; i++ ) intmat_set(m, 1, i, k[i]);
	for ( i=0; i<3; i++ ) intmat_set(m, 2, i, l[i]);

	free(h);
	free(k);
	free(l);

	add_symop(ops, m);
}


static signed int *v(signed int h, signed int k, signed int i, signed int l)
{
	signed int *vec = malloc(3*sizeof(signed int));
	if ( vec == NULL ) return NULL;
	/* Convert back to 3-index form now */
	vec[0] = h-i;  vec[1] = k-i;  vec[2] = l;
	return vec;
}


/**
 * num_equivs:
 * @ops: A %SymOpList
 * @m: A %SymOpMask, which has been shown to special_position()
 *
 * Returns: the number of equivalent reflections for a general reflection
 * in point group "ops", which were not flagged by your call to
 * special_position().
 **/
int num_equivs(const SymOpList *ops, const SymOpMask *m)
{
	int n = num_ops(ops);
	int i;
	int c;

	if ( m == NULL ) return n;

	c = 0;
	for ( i=0; i<n; i++ ) {
		if ( m->mask[i] ) c++;
	}

	return c;
}


static void add_identity(SymOpList *s)
{
	int i, ni;
	int found;

	found = 0;
	ni = num_ops(s);
	for ( i=0; i<ni; i++ ) {
		if ( intmat_is_identity(s->ops[i]) ) {
			found = 1;
			break;
		}
	}
	if ( !found ) {
		add_symop_v(s, v(1,0,0,0), v(0,1,0,0), v(0,0,0,1));  /* I */
	}
}


/* Fill in the other operations for a point group starting from its
 * generators */
static void expand_ops(SymOpList *s)
{
	int added;

	add_identity(s);

	do {

		int i, ni;

		added = 0;

		ni = num_ops(s);
		for ( i=0; i<ni; i++ ) {

			int j;
			IntegerMatrix *opi = s->ops[i];

			/* Apply op 'i' to all the current ops in the list */
			for ( j=0; j<ni; j++ ) {

				IntegerMatrix *opj = s->ops[j];
				IntegerMatrix *m;
				int k, nk;
				int found;

				m = intmat_intmat_mult(opi, opj);
				assert(m != NULL);

				nk = num_ops(s);
				found = 0;
				for ( k=0; k<nk; k++ ) {
					if ( intmat_equals(m, s->ops[k]) ) {
						found = 1;
						intmat_free(m);
						break;
					}
				}

				if ( !found ) {
					add_symop(s, m);
					added++;
				}

			}

		}

	} while ( added );
}


/* Transform all the operations in a SymOpList by a given matrix.
 * The matrix must have a determinant of +/- 1 (otherwise its inverse would
 * not also be an integer matrix). */
static void transform_ops(SymOpList *s, IntegerMatrix *t)
{
	int n, i;
	IntegerMatrix *inv;
	signed int det;

	det = intmat_det(t);
	if ( det == -1 ) {
		ERROR("Warning: mirrored SymOpList.\n");
	} else if ( det != 1 ) {
		ERROR("Invalid transformation for SymOpList.\n");
		return;
	}

	inv = intmat_inverse(t);
	if ( inv == NULL ) {
		ERROR("Failed to invert matrix.\n");
		return;
	}

	n = num_ops(s);
	for ( i=0; i<n; i++ ) {

		IntegerMatrix *r, *f;

		r = intmat_intmat_mult(s->ops[i], t);
		if ( r == NULL ) {
			ERROR("Matrix multiplication failed.\n");
			return;
		}

		f = intmat_intmat_mult(inv, r);
		if ( f == NULL ) {
			ERROR("Matrix multiplication failed.\n");
			return;
		}
		intmat_free(r);

		intmat_free(s->ops[i]);
		s->ops[i] = intmat_copy(f);
		intmat_free(f);

	}

	intmat_free(inv);
}


/********************************* Triclinic **********************************/

static SymOpList *make_1bar()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1));  /* -I */
	new->name = strdup("-1");
	expand_ops(new);
	return new;
}


static SymOpList *make_1()
{
	SymOpList *new = new_symoplist();
	new->name = strdup("1");
	expand_ops(new);
	return new;
}


/********************************* Monoclinic *********************************/

static SymOpList *make_2m()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1)); /* 2 // l */
	add_symop_v(new, v(1,0,0,0), v(0,1,0,0), v(0,0,0,-1));  /* m -| l */
	new->name = strdup("2/m");
	expand_ops(new);
	return new;
}


static SymOpList *make_2()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1)); /* 2 // l */
	new->name = strdup("2");
	expand_ops(new);
	return new;
}


static SymOpList *make_m()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(1,0,0,0), v(0,1,0,0), v(0,0,0,-1));  /* m -| l */
	new->name = strdup("m");
	expand_ops(new);
	return new;
}


/******************************** Orthorhombic ********************************/

static SymOpList *make_mmm()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1)); /* 2 // l */
	add_symop_v(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1)); /* 2 // k */
	add_symop_v(new, v(1,0,0,0), v(0,-1,0,0), v(0,0,0,1));  /* m -| k */
	new->name = strdup("mmm");
	expand_ops(new);
	return new;
}


static SymOpList *make_222()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1)); /* 2 // l */
	add_symop_v(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1)); /* 2 // k */
	new->name = strdup("222");
	expand_ops(new);
	return new;
}


static SymOpList *make_mm2()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1)); /* 2 // l */
	add_symop_v(new, v(1,0,0,0), v(0,-1,0,0), v(0,0,0,1));  /* m -| k */
	new->name = strdup("mm2");
	expand_ops(new);
	return new;
}


/********************************* Tetragonal *********************************/

static SymOpList *make_4m()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,-1,0,0), v(1,0,0,0), v(0,0,0,1)); /* 4 // l */
	add_symop_v(new, v(1,0,0,0), v(0,1,0,0), v(0,0,0,-1)); /* m -| l */
	new->name = strdup("4/m");
	expand_ops(new);
	return new;
}


static SymOpList *make_4()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,-1,0,0), v(1,0,0,0), v(0,0,0,1)); /* 4 // l */
	new->name = strdup("4");
	expand_ops(new);
	return new;
}


static SymOpList *make_4mm()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,-1,0,0), v(1,0,0,0), v(0,0,0,1)); /* 4 // l */
	add_symop_v(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,1)); /* m -| l */
	new->name = strdup("4mm");
	expand_ops(new);
	return new;
}


static SymOpList *make_422()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,-1,0,0), v(1,0,0,0), v(0,0,0,1));  /* 4 // l */
	add_symop_v(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1)); /* 2 // k */
	new->name = strdup("422");
	expand_ops(new);
	return new;
}


static SymOpList *make_4bar()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,1,0,0), v(-1,0,0,0), v(0,0,0,-1)); /* -4 // l */
	new->name = strdup("-4");
	expand_ops(new);
	return new;
}


static SymOpList *make_4bar2m()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,1,0,0), v(-1,0,0,0), v(0,0,0,-1)); /* -4 // l */
	add_symop_v(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1)); /* 2 // k */
	new->name = strdup("-42m");
	expand_ops(new);
	return new;
}


static SymOpList *make_4barm2()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,1,0,0), v(-1,0,0,0), v(0,0,0,-1)); /* -4 // l */
	add_symop_v(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,-1)); /* 2 // h+k */
	new->name = strdup("-4m2");
	expand_ops(new);
	return new;
}


static SymOpList *make_4mmm()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,-1,0,0), v(1,0,0,0), v(0,0,0,1)); /* 4 // l */
	add_symop_v(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,1)); /* m -| k */
	add_symop_v(new, v(1,0,0,0), v(0,1,0,0), v(0,0,0,-1)); /* m -| l */
	new->name = strdup("4/mmm");
	expand_ops(new);
	return new;
}


/************************** Trigonal (Rhombohedral) ***************************/

static SymOpList *make_3_R()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,0,0,1), v(1,0,0,0), v(0,1,0,0)); /* 3 // h+k+l */
	new->name = strdup("3_R");
	expand_ops(new);
	return new;
}


static SymOpList *make_3bar_R()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,0,0,1), v(1,0,0,0), v(0,1,0,0)); /* -3 // h+k+l */
	add_symop_v(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1)); /* -I */
	new->name = strdup("-3_R");
	expand_ops(new);
	return new;
}


static SymOpList *make_32_R()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,0,0,1), v(1,0,0,0), v(0,1,0,0)); /* 3 // h+k+l */
	add_symop_v(new, v(0,-1,0,0), v(-1,0,0,0), v(0,0,0,-1)); /* 2 -| 3 */
	new->name = strdup("32_R");
	expand_ops(new);
	return new;
}


static SymOpList *make_3m_R()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,0,0,1), v(1,0,0,0), v(0,1,0,0)); /* 3 // h+k+l */
	add_symop_v(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,1)); /* m */
	new->name = strdup("3m_R");
	expand_ops(new);
	return new;
}


static SymOpList *make_3barm_R()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,0,0,1), v(1,0,0,0), v(0,1,0,0)); /* -3 // h+k+l */
	add_symop_v(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1)); /* -I */
	add_symop_v(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,1));    /* m */
	new->name = strdup("-3m_R");
	expand_ops(new);
	return new;
}


/*************************** Trigonal (Hexagonal) *****************************/

static SymOpList *make_3_H()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,1)); /* 3 // l */
	new->name = strdup("3_H");
	expand_ops(new);
	return new;
}


static SymOpList *make_3bar_H()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,1));    /* 3 // l */
	add_symop_v(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1)); /* -I */
	new->name = strdup("-3_H");
	expand_ops(new);
	return new;
}


static SymOpList *make_321_H()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,1));  /* 3 // l */
	add_symop_v(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,-1)); /* 2 // h */
	new->name = strdup("321_H");
	expand_ops(new);
	return new;
}


static SymOpList *make_312_H()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,1));    /* 3 // l */
	add_symop_v(new, v(0,-1,0,0), v(-1,0,0,0), v(0,0,0,-1)); /* 2 // h+k */
	new->name = strdup("312_H");
	expand_ops(new);
	return new;
}


static SymOpList *make_3m1_H()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,1)); /* 3 // l */
	add_symop_v(new, v(0,-1,0,0), v(-1,0,0,0), v(0,0,0,1)); /* m -| i */
	new->name = strdup("3m1_H");
	expand_ops(new);
	return new;
}


static SymOpList *make_31m_H()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,1)); /* 3 // l */
	add_symop_v(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,1)); /* m -| (k+i) */
	new->name = strdup("31m_H");
	expand_ops(new);
	return new;
}


static SymOpList *make_3barm1_H()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,1));    /* 3 // l */
	add_symop_v(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1)); /* -I */
	add_symop_v(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,-1));   /* 2 // h */
	new->name = strdup("-3m1_H");
	expand_ops(new);
	return new;
}


static SymOpList *make_3bar1m_H()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,1));    /* 3 // l */
	add_symop_v(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1)); /* -I */
	add_symop_v(new, v(0,-1,0,0), v(-1,0,0,0), v(0,0,0,-1)); /* 2 // h+k */
	new->name = strdup("-31m_H");
	expand_ops(new);
	return new;
}


/********************************** Hexgonal **********************************/

static SymOpList *make_6()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,0,-1,0), v(-1,0,0,0), v(0,0,0,1)); /* 6 // l */
	new->name = strdup("6");
	expand_ops(new);
	return new;
}


static SymOpList *make_6bar()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,-1)); /* -6 // l */
	new->name = strdup("-6");
	expand_ops(new);
	return new;
}


static SymOpList *make_6m()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,0,-1,0), v(-1,0,0,0), v(0,0,0,1)); /* 6 // l */
	add_symop_v(new, v(1,0,0,0), v(0,1,0,0), v(0,0,0,-1));  /* m -| l */
	new->name = strdup("6/m");
	expand_ops(new);
	return new;
}


static SymOpList *make_622()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,0,-1,0), v(-1,0,0,0), v(0,0,0,1)); /* 6 // l */
	add_symop_v(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,-1));   /* 2 // h */
	new->name = strdup("622");
	expand_ops(new);
	return new;
}


static SymOpList *make_6mm()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,0,-1,0), v(-1,0,0,0), v(0,0,0,1)); /* 6 // l */
	add_symop_v(new, v(0,-1,0,0), v(-1,0,0,0), v(0,0,0,1)); /* m -| i */
	new->name = strdup("6mm");
	expand_ops(new);
	return new;
}


static SymOpList *make_6barm2()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,-1)); /* -6 // l */
	add_symop_v(new, v(0,-1,0,0), v(-1,0,0,0), v(0,0,0,1)); /* m -| i */
	new->name = strdup("-6m2");
	expand_ops(new);
	return new;
}


static SymOpList *make_6bar2m()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,-1)); /* -6 // l */
	add_symop_v(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,1));  /* m -| (k+i) */
	new->name = strdup("-62m");
	expand_ops(new);
	return new;
}


static SymOpList *make_6mmm()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,-1)); /* -6 // l */
	add_symop_v(new, v(0,-1,0,0), v(-1,0,0,0), v(0,0,0,1)); /* m -| i */
	add_symop_v(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1)); /* -I */
	new->name = strdup("6/mmm");
	expand_ops(new);
	return new;
}


/************************************ Cubic ***********************************/

static SymOpList *make_23()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1)); /* 2 // l */
	add_symop_v(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1)); /* 2 // k */
	add_symop_v(new, v(0,1,0,0), v(0,0,0,1), v(1,0,0,0)); /* 3 // h+k+l */
	new->name = strdup("23");
	expand_ops(new);
	return new;
}


static SymOpList *make_m3bar()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1)); /* 2 // l */
	add_symop_v(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1)); /* 2 // k */
	add_symop_v(new, v(0,1,0,0), v(0,0,0,1), v(1,0,0,0)); /* 3 // h+k+l */
	add_symop_v(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1)); /* -I */
	new->name = strdup("m-3");
	expand_ops(new);
	return new;
}


static SymOpList *make_432()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,-1,0,0), v(1,0,0,0), v(0,0,0,1)); /* 4 // l */
	add_symop_v(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1));/* 2 // k */
	add_symop_v(new, v(0,1,0,0), v(0,0,0,1), v(1,0,0,0));  /* 3 // h+k+l */
	new->name = strdup("432");
	expand_ops(new);
	return new;
}


static SymOpList *make_4bar3m()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,1,0,0), v(-1,0,0,0), v(0,0,0,-1)); /* -4 // l */
	add_symop_v(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1)); /* 2 // k */
	add_symop_v(new, v(0,1,0,0), v(0,0,0,1), v(1,0,0,0));   /* 3 // h+k+l */
	new->name = strdup("-43m");
	expand_ops(new);
	return new;
}


static SymOpList *make_m3barm()
{
	SymOpList *new = new_symoplist();
	add_symop_v(new, v(0,-1,0,0), v(1,0,0,0), v(0,0,0,1)); /* 4 // l */
	add_symop_v(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1));/* 2 // k */
	add_symop_v(new, v(0,1,0,0), v(0,0,0,1), v(1,0,0,0));  /* 3 // h+k+l */
	add_symop_v(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1)); /* -I */
	new->name = strdup("m-3m");
	expand_ops(new);
	return new;
}


static SymOpList *getpg_uac(const char *sym)
{
	/* Triclinic */
	if ( strcmp(sym, "-1") == 0 ) return make_1bar();
	if ( strcmp(sym, "1") == 0 ) return make_1();

	/* Monoclinic */
	if ( strcmp(sym, "2/m") == 0 ) return make_2m();
	if ( strcmp(sym, "2") == 0 ) return make_2();
	if ( strcmp(sym, "m") == 0 ) return make_m();

	/* Orthorhombic */
	if ( strcmp(sym, "mmm") == 0 ) return make_mmm();
	if ( strcmp(sym, "222") == 0 ) return make_222();
	if ( strcmp(sym, "mm2") == 0 ) return make_mm2();

	/* Tetragonal */
	if ( strcmp(sym, "4/m") == 0 ) return make_4m();
	if ( strcmp(sym, "4") == 0 ) return make_4();
	if ( strcmp(sym, "-4") == 0 ) return make_4bar();
	if ( strcmp(sym, "4/mmm") == 0 ) return make_4mmm();
	if ( strcmp(sym, "422") == 0 ) return make_422();
	if ( strcmp(sym, "-42m") == 0 ) return make_4bar2m();
	if ( strcmp(sym, "-4m2") == 0 ) return make_4barm2();
	if ( strcmp(sym, "4mm") == 0 ) return make_4mm();

	/* Trigonal (rhombohedral) */
	if ( strcmp(sym, "3_R") == 0 ) return make_3_R();
	if ( strcmp(sym, "-3_R") == 0 ) return make_3bar_R();
	if ( strcmp(sym, "32_R") == 0 ) return make_32_R();
	if ( strcmp(sym, "3m_R") == 0 ) return make_3m_R();
	if ( strcmp(sym, "-3m_R") == 0 ) return make_3barm_R();

	/* Trigonal (hexagonal) */
	if ( strcmp(sym, "3_H") == 0 ) return make_3_H();
	if ( strcmp(sym, "-3_H") == 0 ) return make_3bar_H();
	if ( strcmp(sym, "321_H") == 0 ) return make_321_H();
	if ( strcmp(sym, "312_H") == 0 ) return make_312_H();
	if ( strcmp(sym, "3m1_H") == 0 ) return make_3m1_H();
	if ( strcmp(sym, "31m_H") == 0 ) return make_31m_H();
	if ( strcmp(sym, "-3m1_H") == 0 ) return make_3barm1_H();
	if ( strcmp(sym, "-31m_H") == 0 ) return make_3bar1m_H();

	/* Hexagonal */
	if ( strcmp(sym, "6/m") == 0 ) return make_6m();
	if ( strcmp(sym, "6") == 0 ) return make_6();
	if ( strcmp(sym, "-6") == 0 ) return make_6bar();
	if ( strcmp(sym, "6/mmm") == 0 ) return make_6mmm();
	if ( strcmp(sym, "622") == 0 ) return make_622();
	if ( strcmp(sym, "-62m") == 0 ) return make_6bar2m();
	if ( strcmp(sym, "-6m2") == 0 ) return make_6barm2();
	if ( strcmp(sym, "6mm") == 0 ) return make_6mm();

	/* Cubic */
	if ( strcmp(sym, "23") == 0 ) return make_23();
	if ( strcmp(sym, "m-3") == 0 ) return make_m3bar();
	if ( strcmp(sym, "432") == 0 ) return make_432();
	if ( strcmp(sym, "-43m") == 0 ) return make_4bar3m();
	if ( strcmp(sym, "m-3m") == 0 ) return make_m3barm();

	ERROR("Unknown point group '%s'\n", sym);
	return NULL;
}


static int char_count(const char *a, char b)
{
	size_t i;
	int n;

	i = 0;  n = 0;
	do {
		if ( a[i] == b ) n++;
		if ( a[i] == '\0' ) return n;
		i++;
	} while ( 1 );
}


static SymOpList *getpg_arbitrary_ua(const char *sym, size_t s)
{
	char ua;
	char *pg_type;
	SymOpList *pg;
	IntegerMatrix *t;

	if ( strncmp(sym+s, "ua", 2) == 0 ) {
		ua = sym[s+2];
	} else {
		ERROR("Unrecognised point group '%s'\n", sym);
		return NULL;
	}

	pg_type = strndup(sym, s-1);
	if ( pg_type == NULL ) {
		ERROR("Couldn't allocate string.\n");
		return NULL;
	}

	pg = getpg_uac(pg_type);
	if ( pg == NULL ) {
		ERROR("Unrecognised point group type '%s'\n",
		      pg_type);
		return NULL;
	}
	free(pg_type);

	t = intmat_new(3, 3);
	if ( t == NULL ) return NULL;