Hi Robbie,

I contacted Charles Ballard a few weeks ago on the very same issue and
he has efficiently rolled out a fix for it in the latest CCP4 update
:-)

Update ID: 6.4.0-001
Date: 25.10.13

Cheers,
Jon

On 4 December 2013 09:23, Robbie Joosten <[log in to unmask]> wrote:
> Hi Ian,
>
>
>> This probably doesn't affect the MTZ file that comes out of it
>> but it will affect the statistics of the twinning tests.
>
> I'd be interested to see if this solves the cases where ctruncate goes into
> an infinite loop during the twinning tests.
>
> Cheers,
> Robbie
>
> ________________________________
> Date: Sat, 30 Nov 2013 15:37:44 +0000
> From: [log in to unmask]
> Subject: [ccp4bb] TRUNCATE & CTRUNCATE issues.
> To: [log in to unmask]
>
>
>
> Hello All, I'd like to get your views on some changes I'm proposing to make
> to the TRUNCATE source code.  IMO there are some issues with the way
> TRUNCATE does its statistical analyses which need to be fixed.  This
> probably doesn't affect the MTZ file that comes out of it but it will affect
> the statistics of the twinning tests.  I've been meaning to do this for some
> time but it involves some fairly radical changes; now I've finally decided
> to bite the bullet.
>
> One problem is that there seems to be some confusion in the source code
> comments concerning the meaning of what I call the "symmetry enhancement
> factor" (aka "epsilon" or "e" below). This is the point-group dependent
> factor by which the mean intensities of special rows and zones are enhanced
> by symmetry; for example in PG121 the mean I of the 0k0 reflections is
> enhanced by a factor of 2; in PG6 & PG622, it's a factor of 6 for the 00l
> reflections, and so on.
>
> For space groups with screw axes (or glide planes in enantiomorphic SGs,
> i.e. the crystal contains both the asymmetric unit and its mirror image),
> the mean is multiplied by e only if the systematic absences are omitted from
> the sums; if you include them then there is no overall enhancement.  So
> rotation and screw axes (or mirrors and glide planes) need to be treated
> differently, yet the CCP4 library code for this (s/rs EPSLN & EPSLON) treat
> them identically (more on this below).  TRUNCATE calls epsilon alternatively
> a "multiplicity" or "weight" which further adds to the confusion since
> multiplicity (m) usually means something completely different (it's the
> number of times a symmetry-equivalent reflection occurs in a full hemisphere
> of data, so for PG222 axial reflections (h00 etc.) m=1, for zero layers (hk0
> etc.) m=2, and for general hkl m=4: in contrast e = 2, 1 and 1 resp.).
>
> This snippet of code from TRUNCATE which accumulates sums in bins according
> to d*^2 illustrates the problem:
>
>       CALL EPSLON(INHKL,WEIGHT,ISYSAB)
> C
>       FF(NT) = FF(NT) + F/WEIGHT
>       SD(NT) = SD(NT) + SS/WEIGHT
>       N(NT) = N(NT) + 1
>       AMULT = NSYM/WEIGHT
> C
> ...
> C     Accumulate sums for Wilson plot.
> C
>       SN(NT) = SN(NT) + AMULT
>       SW(NT) = SW(NT) + AMULT*FFSCAT
>       SR(NT) = SR(NT) + AMULT*Q
>       SI(NT) = SI(NT) + AMULT*F
>
> Here WEIGHT = epsilon (the ISYSAB flag is ignored throughout), so all sums
> are being accumulated with the terms multiplied by 1/e (NSYM is the no. of
> asymmetric units so is constant and doesn't affect the results).  However
> IMO this factor should be applied only to individual intensities or their
> SDs (note that in the above code and that below, F is actually the
> uncorrected intensity just to further confuse you!).  The problem here is
> that for pure rotation axes the law of conservation of energy requires that
> the overall mean I is unchanged.  In any interference phenomenon energy can
> neither be created nor destroyed, merely transferred from one place to
> another.  So what happens is that the enhanced intensity of the axial
> reflections has been transferred from neighbouring reflections which have
> their intensities diminished in total by the same amount.  In fact an
> oscillating Bessel function centred on the axis is superimposed on the
> intensities so you get cylindical zones of alternating enhanced and
> diminishing intensities with the magnitude of the oscillations dying away as
> you go further from the axis.  However the energy conservation law requires
> that the net overall average I must be unchanged by the presence of the
> axis.
>
> This implies that the 1/e correction factor SHOULD NOT be included for pure
> rotation axes (and mirror planes) when summing for the mean I.  However the
> sums should be performed over a complete hemisphere given only one symmetry
> equivalent per reflection, which means that the multiplicity (m) factor
> SHOULD be included.  This is the direct opposite of what the code above is
> doing (i.e. it includes e but not m!).  The Fortran code would look like:
>
>        SI(NT) = SI(NT) + M*FI     # Total energy is conserved!
>        SN(NT) = SN(NT) + M       # Count reflections in hemisphere.
>
> Note that the statistically valid procedure will be different if one is say
> summing Is for a likelihood function, since this requires that the terms are
> statistically independent so one would then add only one term per
> equivalent, not a complete hemisphere, i.e. in that case the multiplicity
> factors should be omitted.
>
> For screw axes (and glide planes) the situation is different: there is no
> Bessel function and only the axial reflections are affected, so therefore it
> requires different handling in the code (as I said above, this is not
> happening!).  Now, since systematic absences are normally not present in the
> data, the mean intensity of the remaining reflections is enhanced by the e
> factor, purely by the action of omitting the systematic absences of zero I.
> This implies that we need to simulate the presence of the systematic
> absences when taking the mean.  So for example in the PG6 case we would have
> to sum the intensities of the 00l, l=6n reflections WITHOUT correction, but
> then count each 00l reflection as though it were 6 reflections (i.e. also
> counting the omitted sys. abs. in the average), so the code would now look
> like:
>
>        SI(NT) = SI(NT) + M*FI     # Total energy is still conserved!
>        SN(NT) = SN(NT) +M*E    # but also count the sys. abs. that were
> omitted.
>
> This implies that at least for the reflection counts, the e correction
> factor SHOULD be included for screw axes (and glide planes), and for the
> same reason as above the m factor SHOULD also be included.
>
> The differences between rotation and screw axes (or between mirrors and
> glides) arise because Wilson's assumption of uniform random distribution of
> atoms breaks down in the former case: an atom cannot approach a rotation
> axis or mirror plane closer than its VDW radius, so this excluded zone along
> the axis or plane causes an interference effect.  In both cases the main
> effect is actually that in projection along the axis the atomic positions
> are not random: they are correlated by an apparent inversion centre (it's
> actually another interference effect this time from pairs of atoms lying in
> the same plane of reflection and related by symmetry).
>
> All the above is relevant ONLY to taking the average intensity.  For other
> purposes the correct procedure is likely to differ.  For example if one is
> interested in the individual normalised structure amplitudes for direct
> methods (i.e. not just the mean), then the sqrt(1/e) factor clearly SHOULD
> be applied to the individual amplitudes.  Also if one is calculating higher
> moments of Z (= normalised I) then the deviations will not cancel (there's
> nothing in the energy conservation law that says that energy^n is conserved
> if n is not 0 or 1).  In the rotation axis case each large on-axis positive
> deviation tends to be offset by several small off-axis deviations, so in
> that case the optimal procedure would appear to be to multiply the on-axis
> Is by m/e before summing for the higher moments (say n >= 2), e.g.:
>
>       SM(N) = SM(N) + M*(FI/E)**N
>       NM(N) = NM(N) + M
>
> For the screw-axis case with the sys. abs. omitted the situation is again
> different and requires different code; IMO we should distribute an amount
> I/e from each axial reflection equally among the omitted sys. abs. and then
> include them as though they were all present:
>
>       SM(N) = SM(N) + M*E*(FI/E)**N
>       NM(N) = NM(N) + M*E
>
> One further issue needs to be addressed (TRUNCATE is not short of issues -
> the ones I mention here are only a fraction!).  The code for calculating the
> moments in TRUNCATE is:
>
> C---- Sums for moments of I
>          F = FFA(1)/WEIGHT
>          IF(F.GT.0.0) SMMEMM(1,NT,ICEN) = SMMEMM(1,NT,ICEN) + sqrt(F)
>          IF(F.GT.0.0) SMMEMM(3,NT,ICEN) = SMMEMM(3,NT,ICEN) + F*sqrt(F)
>          SMMIMM(1,NT,ICEN) = SMMIMM(1,NT,ICEN) + F
>          SMMIMM(2,NT,ICEN) = SMMIMM(2,NT,ICEN) + F**2
>          SMMIMM(3,NT,ICEN) = SMMIMM(3,NT,ICEN) + F**3
>          SMMIMM(4,NT,ICEN) = SMMIMM(4,NT,ICEN) + F**4
>          NMMNUM(NT,ICEN) = NMMNUM(NT,ICEN) + 1
>
> Notwithstanding that the WEIGHT factor (aka epsilon) probably should not be
> applied to the lower moments, and the multiplicity factor probably should be
> applied everywhere, there is still a problem here.  F is the uncorrected
> intensity but IMO we should be using the corrected intensity (i.e. by F &
> W's Bayesian procedure): why else are we calculating the correction if not
> to apply it?  The uncorrected intensities may be negative which adds
> spurious noise to the moments (a negative intensity makes the same
> contribution to an even moment as an equal positive intensity: this cannot
> possibly be correct).  However using the corrected Is brings another
> problem: when calculating the moments we should be using the expected values
> based on the posterior probability distribution of the Is.  This means
> integrating the moments over all expected Is > 0, multiplied by the
> probability density.  This is non-trivial, however there is a way using the
> parabolic cylinder functions U(a,x)
> (http://en.wikipedia.org/wiki/Parabolic_cylinder_functions).  The nice thing
> about this method is one can generalise it to calculate the
> Bayesian-corrected F and I as well as arbitrary moments and get rid of F &
> W's cubic spline interpolation code with the slightly worrying warning (see
> NEGIAS & NEGICS s/rs in TRUNCATE):
>
> C---- Accuracy - better than 5 percent  (i think).
>
> With the PCFs the code becomes both much simpler and much more accurate (to
> REAL*4 precision, so ~ 0.0001% accuracy!).
>
> I should also say that I believe that CTRUNCATE is not immune from these
> issues: it seems to be a straight translation at least in part from Fortran
> to C++ (though I'm far from being a C++ expert so I'll leave fixing
> CTRUNCATE to those who know what they're doing!).
>
> Sorry about the length of this: at least you can't say I didn't consult you!
>
> Cheers
>
> -- Ian
>



-- 
Dr Jon Agirre
York Structural Biology Laboratory / Department of Chemistry
University of York, Heslington, YO10 5DD, York, UK
http://www.york.ac.uk/chemistry/research/ysbl/people/research/jagirre/
+44 (0) 1904 32 8253