1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
|
/*
* refine-rns.c
*
* Refinement by Random Neighbour Search
*
* (c) 2006-2009 Thomas White <taw27@cam.ac.uk>
*
* synth2d - Two-Dimensional Crystallographic Fourier Synthesis
*
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <stdio.h>
#include <stdlib.h>
#include "refine.h"
#include "model.h"
#include "reflist.h"
#include "statistics.h"
void refine_neighboursearch(AtomicModel *model, ReflectionList *reflections,
RefinementSpec spec, double shift)
{
unsigned int iterations_without_change;
double min_r;
ReflectionList *model_reflections;
model_reflections = model_calculate_f(reflections, model, 69);
if ( spec & REFINE_SPEC_INTENSITIES ) {
/* stat_r2(Fobs, Fcalc) */
min_r = stat_r2(reflections, model_reflections);
} else {
/* stat_r(Fobs, Fcalc) */
min_r = stat_r(reflections, model_reflections);
}
iterations_without_change = 0;
while ( (model->refine_window->run_semaphore)
&& (iterations_without_change < MAX_REFINEMENT_ITERATIONS) ) {
unsigned int i;
double x, y, z;
double new_r;
AtomicModel *copy;
/* Copy the structure */
copy = model_copy(model);
/* 'Jiggle' the atoms */
for ( i=1; i<model->n_atoms; i++ ) {
if ( model->atoms[i].refine ) {
x = model->atoms[i].x;
y = model->atoms[i].y;
z = model->atoms[i].z;
if ( spec & REFINE_SPEC_X ) {
if ( random() > (0.9*RAND_MAX) )
x += shift;
if ( random() > (0.9*RAND_MAX) )
x -= shift;
if ( x<0 ) x = 1+x;
if ( x>1 ) x = x-1;
model->atoms[i].x = x;
}
if ( spec & REFINE_SPEC_Y ) {
if ( random() > (0.9*RAND_MAX) )
y += shift;
if ( random() > (0.9*RAND_MAX) )
y -= shift;
if ( y<0 ) y = 1+y;
if ( y>1 ) y = y-1;
model->atoms[i].y = y;
}
if ( spec & REFINE_SPEC_Z ) {
if ( random() > (0.9*RAND_MAX) )
z += shift;
if ( random() > (0.9*RAND_MAX) )
z -= shift;
if ( z<0 ) z = 1+z;
if ( z>1 ) z = z-1;
model->atoms[i].z = z;
}
}
}
/* Recalculate structure factors and check for improvement */
model_reflections = model_calculate_f(reflections, model, 69);
if ( spec & REFINE_SPEC_INTENSITIES ) {
/* stat_r(Fobs, Fcalc) */
new_r = stat_r2(reflections, model_reflections);
} else {
/* stat_r(Fobs, Fcalc) */
new_r = stat_r(reflections, model_reflections);
}
if ( new_r < min_r ) {
min_r = new_r;
refine_schedule_update(model);
iterations_without_change = 0;
model_free(copy);
} else {
model_move(model, copy);
model_free(copy);
iterations_without_change++;
}
free(model_reflections);
}
}
|
