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 package ffx.xray;
39
40 import ffx.algorithms.AlgorithmListener;
41 import ffx.algorithms.Terminatable;
42 import ffx.numerics.optimization.LBFGS;
43 import ffx.numerics.optimization.LineSearch.LineSearchResult;
44 import ffx.numerics.optimization.OptimizationListener;
45 import ffx.potential.MolecularAssembly;
46 import ffx.potential.bonded.Atom;
47 import ffx.potential.bonded.Molecule;
48 import ffx.potential.bonded.Residue;
49 import ffx.realspace.RealSpaceData;
50
51 import java.util.List;
52 import java.util.logging.Level;
53 import java.util.logging.Logger;
54
55 import static java.lang.String.format;
56
57
58
59
60
61
62
63
64 public class RefinementMinimize implements OptimizationListener, Terminatable {
65
66 private static final Logger logger = Logger.getLogger(RefinementMinimize.class.getName());
67 private static final double toSeconds = 1.0e-9;
68 public final RefinementEnergy refinementEnergy;
69 private final AlgorithmListener listener;
70 private final DataContainer dataContainer;
71 private final Atom[] activeAtomArray;
72 private final int nXYZ;
73 private final int nB;
74 private final int nOcc;
75 private final int n;
76 private final double[] x;
77 private final double[] grad;
78 private final double[] scaling;
79 private final RefinementMode refinementMode;
80 private boolean done = false;
81 private boolean terminate = false;
82 private long time;
83 private double grms;
84 private int nSteps;
85 private double eps = 0.1;
86
87
88
89
90
91
92
93 public RefinementMinimize(DataContainer data) {
94 this(data, RefinementMode.COORDINATES_AND_BFACTORS, null);
95 }
96
97
98
99
100
101
102
103
104
105 public RefinementMinimize(DataContainer data, RefinementMode refinementmode) {
106 this(data, refinementmode, null);
107 }
108
109
110
111
112
113
114
115
116
117
118 public RefinementMinimize(
119 DataContainer data, RefinementMode refinementMode, AlgorithmListener listener) {
120 dataContainer = data;
121 this.listener = listener;
122 this.refinementMode = refinementMode;
123 refinementEnergy = new RefinementEnergy(data, refinementMode, null);
124 nXYZ = refinementEnergy.nXYZ;
125 nB = refinementEnergy.nBFactor;
126 nOcc = refinementEnergy.nOccupancy;
127 n = refinementEnergy.getNumberOfVariables();
128
129 Atom[] atomArray = data.getAtomArray();
130 RefinementModel refinementModel = data.getRefinementModel();
131 int nAtoms = atomArray.length;
132
133
134 int count = 0;
135 for (Atom a : atomArray) {
136 if (a.isActive()) {
137 count++;
138 }
139 }
140 int nActive = count;
141 activeAtomArray = new Atom[count];
142 count = 0;
143 for (Atom a : atomArray) {
144 if (a.isActive()) {
145 activeAtomArray[count++] = a;
146 }
147 }
148
149 x = new double[n];
150 grad = new double[n];
151 scaling = new double[n];
152
153 refinementEnergy.getCoordinates(x);
154 refinementEnergy.setScaling(scaling);
155
156 double xyzscale = 1.0;
157 double bisoscale = 1.0;
158 double anisouscale = 50.0;
159 double occscale = 15.0;
160
161 if (refinementMode == RefinementMode.COORDINATES_AND_BFACTORS
162 || refinementMode == RefinementMode.COORDINATES_AND_BFACTORS_AND_OCCUPANCIES) {
163 bisoscale = 0.4;
164 anisouscale = 40.0;
165 }
166
167 if (refinementMode == RefinementMode.COORDINATES_AND_OCCUPANCIES
168 || refinementMode == RefinementMode.COORDINATES_AND_BFACTORS_AND_OCCUPANCIES) {
169 occscale = 10.0;
170 }
171
172
173 if (refinementMode == RefinementMode.COORDINATES
174 || refinementMode == RefinementMode.COORDINATES_AND_BFACTORS
175 || refinementMode == RefinementMode.COORDINATES_AND_OCCUPANCIES
176 || refinementMode == RefinementMode.COORDINATES_AND_BFACTORS_AND_OCCUPANCIES) {
177 for (int i = 0; i < n; i++) {
178 scaling[i] = xyzscale;
179 }
180 }
181
182 if (refinementMode == RefinementMode.BFACTORS
183 || refinementMode == RefinementMode.BFACTORS_AND_OCCUPANCIES
184 || refinementMode == RefinementMode.COORDINATES_AND_BFACTORS
185 || refinementMode == RefinementMode.COORDINATES_AND_BFACTORS_AND_OCCUPANCIES) {
186 int i = nXYZ;
187 int resnum = -1;
188 if (data instanceof DiffractionData) {
189 DiffractionData diffractionData = (DiffractionData) data;
190 int nres = diffractionData.getnResidueBFactor() + 1;
191 for (Atom a : activeAtomArray) {
192
193 if (a.getAtomicNumber() == 1) {
194 continue;
195 }
196 if (a.getAnisou(null) != null) {
197 for (int j = 0; j < 6; j++) {
198 scaling[i + j] = anisouscale;
199 }
200 i += 6;
201 } else if (diffractionData.isResidueBFactor()) {
202 if (resnum != a.getResidueNumber()) {
203 if (nres >= diffractionData.getnResidueBFactor()) {
204 if (resnum > -1 && i < nXYZ + nB - 1) {
205 i++;
206 }
207 if (i < nXYZ + nB) {
208 scaling[i] = bisoscale;
209 }
210 nres = 1;
211 } else {
212 nres++;
213 }
214 resnum = a.getResidueNumber();
215 }
216 } else {
217 scaling[i] = bisoscale;
218 i++;
219 }
220 }
221 } else {
222 logger.severe(" B refinement not supported for this data type!");
223 }
224 }
225
226 if (refinementMode == RefinementMode.OCCUPANCIES
227 || refinementMode == RefinementMode.BFACTORS_AND_OCCUPANCIES
228 || refinementMode == RefinementMode.COORDINATES_AND_OCCUPANCIES
229 || refinementMode == RefinementMode.COORDINATES_AND_BFACTORS_AND_OCCUPANCIES) {
230 if (nAtoms != nActive) {
231 logger.severe(" Occupancy refinement is not supported for inactive atoms.");
232 }
233 if (data instanceof DiffractionData) {
234
235 int i = nXYZ + nB;
236 for (List<Residue> list : refinementModel.getAltResidues()) {
237 for (int j = 0; j < list.size(); j++) {
238 scaling[i] = occscale;
239 i++;
240 }
241 }
242 for (List<Molecule> list : refinementModel.getAltMolecules()) {
243 for (int j = 0; j < list.size(); j++) {
244 scaling[i] = occscale;
245 i++;
246 }
247 }
248 } else {
249 logger.severe(" Occupancy refinement not supported for this data type!");
250 }
251 }
252 }
253
254
255
256
257
258
259
260 public static RefinementMode parseMode(String str) {
261 try {
262 return RefinementMode.valueOf(str.toUpperCase());
263 } catch (Exception e) {
264 logger.info(
265 format(
266 " Could not parse %s as a refinement mode; defaulting to coordinates.", str));
267 return RefinementMode.COORDINATES;
268 }
269 }
270
271
272
273
274
275
276 public Double getEps() {
277 boolean hasaniso = false;
278
279 for (Atom a : activeAtomArray) {
280
281 if (a.getAtomicNumber() == 1) {
282 continue;
283 }
284 if (a.getAnisou(null) != null) {
285 hasaniso = true;
286 break;
287 }
288 }
289
290 switch (refinementMode) {
291 case COORDINATES:
292 eps = 0.4;
293 break;
294 case BFACTORS, BFACTORS_AND_OCCUPANCIES:
295 if (hasaniso) {
296 eps = 20.0;
297 } else {
298 eps = 0.01;
299 }
300 break;
301 case COORDINATES_AND_BFACTORS, COORDINATES_AND_BFACTORS_AND_OCCUPANCIES:
302 if (hasaniso) {
303 eps = 20.0;
304 } else {
305 eps = 0.2;
306 }
307 break;
308 case OCCUPANCIES:
309 eps = 0.1;
310 break;
311 case COORDINATES_AND_OCCUPANCIES:
312 eps = 0.2;
313 break;
314 }
315
316 return eps;
317 }
318
319
320
321
322
323
324 public int getNB() {
325 return nB;
326 }
327
328
329
330
331
332
333 public int getNOcc() {
334 return nOcc;
335 }
336
337
338
339
340
341
342 public int getNXYZ() {
343 return nXYZ;
344 }
345
346
347
348
349
350
351 public RefinementEnergy minimize() {
352 return minimize(1.0);
353 }
354
355
356
357
358
359
360
361 public RefinementEnergy minimize(double eps) {
362 return minimize(7, eps, Integer.MAX_VALUE - 2);
363 }
364
365
366
367
368
369
370
371 public RefinementEnergy minimize(int maxiter) {
372 return minimize(7, 1.0, maxiter);
373 }
374
375
376
377
378
379
380
381
382 public RefinementEnergy minimize(double eps, int maxiter) {
383 return minimize(7, eps, maxiter);
384 }
385
386
387
388
389
390
391
392
393
394 public RefinementEnergy minimize(int m, double eps, int maxiter) {
395 String typestring;
396 if (dataContainer instanceof DiffractionData) {
397 typestring = "X-ray";
398 } else if (dataContainer instanceof RealSpaceData) {
399 typestring = "Real Space";
400 } else {
401 typestring = "null";
402 }
403
404 logger.info(" Beginning " + typestring + " Refinement");
405 switch (refinementMode) {
406 case COORDINATES:
407 logger.info(" Mode: Coordinates");
408 break;
409 case BFACTORS:
410 logger.info(" Mode: B-Factors");
411 break;
412 case COORDINATES_AND_BFACTORS:
413 logger.info(" Mode: Coordinates and B-Factors");
414 break;
415 case OCCUPANCIES:
416 logger.info(" Mode: Occupancies");
417 break;
418 case COORDINATES_AND_OCCUPANCIES:
419 logger.info(" Mode: Coordinates and Occupancies");
420 break;
421 case BFACTORS_AND_OCCUPANCIES:
422 logger.info(" Mode: B-Factors and Occupancies");
423 break;
424 case COORDINATES_AND_BFACTORS_AND_OCCUPANCIES:
425 logger.info(" Mode: Coordinates, B-Factors and Occupancies");
426 break;
427 }
428 logger.info(" Number of Parameters: " + n);
429
430 refinementEnergy.getCoordinates(x);
431
432
433 for (int i = 0; i < n; i++) {
434 x[i] *= scaling[i];
435 }
436
437 long mtime = -System.nanoTime();
438 time = -System.nanoTime();
439 done = false;
440 int status;
441 double e = refinementEnergy.energyAndGradient(x, grad);
442 status = LBFGS.minimize(n, m, x, e, grad, eps, maxiter, refinementEnergy, this);
443 done = true;
444 switch (status) {
445 case 0:
446 logger.info(format("\n Optimization achieved convergence criteria: %8.5f", grms));
447 break;
448 case 1:
449 logger.info(format("\n Optimization terminated at step %d.", nSteps));
450 break;
451 default:
452 logger.warning("\n Optimization failed.");
453 }
454
455 if (logger.isLoggable(Level.INFO)) {
456 StringBuilder sb = new StringBuilder();
457 mtime += System.nanoTime();
458 sb.append(format(" Optimization time: %g (sec)", mtime * toSeconds));
459 logger.info(sb.toString());
460 }
461
462 return refinementEnergy;
463 }
464
465
466 @Override
467 public boolean optimizationUpdate(
468 int iter,
469 int nBFGS,
470 int nfun,
471 double grms,
472 double xrms,
473 double f,
474 double df,
475 double angle,
476 LineSearchResult info) {
477 long currentTime = System.nanoTime();
478 Double seconds = (currentTime - time) * 1.0e-9;
479 time = currentTime;
480 this.grms = grms;
481 this.nSteps = iter;
482
483
484 if (listener != null) {
485 MolecularAssembly[] molecularAssembly = dataContainer.getMolecularAssemblies();
486 for (MolecularAssembly ma : molecularAssembly) {
487 listener.algorithmUpdate(ma);
488 }
489 }
490
491 if (iter == 0) {
492 if (nBFGS > 0) {
493 logger.info("\n Limited Memory BFGS Quasi-Newton Optimization: \n");
494 } else {
495 logger.info("\n Steepest Decent Optimization: \n");
496 }
497 logger.info(
498 " Cycle Energy G RMS Delta E Delta X Angle Evals Time "
499 + dataContainer.printOptimizationHeader());
500 }
501 if (info == null) {
502 logger.info(format("%6d %12.3f %10.3f", iter, f, grms));
503 } else if (info == LineSearchResult.Success) {
504 StringBuilder sb = new StringBuilder();
505 sb.append(format("%6d %12.3f %10.3f %10.3f %9.4f %8.2f %6d %8.3f ",
506 iter, f, grms, df, xrms, angle, nfun, seconds));
507 sb.append(dataContainer.printOptimizationUpdate());
508 logger.info(sb.toString());
509 } else {
510 logger.info(
511 format("%6d %12.3f %10.3f %10.3f %9.4f %8.2f %6d %8s",
512 iter, f, grms, df, xrms, angle, nfun, info));
513 }
514 if (terminate) {
515 logger.info(" The optimization recieved a termination request.");
516
517 return false;
518 }
519 return true;
520 }
521
522
523 @Override
524 public void terminate() {
525 terminate = true;
526 while (!done) {
527 synchronized (this) {
528 try {
529 wait(1);
530 } catch (Exception e) {
531 logger.log(Level.WARNING, " Exception terminating minimization.\n", e);
532 }
533 }
534 }
535 }
536
537
538 public enum RefinementMode {
539
540
541 COORDINATES,
542
543 BFACTORS,
544
545 COORDINATES_AND_BFACTORS,
546
547 OCCUPANCIES,
548
549 BFACTORS_AND_OCCUPANCIES,
550
551 COORDINATES_AND_OCCUPANCIES,
552
553 COORDINATES_AND_BFACTORS_AND_OCCUPANCIES
554 }
555 }