Stuff related to calculating anomalous structure factors

4 files changed, 115 insertions, 38 deletions

diff --git a/scripts/eV-to-A b/scripts/eV-to-A
new file mode 100755
index 00000000..941e422b
--- /dev/null
+++ b/scripts/eV-to-A
@@ -0,0 +1,10 @@
+#!/usr/bin/perl
+
+my $eV = $ARGV[0];
+
+printf("%5.2f eV\n", $eV);
+
+my $J = $eV*1.6021773e-19;
+my $lambda = (6.62606896e-34 * 299792458)/$J;
+
+printf("%5.2f eV = %5.5f A\n", $eV, $lambda*1e10);
diff --git a/scripts/gen-sfs-ano b/scripts/gen-sfs-ano
index be4598d6..522e0a6d 100755
--- a/scripts/gen-sfs-ano
+++ b/scripts/gen-sfs-ano
@@ -2,16 +2,27 @@
 
 PDB=$1
 SYMM=$2
-WAVEL=$3
-RESOLUTION=$4
+PG=$3
+WAVEL=$4
+RESOLUTION=$5
+
+function show_help()
+{
+	echo -n "Syntax: $0 <PDB file> <space group> <point group>"
+	echo "<wavelength in Angstroms> [<resolution>]"
+	echo
+	echo "The space group and point group must be consistent, it's just"
+	echo "that I don't know how to convert the space group to a point"
+	echo "group."
+}
 
 if [ "x$PDB" = "x" ]; then
-	echo "Syntax: $0 <PDB file> <space group> <wavelength in Angstroms> [<resolution>]"
+	show_help
 	exit
 fi
 
 if [ "x$SYMM" = "x" ]; then
-	echo "Syntax: $0 <PDB file> <space group> <wavelength in Anstroms> [<resolution>]"
+	show_help
 	exit
 fi
 
@@ -25,17 +36,17 @@ ano_sfall.com ${PDB} ${RESOLUTION}A wave=${WAVEL}
 if [ $? -ne 0 ]; then exit 1; fi
 
 echo "Converting structure factors to text..."
-mtz2various hklin ideal_ano.mtz hklout ${PDB}-temp.hkl >> gen-sfs.html <<EOF
+mtz2various hklin ideal_ano.mtz hklout ${PDB}-temp1.hkl >> gen-sfs.html <<EOF
 LABIN H=H K=K L=L DUM1=Fplus DUM2=Fminus
 OUTPUT USER '(3I4,2F9.1)'
 EOF
 if [ $? -ne 0 ]; then exit 1; fi
 rm -f ${PDB}.mtz
-perl < ${PDB}-temp.hkl > ${PDB}.hkl << WIBBLE
+perl > ${PDB}-temp2.hkl << WIBBLE
 use strict;
 
 my \$line;
-open(FILE, "${PDB}-temp.hkl");
+open(FILE, "${PDB}-temp1.hkl");
 
 printf("  h   k   l          I    phase   sigma(I)  1/d(nm^-1)  ".
        "counts  fs/px  ss/px\\n");
