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.realspace;
39
40 import static ffx.numerics.math.ScalarMath.mod;
41 import static java.lang.String.format;
42 import static java.util.Arrays.fill;
43 import static org.apache.commons.math3.util.FastMath.floor;
44
45 import edu.rit.pj.IntegerForLoop;
46 import edu.rit.pj.ParallelRegion;
47 import edu.rit.pj.ParallelTeam;
48 import edu.rit.pj.reduction.SharedDouble;
49 import ffx.crystal.Crystal;
50 import ffx.numerics.atomic.AtomicDoubleArray.AtomicDoubleArrayImpl;
51 import ffx.numerics.atomic.AtomicDoubleArray3D;
52 import ffx.numerics.spline.TriCubicSpline;
53 import ffx.potential.MolecularAssembly;
54 import ffx.potential.Utilities;
55 import ffx.potential.bonded.Atom;
56 import ffx.potential.bonded.Molecule;
57 import ffx.potential.bonded.Residue;
58 import ffx.realspace.parsers.RealSpaceFile;
59 import ffx.xray.DataContainer;
60 import ffx.xray.DiffractionData;
61 import ffx.xray.RefinementMinimize.RefinementMode;
62 import ffx.xray.RefinementModel;
63 import java.util.List;
64 import java.util.logging.Level;
65 import java.util.logging.Logger;
66 import org.apache.commons.configuration2.CompositeConfiguration;
67
68
69
70
71
72
73
74 public class RealSpaceData implements DataContainer {
75
76 private static final Logger logger = Logger.getLogger(RealSpaceData.class.getName());
77
78 private final ParallelTeam parallelTeam;
79
80 private final RealSpaceRegion realSpaceRegion;
81
82 private final RealSpaceFile[] realSpaceFile;
83
84 private final int nRealSpaceData;
85
86 private final RefinementModel refinementModel;
87
88 private final boolean lambdaTerm;
89
90 private final MolecularAssembly[] molecularAssemblies;
91
92 private Crystal[] crystal;
93
94 private RealSpaceRefinementData[] refinementData;
95
96 private double realSpaceEnergy = 0.0;
97
98 private double realSpacedUdL = 0.0;
99
100 private double[] realSpaceGradient;
101
102 private double[] realSpacedUdXdL;
103
104 private double xweight;
105
106 private double lambda = 1.0;
107
108
109
110
111
112
113
114
115
116
117
118 public RealSpaceData(
119 MolecularAssembly molecularAssembly,
120 CompositeConfiguration properties,
121 ParallelTeam parallelTeam,
122 DiffractionData diffractionData) {
123 this(new MolecularAssembly[] {molecularAssembly}, properties, parallelTeam, diffractionData);
124 }
125
126
127
128
129
130
131
132
133
134
135
136 public RealSpaceData(
137 MolecularAssembly[] molecularAssemblies,
138 CompositeConfiguration properties,
139 ParallelTeam parallelTeam,
140 DiffractionData diffractionData) {
141 this.molecularAssemblies = molecularAssemblies;
142 this.parallelTeam = parallelTeam;
143 this.realSpaceFile = null;
144 this.nRealSpaceData = 1;
145 crystal = new Crystal[nRealSpaceData];
146 crystal[0] = diffractionData.getCrystal()[0];
147 refinementData = new RealSpaceRefinementData[nRealSpaceData];
148 refinementData[0] = new RealSpaceRefinementData();
149 refinementData[0].setPeriodic(true);
150
151 xweight = properties.getDouble("xweight", 1.0);
152 lambdaTerm = properties.getBoolean("lambdaterm", false);
153
154 if (logger.isLoggable(Level.INFO)) {
155 StringBuilder sb = new StringBuilder(" Real Space Refinement Settings\n");
156 sb.append(format(" Refinement weight (xweight): %5.3f", xweight));
157 logger.info(sb.toString());
158 }
159
160 if (!diffractionData.getScaled()[0]) {
161 diffractionData.scaleBulkFit();
162 diffractionData.printStats();
163 }
164
165
166 diffractionData.getCrystalReciprocalSpacesFc()[0].computeAtomicGradients(
167 diffractionData.getRefinementData()[0].foFc1,
168 diffractionData.getRefinementData()[0].freeR,
169 diffractionData.getRefinementData()[0].rFreeFlag,
170 RefinementMode.COORDINATES);
171
172 refinementData[0].setOrigin(0, 0, 0);
173 int extx = (int) diffractionData.getCrystalReciprocalSpacesFc()[0].getXDim();
174 int exty = (int) diffractionData.getCrystalReciprocalSpacesFc()[0].getYDim();
175 int extz = (int) diffractionData.getCrystalReciprocalSpacesFc()[0].getZDim();
176 refinementData[0].setExtent(extx, exty, extz);
177 refinementData[0].setNI(extx, exty, extz);
178 refinementData[0].setData(new double[extx * exty * extz]);
179 for (int k = 0; k < extz; k++) {
180 for (int j = 0; j < exty; j++) {
181 for (int i = 0; i < extx; i++) {
182 int index1 = (i + extx * (j + exty * k));
183 int index2 = 2 * index1;
184 refinementData[0].getData()[index1] =
185 diffractionData.getCrystalReciprocalSpacesFc()[0].getDensityGrid()[index2];
186 }
187 }
188 }
189
190
191 refinementModel = new RefinementModel(molecularAssemblies);
192
193
194 int nAtoms = refinementModel.getTotalAtomArray().length;
195 realSpaceRegion =
196 new RealSpaceRegion(parallelTeam.getThreadCount(), nAtoms, refinementData.length);
197 }
198
199
200
201
202
203
204
205
206
207
208 public RealSpaceData(
209 MolecularAssembly assembly, CompositeConfiguration properties, ParallelTeam parallelTeam) {
210 this(new MolecularAssembly[] {assembly}, properties, parallelTeam, new RealSpaceFile(assembly));
211 }
212
213
214
215
216
217
218
219
220
221
222 public RealSpaceData(
223 MolecularAssembly assembly,
224 CompositeConfiguration properties,
225 ParallelTeam parallelTeam,
226 RealSpaceFile... datafile) {
227 this(new MolecularAssembly[] {assembly}, properties, parallelTeam, datafile);
228 }
229
230
231
232
233
234
235
236
237
238
239
240 public RealSpaceData(
241 MolecularAssembly[] molecularAssemblies,
242 CompositeConfiguration properties,
243 ParallelTeam parallelTeam,
244 RealSpaceFile... dataFile) {
245
246 this.molecularAssemblies = molecularAssemblies;
247 this.parallelTeam = parallelTeam;
248 this.realSpaceFile = dataFile;
249 this.nRealSpaceData = dataFile.length;
250 crystal = new Crystal[nRealSpaceData];
251 refinementData = new RealSpaceRefinementData[nRealSpaceData];
252
253 xweight = properties.getDouble("xweight", 1.0);
254 lambdaTerm = properties.getBoolean("lambdaterm", false);
255
256 for (int i = 0; i < nRealSpaceData; i++) {
257 crystal[i] =
258 dataFile[i].getRealSpaceFileFilter().getCrystal(dataFile[i].getFilename(), properties);
259 if (crystal[i] == null) {
260 logger.severe(" CCP4 map file does not contain full crystal information!");
261 }
262 }
263
264 for (int i = 0; i < nRealSpaceData; i++) {
265 refinementData[i] = new RealSpaceRefinementData();
266 dataFile[i]
267 .getRealSpaceFileFilter()
268 .readFile(dataFile[i].getFilename(), refinementData[i], properties);
269
270 if (refinementData[i].getOrigin()[0] == 0
271 && refinementData[i].getOrigin()[1] == 0
272 && refinementData[i].getOrigin()[2] == 0
273 && refinementData[i].getExtent()[0] == refinementData[i].getNi()[0]
274 && refinementData[i].getExtent()[1] == refinementData[i].getNi()[1]
275 && refinementData[i].getExtent()[2] == refinementData[i].getNi()[2]) {
276 refinementData[i].setPeriodic(true);
277 }
278 }
279
280 if (logger.isLoggable(Level.INFO)) {
281 StringBuilder sb = new StringBuilder(" Real Space Refinement Settings\n");
282 sb.append(format(" Refinement weight (xweight): %5.3f", xweight));
283 logger.info(sb.toString());
284 }
285
286
287 refinementModel = new RefinementModel(molecularAssemblies);
288
289
290 int nAtoms = refinementModel.getTotalAtomArray().length;
291 realSpaceRegion =
292 new RealSpaceRegion(parallelTeam.getThreadCount(), nAtoms, refinementData.length);
293 }
294
295
296
297
298
299
300 public boolean destroy() {
301 try {
302 boolean underlyingShutdown = true;
303 for (MolecularAssembly assembly : molecularAssemblies) {
304
305 boolean thisShutdown = assembly.destroy();
306 underlyingShutdown = underlyingShutdown && thisShutdown;
307 }
308 parallelTeam.shutdown();
309 return underlyingShutdown;
310 } catch (Exception ex) {
311 logger.warning(format(" Exception in shutting down a RealSpaceData: %s", ex));
312 logger.info(Utilities.stackTraceToString(ex));
313 return false;
314 }
315 }
316
317
318 @Override
319 public Atom[] getActiveAtomArray() {
320 return refinementModel.getActiveAtomArray();
321 }
322
323
324 @Override
325 public List<List<Molecule>> getAltMolecules() {
326 return refinementModel.getAltMolecules();
327 }
328
329
330 @Override
331 public List<List<Residue>> getAltResidues() {
332 return refinementModel.getAltResidues();
333 }
334
335
336 @Override
337 public Atom[] getAtomArray() {
338 return refinementModel.getTotalAtomArray();
339 }
340
341
342
343
344
345
346 public Crystal[] getCrystal() {
347 return crystal;
348 }
349
350
351
352
353
354
355 public void setCrystal(Crystal[] crystal) {
356 this.crystal = crystal;
357 }
358
359
360
361
362
363
364 public double getLambda() {
365 return lambda;
366 }
367
368
369
370
371
372
373 public void setLambda(double lambda) {
374 this.lambda = lambda;
375 }
376
377
378 @Override
379 public MolecularAssembly[] getMolecularAssemblies() {
380 return molecularAssemblies;
381 }
382
383
384
385
386
387
388
389 public double[] getRealSpaceGradient(double[] gradient) {
390 int nAtoms = refinementModel.getTotalAtomArray().length;
391 int nActiveAtoms = 0;
392 for (int i = 0; i < nAtoms; i++) {
393 if (refinementModel.getTotalAtomArray()[i].isActive()) {
394 nActiveAtoms++;
395 }
396 }
397
398 if (gradient == null || gradient.length < nActiveAtoms * 3) {
399 gradient = new double[nActiveAtoms * 3];
400 }
401 for (int i = 0; i < nActiveAtoms * 3; i++) {
402 gradient[i] += realSpaceGradient[i];
403 }
404 return gradient;
405 }
406
407
408
409
410
411
412 public RealSpaceRefinementData[] getRefinementData() {
413 return refinementData;
414 }
415
416
417
418
419
420
421 public void setRefinementData(RealSpaceRefinementData[] refinementData) {
422 this.refinementData = refinementData;
423 }
424
425
426 @Override
427 public RefinementModel getRefinementModel() {
428 return refinementModel;
429 }
430
431
432 @Override
433 public double getWeight() {
434 return xweight;
435 }
436
437
438 @Override
439 public void setWeight(double weight) {
440 this.xweight = weight;
441 }
442
443
444 @Override
445 public String printEnergyUpdate() {
446 StringBuilder sb = new StringBuilder();
447 for (int i = 0; i < nRealSpaceData; i++) {
448 sb.append(
449 format(
450 " dataset %d (weight: %5.1f): chemical energy: %8.2f density score: %8.2f\n",
451 i + 1,
452 realSpaceFile[i].getWeight(),
453 molecularAssemblies[0].getPotentialEnergy().getTotalEnergy(),
454 realSpaceFile[i].getWeight() * refinementData[i].getDensityScore()));
455 }
456 return sb.toString();
457 }
458
459
460 @Override
461 public String printOptimizationHeader() {
462 return "Density score";
463 }
464
465
466 @Override
467 public String printOptimizationUpdate() {
468 StringBuilder sb = new StringBuilder();
469 for (int i = 0; i < nRealSpaceData; i++) {
470 sb.append(format("%6.2f ", refinementData[i].getDensityScore()));
471 }
472 return sb.toString();
473 }
474
475
476
477
478
479
480 double computeRealSpaceTarget() {
481
482 long time = -System.nanoTime();
483
484 realSpaceEnergy = 0.0;
485
486 realSpacedUdL = 0.0;
487
488 int nActive = 0;
489 int nAtoms = refinementModel.getTotalAtomArray().length;
490 for (int i = 0; i < nAtoms; i++) {
491 if (refinementModel.getTotalAtomArray()[i].isActive()) {
492 nActive++;
493 }
494 }
495
496 int nGrad = nActive * 3;
497 if (realSpaceGradient == null || realSpaceGradient.length < nGrad) {
498 realSpaceGradient = new double[nGrad];
499 } else {
500 fill(realSpaceGradient, 0.0);
501 }
502
503
504 if (realSpacedUdXdL == null || realSpacedUdXdL.length < nGrad) {
505 realSpacedUdXdL = new double[nGrad];
506 } else {
507 fill(realSpacedUdXdL, 0.0);
508 }
509
510 try {
511 parallelTeam.execute(realSpaceRegion);
512 } catch (Exception e) {
513 String message = " Exception computing real space energy";
514 logger.log(Level.SEVERE, message, e);
515 }
516
517 time += System.nanoTime();
518 if (logger.isLoggable(Level.FINE)) {
519 logger.fine(format(" Real space energy time: %16.8f (sec).", time * 1.0e-9));
520 }
521
522 return realSpaceEnergy;
523 }
524
525
526
527
528
529
530 double getdEdL() {
531 return realSpacedUdL;
532 }
533
534
535
536
537
538
539
540 double[] getdEdXdL(double[] gradient) {
541 int nAtoms = refinementModel.getTotalAtomArray().length;
542 int nActiveAtoms = 0;
543 for (int i = 0; i < nAtoms; i++) {
544 if (refinementModel.getTotalAtomArray()[i].isActive()) {
545 nActiveAtoms++;
546 }
547 }
548
549 if (gradient == null || gradient.length < nActiveAtoms * 3) {
550 gradient = new double[nActiveAtoms * 3];
551 }
552
553 for (int i = 0; i < nActiveAtoms * 3; i++) {
554 gradient[i] += realSpacedUdXdL[i];
555 }
556 return gradient;
557 }
558
559
560
561
562
563
564 private RealSpaceFile[] getRealSpaceFile() {
565 return realSpaceFile;
566 }
567
568
569
570
571
572
573 private int getnRealSpaceData() {
574 return nRealSpaceData;
575 }
576
577 private class RealSpaceRegion extends ParallelRegion {
578
579 private final AtomicDoubleArray3D gradient;
580 private final AtomicDoubleArray3D lambdaGrad;
581 private final InitializationLoop[] initializationLoops;
582 private final RealSpaceLoop[] realSpaceLoops;
583 private final SharedDouble[] sharedTarget;
584 private final SharedDouble shareddUdL;
585 private final int nAtoms;
586 private final int nData;
587
588 RealSpaceRegion(int nThreads, int nAtoms, int nData) {
589 this.nAtoms = nAtoms;
590 this.nData = nData;
591 initializationLoops = new InitializationLoop[nThreads];
592 realSpaceLoops = new RealSpaceLoop[nThreads];
593 sharedTarget = new SharedDouble[nData];
594 for (int i = 0; i < nData; i++) {
595 sharedTarget[i] = new SharedDouble();
596 }
597 shareddUdL = new SharedDouble();
598 gradient = new AtomicDoubleArray3D(AtomicDoubleArrayImpl.MULTI, nThreads, nAtoms);
599 lambdaGrad = new AtomicDoubleArray3D(AtomicDoubleArrayImpl.MULTI, nThreads, nAtoms);
600 }
601
602 @Override
603 public void finish() {
604
605 realSpaceEnergy = 0;
606 for (int i = 0; i < nData; i++) {
607 refinementData[i].setDensityScore(sharedTarget[i].get());
608 realSpaceEnergy += getRefinementData()[i].getDensityScore();
609 }
610 realSpacedUdL = shareddUdL.get();
611 int index = 0;
612 for (int i = 0; i < nAtoms; i++) {
613 Atom atom = refinementModel.getTotalAtomArray()[i];
614 if (atom.isActive()) {
615 int ii = index * 3;
616 double gx = gradient.getX(i);
617 double gy = gradient.getY(i);
618 double gz = gradient.getZ(i);
619 realSpaceGradient[ii] = gx;
620 realSpaceGradient[ii + 1] = gy;
621 realSpaceGradient[ii + 2] = gz;
622 atom.setXYZGradient(gx, gy, gz);
623 gx = lambdaGrad.getX(i);
624 gy = lambdaGrad.getY(i);
625 gz = lambdaGrad.getZ(i);
626 realSpacedUdXdL[ii] = gx;
627 realSpacedUdXdL[ii + 1] = gy;
628 realSpacedUdXdL[ii + 2] = gz;
629 atom.setLambdaXYZGradient(gx, gy, gz);
630 index++;
631 }
632 }
633 }
634
635 @Override
636 public void run() throws Exception {
637 int threadID = getThreadIndex();
638
639 if (initializationLoops[threadID] == null) {
640 initializationLoops[threadID] = new InitializationLoop();
641 }
642 execute(0, nAtoms - 1, initializationLoops[threadID]);
643
644 if (realSpaceLoops[threadID] == null) {
645 realSpaceLoops[threadID] = new RealSpaceLoop();
646 }
647 execute(0, nAtoms - 1, realSpaceLoops[threadID]);
648 }
649
650 @Override
651 public void start() {
652 for (int i = 0; i < nData; i++) {
653 sharedTarget[i].set(0.0);
654 }
655 shareddUdL.set(0.0);
656 gradient.alloc(nAtoms);
657 lambdaGrad.alloc(nAtoms);
658 }
659
660
661 private class InitializationLoop extends IntegerForLoop {
662
663 @Override
664 public void run(int lb, int ub) {
665 int threadID = getThreadIndex();
666 gradient.reset(threadID, lb, ub);
667 lambdaGrad.reset(threadID, lb, ub);
668 for (int i = lb; i <= ub; i++) {
669 Atom a = refinementModel.getTotalAtomArray()[i];
670 a.setXYZGradient(0.0, 0.0, 0.0);
671 a.setLambdaXYZGradient(0.0, 0.0, 0.0);
672 }
673 }
674 }
675
676 private class RealSpaceLoop extends IntegerForLoop {
677
678 double[] target = new double[nData];
679 double localdUdL;
680
681 @Override
682 public void finish() {
683 for (int i = 0; i < nData; i++) {
684 sharedTarget[i].addAndGet(target[i]);
685 }
686 shareddUdL.addAndGet(localdUdL);
687 }
688
689 @Override
690 public void run(int first, int last) throws Exception {
691
692 int threadID = getThreadIndex();
693 double[] xyz = new double[3];
694 double[] uvw = new double[3];
695 double[] grad = new double[3];
696 double[][][] scalar = new double[4][4][4];
697 TriCubicSpline spline = new TriCubicSpline();
698
699 for (int i = 0; i < getnRealSpaceData(); i++) {
700
701
702 int extX = getRefinementData()[i].getExtent()[0];
703 int extY = getRefinementData()[i].getExtent()[1];
704 int extZ = getRefinementData()[i].getExtent()[2];
705 int nX = getRefinementData()[i].getNi()[0];
706 int nY = getRefinementData()[i].getNi()[1];
707 int nZ = getRefinementData()[i].getNi()[2];
708 int originX = getRefinementData()[i].getOrigin()[0];
709 int originY = getRefinementData()[i].getOrigin()[1];
710 int originZ = getRefinementData()[i].getOrigin()[2];
711
712 for (int ia = first; ia <= last; ia++) {
713 Atom a = refinementModel.getTotalAtomArray()[ia];
714
715
716 if (!a.getUse()) {
717 continue;
718 }
719
720 double lambdai = 1.0;
721 double dUdL = 0.0;
722 if (lambdaTerm && a.applyLambda()) {
723 lambdai = getLambda();
724 dUdL = 1.0;
725 }
726 a.getXYZ(xyz);
727 getCrystal()[i].toFractionalCoordinates(xyz, uvw);
728
729
730 final double frx = nX * uvw[0];
731 final int ifrx = ((int) floor(frx)) - originX;
732 final double dfrx = frx - floor(frx);
733
734 final double fry = nY * uvw[1];
735 final int ifry = ((int) floor(fry)) - originY;
736 final double dfry = fry - floor(fry);
737
738 final double frz = nZ * uvw[2];
739 final int ifrz = ((int) floor(frz)) - originZ;
740 final double dfrz = frz - floor(frz);
741
742 if (!refinementData[i].isPeriodic()) {
743 if (ifrx - 1 < 0
744 || ifrx + 2 > extX
745 || ifry - 1 < 0
746 || ifry + 2 > extY
747 || ifrz - 1 < 0
748 || ifrz + 2 > extZ) {
749 String message = format(" Atom %s is outside the density will be ignored.", a);
750 logger.warning(message);
751 continue;
752 }
753 }
754
755
756 if (getRefinementData()[i].isPeriodic()) {
757 for (int ui = ifrx - 1; ui < ifrx + 3; ui++) {
758 int uii = ui - (ifrx - 1);
759 int pui = mod(ui, extX);
760 for (int vi = ifry - 1; vi < ifry + 3; vi++) {
761 int vii = vi - (ifry - 1);
762 int pvi = mod(vi, extY);
763 for (int wi = ifrz - 1; wi < ifrz + 3; wi++) {
764 int wii = wi - (ifrz - 1);
765 int pwi = mod(wi, extZ);
766 scalar[uii][vii][wii] = getRefinementData()[i].getDataIndex(pui, pvi, pwi);
767 }
768 }
769 }
770 } else {
771 for (int ui = ifrx - 1; ui < ifrx + 3; ui++) {
772 int uii = ui - (ifrx - 1);
773 for (int vi = ifry - 1; vi < ifry + 3; vi++) {
774 int vii = vi - (ifry - 1);
775 for (int wi = ifrz - 1; wi < ifrz + 3; wi++) {
776 int wii = wi - (ifrz - 1);
777 scalar[uii][vii][wii] = getRefinementData()[i].getDataIndex(ui, vi, wi);
778 }
779 }
780 }
781 }
782
783
784 double scale;
785 double scaledUdL;
786 double atomicWeight = a.getAtomType().atomicWeight;
787 if (getRealSpaceFile() != null) {
788 atomicWeight *= getRealSpaceFile()[i].getWeight();
789 }
790 scale = -1.0 * lambdai * atomicWeight;
791 scaledUdL = -1.0 * dUdL * atomicWeight;
792
793 double val = spline.spline(dfrx, dfry, dfrz, scalar, grad);
794 target[i] += scale * val;
795 localdUdL += scaledUdL * val;
796
797 if (a.isActive()) {
798 grad[0] = grad[0] * nX;
799 grad[1] = grad[1] * nY;
800 grad[2] = grad[2] * nZ;
801
802 xyz[0] =
803 grad[0] * getCrystal()[i].A00
804 + grad[1] * getCrystal()[i].A01
805 + grad[2] * getCrystal()[i].A02;
806 xyz[1] =
807 grad[0] * getCrystal()[i].A10
808 + grad[1] * getCrystal()[i].A11
809 + grad[2] * getCrystal()[i].A12;
810 xyz[2] =
811 grad[0] * getCrystal()[i].A20
812 + grad[1] * getCrystal()[i].A21
813 + grad[2] * getCrystal()[i].A22;
814 gradient.add(threadID, ia, scale * xyz[0], scale * xyz[1], scale * xyz[2]);
815 lambdaGrad.add(
816 threadID, ia, scaledUdL * xyz[0], scaledUdL * xyz[1], scaledUdL * xyz[2]);
817 }
818 }
819 }
820 gradient.reduce(first, last);
821 lambdaGrad.reduce(first, last);
822 }
823
824 @Override
825 public void start() {
826 for (int i = 0; i < nData; i++) {
827 target[i] = 0;
828 }
829 localdUdL = 0;
830 }
831 }
832 }
833 }