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