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 edu.rit.pj.IntegerForLoop;
41 import edu.rit.pj.IntegerSchedule;
42 import edu.rit.pj.ParallelRegion;
43 import edu.rit.pj.ParallelTeam;
44 import edu.rit.pj.reduction.SharedIntegerArray;
45 import ffx.crystal.Crystal;
46 import ffx.crystal.HKL;
47 import ffx.crystal.ReflectionList;
48 import ffx.crystal.Resolution;
49 import ffx.crystal.SymOp;
50 import ffx.numerics.fft.Complex;
51 import ffx.numerics.fft.Complex3DParallel;
52 import ffx.numerics.math.ComplexNumber;
53 import ffx.potential.bonded.Atom;
54 import ffx.potential.nonbonded.RowLoop;
55 import ffx.potential.nonbonded.RowRegion;
56 import ffx.potential.nonbonded.SliceLoop;
57 import ffx.potential.nonbonded.SliceRegion;
58 import ffx.potential.nonbonded.SpatialDensityLoop;
59 import ffx.potential.nonbonded.SpatialDensityRegion;
60 import ffx.xray.RefinementMinimize.RefinementMode;
61
62 import java.util.ArrayList;
63 import java.util.List;
64 import java.util.logging.Level;
65 import java.util.logging.Logger;
66
67 import static ffx.crystal.SymOp.applyTransSymRot;
68 import static ffx.numerics.fft.Complex3D.interleavedIndex;
69 import static ffx.numerics.math.ScalarMath.mod;
70 import static java.lang.String.format;
71 import static java.lang.System.arraycopy;
72 import static java.util.Arrays.fill;
73 import static org.apache.commons.math3.util.FastMath.abs;
74 import static org.apache.commons.math3.util.FastMath.exp;
75 import static org.apache.commons.math3.util.FastMath.floor;
76 import static org.apache.commons.math3.util.FastMath.max;
77 import static org.apache.commons.math3.util.FastMath.min;
78 import static org.apache.commons.math3.util.FastMath.pow;
79 import static org.apache.commons.math3.util.FastMath.sqrt;
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103 public class CrystalReciprocalSpace {
104
105 private static final Logger logger = Logger.getLogger(CrystalReciprocalSpace.class.getName());
106 private static final double toSeconds = 1.0e-9;
107 private final boolean neutron;
108 private final double bAdd;
109 private final double fftScale;
110 private final double[] densityGrid;
111 private final double[] solventGrid;
112 private final int nSymm;
113 private final int bulkNSymm;
114 private final int nAtoms;
115 private final int fftX, fftY, fftZ;
116 private final double ifftX, ifftY, ifftZ;
117 private final int halfFFTX;
118 private final int complexFFT3DSpace;
119 private final int threadCount;
120 private final int bufferSize = 3;
121
122 private final SharedIntegerArray optSliceWeightAtomic;
123 private final SharedIntegerArray optSliceWeightSolvent;
124 private final SharedIntegerArray optSliceWeightBulkSolvent;
125 private final int[] previousOptSliceWeightAtomic;
126 private final int[] previousOptSliceWeightSolvent;
127 private final int[] previousOptSliceWeightBulkSolvent;
128 private final SliceSchedule atomicSliceSchedule;
129 private final SliceSchedule solventSliceSchedule;
130 private final SliceSchedule bulkSolventSliceSchedule;
131
132 private final SharedIntegerArray optRowWeightAtomic;
133 private final SharedIntegerArray optRowWeightSolvent;
134 private final SharedIntegerArray optRowWeightBulkSolvent;
135 private final int[] previousOptRowWeightAtomic;
136 private final int[] previousOptRowWeightSolvent;
137 private final int[] previousOptRowWeightBulkSolvent;
138 private final RowSchedule atomicRowSchedule;
139 private final RowSchedule solventRowSchedule;
140 private final RowSchedule bulkSolventRowSchedule;
141
142 private final ParallelTeam parallelTeam;
143 private final Atom[] atoms;
144 private final Crystal crystal;
145 private final ReflectionList reflectionList;
146 private final FormFactor[][] atomFormFactors;
147 private final FormFactor[][] solventFormFactors;
148 private final GridMethod gridMethod;
149 private final SharedIntegerArray optAtomicGradientWeight;
150 private final int[] previousOptAtomicGradientWeight;
151 private final GradientSchedule atomicGradientSchedule;
152
153
154
155
156 private final SpatialDensityRegion atomicDensityRegion;
157
158 private final SpatialDensityRegion solventDensityRegion;
159
160
161
162
163 private final SliceRegion atomicSliceRegion;
164
165 private final AtomicSliceLoop[] atomicSliceLoops;
166 private final SliceRegion solventSliceRegion;
167 private final SolventSliceLoop[] solventSliceLoops;
168 private final BulkSolventSliceLoop[] bulkSolventSliceLoops;
169
170
171
172
173 private final RowRegion atomicRowRegion;
174
175 private final AtomicRowLoop[] atomicRowLoops;
176 private final RowRegion solventRowRegion;
177 private final SolventRowLoop[] solventRowLoops;
178 private final BulkSolventRowLoop[] bulkSolventRowLoops;
179
180 private final InitRegion initRegion;
181
182 private final ExtractRegion extractRegion;
183 private final AtomicScaleRegion atomicScaleRegion;
184 private final AtomicGradientRegion atomicGradientRegion;
185
186 private final BabinetRegion babinetRegion;
187
188 private final SolventGridRegion solventGridRegion;
189 private final SolventScaleRegion solventScaleRegion;
190 private final SolventGradientRegion solventGradientRegion;
191 private final Complex3DParallel complexFFT3D;
192 protected SolventModel solventModel;
193 boolean lambdaTerm = false;
194 double solventA;
195 double solventB;
196 private final boolean solvent;
197 private boolean useThreeGaussians = true;
198
199 private final double[][][] coordinates;
200 private double weight = 1.0;
201 private final int aRadGrid;
202
203 private boolean nativeEnvironmentApproximation = false;
204
205
206
207
208
209
210
211
212
213
214
215 public CrystalReciprocalSpace(
216 ReflectionList reflectionList,
217 Atom[] atoms,
218 ParallelTeam fftTeam,
219 ParallelTeam parallelTeam) {
220 this(
221 reflectionList,
222 atoms,
223 fftTeam,
224 parallelTeam,
225 false,
226 false,
227 SolventModel.POLYNOMIAL,
228 GridMethod.SLICE);
229 }
230
231
232
233
234
235
236
237
238
239
240
241
242 public CrystalReciprocalSpace(
243 ReflectionList reflectionList,
244 Atom[] atoms,
245 ParallelTeam fftTeam,
246 ParallelTeam parallelTeam,
247 boolean solventMask) {
248 this(
249 reflectionList,
250 atoms,
251 fftTeam,
252 parallelTeam,
253 solventMask,
254 false,
255 SolventModel.POLYNOMIAL,
256 GridMethod.SLICE);
257 }
258
259
260
261
262
263
264
265
266
267
268
269
270 public CrystalReciprocalSpace(
271 ReflectionList reflectionList,
272 Atom[] atoms,
273 ParallelTeam fftTeam,
274 ParallelTeam parallelTeam,
275 boolean solventMask,
276 boolean neutron) {
277 this(
278 reflectionList,
279 atoms,
280 fftTeam,
281 parallelTeam,
282 solventMask,
283 neutron,
284 SolventModel.POLYNOMIAL,
285 GridMethod.SLICE);
286 }
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302 public CrystalReciprocalSpace(
303 ReflectionList reflectionlist,
304 Atom[] atoms,
305 ParallelTeam fftTeam,
306 ParallelTeam parallelTeam,
307 boolean solventMask,
308 boolean neutron,
309 SolventModel solventModel,
310 GridMethod gridMethod) {
311 this.reflectionList = reflectionlist;
312 this.atoms = atoms;
313 this.parallelTeam = parallelTeam;
314 this.solvent = solventMask;
315 this.neutron = neutron;
316 this.solventModel = solventModel;
317 this.gridMethod = gridMethod;
318
319 nAtoms = atoms.length;
320 crystal = reflectionlist.crystal;
321 nSymm = 1;
322 threadCount = parallelTeam.getThreadCount();
323
324
325 bulkNSymm = crystal.spaceGroup.symOps.size();
326 Resolution resolution = reflectionlist.resolution;
327 double gridFactor = resolution.samplingLimit() / 2.0;
328 double gridStep = resolution.resolutionLimit() * gridFactor;
329
330 int nX = (int) Math.floor(crystal.a / gridStep) + 1;
331 int nY = (int) Math.floor(crystal.b / gridStep) + 1;
332 int nZ = (int) Math.floor(crystal.c / gridStep) + 1;
333 if (nX % 2 != 0) {
334 nX += 1;
335 }
336
337 while (!Complex.preferredDimension(nX)) {
338 nX += 2;
339 }
340 while (!Complex.preferredDimension(nY)) {
341 nY += 1;
342 }
343 while (!Complex.preferredDimension(nZ)) {
344 nZ += 1;
345 }
346
347 fftX = nX;
348 fftY = nY;
349 fftZ = nZ;
350 halfFFTX = nX / 2;
351 ifftX = 1.0 / (double) fftX;
352 ifftY = 1.0 / (double) fftY;
353 ifftZ = 1.0 / (double) fftZ;
354 fftScale = crystal.volume;
355 complexFFT3DSpace = fftX * fftY * fftZ * 2;
356 densityGrid = new double[complexFFT3DSpace];
357
358
359 optSliceWeightAtomic = new SharedIntegerArray(fftZ);
360 optSliceWeightSolvent = new SharedIntegerArray(fftZ);
361 optSliceWeightBulkSolvent = new SharedIntegerArray(fftZ);
362 previousOptSliceWeightAtomic = new int[fftZ];
363 previousOptSliceWeightSolvent = new int[fftZ];
364 previousOptSliceWeightBulkSolvent = new int[fftZ];
365
366 atomicSliceSchedule = new SliceSchedule(threadCount, fftZ);
367 solventSliceSchedule = new SliceSchedule(threadCount, fftZ);
368 bulkSolventSliceSchedule = new SliceSchedule(threadCount, fftZ);
369
370
371 optRowWeightAtomic = new SharedIntegerArray(fftZ * fftY);
372 optRowWeightSolvent = new SharedIntegerArray(fftZ * fftY);
373 optRowWeightBulkSolvent = new SharedIntegerArray(fftZ * fftY);
374 previousOptRowWeightAtomic = new int[fftZ * fftY];
375 previousOptRowWeightSolvent = new int[fftZ * fftY];
376 previousOptRowWeightBulkSolvent = new int[fftZ * fftY];
377 atomicRowSchedule = new RowSchedule(threadCount, fftZ, fftY);
378 solventRowSchedule = new RowSchedule(threadCount, fftZ, fftY);
379 bulkSolventRowSchedule = new RowSchedule(threadCount, fftZ, fftY);
380
381 optAtomicGradientWeight = new SharedIntegerArray(nAtoms);
382 previousOptAtomicGradientWeight = new int[nAtoms];
383 atomicGradientSchedule = new GradientSchedule(threadCount, nAtoms);
384
385 if (solvent) {
386 bAdd = 0.0;
387 atomFormFactors = null;
388 double vdwr;
389 switch (solventModel) {
390 case BINARY:
391 solventA = 1.0;
392 solventB = 1.0;
393 solventGrid = new double[complexFFT3DSpace];
394 solventFormFactors = new SolventBinaryFormFactor[bulkNSymm][nAtoms];
395 for (int iSymm = 0; iSymm < bulkNSymm; iSymm++) {
396 for (int i = 0; i < nAtoms; i++) {
397 vdwr = atoms[i].getVDWType().radius * 0.5;
398 solventFormFactors[iSymm][i] = new SolventBinaryFormFactor(atoms[i], vdwr + solventA);
399 }
400 }
401 break;
402 case GAUSSIAN:
403 solventA = 11.5;
404 solventB = 0.55;
405 solventGrid = new double[complexFFT3DSpace];
406 solventFormFactors = new SolventGaussFormFactor[bulkNSymm][nAtoms];
407 for (int iSymm = 0; iSymm < bulkNSymm; iSymm++) {
408 for (int i = 0; i < nAtoms; i++) {
409 vdwr = atoms[i].getVDWType().radius * 0.5;
410 solventFormFactors[iSymm][i] = new SolventGaussFormFactor(atoms[i], vdwr * solventB);
411 }
412 }
413 break;
414 case POLYNOMIAL:
415 solventA = 0.0;
416 solventB = 0.8;
417 solventGrid = new double[complexFFT3DSpace];
418 solventFormFactors = new SolventPolyFormFactor[bulkNSymm][nAtoms];
419 for (int iSymm = 0; iSymm < bulkNSymm; iSymm++) {
420 for (int i = 0; i < nAtoms; i++) {
421 vdwr = atoms[i].getVDWType().radius * 0.5;
422 solventFormFactors[iSymm][i] =
423 new SolventPolyFormFactor(atoms[i], vdwr + solventA, solventB);
424 }
425 }
426 break;
427 case NONE:
428 default:
429 solventGrid = null;
430 solventFormFactors = null;
431 break;
432 }
433 } else {
434 bAdd = 2.0;
435 if (neutron) {
436 atomFormFactors = new NeutronFormFactor[bulkNSymm][nAtoms];
437 for (int iSymm = 0; iSymm < bulkNSymm; iSymm++) {
438 for (int i = 0; i < nAtoms; i++) {
439 atomFormFactors[iSymm][i] = new NeutronFormFactor(atoms[i], bAdd);
440 }
441 }
442 } else {
443 atomFormFactors = new XRayFormFactor[bulkNSymm][nAtoms];
444 for (int iSymm = 0; iSymm < bulkNSymm; iSymm++) {
445 for (int i = 0; i < nAtoms; i++) {
446 atomFormFactors[iSymm][i] = new XRayFormFactor(atoms[i], useThreeGaussians, bAdd);
447 }
448 }
449 }
450 solventGrid = null;
451 solventFormFactors = null;
452 }
453
454
455 double aRad = getaRad(atoms, solventModel);
456 aRadGrid = (int) Math.floor(aRad * nX / crystal.a) + 1;
457
458
459 coordinates = new double[bulkNSymm][3][nAtoms];
460 List<SymOp> symops = crystal.spaceGroup.symOps;
461 double[] xyz = new double[3];
462 for (int i = 0; i < bulkNSymm; i++) {
463 double[] x = coordinates[i][0];
464 double[] y = coordinates[i][1];
465 double[] z = coordinates[i][2];
466 for (int j = 0; j < nAtoms; j++) {
467 Atom aj = atoms[j];
468 crystal.applySymOp(aj.getXYZ(null), xyz, symops.get(i));
469 x[j] = xyz[0];
470 y[j] = xyz[1];
471 z[j] = xyz[2];
472 }
473 }
474
475 if (logger.isLoggable(Level.INFO)) {
476 StringBuilder sb = new StringBuilder();
477 if (solvent) {
478 sb.append(" Bulk Solvent Grid\n");
479 sb.append(format(" Bulk solvent model type: %8s\n", solventModel));
480 } else {
481 sb.append(" Atomic Scattering Grid\n");
482 }
483 sb.append(format(" Form factor grid points: %8d\n", aRadGrid * 2));
484 sb.append(format(" Form factor grid diameter: %8.3f\n", aRad * 2));
485 sb.append(format(" Grid density: %8.3f\n", 1.0 / gridStep));
486 sb.append(format(" Grid dimensions: (%3d,%3d,%3d)\n", fftX, fftY, fftZ));
487 sb.append(format(" Grid method: %8s\n", gridMethod));
488 logger.info(sb.toString());
489 }
490
491 if (solvent) {
492 int minWork = nSymm;
493 if (solventModel != SolventModel.NONE) {
494 solventGradientRegion =
495 new SolventGradientRegion(RefinementMode.COORDINATES_AND_BFACTORS_AND_OCCUPANCIES);
496
497 atomicGradientRegion = null;
498 atomicSliceLoops = null;
499 atomicRowLoops = null;
500 switch (gridMethod) {
501 case SPATIAL:
502 atomicDensityRegion =
503 new SpatialDensityRegion(
504 fftX,
505 fftY,
506 fftZ,
507 densityGrid,
508 (aRadGrid + 2) * 2,
509 nSymm,
510 minWork,
511 threadCount,
512 crystal,
513 atoms,
514 coordinates);
515 solventDensityRegion =
516 new BulkSolventDensityRegion(
517 fftX,
518 fftY,
519 fftZ,
520 solventGrid,
521 (aRadGrid + 2) * 2,
522 bulkNSymm,
523 minWork,
524 threadCount,
525 crystal,
526 atoms,
527 coordinates,
528 4.0,
529 parallelTeam);
530 if (solventModel == SolventModel.GAUSSIAN) {
531 atomicDensityRegion.setInitValue(0.0);
532 } else {
533 atomicDensityRegion.setInitValue(1.0);
534 }
535 SolventDensityLoop[] solventDensityLoops = new SolventDensityLoop[threadCount];
536 SolventDensityLoop[] bulkSolventDensityLoops = new SolventDensityLoop[threadCount];
537 for (int i = 0; i < threadCount; i++) {
538 solventDensityLoops[i] = new SolventDensityLoop(atomicDensityRegion);
539 bulkSolventDensityLoops[i] = new SolventDensityLoop(solventDensityRegion);
540 }
541 atomicDensityRegion.setDensityLoop(solventDensityLoops);
542 solventDensityRegion.setDensityLoop(bulkSolventDensityLoops);
543
544 atomicSliceRegion = null;
545 solventSliceRegion = null;
546 bulkSolventSliceLoops = null;
547 solventSliceLoops = null;
548
549 atomicRowRegion = null;
550 solventRowRegion = null;
551 bulkSolventRowLoops = null;
552 solventRowLoops = null;
553 break;
554 case ROW:
555 atomicRowRegion =
556 new RowRegion(
557 fftX, fftY, fftZ, densityGrid, nSymm, threadCount, atoms, coordinates);
558 solventRowRegion =
559 new BulkSolventRowRegion(
560 fftX,
561 fftY,
562 fftZ,
563 solventGrid,
564 bulkNSymm,
565 threadCount,
566 crystal,
567 atoms,
568 coordinates,
569 4.0,
570 parallelTeam);
571 if (solventModel == SolventModel.GAUSSIAN) {
572 atomicRowRegion.setInitValue(0.0);
573 } else {
574 atomicRowRegion.setInitValue(1.0);
575 }
576 solventRowLoops = new SolventRowLoop[threadCount];
577 bulkSolventRowLoops = new BulkSolventRowLoop[threadCount];
578 for (int i = 0; i < threadCount; i++) {
579 solventRowLoops[i] = new SolventRowLoop(atomicRowRegion);
580 bulkSolventRowLoops[i] = new BulkSolventRowLoop(solventRowRegion);
581 }
582 atomicRowRegion.setDensityLoop(solventRowLoops);
583 solventRowRegion.setDensityLoop(bulkSolventRowLoops);
584 atomicDensityRegion = null;
585 solventDensityRegion = null;
586 atomicSliceRegion = null;
587 solventSliceRegion = null;
588 bulkSolventSliceLoops = null;
589 solventSliceLoops = null;
590 break;
591 case SLICE:
592 default:
593 atomicSliceRegion =
594 new SliceRegion(
595 fftX, fftY, fftZ, densityGrid, nSymm, threadCount, atoms, coordinates);
596 solventSliceRegion =
597 new BulkSolventSliceRegion(
598 fftX,
599 fftY,
600 fftZ,
601 solventGrid,
602 bulkNSymm,
603 threadCount,
604 crystal,
605 atoms,
606 coordinates,
607 4.0,
608 parallelTeam);
609 if (solventModel == SolventModel.GAUSSIAN) {
610 atomicSliceRegion.setInitValue(0.0);
611 } else {
612 atomicSliceRegion.setInitValue(1.0);
613 }
614 solventSliceLoops = new SolventSliceLoop[threadCount];
615 bulkSolventSliceLoops = new BulkSolventSliceLoop[threadCount];
616 for (int i = 0; i < threadCount; i++) {
617 solventSliceLoops[i] = new SolventSliceLoop(atomicSliceRegion);
618 bulkSolventSliceLoops[i] = new BulkSolventSliceLoop(solventSliceRegion);
619 }
620 atomicSliceRegion.setDensityLoop(solventSliceLoops);
621 solventSliceRegion.setDensityLoop(bulkSolventSliceLoops);
622
623 atomicDensityRegion = null;
624 solventDensityRegion = null;
625 atomicRowRegion = null;
626 solventRowRegion = null;
627 bulkSolventRowLoops = null;
628 solventRowLoops = null;
629 }
630 } else {
631 atomicDensityRegion = null;
632 solventDensityRegion = null;
633 atomicGradientRegion = null;
634 solventGradientRegion = null;
635 atomicSliceRegion = null;
636 atomicSliceLoops = null;
637 solventSliceRegion = null;
638 bulkSolventSliceLoops = null;
639 solventSliceLoops = null;
640 atomicRowRegion = null;
641 atomicRowLoops = null;
642 solventRowRegion = null;
643 bulkSolventRowLoops = null;
644 solventRowLoops = null;
645 }
646 } else {
647 atomicGradientRegion =
648 new AtomicGradientRegion(RefinementMode.COORDINATES_AND_BFACTORS_AND_OCCUPANCIES);
649
650
651 solventDensityRegion = null;
652 solventSliceRegion = null;
653 solventSliceLoops = null;
654 bulkSolventSliceLoops = null;
655 solventRowRegion = null;
656 solventRowLoops = null;
657 bulkSolventRowLoops = null;
658 solventGradientRegion = null;
659
660
661
662
663
664 switch (gridMethod) {
665 case SPATIAL:
666 int minWork = nSymm;
667 atomicDensityRegion =
668 new SpatialDensityRegion(
669 fftX,
670 fftY,
671 fftZ,
672 densityGrid,
673 (aRadGrid + 2) * 2,
674 nSymm,
675 minWork,
676 threadCount,
677 crystal,
678 atoms,
679 coordinates);
680 atomicDensityRegion.setInitValue(0.0);
681 AtomicDensityLoop[] atomicDensityLoops = new AtomicDensityLoop[threadCount];
682 for (int i = 0; i < threadCount; i++) {
683 atomicDensityLoops[i] = new AtomicDensityLoop(atomicDensityRegion);
684 }
685 atomicDensityRegion.setDensityLoop(atomicDensityLoops);
686 atomicSliceRegion = null;
687 atomicSliceLoops = null;
688 atomicRowRegion = null;
689 atomicRowLoops = null;
690 break;
691
692 case ROW:
693 atomicRowRegion =
694 new RowRegion(fftX, fftY, fftZ, densityGrid, nSymm, threadCount, atoms, coordinates);
695 atomicRowRegion.setInitValue(0.0);
696 atomicRowLoops = new AtomicRowLoop[threadCount];
697 for (int i = 0; i < threadCount; i++) {
698 atomicRowLoops[i] = new AtomicRowLoop(atomicRowRegion);
699 }
700 atomicRowRegion.setDensityLoop(atomicRowLoops);
701 atomicDensityRegion = null;
702 atomicSliceRegion = null;
703 atomicSliceLoops = null;
704 break;
705
706 case SLICE:
707 default:
708 atomicSliceRegion =
709 new SliceRegion(
710 fftX, fftY, fftZ, densityGrid, nSymm, threadCount, atoms, coordinates);
711 atomicSliceRegion.setInitValue(0.0);
712 atomicSliceLoops = new AtomicSliceLoop[threadCount];
713 for (int i = 0; i < threadCount; i++) {
714 atomicSliceLoops[i] = new AtomicSliceLoop(atomicSliceRegion);
715 }
716 atomicSliceRegion.setDensityLoop(atomicSliceLoops);
717 atomicDensityRegion = null;
718 atomicRowRegion = null;
719 atomicRowLoops = null;
720 }
721 }
722
723 extractRegion = new ExtractRegion(threadCount);
724 babinetRegion = new BabinetRegion(threadCount);
725 solventGridRegion = new SolventGridRegion(threadCount);
726 initRegion = new InitRegion(threadCount);
727 solventScaleRegion = new SolventScaleRegion(threadCount);
728 atomicScaleRegion = new AtomicScaleRegion(threadCount);
729 complexFFT3D = new Complex3DParallel(fftX, fftY, fftZ, fftTeam);
730 }
731
732 private double getaRad(Atom[] atoms, SolventModel solventModel) {
733 double aRad = -1.0;
734 for (Atom a : atoms) {
735 double vdwr = a.getVDWType().radius * 0.5;
736 if (!solvent) {
737 aRad = Math.max(aRad, a.getFormFactorWidth());
738 } else {
739 switch (solventModel) {
740 case BINARY:
741 aRad = Math.max(aRad, vdwr + solventA + 0.2);
742 break;
743 case GAUSSIAN:
744 aRad = Math.max(aRad, vdwr * solventB + 2.0);
745 break;
746 case POLYNOMIAL:
747 aRad = Math.max(aRad, vdwr + solventB + 0.2);
748 break;
749 case NONE:
750 default:
751 aRad = 0.0;
752 }
753 }
754 }
755 return aRad;
756 }
757
758
759
760
761
762
763
764 public static SolventModel parseSolventModel(String str) {
765 try {
766 return SolventModel.valueOf(str.toUpperCase().replaceAll("\\s+", ""));
767 } catch (Exception e) {
768 logger.info(format(" Could not parse %s as a Solvent Model; defaulting to Polynomial.", str));
769 return SolventModel.POLYNOMIAL;
770 }
771 }
772
773
774
775
776
777
778
779
780
781
782
783 public void computeAtomicGradients(
784 double[][] hklData, int[] freer, int flag, RefinementMode refinementMode) {
785 computeAtomicGradients(hklData, freer, flag, refinementMode, false);
786 }
787
788
789
790
791
792
793
794 public void computeSolventDensity(double[][] hklData) {
795 computeSolventDensity(hklData, false);
796 }
797
798
799
800
801
802
803
804
805 public void densityNorm(double[] data, double[] meansd, boolean norm) {
806 double mean, sd;
807
808 mean = sd = 0.0;
809 int n = 0;
810 for (int k = 0; k < fftZ; k++) {
811 for (int j = 0; j < fftY; j++) {
812 for (int i = 0; i < fftX; i++) {
813 int index = interleavedIndex(i, j, k, fftX, fftY);
814 n++;
815 mean += (data[index] - mean) / n;
816 }
817 }
818 }
819
820 n = 0;
821 for (int k = 0; k < fftZ; k++) {
822 for (int j = 0; j < fftY; j++) {
823 for (int i = 0; i < fftX; i++) {
824 int index = interleavedIndex(i, j, k, fftX, fftY);
825 sd += pow(data[index] - mean, 2);
826 n++;
827 }
828 }
829 }
830 sd = sqrt(sd / n);
831
832 if (meansd != null) {
833 meansd[0] = mean;
834 meansd[1] = sd;
835 }
836
837 if (norm) {
838 double isd = 1.0 / sd;
839 for (int k = 0; k < fftZ; k++) {
840 for (int j = 0; j < fftY; j++) {
841 for (int i = 0; i < fftX; i++) {
842 int index = interleavedIndex(i, j, k, fftX, fftY);
843 data[index] = (data[index] - mean) * isd;
844 }
845 }
846 }
847 }
848 }
849
850
851
852
853
854
855 public double[] getDensityGrid() {
856 return densityGrid;
857 }
858
859
860
861
862
863
864 public double getWeight() {
865 return weight;
866 }
867
868
869
870
871
872
873 public void setWeight(double weight) {
874 this.weight = weight;
875 }
876
877
878
879
880
881
882 public double getXDim() {
883 return fftX;
884 }
885
886
887
888
889
890
891 public double getYDim() {
892 return fftY;
893 }
894
895
896
897
898
899
900 public double getZDim() {
901 return fftZ;
902 }
903
904
905
906
907
908
909 public void setCoordinates(double[] coords) {
910 assert (coords != null);
911 List<SymOp> symops = crystal.spaceGroup.symOps;
912 double[] xyz = new double[3];
913 double[] symXYZ = new double[3];
914 int index = 0;
915 for (int i = 0; i < nAtoms; i++) {
916 if (!atoms[i].isActive()) {
917
918 continue;
919 }
920 xyz[0] = coords[index++];
921 xyz[1] = coords[index++];
922 xyz[2] = coords[index++];
923 for (int j = 0; j < nSymm; j++) {
924 crystal.applySymOp(xyz, symXYZ, symops.get(j));
925 coordinates[j][0][i] = symXYZ[0];
926 coordinates[j][1][i] = symXYZ[1];
927 coordinates[j][2][i] = symXYZ[2];
928 }
929 }
930 }
931
932
933
934
935
936
937 void setNativeEnvironmentApproximation(boolean nativeEnvironmentApproximation) {
938 this.nativeEnvironmentApproximation = nativeEnvironmentApproximation;
939 }
940
941
942
943
944
945
946 double[] getSolventGrid() {
947 return solventGrid;
948 }
949
950
951
952
953
954
955 void setLambdaTerm(boolean lambdaTerm) {
956 this.lambdaTerm = lambdaTerm;
957 }
958
959
960
961
962
963
964
965 void setSolventAB(double a, double b) {
966 double vdwr;
967 this.solventA = a;
968 this.solventB = b;
969 if (!solvent) {
970 return;
971 }
972 switch (solventModel) {
973 case BINARY:
974 for (int iSymm = 0; iSymm < bulkNSymm; iSymm++) {
975 for (int i = 0; i < nAtoms; i++) {
976 vdwr = atoms[i].getVDWType().radius * 0.5;
977 solventFormFactors[iSymm][i] = new SolventBinaryFormFactor(atoms[i], vdwr + solventA);
978 }
979 }
980 break;
981 case GAUSSIAN:
982 for (int iSymm = 0; iSymm < bulkNSymm; iSymm++) {
983 for (int i = 0; i < nAtoms; i++) {
984 vdwr = atoms[i].getVDWType().radius * 0.5;
985 solventFormFactors[iSymm][i] = new SolventGaussFormFactor(atoms[i], vdwr * solventB);
986 }
987 }
988 break;
989 case POLYNOMIAL:
990 for (int iSymm = 0; iSymm < bulkNSymm; iSymm++) {
991 for (int i = 0; i < nAtoms; i++) {
992 vdwr = atoms[i].getVDWType().radius * 0.5;
993 solventFormFactors[iSymm][i] =
994 new SolventPolyFormFactor(atoms[i], vdwr + solventA, solventB);
995 }
996 }
997 break;
998 }
999 }
1000
1001
1002
1003
1004
1005
1006 void setUse3G(boolean useThreeGaussians) {
1007 this.useThreeGaussians = useThreeGaussians;
1008 if (solvent || neutron) {
1009 return;
1010 }
1011 for (int iSymm = 0; iSymm < bulkNSymm; iSymm++) {
1012 for (int i = 0; i < nAtoms; i++) {
1013 atomFormFactors[iSymm][i] = new XRayFormFactor(atoms[i], useThreeGaussians, bAdd);
1014 }
1015 }
1016 }
1017
1018
1019
1020
1021
1022
1023
1024 void deltaX(int n, double delta) {
1025 List<SymOp> symops = crystal.spaceGroup.symOps;
1026 double[] xyz = new double[3];
1027 double[] symxyz = new double[3];
1028 atoms[n].getXYZ(xyz);
1029 xyz[0] += delta;
1030 for (int i = 0; i < nSymm; i++) {
1031 crystal.applySymOp(xyz, symxyz, symops.get(i));
1032 coordinates[i][0][n] = symxyz[0];
1033 }
1034 }
1035
1036
1037
1038
1039
1040
1041
1042 void deltaY(int n, double delta) {
1043 List<SymOp> symops = crystal.spaceGroup.symOps;
1044 double[] xyz = new double[3];
1045 double[] symxyz = new double[3];
1046 atoms[n].getXYZ(xyz);
1047 xyz[1] += delta;
1048 for (int i = 0; i < nSymm; i++) {
1049 crystal.applySymOp(xyz, symxyz, symops.get(i));
1050 coordinates[i][1][n] = symxyz[1];
1051 }
1052 }
1053
1054
1055
1056
1057
1058
1059
1060 void deltaZ(int n, double delta) {
1061 List<SymOp> symops = crystal.spaceGroup.symOps;
1062 double[] xyz = new double[3];
1063 double[] symxyz = new double[3];
1064 atoms[n].getXYZ(xyz);
1065 xyz[2] += delta;
1066 for (int i = 0; i < nSymm; i++) {
1067 crystal.applySymOp(xyz, symxyz, symops.get(i));
1068 coordinates[i][2][n] = symxyz[2];
1069 }
1070 }
1071
1072
1073
1074
1075
1076
1077
1078 void computeDensity(double[][] hklData) {
1079 computeDensity(hklData, false);
1080 }
1081
1082
1083
1084
1085
1086
1087
1088
1089 void computeDensity(double[][] hklData, boolean print) {
1090 if (solvent) {
1091 if (solventModel != SolventModel.NONE) {
1092 computeSolventDensity(hklData, print);
1093 }
1094 } else {
1095 computeAtomicDensity(hklData, print);
1096 }
1097 }
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110 void computeAtomicGradients(
1111 double[][] hklData, int[] freer, int flag, RefinementMode refinementMode, boolean print) {
1112
1113 if (solvent && solventModel == SolventModel.NONE) {
1114 return;
1115 }
1116
1117
1118 for (int i = 0; i < complexFFT3DSpace; i++) {
1119 densityGrid[i] = 0.0;
1120 }
1121
1122 int nfree = 0;
1123 StringBuilder sb = new StringBuilder();
1124 long symtime = -System.nanoTime();
1125 int nsym = crystal.spaceGroup.symOps.size();
1126
1127 List<SymOp> symops = crystal.spaceGroup.symOps;
1128 ComplexNumber c = new ComplexNumber();
1129 ComplexNumber cj = new ComplexNumber();
1130 HKL ij = new HKL();
1131 for (HKL ih : reflectionList.hklList) {
1132 double[] fc = hklData[ih.getIndex()];
1133 if (Double.isNaN(fc[0])) {
1134 continue;
1135 }
1136
1137
1138 if (freer != null) {
1139 if (freer[ih.getIndex()] == flag) {
1140 nfree++;
1141 continue;
1142 }
1143 }
1144
1145 c.re(fc[0]);
1146 c.im(fc[1]);
1147
1148 c.timesIP(2.0 / fftScale);
1149
1150
1151 for (int j = 0; j < nsym; j++) {
1152 cj.copy(c);
1153 SymOp symOp = symops.get(j);
1154 applyTransSymRot(ih, ij, symOp);
1155 double shift = symOp.symPhaseShift(ih);
1156
1157 int h = mod(ij.getH(), fftX);
1158 int k = mod(ij.getK(), fftY);
1159 int l = mod(ij.getL(), fftZ);
1160
1161 if (h < halfFFTX + 1) {
1162 final int ii = interleavedIndex(h, k, l, fftX, fftY);
1163 cj.phaseShiftIP(shift);
1164 densityGrid[ii] += cj.re();
1165
1166
1167 densityGrid[ii + 1] += (-cj.im());
1168 } else {
1169 h = (fftX - h) % fftX;
1170 k = (fftY - k) % fftY;
1171 l = (fftZ - l) % fftZ;
1172 final int ii = interleavedIndex(h, k, l, fftX, fftY);
1173 cj.phaseShiftIP(shift);
1174 densityGrid[ii] += cj.re();
1175 densityGrid[ii + 1] += cj.im();
1176 }
1177 }
1178 }
1179 symtime += System.nanoTime();
1180
1181 long startTime = System.nanoTime();
1182 complexFFT3D.ifft(densityGrid);
1183 long fftTime = System.nanoTime() - startTime;
1184
1185
1186
1187
1188
1189 startTime = System.nanoTime();
1190 long permanentDensityTime = 0;
1191 try {
1192 if (solvent) {
1193 solventGradientRegion.setRefinementMode(refinementMode);
1194 parallelTeam.execute(solventGradientRegion);
1195 } else {
1196 for (int i = 0; i < nAtoms; i++) {
1197 optAtomicGradientWeight.set(i, 0);
1198 }
1199 atomicGradientSchedule.updateWeights(previousOptAtomicGradientWeight);
1200 atomicGradientRegion.setRefinementMode(refinementMode);
1201 parallelTeam.execute(atomicGradientRegion);
1202 for (int i = 0; i < nAtoms; i++) {
1203 previousOptAtomicGradientWeight[i] = optAtomicGradientWeight.get(i);
1204 }
1205 }
1206 permanentDensityTime = System.nanoTime() - startTime;
1207 } catch (Exception e) {
1208 String message = "Exception computing atomic gradients.";
1209 logger.log(Level.SEVERE, message, e);
1210 }
1211
1212 if (solvent) {
1213 sb.append(format(" Solvent symmetry: %8.3f\n", symtime * toSeconds));
1214 sb.append(format(" Solvent inverse FFT: %8.3f\n", fftTime * toSeconds));
1215 sb.append(format(" Grid solvent gradients: %8.3f\n", permanentDensityTime * toSeconds));
1216 sb.append(format(" %d reflections ignored (cross validation set)\n", nfree));
1217 } else {
1218 sb.append(format(" Atomic symmetry: %8.3f\n", symtime * toSeconds));
1219 sb.append(format(" Atomic inverse FFT: %8.3f\n", fftTime * toSeconds));
1220 sb.append(format(" Grid atomic gradients: %8.3f\n", permanentDensityTime * toSeconds));
1221 sb.append(format(" %d reflections ignored (cross validation set)\n", nfree));
1222 }
1223
1224 if (logger.isLoggable(Level.INFO) && print) {
1225 logger.info(sb.toString());
1226 }
1227 }
1228
1229
1230
1231
1232
1233
1234
1235 void computeAtomicDensity(double[][] hklData) {
1236 computeAtomicDensity(hklData, false);
1237 }
1238
1239
1240
1241
1242
1243
1244
1245
1246 private void computeAtomicDensity(double[][] hklData, boolean print) {
1247
1248 long initTime = -System.nanoTime();
1249 try {
1250 initRegion.setHKL(hklData);
1251 parallelTeam.execute(initRegion);
1252 } catch (Exception e) {
1253 String message = "Fatal exception zeroing out reflection data.";
1254 logger.log(Level.SEVERE, message, e);
1255 }
1256 initTime += System.nanoTime();
1257
1258
1259 long atomicGridTime = -System.nanoTime();
1260 try {
1261 switch (gridMethod) {
1262 case SPATIAL:
1263 atomicDensityRegion.assignAtomsToCells();
1264 parallelTeam.execute(atomicDensityRegion);
1265 break;
1266
1267 case ROW:
1268 for (int i = 0; i < fftZ * fftY; i++) {
1269 optRowWeightAtomic.set(i, 0);
1270 }
1271 atomicRowSchedule.updateWeights(previousOptRowWeightAtomic);
1272 parallelTeam.execute(atomicRowRegion);
1273 for (int i = 0; i < fftZ * fftY; i++) {
1274 previousOptRowWeightAtomic[i] = optRowWeightAtomic.get(i);
1275 }
1276 break;
1277 case SLICE:
1278 default:
1279 for (int i = 0; i < fftZ; i++) {
1280 optSliceWeightAtomic.set(i, 0);
1281 }
1282 atomicSliceSchedule.updateWeights(previousOptSliceWeightAtomic);
1283 parallelTeam.execute(atomicSliceRegion);
1284 for (int i = 0; i < fftZ; i++) {
1285 previousOptSliceWeightAtomic[i] = optSliceWeightAtomic.get(i);
1286 }
1287 break;
1288 }
1289 } catch (Exception e) {
1290 String message = "Fatal exception evaluating atomic electron density.";
1291 logger.log(Level.SEVERE, message, e);
1292 }
1293 atomicGridTime += System.nanoTime();
1294
1295
1296 long startTime = System.nanoTime();
1297 complexFFT3D.fft(densityGrid);
1298 long fftTime = System.nanoTime() - startTime;
1299
1300
1301 long symTime = -System.nanoTime();
1302 try {
1303 extractRegion.setHKL(hklData);
1304 parallelTeam.execute(extractRegion);
1305 double scale = (fftScale) / (fftX * fftY * fftZ);
1306 atomicScaleRegion.setHKL(hklData);
1307 atomicScaleRegion.setScale(scale);
1308 parallelTeam.execute(atomicScaleRegion);
1309 } catch (Exception e) {
1310 String message = "Fatal exception extracting structure factors.";
1311 logger.log(Level.SEVERE, message, e);
1312 }
1313 symTime += System.nanoTime();
1314
1315 if (logger.isLoggable(Level.INFO) && print) {
1316 StringBuilder sb = new StringBuilder();
1317 sb.append(format("\n Fc Initialization: %8.4f\n", initTime * toSeconds));
1318 sb.append(format(" Atomic grid density: %8.4f\n", atomicGridTime * toSeconds));
1319 sb.append(format(" Atomic FFT: %8.4f\n", fftTime * toSeconds));
1320 sb.append(format(" Atomic symmetry & scaling: %8.4f\n", symTime * toSeconds));
1321 logger.info(sb.toString());
1322
1323 StringBuilder ASB = new StringBuilder();
1324 switch (gridMethod) {
1325 case ROW:
1326 double atomicRowTotal = 0;
1327 double atomicRowMin = atomicRowLoops[0].getThreadTime();
1328 double atomicRowMax = 0;
1329 double atomicRowWeightTotal = 0;
1330
1331 for (int i = 0; i < threadCount; i++) {
1332 atomicRowTotal += atomicRowLoops[i].getThreadTime();
1333 atomicRowMax = max(atomicRowLoops[i].getThreadTime(), atomicRowMax);
1334 atomicRowMin = min(atomicRowLoops[i].getThreadTime(), atomicRowMin);
1335 atomicRowWeightTotal += atomicRowLoops[i].getThreadWeight();
1336 }
1337
1338
1339 ASB.append(
1340 format(
1341 "\n RowLoop (Atomic): %7.5f (sec) | Total Weight: %7.0f\n",
1342 atomicGridTime * toSeconds, atomicRowWeightTotal));
1343 ASB.append(
1344 " Thread LoopTime Balance(%) Normalized | Rows Weight Balance(%) Normalized\n");
1345
1346
1347 if (atomicRowWeightTotal == 0) {
1348 atomicRowWeightTotal = 1;
1349 }
1350
1351 for (int i = 0; i < threadCount; i++) {
1352 ASB.append(
1353 format(
1354 " %3d %7.5f %7.1f %7.1f | %7d %7d %7.1f %7.1f\n",
1355 i, atomicRowLoops[i].getThreadTime() * toSeconds,
1356 ((atomicRowLoops[i].getThreadTime()) / (atomicRowTotal)) * 100.00,
1357 ((atomicRowLoops[i].getThreadTime()) / (atomicRowTotal)) * (100.00
1358 * threadCount),
1359 atomicRowLoops[i].getNumberofSlices(),
1360 atomicRowLoops[i].getThreadWeight(),
1361 100.00 * (atomicRowLoops[i].getThreadWeight() / atomicRowWeightTotal),
1362 (100.00 * threadCount) * (atomicRowLoops[i].getThreadWeight()
1363 / atomicRowWeightTotal)));
1364 }
1365 ASB.append(format(" Min %7.5f\n", atomicRowMin * toSeconds));
1366 ASB.append(format(" Max %7.5f\n", atomicRowMax * toSeconds));
1367 ASB.append(format(" Delta %7.5f\n", (atomicRowMax - atomicRowMin) * toSeconds));
1368 logger.info(ASB.toString());
1369
1370 break;
1371 case SPATIAL:
1372 break;
1373 case SLICE:
1374 double atomicSliceTotal = 0;
1375 double atomicSliceMin = atomicSliceLoops[0].getThreadTime();
1376 double atomicSliceMax = 0;
1377 double atomicSliceWeightTotal = 0;
1378
1379 for (int i = 0; i < threadCount; i++) {
1380 atomicSliceTotal += atomicSliceLoops[i].getThreadTime();
1381 atomicSliceMax = max(atomicSliceLoops[i].getThreadTime(), atomicSliceMax);
1382 atomicSliceMin = min(atomicSliceLoops[i].getThreadTime(), atomicSliceMin);
1383 atomicSliceWeightTotal += atomicSliceLoops[i].getThreadWeight();
1384 }
1385
1386
1387 ASB.append(
1388 format(
1389 "\n SliceLoop (Atomic): %7.5f (sec) | Total Weight: %7.0f\n",
1390 atomicGridTime * toSeconds, atomicSliceWeightTotal));
1391 ASB.append(
1392 " Thread LoopTime Balance(%) Normalized | Slices Weight Balance(%) Normalized\n");
1393
1394
1395 if (atomicSliceWeightTotal == 0) {
1396 atomicSliceWeightTotal = 1;
1397 }
1398
1399 for (int i = 0; i < threadCount; i++) {
1400 ASB.append(
1401 format(
1402 " %3d %7.5f %7.1f %7.1f | %7d %7d %7.1f %7.1f\n",
1403 i,
1404 atomicSliceLoops[i].getThreadTime() * toSeconds,
1405 ((atomicSliceLoops[i].getThreadTime()) / (atomicSliceTotal)) * 100.00,
1406 ((atomicSliceLoops[i].getThreadTime()) / (atomicSliceTotal)) * (100.00
1407 * threadCount),
1408 atomicSliceLoops[i].getNumberofSlices(),
1409 atomicSliceLoops[i].getThreadWeight(),
1410 100.00 * (atomicSliceLoops[i].getThreadWeight() / atomicSliceWeightTotal),
1411 (100.00 * threadCount) * (atomicSliceLoops[i].getThreadWeight()
1412 / atomicSliceWeightTotal)));
1413 }
1414 ASB.append(format(" Min %7.5f\n", atomicSliceMin * toSeconds));
1415 ASB.append(format(" Max %7.5f\n", atomicSliceMax * toSeconds));
1416 ASB.append(format(" Delta %7.5f\n", (atomicSliceMax - atomicSliceMin) * toSeconds));
1417 logger.info(ASB.toString());
1418 break;
1419 default:
1420 }
1421 }
1422 }
1423
1424
1425
1426
1427
1428
1429
1430
1431 private void computeSolventDensity(double[][] hklData, boolean print) {
1432
1433 long initTime = -System.nanoTime();
1434 try {
1435 initRegion.setHKL(hklData);
1436 parallelTeam.execute(initRegion);
1437 } catch (Exception e) {
1438 String message = "Fatal exception initializing structure factors.";
1439 logger.log(Level.SEVERE, message, e);
1440 }
1441 initTime += System.nanoTime();
1442
1443 long solventGridTime = -System.nanoTime();
1444 try {
1445 switch (gridMethod) {
1446 case SPATIAL:
1447 atomicDensityRegion.assignAtomsToCells();
1448 parallelTeam.execute(atomicDensityRegion);
1449 break;
1450 case ROW:
1451 for (int i = 0; i < fftZ * fftY; i++) {
1452 optRowWeightSolvent.set(i, 0);
1453 }
1454 solventRowSchedule.updateWeights(previousOptRowWeightSolvent);
1455 parallelTeam.execute(atomicRowRegion);
1456 for (int i = 0; i < fftZ * fftY; i++) {
1457 previousOptRowWeightSolvent[i] = optRowWeightSolvent.get(i);
1458 }
1459 break;
1460 case SLICE:
1461 default:
1462 for (int i = 0; i < fftZ; i++) {
1463 optSliceWeightSolvent.set(i, 0);
1464 }
1465 solventSliceSchedule.updateWeights(previousOptSliceWeightSolvent);
1466 parallelTeam.execute(atomicSliceRegion);
1467 for (int i = 0; i < fftZ; i++) {
1468 previousOptSliceWeightSolvent[i] = optSliceWeightSolvent.get(i);
1469 }
1470 break;
1471 }
1472 } catch (Exception e) {
1473 String message = "Fatal exception evaluating solvent electron density.";
1474 logger.log(Level.SEVERE, message, e);
1475 }
1476 solventGridTime += System.nanoTime();
1477
1478 long expTime = -System.nanoTime();
1479 if (solventModel == SolventModel.BINARY) {
1480
1481
1482
1483
1484 int nmap = densityGrid.length;
1485 arraycopy(densityGrid, 0, solventGrid, 0, nmap);
1486
1487 double[] xyz = {solventB, solventB, solventB};
1488 double[] uvw = new double[3];
1489 crystal.toFractionalCoordinates(xyz, uvw);
1490 final double frx = fftX * uvw[0];
1491 final int ifrx = (int) frx;
1492 final double fry = fftY * uvw[1];
1493 final int ifry = (int) fry;
1494 final double frz = fftZ * uvw[2];
1495 final int ifrz = (int) frz;
1496 for (int k = 0; k < fftZ; k++) {
1497 for (int j = 0; j < fftY; j++) {
1498 for (int i = 0; i < fftX; i++) {
1499 final int ii = interleavedIndex(i, j, k, fftX, fftY);
1500 if (densityGrid[ii] == 1.0) {
1501 continue;
1502 }
1503 outerLoop:
1504 for (int iz = k - ifrz; iz <= k + ifrz; iz++) {
1505 int giz = mod(iz, fftZ);
1506 for (int iy = j - ifry; iy <= j + ifry; iy++) {
1507 int giy = mod(iy, fftY);
1508 for (int ix = i - ifrx; ix <= i + ifrx; ix++) {
1509 int gix = mod(ix, fftX);
1510 final int jj = interleavedIndex(gix, giy, giz, fftX, fftY);
1511 if (densityGrid[jj] == 1.0) {
1512 solventGrid[ii] = 1.0;
1513 break outerLoop;
1514 }
1515 }
1516 }
1517 }
1518 }
1519 }
1520 }
1521
1522
1523 arraycopy(solventGrid, 0, densityGrid, 0, nmap);
1524 }
1525
1526
1527 ArrayCopyRegion arrayCopyRegion = new ArrayCopyRegion();
1528 try {
1529 parallelTeam.execute(arrayCopyRegion);
1530 } catch (Exception e) {
1531 logger.info(e.toString());
1532 }
1533
1534
1535 try {
1536 parallelTeam.execute(babinetRegion);
1537 } catch (Exception e) {
1538 String message = "Fatal exception Babinet Principle.";
1539 logger.log(Level.SEVERE, message, e);
1540 }
1541 expTime += System.nanoTime();
1542
1543 long fftTime = -System.nanoTime();
1544 complexFFT3D.fft(densityGrid);
1545 fftTime += System.nanoTime();
1546
1547
1548 long symTime = -System.nanoTime();
1549 try {
1550 extractRegion.setHKL(hklData);
1551 parallelTeam.execute(extractRegion);
1552 final double scale = (fftScale) / (fftX * fftY * fftZ);
1553 solventScaleRegion.setHKL(hklData);
1554 solventScaleRegion.setScale(scale);
1555 parallelTeam.execute(solventScaleRegion);
1556 } catch (Exception e) {
1557 String message = "Fatal exception extracting structure factors.";
1558 logger.log(Level.SEVERE, message, e);
1559 }
1560 symTime += System.nanoTime();
1561
1562
1563 long solventDensityTime = -System.nanoTime();
1564 try {
1565 switch (gridMethod) {
1566 case SPATIAL:
1567 solventDensityRegion.assignAtomsToCells();
1568 parallelTeam.execute(solventDensityRegion);
1569 break;
1570
1571 case ROW:
1572 for (int i = 0; i < (fftZ * fftY); i++) {
1573 optRowWeightBulkSolvent.set(i, 0);
1574 }
1575 bulkSolventRowSchedule.updateWeights(previousOptRowWeightBulkSolvent);
1576 parallelTeam.execute(solventRowRegion);
1577 for (int i = 0; i < fftZ * fftY; i++) {
1578 previousOptRowWeightBulkSolvent[i] = optRowWeightBulkSolvent.get(i);
1579 }
1580 break;
1581 case SLICE:
1582 default:
1583 for (int i = 0; i < (fftZ); i++) {
1584 optSliceWeightBulkSolvent.set(i, 0);
1585 }
1586
1587 bulkSolventSliceSchedule.updateWeights(previousOptSliceWeightBulkSolvent);
1588 parallelTeam.execute(solventSliceRegion);
1589 for (int i = 0; i < fftZ; i++) {
1590 previousOptSliceWeightBulkSolvent[i] = optSliceWeightBulkSolvent.get(i);
1591 }
1592 break;
1593 }
1594 } catch (Exception e) {
1595 String message = "Fatal exception evaluating solvent electron density.";
1596 logger.log(Level.SEVERE, message, e);
1597 }
1598 solventDensityTime += System.nanoTime();
1599
1600 long solventExpTime = 0;
1601 if (solventModel == SolventModel.GAUSSIAN) {
1602 solventExpTime = -System.nanoTime();
1603 try {
1604 parallelTeam.execute(solventGridRegion);
1605 } catch (Exception e) {
1606 String message = "Fatal exception solvent grid region.";
1607 logger.log(Level.SEVERE, message, e);
1608 }
1609 solventExpTime += System.nanoTime();
1610 }
1611
1612 if (logger.isLoggable(Level.INFO) && print) {
1613 StringBuilder sb = new StringBuilder();
1614 sb.append(format("\n Fc Initialization: %8.4f\n", initTime * toSeconds));
1615 sb.append(
1616 format(" Solvent grid density: %8.4f\n", solventGridTime * toSeconds));
1617 sb.append(format(" Bulk solvent exponentiation: %8.4f\n", expTime * toSeconds));
1618 sb.append(format(" Solvent FFT: %8.4f\n", fftTime * toSeconds));
1619 sb.append(format(" Solvent symmetry & scaling: %8.4f\n", symTime * toSeconds));
1620 sb.append(
1621 format(
1622 " Solvent grid expansion: %8.4f\n", solventDensityTime * toSeconds));
1623 if (solventModel == SolventModel.GAUSSIAN) {
1624 sb.append(format(" Gaussian grid exponentiation: %8.4f\n", solventExpTime * toSeconds));
1625 }
1626 logger.info(sb.toString());
1627
1628 StringBuilder SSB = new StringBuilder();
1629 StringBuilder BSSB = new StringBuilder();
1630 switch (gridMethod) {
1631 case ROW:
1632 double solventRowTotal = 0;
1633 double solventRowMin = solventRowLoops[0].getThreadTime();
1634 double solventRowMax = 0;
1635 double solventRowWeightTotal = 0;
1636
1637 double bulkSolventRowTotal = 0;
1638 double bulkSolventRowMin = bulkSolventRowLoops[0].getTimePerThread();
1639 double bulkSolventRowMax = 0;
1640 double bulkSolventRowWeightTotal = 0;
1641
1642 for (int i = 0; i < threadCount; i++) {
1643 solventRowTotal += solventRowLoops[i].getThreadTime();
1644 solventRowMax = max(solventRowLoops[i].getThreadTime(), solventRowMax);
1645 solventRowMin = min(solventRowLoops[i].getThreadTime(), solventRowMin);
1646 solventRowWeightTotal += solventRowLoops[i].getThreadWeight();
1647
1648 bulkSolventRowTotal += bulkSolventRowLoops[i].getTimePerThread();
1649 bulkSolventRowMax = max(bulkSolventRowLoops[i].getTimePerThread(), bulkSolventRowMax);
1650 bulkSolventRowMin = min(bulkSolventRowLoops[i].getTimePerThread(), bulkSolventRowMin);
1651 bulkSolventRowWeightTotal += bulkSolventRowLoops[i].getWeightPerThread()[i];
1652 }
1653
1654
1655 SSB.append(
1656 format(
1657 "\n RowLoop (Solvent): %7.5f (sec) | Total Weight: %7.0f\n",
1658 solventGridTime * toSeconds, solventRowWeightTotal));
1659 SSB.append(
1660 " Thread LoopTime Balance(%) Normalized | Rows Weight Balance(%) Normalized\n");
1661
1662
1663 if (solventRowWeightTotal == 0) {
1664 solventRowWeightTotal = 1;
1665 }
1666
1667 for (int i = 0; i < threadCount; i++) {
1668 SSB.append(
1669 format(
1670 " %3d %7.5f %7.1f %7.1f | %7d %7d %7.1f %7.1f\n",
1671 i, solventRowLoops[i].getThreadTime() * toSeconds,
1672 ((solventRowLoops[i].getThreadTime()) / (solventRowTotal)) * 100.00,
1673 ((solventRowLoops[i].getThreadTime()) / (solventRowTotal))
1674 * (100.00 * threadCount),
1675 solventRowLoops[i].getNumberofSlices(),
1676 solventRowLoops[i].getThreadWeight(),
1677 100.00 * (solventRowLoops[i].getThreadWeight() / solventRowWeightTotal),
1678 (100.00 * threadCount)
1679 * (solventRowLoops[i].getThreadWeight() / solventRowWeightTotal)));
1680 }
1681 SSB.append(format(" Min %7.5f\n", solventRowMin * toSeconds));
1682 SSB.append(format(" Max %7.5f\n", solventRowMax * toSeconds));
1683 SSB.append(format(" Delta %7.5f\n", (solventRowMax - solventRowMin) * toSeconds));
1684 logger.info(SSB.toString());
1685
1686
1687 BSSB.append(
1688 format(
1689 "\n RowLoop (Bulk Solvent): %7.5f (sec) | Total Weight: %7.0f\n",
1690 bulkSolventRowTotal * toSeconds, bulkSolventRowWeightTotal));
1691 BSSB.append(
1692 " Thread LoopTime Balance(%) Normalized | Rows Weight Balance(%) Normalized\n");
1693
1694
1695 if (bulkSolventRowWeightTotal == 0) {
1696 bulkSolventRowWeightTotal = 1;
1697 }
1698
1699 for (int i = 0; i < threadCount; i++) {
1700 BSSB.append(
1701 format(
1702 " %3d %7.5f %7.1f %7.1f | %7d %7d %7.1f %7.1f\n",
1703 i, bulkSolventRowLoops[i].getTimePerThread() * toSeconds,
1704 ((bulkSolventRowLoops[i].getTimePerThread()) / (bulkSolventRowTotal)) * 100.00,
1705 ((bulkSolventRowLoops[i].getTimePerThread()) / (bulkSolventRowTotal))
1706 * (100.00 * threadCount),
1707 bulkSolventRowLoops[i].getBoundsPerThread()[i],
1708 bulkSolventRowLoops[i].getWeightPerThread()[i],
1709 100.00
1710 * (bulkSolventRowLoops[i].getWeightPerThread()[i]
1711 / bulkSolventRowWeightTotal),
1712 (100.00 * threadCount)
1713 * (bulkSolventRowLoops[i].getWeightPerThread()[i]
1714 / bulkSolventRowWeightTotal)));
1715 }
1716 BSSB.append(format(" Min %7.5f\n", bulkSolventRowMin * toSeconds));
1717 BSSB.append(format(" Max %7.5f\n", bulkSolventRowMax * toSeconds));
1718 BSSB.append(
1719 format(" Delta %7.5f\n", (bulkSolventRowMax - bulkSolventRowMin) * toSeconds));
1720 logger.info(BSSB.toString());
1721 break;
1722 case SPATIAL:
1723 break;
1724 case SLICE:
1725 double solventSliceTotal = 0;
1726 double solventSliceMin = solventSliceLoops[0].getThreadTime();
1727 double solventSliceMax = 0;
1728 double solventSliceWeightTotal = 0;
1729
1730 double bulkSolventSliceTotal = 0;
1731 double bulkSolventSliceMin = bulkSolventSliceLoops[0].getTimePerThread();
1732 double bulkSolventSliceMax = 0;
1733 double bulkSolventSliceWeightTotal = 0;
1734
1735 for (int i = 0; i < threadCount; i++) {
1736 solventSliceTotal += solventSliceLoops[i].getThreadTime();
1737 solventSliceMax = max(solventSliceLoops[i].getThreadTime(), solventSliceMax);
1738 solventSliceMin = min(solventSliceLoops[i].getThreadTime(), solventSliceMin);
1739 solventSliceWeightTotal += solventSliceLoops[i].getThreadWeight();
1740
1741 bulkSolventSliceTotal += bulkSolventSliceLoops[i].getTimePerThread();
1742 bulkSolventSliceMax =
1743 max(bulkSolventSliceLoops[i].getTimePerThread(), bulkSolventSliceMax);
1744 bulkSolventSliceMin =
1745 min(bulkSolventSliceLoops[i].getTimePerThread(), bulkSolventSliceMin);
1746 bulkSolventSliceWeightTotal += bulkSolventSliceLoops[i].getWeightPerThread()[i];
1747 }
1748
1749
1750 SSB.append(
1751 format(
1752 "\n SliceLoop (Solvent): %7.5f (sec) | Total Weight: %7.0f\n",
1753 solventGridTime * toSeconds, solventSliceWeightTotal));
1754 SSB.append(
1755 " Thread LoopTime Balance(%) Normalized | Slices Weight Balance(%) Normalized\n");
1756
1757
1758 if (solventSliceWeightTotal == 0) {
1759 solventSliceWeightTotal = 1;
1760 }
1761
1762 for (int i = 0; i < threadCount; i++) {
1763 SSB.append(
1764 format(
1765 " %3d %7.5f %7.1f %7.1f | %7d %7d %7.1f %7.1f\n",
1766 i, solventSliceLoops[i].getThreadTime() * toSeconds,
1767 ((solventSliceLoops[i].getThreadTime()) / (solventSliceTotal)) * 100.00,
1768 ((solventSliceLoops[i].getThreadTime()) / (solventSliceTotal))
1769 * (100.00 * threadCount),
1770 solventSliceLoops[i].getNumberofSlices(),
1771 solventSliceLoops[i].getThreadWeight(),
1772 100.00 * (solventSliceLoops[i].getThreadWeight() / solventSliceWeightTotal),
1773 (100.00 * threadCount)
1774 * (solventSliceLoops[i].getThreadWeight() / solventSliceWeightTotal)));
1775 }
1776 SSB.append(format(" Min %7.5f\n", solventSliceMin * toSeconds));
1777 SSB.append(format(" Max %7.5f\n", solventSliceMax * toSeconds));
1778 SSB.append(
1779 format(" Delta %7.5f\n", (solventSliceMax - solventSliceMin) * toSeconds));
1780 logger.info(SSB.toString());
1781
1782
1783 BSSB.append(
1784 format(
1785 "\n SliceLoop (Bulk Solvent): %7.5f (sec) | Total Weight: %7.0f\n",
1786 bulkSolventSliceTotal * toSeconds, bulkSolventSliceWeightTotal));
1787 BSSB.append(
1788 " Thread LoopTime Balance(%) Normalized | Slices Weight Balance(%) Normalized\n");
1789
1790
1791 if (bulkSolventSliceWeightTotal == 0) {
1792 bulkSolventSliceWeightTotal = 1;
1793 }
1794
1795 for (int i = 0; i < threadCount; i++) {
1796 BSSB.append(
1797 format(
1798 " %3d %7.5f %7.1f %7.1f | %7d %7d %7.1f %7.1f\n",
1799 i, bulkSolventSliceLoops[i].getTimePerThread() * toSeconds,
1800 ((bulkSolventSliceLoops[i].getTimePerThread()) / (bulkSolventSliceTotal))
1801 * 100.00,
1802 ((bulkSolventSliceLoops[i].getTimePerThread()) / (bulkSolventSliceTotal))
1803 * (100.00 * threadCount),
1804 bulkSolventSliceLoops[i].getBoundsPerThread()[i],
1805 bulkSolventSliceLoops[i].getWeightPerThread()[i],
1806 100.00
1807 * (bulkSolventSliceLoops[i].getWeightPerThread()[i]
1808 / bulkSolventSliceWeightTotal),
1809 (100.00 * threadCount)
1810 * (bulkSolventSliceLoops[i].getWeightPerThread()[i]
1811 / bulkSolventSliceWeightTotal)));
1812 }
1813 BSSB.append(format(" Min %7.5f\n", bulkSolventSliceMin * toSeconds));
1814 BSSB.append(format(" Max %7.5f\n", bulkSolventSliceMax * toSeconds));
1815 BSSB.append(
1816 format(
1817 " Delta %7.5f\n", (bulkSolventSliceMax - bulkSolventSliceMin) * toSeconds));
1818 logger.info(BSSB.toString());
1819 break;
1820 default:
1821 }
1822 }
1823 }
1824
1825
1826 public enum SolventModel {
1827
1828
1829 NONE,
1830
1831 BINARY,
1832
1833 GAUSSIAN,
1834
1835 POLYNOMIAL
1836 }
1837
1838 public enum GridMethod {
1839 SPATIAL,
1840 SLICE,
1841 ROW
1842 }
1843
1844 private class AtomicDensityLoop extends SpatialDensityLoop {
1845
1846 final double[] xyz = new double[3];
1847 final double[] uvw = new double[3];
1848 final double[] xc = new double[3];
1849 final double[] xf = new double[3];
1850 final double[] grid;
1851
1852 AtomicDensityLoop(SpatialDensityRegion region) {
1853 super(region, region.getNsymm(), region.actualCount);
1854 grid = region.getGrid();
1855 }
1856
1857 @Override
1858 public void gridDensity(int iSymm, int n) {
1859 if (!atoms[n].getUse() && !nativeEnvironmentApproximation) {
1860 return;
1861 }
1862
1863 if (lambdaTerm && atoms[n].applyLambda()) {
1864 return;
1865 }
1866
1867 xyz[0] = coordinates[iSymm][0][n];
1868 xyz[1] = coordinates[iSymm][1][n];
1869 xyz[2] = coordinates[iSymm][2][n];
1870 FormFactor atomff = atomFormFactors[iSymm][n];
1871 crystal.toFractionalCoordinates(xyz, uvw);
1872 final int frad =
1873 Math.min(
1874 aRadGrid, (int) Math.floor(atoms[n].getFormFactorWidth() * fftX / crystal.a) + 1);
1875
1876
1877 final double frx = fftX * uvw[0];
1878 final int ifrx = (int) frx;
1879 final int ifrxu = ifrx + frad;
1880
1881 final double fry = fftY * uvw[1];
1882 final int ifry = (int) fry;
1883 final int ifryu = ifry + frad;
1884
1885 final double frz = fftZ * uvw[2];
1886 final int ifrz = (int) frz;
1887 final int ifrzu = ifrz + frad;
1888
1889 for (int iz = ifrz - frad; iz <= ifrzu; iz++) {
1890 int giz = mod(iz, fftZ);
1891 xf[2] = iz * ifftZ;
1892 for (int iy = ifry - frad; iy <= ifryu; iy++) {
1893 int giy = mod(iy, fftY);
1894 xf[1] = iy * ifftY;
1895 for (int ix = ifrx - frad; ix <= ifrxu; ix++) {
1896 int gix = mod(ix, fftX);
1897 xf[0] = ix * ifftX;
1898 crystal.toCartesianCoordinates(xf, xc);
1899 final int ii = interleavedIndex(gix, giy, giz, fftX, fftY);
1900 grid[ii] = atomff.rho(grid[ii], 1.0, xc);
1901 }
1902 }
1903 }
1904 }
1905 }
1906
1907 private class SolventDensityLoop extends SpatialDensityLoop {
1908
1909 final double[] xyz = new double[3];
1910 final double[] uvw = new double[3];
1911 final double[] xc = new double[3];
1912 final double[] xf = new double[3];
1913 final double[] grid;
1914
1915 SolventDensityLoop(SpatialDensityRegion region) {
1916 super(region, region.getNsymm(), region.actualCount);
1917 grid = region.getGrid();
1918 }
1919
1920 @Override
1921 public void gridDensity(int iSymm, int n) {
1922 if (!atoms[n].getUse() && !nativeEnvironmentApproximation) {
1923 return;
1924 }
1925
1926 if (lambdaTerm && atoms[n].applyLambda()) {
1927 return;
1928 }
1929
1930 xyz[0] = coordinates[iSymm][0][n];
1931 xyz[1] = coordinates[iSymm][1][n];
1932 xyz[2] = coordinates[iSymm][2][n];
1933 FormFactor solventff = solventFormFactors[iSymm][n];
1934 crystal.toFractionalCoordinates(xyz, uvw);
1935 double vdwr = atoms[n].getVDWType().radius * 0.5;
1936 int frad = aRadGrid;
1937 switch (solventModel) {
1938 case BINARY:
1939 frad =
1940 Math.min(aRadGrid, (int) Math.floor((vdwr + solventA + 0.2) * fftX / crystal.a) + 1);
1941 break;
1942 case GAUSSIAN:
1943 frad =
1944 Math.min(aRadGrid, (int) Math.floor((vdwr * solventB + 2.0) * fftX / crystal.a) + 1);
1945 break;
1946 case POLYNOMIAL:
1947 frad =
1948 Math.min(aRadGrid, (int) Math.floor((vdwr + solventB + 0.2) * fftX / crystal.a) + 1);
1949 break;
1950 case NONE:
1951 default:
1952 return;
1953 }
1954
1955
1956 final double frx = fftX * uvw[0];
1957 final int ifrx = (int) frx;
1958 final int ifrxu = ifrx + frad;
1959
1960 final double fry = fftY * uvw[1];
1961 final int ifry = (int) fry;
1962 final int ifryu = ifry + frad;
1963
1964 final double frz = fftZ * uvw[2];
1965 final int ifrz = (int) frz;
1966 final int ifrzu = ifrz + frad;
1967
1968 for (int iz = ifrz - frad; iz <= ifrzu; iz++) {
1969 int giz = mod(iz, fftZ);
1970 xf[2] = iz * ifftZ;
1971 for (int iy = ifry - frad; iy <= ifryu; iy++) {
1972 int giy = mod(iy, fftY);
1973 xf[1] = iy * ifftY;
1974 for (int ix = ifrx - frad; ix <= ifrxu; ix++) {
1975 int gix = mod(ix, fftX);
1976 xf[0] = ix * ifftX;
1977 crystal.toCartesianCoordinates(xf, xc);
1978 final int ii = interleavedIndex(gix, giy, giz, fftX, fftY);
1979 grid[ii] = solventff.rho(grid[ii], 1.0, xc);
1980 }
1981 }
1982 }
1983 }
1984 }
1985
1986 private class AtomicRowLoop extends RowLoop {
1987
1988 final double[] xyz = new double[3];
1989 final double[] uvw = new double[3];
1990 final double[] xc = new double[3];
1991 final double[] xf = new double[3];
1992 final double[] grid;
1993 final int[] optLocal;
1994 long timer;
1995 int previousUB, previousLB;
1996 int actualWeight;
1997
1998 AtomicRowLoop(RowRegion region) {
1999 super(region.getNatoms(), region.getNsymm(), region);
2000 grid = region.getGrid();
2001 optLocal = new int[fftZ * fftY];
2002 }
2003
2004 public void buildList(int iSymm, int iAtom, int lb, int ub) {
2005 if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2006 return;
2007 }
2008 xyz[0] = coordinates[iSymm][0][iAtom];
2009 xyz[1] = coordinates[iSymm][1][iAtom];
2010 xyz[2] = coordinates[iSymm][2][iAtom];
2011 crystal.toFractionalCoordinates(xyz, uvw);
2012 final int frad =
2013 min(aRadGrid, (int) floor(atoms[iAtom].getFormFactorWidth() * fftX / crystal.a) + 1);
2014
2015 final double frz = fftZ * uvw[2];
2016 final int ifrz = (int) frz;
2017 final int ifrzu = ifrz + frad;
2018 final int ifrzl = ifrz - frad;
2019
2020 final double fry = fftY * uvw[1];
2021 final int ifry = (int) fry;
2022 final int ifryu = ifry + frad;
2023 final int ifryl = ifry - frad;
2024
2025
2026
2027
2028 int buff = bufferSize;
2029
2030 int lbZ = rowRegion.zFromRowIndex(lb);
2031 int ubZ = rowRegion.zFromRowIndex(ub);
2032
2033 for (int iz = ifrzl - buff; iz <= ifrzu + buff; iz++) {
2034 int giz = mod(iz, fftZ);
2035 if (lbZ > giz || giz > ubZ) {
2036 continue;
2037 }
2038 int rowLB = rowRegion.rowIndexForYZ(mod(ifryl - buff, fftY), giz);
2039 int rowUB = rowRegion.rowIndexForYZ(mod(ifryu + buff, fftY), giz);
2040 if (lb >= rowLB || rowUB <= ub) {
2041 buildListA.add(iAtom);
2042 buildListS.add(iSymm);
2043 break;
2044 }
2045 }
2046 }
2047
2048 @Override
2049 public boolean checkList(int[][][] zyAtListBuild, int buff) {
2050 for (int iSymm = 0; iSymm < nSymm; iSymm++) {
2051 for (int iAtom = 0; iAtom < nAtoms; iAtom++) {
2052 if (rowRegion.select[iSymm][iAtom]) {
2053 if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2054 continue;
2055 }
2056 xyz[0] = coordinates[iSymm][0][iAtom];
2057 xyz[1] = coordinates[iSymm][1][iAtom];
2058 xyz[2] = coordinates[iSymm][2][iAtom];
2059 crystal.toFractionalCoordinates(xyz, uvw);
2060 final double frz = fftZ * uvw[2];
2061 final int ifrz = (int) frz;
2062 final int previousZ = zyAtListBuild[iSymm][iAtom][0];
2063
2064 final double fry = fftY * uvw[1];
2065 final int ifry = (int) fry;
2066 final int previousY = zyAtListBuild[iSymm][iAtom][1];
2067 if (abs(ifrz - previousZ) >= buff || abs(ifry - previousY) >= buff) {
2068 return true;
2069 }
2070 }
2071 }
2072 }
2073 return false;
2074 }
2075
2076 @Override
2077 public void finish() {
2078 for (int i = 0; i < fftZ * fftY; i++) {
2079 optRowWeightAtomic.addAndGet(i, optLocal[i]);
2080 }
2081 timer += System.nanoTime();
2082 }
2083
2084 @Override
2085 public void gridDensity(int iSymm, int iAtom, int lb, int ub) {
2086 if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2087 return;
2088 }
2089
2090 if (lambdaTerm && atoms[iAtom].applyLambda()) {
2091 return;
2092 }
2093
2094 int lbZ = rowRegion.zFromRowIndex(lb);
2095 int ubZ = rowRegion.zFromRowIndex(ub);
2096
2097 xyz[0] = coordinates[iSymm][0][iAtom];
2098 xyz[1] = coordinates[iSymm][1][iAtom];
2099 xyz[2] = coordinates[iSymm][2][iAtom];
2100 FormFactor atomff = atomFormFactors[iSymm][iAtom];
2101 crystal.toFractionalCoordinates(xyz, uvw);
2102 final int frad =
2103 min(aRadGrid, (int) floor(atoms[iAtom].getFormFactorWidth() * fftX / crystal.a) + 1);
2104
2105
2106 final double frx = fftX * uvw[0];
2107 final int ifrx = (int) frx;
2108 final int ifrxu = ifrx + frad;
2109
2110 final double fry = fftY * uvw[1];
2111 final int ifry = (int) fry;
2112 final int ifryu = ifry + frad;
2113
2114 final double frz = fftZ * uvw[2];
2115 final int ifrz = (int) frz;
2116 final int ifrzu = ifrz + frad;
2117
2118 for (int iz = ifrz - frad; iz <= ifrzu; iz++) {
2119 int giz = mod(iz, fftZ);
2120 if (lbZ > giz || giz > ubZ) {
2121 continue;
2122 }
2123 xf[2] = iz * ifftZ;
2124 for (int iy = ifry - frad; iy <= ifryu; iy++) {
2125 int giy = mod(iy, fftY);
2126 int rowIndex = rowRegion.rowIndexForYZ(giy, giz);
2127 if (lb > rowIndex || rowIndex > ub) {
2128 continue;
2129 }
2130 xf[1] = iy * ifftY;
2131 for (int ix = ifrx - frad; ix <= ifrxu; ix++) {
2132 int gix = mod(ix, fftX);
2133 xf[0] = ix * ifftX;
2134 crystal.toCartesianCoordinates(xf, xc);
2135 optLocal[rowIndex]++;
2136 actualWeight++;
2137 final int ii = interleavedIndex(gix, giy, giz, fftX, fftY);
2138 grid[ii] = atomff.rho(grid[ii], 1.0, xc);
2139 }
2140 }
2141 }
2142 }
2143
2144 @Override
2145 public void run(int lb, int ub) throws Exception {
2146 boolean boundsChange = false;
2147 if (previousLB != lb || previousUB != ub) {
2148 boundsChange = true;
2149 }
2150 previousLB = lb;
2151 previousUB = ub;
2152 if (rebuildList || boundsChange) {
2153 buildListA = new ArrayList<>();
2154 buildListS = new ArrayList<>();
2155 for (int iSymm = 0; iSymm < nSymm; iSymm++) {
2156 for (int iAtom = 0; iAtom < nAtoms; iAtom++) {
2157 if (rowRegion.select[iSymm][iAtom]) {
2158 buildList(iSymm, iAtom, lb, ub);
2159 }
2160 }
2161 }
2162 }
2163 for (int i = 0; i < buildListA.size(); i++) {
2164 if (rowRegion.select[buildListS.get(i)][buildListA.get(i)]) {
2165 gridDensity(buildListS.get(i), buildListA.get(i), lb, ub);
2166 }
2167 }
2168 }
2169
2170 @Override
2171 public void saveZYValues(int[][][] zyAtListBuild) {
2172 for (int iSymm = 0; iSymm < nSymm; iSymm++) {
2173 for (int iAtom = 0; iAtom < nAtoms; iAtom++) {
2174 zyAtListBuild[iSymm][iAtom][0] = 0;
2175 zyAtListBuild[iSymm][iAtom][1] = 0;
2176 if (rowRegion.select[iSymm][iAtom]) {
2177 if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2178 continue;
2179 }
2180 xyz[0] = coordinates[iSymm][0][iAtom];
2181 xyz[1] = coordinates[iSymm][1][iAtom];
2182 xyz[2] = coordinates[iSymm][2][iAtom];
2183 crystal.toFractionalCoordinates(xyz, uvw);
2184 final double frz = fftZ * uvw[2];
2185 final int ifrz = (int) frz;
2186
2187 final double fry = fftY * uvw[1];
2188 final int ifry = (int) fry;
2189 zyAtListBuild[iSymm][iAtom][0] = ifrz;
2190 zyAtListBuild[iSymm][iAtom][1] = ifry;
2191 }
2192 }
2193 }
2194 }
2195
2196 @Override
2197 public IntegerSchedule schedule() {
2198 return atomicRowSchedule;
2199 }
2200
2201 @Override
2202 public void start() {
2203 fill(optLocal, 0);
2204 actualWeight = 0;
2205 timer = -System.nanoTime();
2206 }
2207
2208 double getThreadTime() {
2209 return timer;
2210 }
2211
2212 int getNumberofSlices() {
2213 return (previousUB - previousLB + 1);
2214 }
2215
2216 int getThreadWeight() {
2217 return actualWeight;
2218 }
2219 }
2220
2221 private class SolventRowLoop extends RowLoop {
2222
2223 final double[] xyz = new double[3];
2224 final double[] uvw = new double[3];
2225 final double[] xc = new double[3];
2226 final double[] xf = new double[3];
2227 final double[] grid;
2228 final int[] optLocal;
2229 long timer;
2230 int previousUB, previousLB;
2231 int actualWeight;
2232
2233 public SolventRowLoop(RowRegion region) {
2234 super(region.getNatoms(), region.getNsymm(), region);
2235 grid = region.getGrid();
2236 optLocal = new int[fftZ * fftY];
2237 }
2238
2239 public void buildList(int iSymm, int iAtom, int lb, int ub) {
2240 if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2241 return;
2242 }
2243 xyz[0] = coordinates[iSymm][0][iAtom];
2244 xyz[1] = coordinates[iSymm][1][iAtom];
2245 xyz[2] = coordinates[iSymm][2][iAtom];
2246 crystal.toFractionalCoordinates(xyz, uvw);
2247 final int frad =
2248 min(aRadGrid, (int) floor(atoms[iAtom].getFormFactorWidth() * fftX / crystal.a) + 1);
2249
2250 final double frz = fftZ * uvw[2];
2251 final int ifrz = (int) frz;
2252 final int ifrzu = ifrz + frad;
2253 final int ifrzl = ifrz - frad;
2254
2255 final double fry = fftY * uvw[1];
2256 final int ifry = (int) fry;
2257 final int ifryu = ifry + frad;
2258 final int ifryl = ifry - frad;
2259
2260
2261
2262
2263 int buff = bufferSize;
2264
2265 int lbZ = rowRegion.zFromRowIndex(lb);
2266 int ubZ = rowRegion.zFromRowIndex(ub);
2267
2268 for (int iz = ifrzl - buff; iz <= ifrzu + buff; iz++) {
2269 int giz = mod(iz, fftZ);
2270 if (lbZ > giz || giz > ubZ) {
2271 continue;
2272 }
2273 int rowLB = rowRegion.rowIndexForYZ(mod(ifryl - buff, fftY), giz);
2274 int rowUB = rowRegion.rowIndexForYZ(mod(ifryu + buff, fftY), giz);
2275 if (lb >= rowLB || rowUB <= ub) {
2276 buildListA.add(iAtom);
2277 buildListS.add(iSymm);
2278 break;
2279 }
2280 }
2281 }
2282
2283 @Override
2284 public boolean checkList(int[][][] zyAtListBuild, int buff) {
2285 for (int iSymm = 0; iSymm < nSymm; iSymm++) {
2286 for (int iAtom = 0; iAtom < nAtoms; iAtom++) {
2287 if (rowRegion.select[iSymm][iAtom]) {
2288 if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2289 continue;
2290 }
2291 xyz[0] = coordinates[iSymm][0][iAtom];
2292 xyz[1] = coordinates[iSymm][1][iAtom];
2293 xyz[2] = coordinates[iSymm][2][iAtom];
2294 crystal.toFractionalCoordinates(xyz, uvw);
2295 final double frz = fftZ * uvw[2];
2296 final int ifrz = (int) frz;
2297 final int previousZ = zyAtListBuild[iSymm][iAtom][0];
2298
2299 final double fry = fftY * uvw[1];
2300 final int ifry = (int) fry;
2301 final int previousY = zyAtListBuild[iSymm][iAtom][1];
2302 if (abs(ifrz - previousZ) >= buff || abs(ifry - previousY) >= buff) {
2303 return true;
2304 }
2305 }
2306 }
2307 }
2308 return false;
2309 }
2310
2311 @Override
2312 public void finish() {
2313 timer += System.nanoTime();
2314 for (int i = 0; i < fftZ * fftY; i++) {
2315 optRowWeightSolvent.addAndGet(i, optLocal[i]);
2316 }
2317 }
2318
2319 public int getNumberofSlices() {
2320 return (previousUB - previousLB + 1);
2321 }
2322
2323 public double getThreadTime() {
2324 return timer;
2325 }
2326
2327 public int getThreadWeight() {
2328 return actualWeight;
2329 }
2330
2331 @Override
2332 public void gridDensity(int iSymm, int iAtom, int lb, int ub) {
2333
2334 if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2335 return;
2336 }
2337
2338 if (lambdaTerm && atoms[iAtom].applyLambda()) {
2339 return;
2340 }
2341
2342 xyz[0] = coordinates[iSymm][0][iAtom];
2343 xyz[1] = coordinates[iSymm][1][iAtom];
2344 xyz[2] = coordinates[iSymm][2][iAtom];
2345 FormFactor formFactor = solventFormFactors[iSymm][iAtom];
2346 crystal.toFractionalCoordinates(xyz, uvw);
2347 double vdwr = atoms[iAtom].getVDWType().radius * 0.5;
2348 int frad = aRadGrid;
2349 switch (solventModel) {
2350 case BINARY:
2351 frad = min(aRadGrid, (int) floor((vdwr + solventA + 0.2) * fftX / crystal.a) + 1);
2352 break;
2353 case GAUSSIAN:
2354 frad = min(aRadGrid, (int) floor((vdwr * solventB + 2.0) * fftX / crystal.a) + 1);
2355 break;
2356 case POLYNOMIAL:
2357 frad = min(aRadGrid, (int) floor((vdwr + solventB + 0.2) * fftX / crystal.a) + 1);
2358 break;
2359 case NONE:
2360 default:
2361 return;
2362 }
2363
2364 int ubZ = rowRegion.zFromRowIndex(ub);
2365 int lbZ = rowRegion.zFromRowIndex(lb);
2366
2367
2368 final double frx = fftX * uvw[0];
2369 final int ifrx = (int) frx;
2370 final int ifrxu = ifrx + frad;
2371
2372 final double fry = fftY * uvw[1];
2373 final int ifry = (int) fry;
2374 final int ifryu = ifry + frad;
2375
2376 final double frz = fftZ * uvw[2];
2377 final int ifrz = (int) frz;
2378 final int ifrzu = ifrz + frad;
2379
2380 for (int iz = ifrz - frad; iz <= ifrzu; iz++) {
2381 int giz = mod(iz, fftZ);
2382 if (lbZ > giz || giz > ubZ) {
2383 continue;
2384 }
2385 xf[2] = iz * ifftZ;
2386 for (int iy = ifry - frad; iy <= ifryu; iy++) {
2387 int giy = mod(iy, fftY);
2388 int rowIndex = rowRegion.rowIndexForYZ(giy, giz);
2389 if (lb > rowIndex || rowIndex > ub) {
2390 continue;
2391 }
2392 xf[1] = iy * ifftY;
2393 for (int ix = ifrx - frad; ix <= ifrxu; ix++) {
2394 int gix = mod(ix, fftX);
2395 xf[0] = ix * ifftX;
2396 crystal.toCartesianCoordinates(xf, xc);
2397 optLocal[rowIndex]++;
2398 actualWeight++;
2399 final int ii = interleavedIndex(gix, giy, giz, fftX, fftY);
2400 grid[ii] = formFactor.rho(grid[ii], 1.0, xc);
2401 }
2402 }
2403 }
2404 }
2405
2406 @Override
2407 public void run(int lb, int ub) throws Exception {
2408 boolean boundsChange = false;
2409 if (previousLB != lb || previousUB != ub) {
2410 boundsChange = true;
2411 }
2412 previousLB = lb;
2413 previousUB = ub;
2414 if (rebuildList || boundsChange) {
2415 buildListA = new ArrayList<>();
2416 buildListS = new ArrayList<>();
2417 for (int iSymm = 0; iSymm < nSymm; iSymm++) {
2418 for (int iAtom = 0; iAtom < nAtoms; iAtom++) {
2419 if (rowRegion.select[iSymm][iAtom]) {
2420 buildList(iSymm, iAtom, lb, ub);
2421 }
2422 }
2423 }
2424 }
2425 for (int i = 0; i < buildListA.size(); i++) {
2426 if (rowRegion.select[buildListS.get(i)][buildListA.get(i)]) {
2427 gridDensity(buildListS.get(i), buildListA.get(i), lb, ub);
2428 }
2429 }
2430 }
2431
2432 @Override
2433 public void saveZYValues(int[][][] zyAtListBuild) {
2434 for (int iSymm = 0; iSymm < nSymm; iSymm++) {
2435 for (int iAtom = 0; iAtom < nAtoms; iAtom++) {
2436 zyAtListBuild[iSymm][iAtom][0] = 0;
2437 zyAtListBuild[iSymm][iAtom][1] = 0;
2438 if (rowRegion.select[iSymm][iAtom]) {
2439 if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2440 continue;
2441 }
2442 xyz[0] = coordinates[iSymm][0][iAtom];
2443 xyz[1] = coordinates[iSymm][1][iAtom];
2444 xyz[2] = coordinates[iSymm][2][iAtom];
2445 crystal.toFractionalCoordinates(xyz, uvw);
2446 final double frz = fftZ * uvw[2];
2447 final int ifrz = (int) frz;
2448
2449 final double fry = fftY * uvw[1];
2450 final int ifry = (int) fry;
2451 zyAtListBuild[iSymm][iAtom][0] = ifrz;
2452 zyAtListBuild[iSymm][iAtom][1] = ifry;
2453 }
2454 }
2455 }
2456 }
2457
2458 @Override
2459 public IntegerSchedule schedule() {
2460 return solventRowSchedule;
2461 }
2462
2463 @Override
2464 public void start() {
2465 fill(optLocal, 0);
2466 timer = -System.nanoTime();
2467 }
2468 }
2469
2470 private class BulkSolventRowLoop extends RowLoop {
2471
2472 final double[] xyz = new double[3];
2473 final double[] uvw = new double[3];
2474 final double[] xc = new double[3];
2475 final double[] xf = new double[3];
2476 final double[] grid;
2477 final int[] optLocal;
2478 long timer;
2479 int[] threadWeights;
2480 int[] threadBounds;
2481
2482 public BulkSolventRowLoop(RowRegion region) {
2483 super(region.getNatoms(), region.getNsymm(), region);
2484 grid = region.getGrid();
2485 optLocal = new int[fftZ * fftY];
2486 }
2487
2488 @Override
2489 public void finish() {
2490 for (int i = 0; i < fftZ * fftY; i++) {
2491 optRowWeightBulkSolvent.addAndGet(i, optLocal[i]);
2492 }
2493 timer += System.nanoTime();
2494 threadBounds = bulkSolventRowSchedule.getLowerBounds().clone();
2495 threadWeights = bulkSolventRowSchedule.getThreadWeights().clone();
2496 for (int i = threadCount - 1; i > 0; i--) {
2497 threadWeights[i] -= threadWeights[i - 1];
2498 }
2499 }
2500
2501 public int[] getBoundsPerThread() {
2502 return threadBounds;
2503 }
2504
2505 public double getTimePerThread() {
2506 return timer;
2507 }
2508
2509 public int[] getWeightPerThread() {
2510 return optLocal;
2511 }
2512
2513 @Override
2514 public void gridDensity(int iSymm, int iAtom, int lb, int ub) {
2515 if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2516 return;
2517 }
2518
2519 if (lambdaTerm && atoms[iAtom].applyLambda()) {
2520 return;
2521 }
2522
2523 xyz[0] = coordinates[iSymm][0][iAtom];
2524 xyz[1] = coordinates[iSymm][1][iAtom];
2525 xyz[2] = coordinates[iSymm][2][iAtom];
2526 FormFactor formFactor = solventFormFactors[iSymm][iAtom];
2527 crystal.toFractionalCoordinates(xyz, uvw);
2528 double vdwr = atoms[iAtom].getVDWType().radius * 0.5;
2529 int frad = aRadGrid;
2530 switch (solventModel) {
2531 case BINARY:
2532 frad = min(aRadGrid, (int) floor((vdwr + solventA + 0.2) * fftX / crystal.a) + 1);
2533 break;
2534 case GAUSSIAN:
2535 frad = min(aRadGrid, (int) floor((vdwr * solventB + 2.0) * fftX / crystal.a) + 1);
2536 break;
2537 case POLYNOMIAL:
2538 frad = min(aRadGrid, (int) floor((vdwr + solventB + 0.2) * fftX / crystal.a) + 1);
2539 break;
2540 case NONE:
2541 default:
2542 return;
2543 }
2544
2545 int ubZ = rowRegion.zFromRowIndex(ub);
2546 int lbZ = rowRegion.zFromRowIndex(lb);
2547
2548
2549 final double frx = fftX * uvw[0];
2550 final int ifrx = (int) frx;
2551 final int ifrxu = ifrx + frad;
2552
2553 final double fry = fftY * uvw[1];
2554 final int ifry = (int) fry;
2555 final int ifryu = ifry + frad;
2556
2557 final double frz = fftZ * uvw[2];
2558 final int ifrz = (int) frz;
2559 final int ifrzu = ifrz + frad;
2560
2561 for (int iz = ifrz - frad; iz <= ifrzu; iz++) {
2562 int giz = mod(iz, fftZ);
2563 if (lbZ > giz || giz > ubZ) {
2564 continue;
2565 }
2566 xf[2] = iz * ifftZ;
2567 for (int iy = ifry - frad; iy <= ifryu; iy++) {
2568 int giy = mod(iy, fftY);
2569 int rowIndex = rowRegion.rowIndexForYZ(giy, giz);
2570 if (lb > rowIndex || rowIndex > ub) {
2571 continue;
2572 }
2573 xf[1] = iy * ifftY;
2574 for (int ix = ifrx - frad; ix <= ifrxu; ix++) {
2575 int gix = mod(ix, fftX);
2576 xf[0] = ix * ifftX;
2577 crystal.toCartesianCoordinates(xf, xc);
2578 optLocal[rowIndex]++;
2579 final int ii = interleavedIndex(gix, giy, giz, fftX, fftY);
2580 grid[ii] = formFactor.rho(grid[ii], 1.0, xc);
2581 }
2582 }
2583 }
2584 }
2585
2586 @Override
2587 public IntegerSchedule schedule() {
2588 return bulkSolventRowSchedule;
2589 }
2590
2591 @Override
2592 public void start() {
2593 timer = -System.nanoTime();
2594 fill(optLocal, 0);
2595 }
2596 }
2597
2598 private class AtomicSliceLoop extends SliceLoop {
2599
2600 final double[] xyz = new double[3];
2601 final double[] uvw = new double[3];
2602 final double[] xc = new double[3];
2603 final double[] xf = new double[3];
2604 final double[] grid;
2605 final int[] optLocal;
2606 long timer;
2607 int previousUB, previousLB;
2608 int actualWeight;
2609
2610 public AtomicSliceLoop(SliceRegion region) {
2611 super(region.getNatoms(), region.getNsymm(), region);
2612 grid = region.getGrid();
2613 optLocal = new int[fftZ];
2614 }
2615
2616 public void buildList(int iSymm, int iAtom, int lb, int ub) {
2617 if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2618 return;
2619 }
2620 xyz[0] = coordinates[iSymm][0][iAtom];
2621 xyz[1] = coordinates[iSymm][1][iAtom];
2622 xyz[2] = coordinates[iSymm][2][iAtom];
2623 crystal.toFractionalCoordinates(xyz, uvw);
2624 final int frad =
2625 min(aRadGrid, (int) floor(atoms[iAtom].getFormFactorWidth() * fftX / crystal.a) + 1);
2626 final double frz = fftZ * uvw[2];
2627 final int ifrz = (int) frz;
2628 final int ifrzu = ifrz + frad;
2629 final int ifrzl = ifrz - frad;
2630
2631
2632
2633 int buff = bufferSize;
2634 for (int iz = ifrzl - buff; iz <= ifrzu + buff; iz++) {
2635 int giz = mod(iz, fftZ);
2636 if (giz >= lb && giz <= ub) {
2637 buildListA.add(iAtom);
2638 buildListS.add(iSymm);
2639 break;
2640 }
2641 }
2642 }
2643
2644 @Override
2645 public boolean checkList(int[][] zAtListBuild, int buff) {
2646 for (int iSymm = 0; iSymm < nSymm; iSymm++) {
2647 for (int iAtom = 0; iAtom < nAtoms; iAtom++) {
2648 if (sliceRegion.select[iSymm][iAtom]) {
2649 if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2650 continue;
2651 }
2652 xyz[0] = coordinates[iSymm][0][iAtom];
2653 xyz[1] = coordinates[iSymm][1][iAtom];
2654 xyz[2] = coordinates[iSymm][2][iAtom];
2655 crystal.toFractionalCoordinates(xyz, uvw);
2656 final double frz = fftZ * uvw[2];
2657 final int ifrz = (int) frz;
2658 final int previousZ = zAtListBuild[iSymm][iAtom];
2659 if (abs(ifrz - previousZ) >= buff) {
2660 return true;
2661 }
2662 }
2663 }
2664 }
2665 return false;
2666 }
2667
2668 @Override
2669 public void finish() {
2670 for (int i = 0; i < fftZ; i++) {
2671 optSliceWeightAtomic.addAndGet(i, optLocal[i]);
2672 }
2673 timer += System.nanoTime();
2674 }
2675
2676 public int getNumberofSlices() {
2677 return (previousUB - previousLB + 1);
2678 }
2679
2680 public double getThreadTime() {
2681 return timer;
2682 }
2683
2684 public int getThreadWeight() {
2685 return actualWeight;
2686 }
2687
2688 @Override
2689 public void gridDensity(int iSymm, int iAtom, int lb, int ub) {
2690 if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2691 return;
2692 }
2693
2694 if (lambdaTerm && atoms[iAtom].applyLambda()) {
2695 return;
2696 }
2697
2698 xyz[0] = coordinates[iSymm][0][iAtom];
2699 xyz[1] = coordinates[iSymm][1][iAtom];
2700 xyz[2] = coordinates[iSymm][2][iAtom];
2701 FormFactor atomff = atomFormFactors[iSymm][iAtom];
2702 crystal.toFractionalCoordinates(xyz, uvw);
2703 final int frad =
2704 min(aRadGrid, (int) floor(atoms[iAtom].getFormFactorWidth() * fftX / crystal.a) + 1);
2705
2706
2707 final double frx = fftX * uvw[0];
2708 final int ifrx = (int) frx;
2709 final int ifrxu = ifrx + frad;
2710
2711 final double fry = fftY * uvw[1];
2712 final int ifry = (int) fry;
2713 final int ifryu = ifry + frad;
2714
2715 final double frz = fftZ * uvw[2];
2716 final int ifrz = (int) frz;
2717 final int ifrzu = ifrz + frad;
2718
2719 for (int iz = ifrz - frad; iz <= ifrzu; iz++) {
2720 int giz = mod(iz, fftZ);
2721 if (lb > giz || giz > ub) {
2722 continue;
2723 }
2724 xf[2] = iz * ifftZ;
2725 for (int iy = ifry - frad; iy <= ifryu; iy++) {
2726 int giy = mod(iy, fftY);
2727 xf[1] = iy * ifftY;
2728 for (int ix = ifrx - frad; ix <= ifrxu; ix++) {
2729 int gix = mod(ix, fftX);
2730 xf[0] = ix * ifftX;
2731 crystal.toCartesianCoordinates(xf, xc);
2732 optLocal[giz]++;
2733 actualWeight++;
2734 final int ii = interleavedIndex(gix, giy, giz, fftX, fftY);
2735 grid[ii] = atomff.rho(grid[ii], 1.0, xc);
2736 }
2737 }
2738 }
2739 }
2740
2741 @Override
2742 public void run(int lb, int ub) throws Exception {
2743 boolean boundsChange = false;
2744 if (previousLB != lb || previousUB != ub) {
2745 boundsChange = true;
2746 }
2747 previousLB = lb;
2748 previousUB = ub;
2749 if (rebuildList || boundsChange) {
2750 buildListA = new ArrayList<>();
2751 buildListS = new ArrayList<>();
2752 for (int iSymm = 0; iSymm < nSymm; iSymm++) {
2753 for (int iAtom = 0; iAtom < nAtoms; iAtom++) {
2754 if (sliceRegion.select[iSymm][iAtom]) {
2755 buildList(iSymm, iAtom, lb, ub);
2756 }
2757 }
2758 }
2759 }
2760 for (int i = 0; i < buildListA.size(); i++) {
2761 if (sliceRegion.select[buildListS.get(i)][buildListA.get(i)]) {
2762 gridDensity(buildListS.get(i), buildListA.get(i), lb, ub);
2763 }
2764 }
2765 }
2766
2767 @Override
2768 public void saveZValues(int[][] zAtListBuild) {
2769 for (int iSymm = 0; iSymm < nSymm; iSymm++) {
2770 for (int iAtom = 0; iAtom < nAtoms; iAtom++) {
2771 zAtListBuild[iSymm][iAtom] = 0;
2772 if (sliceRegion.select[iSymm][iAtom]) {
2773 if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2774 continue;
2775 }
2776 xyz[0] = coordinates[iSymm][0][iAtom];
2777 xyz[1] = coordinates[iSymm][1][iAtom];
2778 xyz[2] = coordinates[iSymm][2][iAtom];
2779 crystal.toFractionalCoordinates(xyz, uvw);
2780 final double frz = fftZ * uvw[2];
2781 final int ifrz = (int) frz;
2782 zAtListBuild[iSymm][iAtom] = ifrz;
2783 }
2784 }
2785 }
2786 }
2787
2788 @Override
2789 public IntegerSchedule schedule() {
2790 return atomicSliceSchedule;
2791 }
2792
2793 @Override
2794 public void start() {
2795 fill(optLocal, 0);
2796 actualWeight = 0;
2797 timer = -System.nanoTime();
2798 }
2799 }
2800
2801 private class SolventSliceLoop extends SliceLoop {
2802
2803 final double[] xyz = new double[3];
2804 final double[] uvw = new double[3];
2805 final double[] xc = new double[3];
2806 final double[] xf = new double[3];
2807 final double[] grid;
2808 final int[] optLocal;
2809 long timer;
2810 int previousUB, previousLB;
2811 int actualWeight;
2812
2813 public SolventSliceLoop(SliceRegion region) {
2814 super(region.getNatoms(), region.getNsymm(), region);
2815 grid = region.getGrid();
2816 optLocal = new int[fftZ];
2817 }
2818
2819 public void buildList(int iSymm, int iAtom, int lb, int ub) {
2820 if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2821 return;
2822 }
2823 xyz[0] = coordinates[iSymm][0][iAtom];
2824 xyz[1] = coordinates[iSymm][1][iAtom];
2825 xyz[2] = coordinates[iSymm][2][iAtom];
2826 crystal.toFractionalCoordinates(xyz, uvw);
2827 final int frad =
2828 min(aRadGrid, (int) floor(atoms[iAtom].getFormFactorWidth() * fftX / crystal.a) + 1);
2829 final double frz = fftZ * uvw[2];
2830 final int ifrz = (int) frz;
2831 final int ifrzu = ifrz + frad;
2832 final int ifrzl = ifrz - frad;
2833
2834
2835
2836 int buff = bufferSize;
2837 for (int iz = ifrzl - buff; iz <= ifrzu + buff; iz++) {
2838 int giz = mod(iz, fftZ);
2839 if (giz >= lb && giz <= ub) {
2840 buildListA.add(iAtom);
2841 buildListS.add(iSymm);
2842 break;
2843 }
2844 }
2845 }
2846
2847 @Override
2848 public boolean checkList(int[][] zAtListBuild, int buff) {
2849 for (int iSymm = 0; iSymm < nSymm; iSymm++) {
2850 for (int iAtom = 0; iAtom < nAtoms; iAtom++) {
2851 if (sliceRegion.select[iSymm][iAtom]) {
2852 if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2853 continue;
2854 }
2855 xyz[0] = coordinates[iSymm][0][iAtom];
2856 xyz[1] = coordinates[iSymm][1][iAtom];
2857 xyz[2] = coordinates[iSymm][2][iAtom];
2858 crystal.toFractionalCoordinates(xyz, uvw);
2859 final double frz = fftZ * uvw[2];
2860 final int ifrz = (int) frz;
2861 final int previousZ = zAtListBuild[iSymm][iAtom];
2862 if (abs(ifrz - previousZ) >= buff) {
2863 return true;
2864 }
2865 }
2866 }
2867 }
2868 return false;
2869 }
2870
2871 @Override
2872 public void finish() {
2873 timer += System.nanoTime();
2874 for (int i = 0; i < fftZ; i++) {
2875 optSliceWeightSolvent.addAndGet(i, optLocal[i]);
2876 }
2877 }
2878
2879 public int getNumberofSlices() {
2880 return (previousUB - previousLB + 1);
2881 }
2882
2883 public double getThreadTime() {
2884 return timer;
2885 }
2886
2887 public int getThreadWeight() {
2888 return actualWeight;
2889 }
2890
2891 @Override
2892 public void gridDensity(int iSymm, int iAtom, int lb, int ub) {
2893 if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2894 return;
2895 }
2896
2897 if (lambdaTerm && atoms[iAtom].applyLambda()) {
2898 return;
2899 }
2900
2901 xyz[0] = coordinates[iSymm][0][iAtom];
2902 xyz[1] = coordinates[iSymm][1][iAtom];
2903 xyz[2] = coordinates[iSymm][2][iAtom];
2904 FormFactor formFactor = solventFormFactors[iSymm][iAtom];
2905 crystal.toFractionalCoordinates(xyz, uvw);
2906 double vdwr = atoms[iAtom].getVDWType().radius * 0.5;
2907 int frad = aRadGrid;
2908 switch (solventModel) {
2909 case BINARY:
2910 frad = min(aRadGrid, (int) floor((vdwr + solventA + 0.2) * fftX / crystal.a) + 1);
2911 break;
2912 case GAUSSIAN:
2913 frad = min(aRadGrid, (int) floor((vdwr * solventB + 2.0) * fftX / crystal.a) + 1);
2914 break;
2915 case POLYNOMIAL:
2916 frad = min(aRadGrid, (int) floor((vdwr + solventB + 0.2) * fftX / crystal.a) + 1);
2917 break;
2918 case NONE:
2919 default:
2920 return;
2921 }
2922
2923
2924 final double frx = fftX * uvw[0];
2925 final int ifrx = (int) frx;
2926 final int ifrxu = ifrx + frad;
2927
2928 final double fry = fftY * uvw[1];
2929 final int ifry = (int) fry;
2930 final int ifryu = ifry + frad;
2931
2932 final double frz = fftZ * uvw[2];
2933 final int ifrz = (int) frz;
2934 final int ifrzu = ifrz + frad;
2935
2936 for (int iz = ifrz - frad; iz <= ifrzu; iz++) {
2937 int giz = mod(iz, fftZ);
2938 if (lb > giz || giz > ub) {
2939 continue;
2940 }
2941 xf[2] = iz * ifftZ;
2942 for (int iy = ifry - frad; iy <= ifryu; iy++) {
2943 int giy = mod(iy, fftY);
2944 xf[1] = iy * ifftY;
2945 for (int ix = ifrx - frad; ix <= ifrxu; ix++) {
2946 int gix = mod(ix, fftX);
2947 xf[0] = ix * ifftX;
2948 crystal.toCartesianCoordinates(xf, xc);
2949 optLocal[giz]++;
2950 actualWeight++;
2951 final int ii = interleavedIndex(gix, giy, giz, fftX, fftY);
2952 grid[ii] = formFactor.rho(grid[ii], 1.0, xc);
2953 }
2954 }
2955 }
2956 }
2957
2958 @Override
2959 public void run(int lb, int ub) throws Exception {
2960 boolean boundsChange = false;
2961 if (previousLB != lb || previousUB != ub) {
2962 boundsChange = true;
2963 }
2964 previousLB = lb;
2965 previousUB = ub;
2966 if (rebuildList || boundsChange) {
2967 buildListA = new ArrayList<>();
2968 buildListS = new ArrayList<>();
2969 for (int iSymm = 0; iSymm < nSymm; iSymm++) {
2970 for (int iAtom = 0; iAtom < nAtoms; iAtom++) {
2971 if (sliceRegion.select[iSymm][iAtom]) {
2972 buildList(iSymm, iAtom, lb, ub);
2973 }
2974 }
2975 }
2976 }
2977 for (int i = 0; i < buildListA.size(); i++) {
2978 if (sliceRegion.select[buildListS.get(i)][buildListA.get(i)]) {
2979 gridDensity(buildListS.get(i), buildListA.get(i), lb, ub);
2980 }
2981 }
2982 }
2983
2984 @Override
2985 public void saveZValues(int[][] zAtListBuild) {
2986 for (int iSymm = 0; iSymm < nSymm; iSymm++) {
2987 for (int iAtom = 0; iAtom < nAtoms; iAtom++) {
2988 zAtListBuild[iSymm][iAtom] = 0;
2989 if (sliceRegion.select[iSymm][iAtom]) {
2990 if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2991 continue;
2992 }
2993 xyz[0] = coordinates[iSymm][0][iAtom];
2994 xyz[1] = coordinates[iSymm][1][iAtom];
2995 xyz[2] = coordinates[iSymm][2][iAtom];
2996 crystal.toFractionalCoordinates(xyz, uvw);
2997 final double frz = fftZ * uvw[2];
2998 final int ifrz = (int) frz;
2999 zAtListBuild[iSymm][iAtom] = ifrz;
3000 }
3001 }
3002 }
3003 }
3004
3005 @Override
3006 public IntegerSchedule schedule() {
3007 return solventSliceSchedule;
3008 }
3009
3010 @Override
3011 public void start() {
3012 fill(optLocal, 0);
3013 actualWeight = 0;
3014 timer = -System.nanoTime();
3015 }
3016 }
3017
3018 private class BulkSolventSliceLoop extends SliceLoop {
3019
3020 final double[] xyz = new double[3];
3021 final double[] uvw = new double[3];
3022 final double[] xc = new double[3];
3023 final double[] xf = new double[3];
3024 final double[] grid;
3025 final int[] optLocal;
3026 long timer;
3027 int[] threadBounds;
3028 int[] threadWeights;
3029
3030 public BulkSolventSliceLoop(SliceRegion region) {
3031 super(region.getNatoms(), region.getNsymm(), region);
3032 grid = region.getGrid();
3033 optLocal = new int[fftZ];
3034 }
3035
3036 @Override
3037 public void finish() {
3038 timer += System.nanoTime();
3039 for (int i = 0; i < fftZ; i++) {
3040 optSliceWeightBulkSolvent.addAndGet(i, optLocal[i]);
3041 }
3042 threadBounds = bulkSolventSliceSchedule.getLowerBounds().clone();
3043 threadWeights = bulkSolventSliceSchedule.getThreadWeights().clone();
3044 for (int i = threadBounds.length - 1; i > 0; i--) {
3045 threadBounds[i] -= threadBounds[i - 1];
3046 }
3047 }
3048
3049 public int[] getBoundsPerThread() {
3050 return threadBounds;
3051 }
3052
3053 public double getTimePerThread() {
3054 return timer;
3055 }
3056
3057 public int[] getWeightPerThread() {
3058 return threadWeights;
3059 }
3060
3061 @Override
3062 public void gridDensity(int iSymm, int iAtom, int lb, int ub) {
3063 if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
3064 return;
3065 }
3066
3067 if (lambdaTerm && atoms[iAtom].applyLambda()) {
3068 return;
3069 }
3070
3071 xyz[0] = coordinates[iSymm][0][iAtom];
3072 xyz[1] = coordinates[iSymm][1][iAtom];
3073 xyz[2] = coordinates[iSymm][2][iAtom];
3074 FormFactor formFactor = solventFormFactors[iSymm][iAtom];
3075 crystal.toFractionalCoordinates(xyz, uvw);
3076 double vdwr = atoms[iAtom].getVDWType().radius * 0.5;
3077 int frad = aRadGrid;
3078 switch (solventModel) {
3079 case BINARY:
3080 frad = min(aRadGrid, (int) floor((vdwr + solventA + 0.2) * fftX / crystal.a) + 1);
3081 break;
3082 case GAUSSIAN:
3083 frad = min(aRadGrid, (int) floor((vdwr * solventB + 2.0) * fftX / crystal.a) + 1);
3084 break;
3085 case POLYNOMIAL:
3086 frad = min(aRadGrid, (int) floor((vdwr + solventB + 0.2) * fftX / crystal.a) + 1);
3087 break;
3088 case NONE:
3089 default:
3090 return;
3091 }
3092
3093
3094 final double frx = fftX * uvw[0];
3095 final int ifrx = (int) frx;
3096 final int ifrxu = ifrx + frad;
3097
3098 final double fry = fftY * uvw[1];
3099 final int ifry = (int) fry;
3100 final int ifryu = ifry + frad;
3101
3102 final double frz = fftZ * uvw[2];
3103 final int ifrz = (int) frz;
3104 final int ifrzu = ifrz + frad;
3105
3106 for (int iz = ifrz - frad; iz <= ifrzu; iz++) {
3107 int giz = mod(iz, fftZ);
3108 if (lb > giz || giz > ub) {
3109 continue;
3110 }
3111 xf[2] = iz * ifftZ;
3112 for (int iy = ifry - frad; iy <= ifryu; iy++) {
3113 int giy = mod(iy, fftY);
3114 xf[1] = iy * ifftY;
3115 for (int ix = ifrx - frad; ix <= ifrxu; ix++) {
3116 int gix = mod(ix, fftX);
3117 xf[0] = ix * ifftX;
3118 crystal.toCartesianCoordinates(xf, xc);
3119 optLocal[giz]++;
3120 final int ii = interleavedIndex(gix, giy, giz, fftX, fftY);
3121 grid[ii] = formFactor.rho(grid[ii], 1.0, xc);
3122 }
3123 }
3124 }
3125 }
3126
3127 @Override
3128 public IntegerSchedule schedule() {
3129 return bulkSolventSliceSchedule;
3130 }
3131
3132 @Override
3133 public void start() {
3134 fill(optLocal, 0);
3135 timer = -System.nanoTime();
3136 }
3137 }
3138
3139 private class AtomicGradientRegion extends ParallelRegion {
3140
3141 private final AtomicGradientLoop[] atomicGradientLoop;
3142 private RefinementMode refinementmode;
3143
3144 public AtomicGradientRegion(RefinementMode refinementmode) {
3145 this.refinementmode = refinementmode;
3146 atomicGradientLoop = new AtomicGradientLoop[threadCount];
3147 for (int i = 0; i < threadCount; i++) {
3148 atomicGradientLoop[i] = new AtomicGradientLoop();
3149 }
3150 }
3151
3152 public RefinementMode getRefinementMode() {
3153 return refinementmode;
3154 }
3155
3156 public void setRefinementMode(RefinementMode refinementmode) {
3157 this.refinementmode = refinementmode;
3158 }
3159
3160 @Override
3161 public void run() {
3162 try {
3163 execute(0, nAtoms - 1, atomicGradientLoop[getThreadIndex()]);
3164 } catch (Exception e) {
3165 logger.severe(e.toString());
3166 }
3167 }
3168
3169 private class AtomicGradientLoop extends IntegerForLoop {
3170
3171 final double[] xyz = new double[3];
3172 final double[] uvw = new double[3];
3173 final double[] xc = new double[3];
3174 final double[] xf = new double[3];
3175 final int[] optLocal = new int[nAtoms];
3176 long timer;
3177
3178 @Override
3179 public void finish() {
3180 timer += System.nanoTime();
3181 for (int i = 0; i < nAtoms; i++) {
3182 optAtomicGradientWeight.addAndGet(i, optLocal[i]);
3183 }
3184 }
3185
3186 @Override
3187 public void run(final int lb, final int ub) {
3188 for (int n = lb; n <= ub; n++) {
3189 if (!atoms[n].getUse() && !nativeEnvironmentApproximation) {
3190 continue;
3191 }
3192 if (lambdaTerm && atoms[n].applyLambda()) {
3193 continue;
3194 }
3195
3196 xyz[0] = coordinates[0][0][n];
3197 xyz[1] = coordinates[0][1][n];
3198 xyz[2] = coordinates[0][2][n];
3199 FormFactor atomff = atomFormFactors[0][n];
3200 atomff.update(xyz, 0.0);
3201 crystal.toFractionalCoordinates(xyz, uvw);
3202 final int dfrad =
3203 Math.min(
3204 aRadGrid, (int) Math.floor(atoms[n].getFormFactorWidth() * fftX / crystal.a) + 1);
3205
3206
3207 final double frx = fftX * uvw[0];
3208 final int ifrx = (int) frx;
3209
3210 final double fry = fftY * uvw[1];
3211 final int ifry = (int) fry;
3212
3213 final double frz = fftZ * uvw[2];
3214 final int ifrz = (int) frz;
3215
3216 for (int ix = ifrx - dfrad; ix <= ifrx + dfrad; ix++) {
3217 int gix = mod(ix, fftX);
3218 xf[0] = ix * ifftX;
3219 for (int iy = ifry - dfrad; iy <= ifry + dfrad; iy++) {
3220 int giy = mod(iy, fftY);
3221 xf[1] = iy * ifftY;
3222 for (int iz = ifrz - dfrad; iz <= ifrz + dfrad; iz++) {
3223 int giz = mod(iz, fftZ);
3224 xf[2] = iz * ifftZ;
3225 crystal.toCartesianCoordinates(xf, xc);
3226 optLocal[n]++;
3227 final int ii = interleavedIndex(gix, giy, giz, fftX, fftY);
3228 atomff.rhoGrad(xc, weight * getDensityGrid()[ii], refinementmode);
3229 }
3230 }
3231 }
3232 }
3233 }
3234
3235 @Override
3236 public IntegerSchedule schedule() {
3237 return atomicGradientSchedule;
3238 }
3239
3240 @Override
3241 public void start() {
3242 fill(optLocal, 0);
3243 timer = -System.nanoTime();
3244 }
3245 }
3246 }
3247
3248 private class SolventGradientRegion extends ParallelRegion {
3249
3250 private final SolventGradientLoop[] solventGradientLoop;
3251 private RefinementMode refinementmode;
3252
3253 public SolventGradientRegion(RefinementMode refinementmode) {
3254 this.refinementmode = refinementmode;
3255 solventGradientLoop = new SolventGradientLoop[threadCount];
3256 for (int i = 0; i < threadCount; i++) {
3257 solventGradientLoop[i] = new SolventGradientLoop();
3258 }
3259 }
3260
3261 public RefinementMode getRefinementMode() {
3262 return refinementmode;
3263 }
3264
3265 public void setRefinementMode(RefinementMode refinementmode) {
3266 this.refinementmode = refinementmode;
3267 }
3268
3269 @Override
3270 public void run() {
3271 try {
3272 execute(0, nAtoms - 1, solventGradientLoop[getThreadIndex()]);
3273 } catch (Exception e) {
3274 logger.severe(e.toString());
3275 }
3276 }
3277
3278 private class SolventGradientLoop extends IntegerForLoop {
3279
3280 double[] xyz = new double[3];
3281 double[] uvw = new double[3];
3282 double[] xc = new double[3];
3283 double[] xf = new double[3];
3284
3285 @Override
3286 public void run(final int lb, final int ub) {
3287 for (int n = lb; n <= ub; n++) {
3288 if (!atoms[n].getUse() && !nativeEnvironmentApproximation) {
3289 continue;
3290 }
3291 if (lambdaTerm && atoms[n].applyLambda()) {
3292 continue;
3293 }
3294 xyz[0] = coordinates[0][0][n];
3295 xyz[1] = coordinates[0][1][n];
3296 xyz[2] = coordinates[0][2][n];
3297 FormFactor solventff = solventFormFactors[0][n];
3298 solventff.update(xyz);
3299 crystal.toFractionalCoordinates(xyz, uvw);
3300 double vdwr = atoms[n].getVDWType().radius * 0.5;
3301 double dfcmult = 1.0;
3302 int dfrad = aRadGrid;
3303 switch (solventModel) {
3304 case BINARY:
3305 dfrad =
3306 Math.min(
3307 aRadGrid, (int) Math.floor((vdwr + solventA + 0.2) * fftX / crystal.a) + 1);
3308 break;
3309 case GAUSSIAN:
3310 dfrad =
3311 Math.min(
3312 aRadGrid, (int) Math.floor((vdwr * solventB + 2.0) * fftX / crystal.a) + 1);
3313 dfcmult = solventA;
3314 break;
3315 case POLYNOMIAL:
3316 dfrad =
3317 Math.min(
3318 aRadGrid, (int) Math.floor((vdwr + solventB + 0.2) * fftX / crystal.a) + 1);
3319 break;
3320 case NONE:
3321 return;
3322 }
3323
3324
3325 final double frx = fftX * uvw[0];
3326 final int ifrx = (int) frx;
3327 final int ifrxu = ifrx + dfrad;
3328
3329 final double fry = fftY * uvw[1];
3330 final int ifry = (int) fry;
3331 final int ifryu = ifry + dfrad;
3332
3333 final double frz = fftZ * uvw[2];
3334 final int ifrz = (int) frz;
3335 final int ifrzu = ifrz + dfrad;
3336
3337 for (int ix = ifrx - dfrad; ix <= ifrxu; ix++) {
3338 int gix = mod(ix, fftX);
3339 xf[0] = ix * ifftX;
3340 for (int iy = ifry - dfrad; iy <= ifryu; iy++) {
3341 int giy = mod(iy, fftY);
3342 xf[1] = iy * ifftY;
3343 for (int iz = ifrz - dfrad; iz <= ifrzu; iz++) {
3344 int giz = mod(iz, fftZ);
3345 xf[2] = iz * ifftZ;
3346 crystal.toCartesianCoordinates(xf, xc);
3347 final int ii = interleavedIndex(gix, giy, giz, fftX, fftY);
3348 solventff.rhoGrad(
3349 xc,
3350 weight * getDensityGrid()[ii] * dfcmult * getSolventGrid()[ii],
3351 refinementmode);
3352 }
3353 }
3354 }
3355 }
3356 }
3357 }
3358 }
3359
3360 private class BabinetRegion extends ParallelRegion {
3361
3362 BabinetLoop[] babinetLoops;
3363
3364 public BabinetRegion(int nThreads) {
3365 babinetLoops = new BabinetLoop[nThreads];
3366 }
3367
3368 @Override
3369 public void run() throws Exception {
3370 int ti = getThreadIndex();
3371
3372 if (babinetLoops[ti] == null) {
3373 babinetLoops[ti] = new BabinetLoop();
3374 }
3375
3376 try {
3377 execute(0, fftZ - 1, babinetLoops[ti]);
3378 } catch (Exception e) {
3379 logger.info(e.toString());
3380 }
3381 }
3382
3383 private class BabinetLoop extends IntegerForLoop {
3384
3385 @Override
3386 public void run(int lb, int ub) throws Exception {
3387 switch (solventModel) {
3388 case BINARY:
3389 case POLYNOMIAL:
3390 for (int k = lb; k <= ub; k++) {
3391 for (int j = 0; j < fftY; j++) {
3392 for (int i = 0; i < fftX; i++) {
3393 final int ii = interleavedIndex(i, j, k, fftX, fftY);
3394 densityGrid[ii] = 1.0 - getDensityGrid()[ii];
3395 }
3396 }
3397 }
3398 break;
3399 case GAUSSIAN:
3400 for (int k = lb; k <= ub; k++) {
3401 for (int j = 0; j < fftY; j++) {
3402 for (int i = 0; i < fftX; i++) {
3403 final int ii = interleavedIndex(i, j, k, fftX, fftY);
3404 densityGrid[ii] = 1.0 - exp(-solventA * getDensityGrid()[ii]);
3405 }
3406 }
3407 }
3408 break;
3409 default:
3410 }
3411 }
3412 }
3413 }
3414
3415 private class SolventGridRegion extends ParallelRegion {
3416
3417 SolventGridLoop[] solventGridLoops;
3418
3419 public SolventGridRegion(int nThreads) {
3420 solventGridLoops = new SolventGridLoop[nThreads];
3421 }
3422
3423 @Override
3424 public void run() throws Exception {
3425 int ti = getThreadIndex();
3426
3427 if (solventGridLoops[ti] == null) {
3428 solventGridLoops[ti] = new SolventGridLoop();
3429 }
3430
3431 try {
3432 execute(0, fftZ - 1, solventGridLoops[ti]);
3433 } catch (Exception e) {
3434 logger.info(e.toString());
3435 }
3436 }
3437
3438 private class SolventGridLoop extends IntegerForLoop {
3439
3440 @Override
3441 public void run(int lb, int ub) throws Exception {
3442
3443 for (int k = lb; k <= ub; k++) {
3444 for (int j = 0; j < fftY; j++) {
3445 for (int i = 0; i < fftX; i++) {
3446 final int ii = interleavedIndex(i, j, k, fftX, fftY);
3447 solventGrid[ii] = exp(-solventA * getSolventGrid()[ii]);
3448 }
3449 }
3450 }
3451 }
3452 }
3453 }
3454
3455 private class InitRegion extends ParallelRegion {
3456
3457 InitLoop[] initLoops;
3458 FormFactorUpdateLoop[] formFactorUpdateLoops;
3459 int nHKL = reflectionList.hklList.size();
3460 double[][] hkldata = null;
3461
3462 public InitRegion(int nThreads) {
3463 initLoops = new InitLoop[nThreads];
3464 formFactorUpdateLoops = new FormFactorUpdateLoop[nThreads];
3465 }
3466
3467 @Override
3468 public void run() throws Exception {
3469 int ti = getThreadIndex();
3470
3471 if (initLoops[ti] == null) {
3472 initLoops[ti] = new InitLoop();
3473 formFactorUpdateLoops[ti] = new FormFactorUpdateLoop();
3474 }
3475
3476 try {
3477 execute(0, nHKL - 1, initLoops[ti]);
3478 execute(0, nAtoms - 1, formFactorUpdateLoops[ti]);
3479 } catch (Exception e) {
3480 logger.info(e.toString());
3481 }
3482 }
3483
3484 public void setHKL(double[][] hkldata) {
3485 this.hkldata = hkldata;
3486 }
3487
3488 private class InitLoop extends IntegerForLoop {
3489
3490 @Override
3491 public void run(int lb, int ub) throws Exception {
3492 for (int i = lb; i <= ub; i++) {
3493 hkldata[i][0] = 0.0;
3494 hkldata[i][1] = 0.0;
3495 }
3496 }
3497 }
3498
3499 private class FormFactorUpdateLoop extends IntegerForLoop {
3500
3501 private final double[] xyz = new double[3];
3502
3503 @Override
3504 public void run(int lb, int ub) throws Exception {
3505 for (int iSymm = 0; iSymm < bulkNSymm; iSymm++) {
3506 for (int i = lb; i <= ub; i++) {
3507 xyz[0] = coordinates[iSymm][0][i];
3508 xyz[1] = coordinates[iSymm][1][i];
3509 xyz[2] = coordinates[iSymm][2][i];
3510 if (solventFormFactors != null) {
3511 solventFormFactors[iSymm][i].update(xyz);
3512 }
3513 if (atomFormFactors != null) {
3514 atomFormFactors[iSymm][i].update(xyz, bAdd);
3515 }
3516 }
3517 }
3518 }
3519 }
3520 }
3521
3522 private class AtomicScaleRegion extends ParallelRegion {
3523
3524 AtomicScaleLoop[] atomicScaleLoops;
3525 int nHKL = reflectionList.hklList.size();
3526 double[][] hklData = null;
3527 double scale;
3528
3529 public AtomicScaleRegion(int nThreads) {
3530 atomicScaleLoops = new AtomicScaleLoop[nThreads];
3531 }
3532
3533 @Override
3534 public void run() throws Exception {
3535 int ti = getThreadIndex();
3536
3537 if (atomicScaleLoops[ti] == null) {
3538 atomicScaleLoops[ti] = new AtomicScaleLoop();
3539 }
3540
3541 try {
3542 execute(0, nHKL - 1, atomicScaleLoops[ti]);
3543 } catch (Exception e) {
3544 logger.info(e.toString());
3545 }
3546 }
3547
3548 public void setHKL(double[][] hkldata) {
3549 this.hklData = hkldata;
3550 }
3551
3552 public void setScale(double scale) {
3553 this.scale = scale;
3554 }
3555
3556 private class AtomicScaleLoop extends IntegerForLoop {
3557
3558 ComplexNumber c;
3559
3560 public AtomicScaleLoop() {
3561 c = new ComplexNumber();
3562 }
3563
3564 @Override
3565 public void run(int lb, int ub) throws Exception {
3566 for (int i = lb; i <= ub; i++) {
3567 HKL ih = reflectionList.hklList.get(i);
3568 double[] fc = hklData[ih.getIndex()];
3569 c.re(fc[0]);
3570 c.im(fc[1]);
3571
3572 double s = crystal.invressq(ih);
3573 c.timesIP(scale * exp(0.25 * bAdd * s));
3574 c.conjugateIP();
3575 fc[0] = c.re();
3576 fc[1] = c.im();
3577 }
3578 }
3579 }
3580 }
3581
3582 private class SolventScaleRegion extends ParallelRegion {
3583
3584 SolventScaleLoop[] solventScaleLoops;
3585 int nHKL = reflectionList.hklList.size();
3586 double[][] hkldata = null;
3587 double scale;
3588
3589 public SolventScaleRegion(int nThreads) {
3590 solventScaleLoops = new SolventScaleLoop[nThreads];
3591 }
3592
3593 @Override
3594 public void run() throws Exception {
3595 int ti = getThreadIndex();
3596
3597 if (solventScaleLoops[ti] == null) {
3598 solventScaleLoops[ti] = new SolventScaleLoop();
3599 }
3600
3601 try {
3602 execute(0, nHKL - 1, solventScaleLoops[ti]);
3603 } catch (Exception e) {
3604 logger.info(e.toString());
3605 }
3606 }
3607
3608 public void setHKL(double[][] hkldata) {
3609 this.hkldata = hkldata;
3610 }
3611
3612 public void setScale(double scale) {
3613 this.scale = scale;
3614 }
3615
3616 private class SolventScaleLoop extends IntegerForLoop {
3617
3618 ComplexNumber c;
3619
3620 public SolventScaleLoop() {
3621 c = new ComplexNumber();
3622 }
3623
3624 @Override
3625 public void run(int lb, int ub) throws Exception {
3626 for (int i = lb; i <= ub; i++) {
3627 HKL ih = reflectionList.hklList.get(i);
3628 double[] fc = hkldata[ih.getIndex()];
3629 c.re(fc[0]);
3630 c.im(fc[1]);
3631 c.timesIP(scale);
3632 c.conjugateIP();
3633
3634 fc[0] = -c.re();
3635 fc[1] = -c.im();
3636 }
3637 }
3638 }
3639 }
3640
3641 private class ExtractRegion extends ParallelRegion {
3642
3643 ExtractLoop[] extractLoops;
3644 int nHKL = reflectionList.hklList.size();
3645 double[][] hkldata = null;
3646
3647 public ExtractRegion(int nThreads) {
3648 extractLoops = new ExtractLoop[nThreads];
3649 }
3650
3651 @Override
3652 public void run() throws Exception {
3653 int ti = getThreadIndex();
3654
3655 if (extractLoops[ti] == null) {
3656 extractLoops[ti] = new ExtractLoop();
3657 }
3658
3659 try {
3660 execute(0, nHKL - 1, extractLoops[ti]);
3661 } catch (Exception e) {
3662 logger.info(e.toString());
3663 }
3664 }
3665
3666 public void setHKL(double[][] hkldata) {
3667 this.hkldata = hkldata;
3668 }
3669
3670 private class ExtractLoop extends IntegerForLoop {
3671
3672 ComplexNumber c;
3673 int nsym;
3674 List<SymOp> symops;
3675 HKL ij;
3676
3677 public ExtractLoop() {
3678 c = new ComplexNumber();
3679 nsym = crystal.spaceGroup.symOps.size();
3680 symops = crystal.spaceGroup.symOps;
3681 ij = new HKL();
3682 }
3683
3684 @Override
3685 public void run(int lb, int ub) throws Exception {
3686 for (int i = lb; i <= ub; i++) {
3687 HKL ih = reflectionList.hklList.get(i);
3688 double[] fc = hkldata[ih.getIndex()];
3689
3690 for (int j = 0; j < nsym; j++) {
3691 SymOp symOp = symops.get(j);
3692 applyTransSymRot(ih, ij, symOp);
3693 double shift = symOp.symPhaseShift(ih);
3694 int h = mod(ij.getH(), fftX);
3695 int k = mod(ij.getK(), fftY);
3696 int l = mod(ij.getL(), fftZ);
3697 if (h < halfFFTX + 1) {
3698 final int ii = interleavedIndex(h, k, l, fftX, fftY);
3699 c.re(getDensityGrid()[ii]);
3700 c.im(getDensityGrid()[ii + 1]);
3701 } else {
3702 h = (fftX - h) % fftX;
3703 k = (fftY - k) % fftY;
3704 l = (fftZ - l) % fftZ;
3705 final int ii = interleavedIndex(h, k, l, fftX, fftY);
3706 c.re(getDensityGrid()[ii]);
3707 c.im(-getDensityGrid()[ii + 1]);
3708 }
3709 c.phaseShiftIP(shift);
3710 fc[0] += c.re();
3711 fc[1] += c.im();
3712 }
3713 }
3714 }
3715 }
3716 }
3717
3718 private class ArrayCopyRegion extends ParallelRegion {
3719
3720 ArrayCopyLoop[] arrayCopyLoops;
3721
3722 public ArrayCopyRegion() {
3723 arrayCopyLoops = new ArrayCopyLoop[threadCount];
3724 }
3725
3726 @Override
3727 public void run() throws Exception {
3728 int ti = getThreadIndex();
3729
3730 if (arrayCopyLoops[ti] == null) {
3731 arrayCopyLoops[ti] = new ArrayCopyLoop();
3732 }
3733
3734 try {
3735 execute(0, complexFFT3DSpace - 1, arrayCopyLoops[ti]);
3736 } catch (Exception e) {
3737 logger.info(e.toString());
3738 }
3739 }
3740
3741 private class ArrayCopyLoop extends IntegerForLoop {
3742
3743 @Override
3744 public void run(int lb, int ub) throws Exception {
3745 arraycopy(getDensityGrid(), lb, getSolventGrid(), lb, ub - lb);
3746 }
3747 }
3748 }
3749 }