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