Add scripts/gen-sfs for calculating scattering

2 files changed, 85 insertions, 5 deletions

diff --git a/doc/man/pattern_sim.1 b/doc/man/pattern_sim.1
index b9c8fa33..275c25cd 100644
--- a/doc/man/pattern_sim.1
+++ b/doc/man/pattern_sim.1
@@ -26,10 +26,10 @@ The result will be written to an HDF5 file in the current directory with the nam
 
 .SH REFLECTION LISTS
 
-pattern_sim does not know about symmetry, so your input reflection list
-(give with "-i") must be expanded.  You can do this with:
+You'll need to create a file containing the intensities of the reflections.  The normal way to do this is to use CCP4 via the "gen-sfs" script in CrystFEL's script folder.  Run it like this:
 
-$ get_hkl -i myfile.hkl -o output.hkl -y mypointgroup -e 1
+$ gen-sfs mymodel.pdb "P6" 3
 
-Please be sure to read the "Note about Unit Cell Settings" in the documentation
-for indexamajig.
+You need to give the PDB model, the symmetry of the output reflections (use the lowest symmetry space group with the right point group), and optionally the maximum resolution in Angstroms.  If you don't specify the resolution, it'll use 3 Angstroms.
+
+The reflections will be output as "mymodel.pdb.hkl" ready for input to pattern_sim.  You'll need to give the Laue class of the symmetry you gave to gen-sfs, "6/m" in this case, to pattern_sim with the "-y" option.  By default, gen-sfs calculates the values for CuKa radiation (8.3 keV, 1.5 A).  It will not calculate the anomalous contribution to scattering, i.e. the differences in intensities between Bijoet pairs.  Both of these are the default behaviour for "sfall" in CCP4, so read the manual for that for further details.  If you need something different, get the "ano_sfall.com" script from James Holton and turn the output into a CrystFEL reflection list in a similar way.
diff --git a/scripts/gen-sfs b/scripts/gen-sfs
new file mode 100755
index 00000000..bf86fb01
--- /dev/null
+++ b/scripts/gen-sfs
@@ -0,0 +1,80 @@
+#!/bin/sh
+
+PDB=$1
+SYMM=$2
+RESOLUTION=$3
+
+if [ "x$PDB" = "x" ]; then
+	echo "Syntax: $0 <PDB file> <space group> [<resolution>]"
+	exit
+fi
+
+if [ "x$SYMM" = "x" ]; then
+	echo "Syntax: $0 <PDB file> <space group> [<resolution>]"
+	exit
+fi
+
+if [ "x$RESOLUTION" = "x" ]; then
+	echo "Resolution not given.  Using 3 Angstroms."
+	RESOLUTION=3
+fi
+
+echo "Running sfall to calculate structure factors..."
+sfall XYZIN $PDB HKLOUT ${PDB}.mtz > gen-sfs.html << EOF
+MODE SFCALC XYZIN
+RESOLUTION $RESOLUTION 1000
+FORM NGAUSS 5
+SYMM $SYMM
+END
+EOF
+if [ $? -ne 0 ]; then exit 1; fi
+
+echo "Running cad to get the right asymmetric unit..."
+cad HKLIN1 ${PDB}.mtz HKLOUT ${PDB}-sorted.mtz >> gen-sfs.html <<EOF
+TITLE Sorted blah
+LABIN FILE 1 E1=FC E2=PHIC
+CTYPE FILE 1 E1=F  E2=P
+EOF
+if [ $? -ne 0 ]; then exit 1; fi
+
+echo "Converting structure factors to text..."
+mtz2various hklin ${PDB}-sorted.mtz hklout ${PDB}-temp.hkl >> gen-sfs.html <<EOF
+LABIN H=H K=K L=L FC=FC PHIC=PHIC
+OUTPUT USER '(3I4,2F9.1)'
+EOF
+if [ $? -ne 0 ]; then exit 1; fi
+rm -f ${PDB}.mtz
+perl < ${PDB}-temp.hkl > ${PDB}.hkl << WIBBLE
+use strict;
+
+my \$line;
+open(FILE, "${PDB}-temp.hkl");
+
+printf("  h   k   l          I    phase   sigma(I)  1/d(nm^-1)  ".
+       "counts  fs/px  ss/px\\n");
+
+while ( \$line = <FILE> ) {
+
+        if ( \$line =~ /^\s*([\d\-]+)\s+([\d\-]+)\s+([\d\-]+)\s+([\d\-\.]+)\s+([\d\-\.]+)/ ) {
+
+                my \$h = \$1;
+                my \$k = \$2;
+                my \$l = \$3;
+                my \$intensity = \$4*\$4;  # Square to convert F->I
+                my \$phase = \$5;
+
+                printf("%3i %3i %3i %10.2f %8.2f %10.2f  %s %7i %6.1f %6.1f\n",
+                       \$h, \$k, \$l, \$intensity, \$phase, 0.0,
+                       "         -", 1, 0.0, 0.0);
+
+        } else {
+        	printf(STDERR "Couldn't understand line '%s'\n", \$line);
+        }
+
+}
+close(FILE);
+printf("End of reflections\n");
+WIBBLE
+exit
+
+rm -f ${PDB}-temp.hkl