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.potential.nonbonded.pme;
39
40 import edu.rit.pj.IntegerForLoop;
41 import edu.rit.pj.IntegerSchedule;
42 import edu.rit.pj.ParallelRegion;
43 import edu.rit.pj.reduction.SharedDouble;
44 import ffx.potential.bonded.Atom;
45 import java.util.logging.Level;
46 import java.util.logging.Logger;
47
48
49
50
51
52
53
54 public class PolarizationEnergyRegion extends ParallelRegion {
55
56 private static final Logger logger = Logger.getLogger(PolarizationEnergyRegion.class.getName());
57 private final double electric;
58 private final PolarizationEnergyLoop[] polarizationLoop;
59 private final SharedDouble polarizationEnergy = new SharedDouble();
60
61 private Atom[] atoms;
62
63 private double[] polarizability;
64
65 private double[][][] inducedDipole;
66
67 private double[][] directDipoleCR;
68
69 private double polarizationScale;
70
71 public PolarizationEnergyRegion(int nt, double electric) {
72 this.electric = electric;
73 polarizationLoop = new PolarizationEnergyLoop[nt];
74 }
75
76 @Override
77 public void finish() {
78 double energy = polarizationEnergy.get();
79 polarizationEnergy.set(energy * polarizationScale * electric);
80 }
81
82
83
84
85
86
87 public double getPolarizationEnergy() {
88 return polarizationEnergy.get();
89 }
90
91
92
93
94
95
96 public void setPolarizationEnergy(double energy) {
97 polarizationEnergy.set(energy);
98 }
99
100 public void init(
101 Atom[] atoms,
102 double[] polarizability,
103 double[][][] inducedDipole,
104 double[][] directDipoleCR,
105 double polarizationScale) {
106 this.atoms = atoms;
107 this.polarizability = polarizability;
108 this.inducedDipole = inducedDipole;
109 this.directDipoleCR = directDipoleCR;
110 this.polarizationScale = polarizationScale;
111 }
112
113 @Override
114 public void run() throws Exception {
115 try {
116 int ti = getThreadIndex();
117 if (polarizationLoop[ti] == null) {
118 polarizationLoop[ti] = new PolarizationEnergyLoop();
119 }
120 int nAtoms = atoms.length;
121 execute(0, nAtoms - 1, polarizationLoop[ti]);
122 } catch (RuntimeException ex) {
123 logger.warning(
124 "Runtime exception computing polarization energy in thread " + getThreadIndex());
125 throw ex;
126 } catch (Exception e) {
127 String message =
128 "Fatal exception computing polarization energy in thread " + getThreadIndex() + "\n";
129 logger.log(Level.SEVERE, message, e);
130 }
131 }
132
133 @Override
134 public void start() {
135 polarizationEnergy.set(0.0);
136 }
137
138 private class PolarizationEnergyLoop extends IntegerForLoop {
139
140 @Override
141 public void run(int lb, int ub) throws Exception {
142 double energy = 0.0;
143 for (int i = lb; i <= ub; i++) {
144 if (polarizability[i] == 0.0) {
145 continue;
146 }
147 double uix = inducedDipole[0][i][0];
148 double uiy = inducedDipole[0][i][1];
149 double uiz = inducedDipole[0][i][2];
150 double pix = directDipoleCR[i][0];
151 double piy = directDipoleCR[i][1];
152 double piz = directDipoleCR[i][2];
153 energy += (uix * pix + uiy * piy + uiz * piz) / polarizability[i];
154 }
155 polarizationEnergy.addAndGet(-0.5 * energy);
156 }
157
158 @Override
159 public IntegerSchedule schedule() {
160 return IntegerSchedule.fixed();
161 }
162 }
163 }