@@ -51,13 +62,10 @@ while ( \$line = <FILE> ) {
                 my \$iminus = \$5*\$5;
 
                 printf("%3i %3i %3i %10.2f -     0.0  %s %7i %6.1f %6.1f\n",
-                       \$h, \$k, \$l, \$iplus, 0.0,
-                       "         -", 1, 0.0, 0.0);
+                       \$h, \$k, \$l, \$iplus, "         -", 1, 0.0, 0.0);
 
                 printf("%3i %3i %3i %10.2f -     0.0  %s %7i %6.1f %6.1f\n",
-                       -\$h, -\$k, -\$l, \$iminus, 0.0,
-                       "         -", 1, 0.0, 0.0);
-
+                       -\$h, -\$k, -\$l, \$iminus, "         -", 1, 0.0, 0.0);
 
         } else {
         	printf(STDERR "Couldn't understand line '%s'\n", \$line);
@@ -67,6 +75,8 @@ while ( \$line = <FILE> ) {
 close(FILE);
 printf("End of reflections\n");
 WIBBLE
-exit
 
-rm -f ${PDB}-temp.hkl
+get_hkl -i ${PDB}-temp2.hkl -o ${PDB}.hkl -p ${PDB} -y ${PG} --trim-centrics
+
+rm -f ${PDB}-temp1.hkl
+rm -f ${PDB}-temp2.hkl
diff --git a/src/get_hkl.c b/src/get_hkl.c
index cb0446d8..e5e1a180 100644
--- a/src/get_hkl.c
+++ b/src/get_hkl.c
@@ -51,6 +51,11 @@ static void show_help(const char *s)
 "                              point group.\n"
 "  -e, --expand=<sym>         Expand reflections to this point group.\n"
 "\n"
+"Use this option with care, and only if you understand why it might sometimes\n"
+" be necessary:\n"
+"  --trim-centrics            Remove reflections which are duplicated in the\n"
+"                              point group specified with the '-y' option.\n"
+"\n"
 "You can restrict which reflections are written out:\n"
 "  -t, --template=<filename>  Only include reflections mentioned in file.\n"
 "\n"
@@ -294,12 +299,49 @@ static RefList *expand_reflections(RefList *in, const SymOpList *target,
 }
 
 
+static RefList *trim_centrics(RefList *in, const SymOpList *sym)
+{
+	Reflection *refl;
+	RefListIterator *iter;
+	RefList *out;
+
+	out = reflist_new();
+
+	for ( refl = first_refl(in, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) )
+	{
+		signed int h, k, l;
+		signed int ha, ka, la;
+		Reflection *new;
+
+		get_indices(refl, &h, &k, &l);
+
+		/* Put it into the asymmetric unit */
+		get_asymm(sym, h, k, l, &ha, &ka, &la);
+
+		new = find_refl(out, ha, ka, la);
+		if ( new != NULL ) {
+			STATUS("Trimmed %i %i %i\n", h, k, l);
+			continue;
+		}
+
+		/* Add new reflection under asymmetric (unique) indices */
+		new = add_refl(out, ha, ka, la);
+		copy_data(new, refl);
+	}
+
+	return out;
+}
+
+
 int main(int argc, char *argv[])
 {
 	int c;
 	int config_noise = 0;
 	int config_poisson = 0;
 	int config_multi = 0;
+	int config_trimc = 0;
 	char *holo_str = NULL;
 	char *mero_str = NULL;
 	char *expand_str = NULL;
@@ -328,6 +370,7 @@ int main(int argc, char *argv[])
 		{"multiplicity",       0, &config_multi,       1},
 		{"beam",               1, NULL,               'b'},
 		{"pdb",                1, NULL,               'p'},
+		{"trim-centrics",      0, &config_trimc,       1},
 		{0, 0, NULL, 0}
 	};
 
@@ -430,7 +473,10 @@ int main(int argc, char *argv[])
 		return 1;
 	}
 	free(input_file);
-	if ( check_list_symmetry(input, mero) ) {
+
+	STATUS("%i reflections in input.\n", num_reflections(input));
+
+	if ( !config_trimc && check_list_symmetry(input, mero) ) {
 		ERROR("The input reflection list does not appear to"
 		      " have symmetry %s\n", symmetry_name(mero));
 		return 1;
@@ -478,6 +524,19 @@ int main(int argc, char *argv[])
 
 	}
 
+	if ( config_trimc ) {
+
+		RefList *new;
+
+		STATUS("Trimming duplicate reflections in %s\n",
+		       symmetry_name(mero));
+		new = trim_centrics(input, mero);
+		reflist_free(input);
+		input = new;
+		STATUS("%i output reflections\n", num_reflections(input));
+
+	}
+
 	if ( config_multi ) {
 
 		Reflection *refl;
diff --git a/src/render_hkl.c b/src/render_hkl.c
index b37cc00b..4c37c005 100644
--- a/src/render_hkl.c
+++ b/src/render_hkl.c
@@ -559,9 +559,9 @@ int main(int argc, char *argv[])
 	UnitCell *cell;
 	RefList *list;
 	char *infile;
-	int config_zoneaxis = 0;
 	int config_sqrt = 0;
 	int config_colkey = 0;
+	int config_zawhinge = 0;
 	unsigned int nproc = 1;
 	char *pdb = NULL;
 	int r = 0;
@@ -583,7 +583,7 @@ int main(int argc, char *argv[])
 	/* Long options */
 	const struct option longopts[] = {
 		{"help",               0, NULL,               'h'},
-		{"zone-axis",          0, &config_zoneaxis,    1},
+		{"zone-axis",          0, &config_zawhinge,    1},
 		{"output",             1, NULL,               'o'},
 		{"pdb",                1, NULL,               'p'},
 		{"boost",              1, NULL,               'b'},
@@ -667,7 +667,7 @@ int main(int argc, char *argv[])
 
 	}
 
-	if ( config_zoneaxis ) {
+	if ( config_zawhinge ) {
 		ERROR("Friendly warning: The --zone-axis option isn't needed"
 		      " any longer (I ignored it for you).\n");
 	}
@@ -733,28 +733,26 @@ int main(int argc, char *argv[])
 		return render_key(colscale, scale_top);
 	}
 
-	if ( config_zoneaxis ) {
-		if ( (( down == NULL ) && ( right != NULL ))
-		  || (( down != NULL ) && ( right == NULL )) ) {
-			ERROR("Either specify both 'down' and 'right',"
-			      " or neither.\n");
+	if ( (( down == NULL ) && ( right != NULL ))
+	  || (( down != NULL ) && ( right == NULL )) ) {
+		ERROR("Either specify both 'down' and 'right',"
+		      " or neither.\n");
+		return 1;
+	}
+	if ( down != NULL ) {
+		int r;
+		r = sscanf(down, "%i,%i,%i", &dh, &dk, &dl);
+		if ( r != 3 ) {
+			ERROR("Invalid format for 'down'\n");
 			return 1;
 		}
-		if ( down != NULL ) {
-			int r;
-			r = sscanf(down, "%i,%i,%i", &dh, &dk, &dl);
-			if ( r != 3 ) {
-				ERROR("Invalid format for 'down'\n");
-				return 1;
-			}
-		}
-		if ( right != NULL ) {
-			int r;
-			r = sscanf(right, "%i,%i,%i", &rh, &rk, &rl);
-			if ( r != 3 ) {
-				ERROR("Invalid format for 'right'\n");
-				return 1;
-			}
+	}
+	if ( right != NULL ) {
+		int r;
+		r = sscanf(right, "%i,%i,%i", &rh, &rk, &rl);
+		if ( r != 3 ) {
+			ERROR("Invalid format for 'right'\n");
+			return 1;
 		}
 	}