aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2011-11-18 16:36:43 +0100
committerThomas White <taw@physics.org>2012-02-22 15:27:41 +0100
commitbd640fa2bbb6bbb219e8dbcdaf88df093cac5c8a (patch)
tree002c2a6d1d5c010f2c6c74975060a5914eee434b
parentd5ddfec8d79c9fc829a32cdf51082f626d3184bb (diff)
Stuff related to calculating anomalous structure factors
-rwxr-xr-xscripts/eV-to-A10
-rwxr-xr-xscripts/gen-sfs-ano38
-rw-r--r--src/get_hkl.c61
-rw-r--r--src/render_hkl.c44
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;
}
}