View Javadoc
1   package ffx.numerics.multipole;
2   
3   import org.junit.Test;
4   import org.junit.runner.RunWith;
5   import org.junit.runners.Parameterized;
6   
7   import java.util.Arrays;
8   import java.util.Collection;
9   import java.util.logging.Logger;
10  
11  import static ffx.numerics.multipole.MultipoleTensorTest.*;
12  import static java.lang.Math.abs;
13  import static java.lang.Math.sqrt;
14  import static org.junit.Assert.assertEquals;
15  
16  @RunWith(Parameterized.class)
17  public class CombinedTensorTest {
18      /**
19       * Logger for the MultipoleTensor class.
20       */
21      private static final Logger logger = Logger.getLogger(AMOEBAPlusMultipoleTensorGlobalTest.class.getName());
22      private final double tolerance = 1.0e-13;
23      private final double fdTolerance = 1.0e-6;
24      private final double[] r2 = new double[3];
25      private final int order;
26      private final String info;
27      private final Operator operator;
28      private final double ewaldA;
29      private final double thole;
30      private final double aiak;
31      private final double alphaI;
32      private final double alphaK;
33      private double[] ewaldTerms;
34      private double[] dampTerms;
35      private double[] overlapTerms;
36      private double[] tholeDirect;
37      private double[] tholeTerms;
38      private double R;
39  
40      public CombinedTensorTest(
41              String info,
42              int order,
43              Operator operator) {
44          this.info = info;
45          this.order = order;
46          this.operator = operator;
47          this.ewaldA = 0.54459051633620059;
48          this.thole = .7;
49          this.aiak = 1.0 / 0.86460072979796410;
50          this.alphaI = MultipoleTensorTest.chargePenAlphaOx;
51          this.alphaK = MultipoleTensorTest.chargePenAlphaHyd;
52          updateSources(MultipoleTensorTest.xyzEwaldAPlus);
53      }
54  
55      @Parameterized.Parameters
56      public static Collection<Object[]> data() {
57          return Arrays.asList(
58                  new Object[][]{
59                          {"Order 5 Aplus Direct", 5, Operator.AMOEBA_PLUS_OVERLAP_FIELD}
60                  });
61      }
62  
63      private void updateSources(double[] xyz){
64          this.R = sqrt(xyz[0]*xyz[0]
65                  + xyz[1]*xyz[1]
66                  + xyz[2]*xyz[2]);
67          // Use static tensor methods for terms fill with ones to not include coulomb terms already
68          this.ewaldTerms = new double[order+1];
69          Arrays.fill(ewaldTerms, 1.0);
70          EwaldTensorGlobal.fillEwaldSource(order, ewaldA,
71                  EwaldTensorGlobal.initEwaldSource(order, ewaldA, new double[order+1]),
72                  R, ewaldTerms);
73          this.dampTerms = new double[order+1];
74          Arrays.fill(dampTerms, 1.0);
75          AmoebaPlusDampTensorGlobal.dampSource(alphaI, R, dampTerms);
76          this.overlapTerms = new double[order+1];
77          Arrays.fill(overlapTerms, 1.0);
78          AmoebaPlusOverlapTensorGlobal.overlapSource(alphaI, alphaK,
79                  alphaK*alphaK/(alphaK*alphaK - alphaI*alphaI),
80                  alphaI*alphaI/(alphaI*alphaI - alphaK*alphaK), R,
81                  overlapTerms);
82          this.tholeDirect = new double[order+1];
83          Arrays.fill(tholeDirect, 1.0);
84          TholeTensorGlobal.tholeSource(thole, aiak, R, true, tholeDirect);
85          this.tholeTerms = new double[order+1];
86          Arrays.fill(tholeTerms, 1.0);
87          TholeTensorGlobal.tholeSource(thole, aiak, R, false, tholeTerms);
88      }
89  
90      private void buildOverlapSource(CombinedTensorGlobal multipoleTensor){
91          multipoleTensor.resetSource();
92          multipoleTensor.addTermsSeparate(ewaldTerms);
93          multipoleTensor.addTerms(overlapTerms);
94          multipoleTensor.addCoulombMultiplier(-1.0);
95      }
96  
97      @Test
98      public void permanentMultipoleEnergyAndGradTest() {
99          double delta = 1.0e-8;
100         double delta2 = 2.0 * delta;
101         double[] Fi = new double[3];
102         double[] Fk = new double[3];
103         double[] Ti = new double[3];
104         double[] Tk = new double[3];
105 
106         // Oxygen-Hydrogen A+ Real-Space Energy under Ewald
107         double[] uI = new double[]{-9.5291976837210871E-002,  -5.1217206248926346E-002,   9.0954068003896951E-002};
108         double[] uK = new double[]{-1.4736767767345298E-002,   1.3165041117687385E-002,  -7.4174987487104807E-005};
109         PolarizableMultipole mI = new PolarizableMultipole(QlAPlusEwal,uI,uI,8.0);
110         PolarizableMultipole mK = new PolarizableMultipole(QmAplusEwal,uK,uK,1.0);
111         CombinedTensorGlobal multipoleTensor = new CombinedTensorGlobal(order);
112         multipoleTensor.setR(xyzEwaldAPlus);
113 
114         //Amoeba+ = M-pol*(Ewald - Coulomb -CDamp)*M-pol
115         buildOverlapSource(multipoleTensor);
116         multipoleTensor.generateTensor();
117         double eOverlap = multipoleTensor.multipoleEnergyAndGradient(mI, mK, Fi, Fk, Ti, Tk);
118         double overlapAns = 1.4311540110919852E-003;
119         assertEquals("Real-Space PME Mpole-Mpole interaction", overlapAns, eOverlap, tolerance);
120 
121         // Analytic gradient.
122         double aX = Fk[0];
123         double aY = Fk[1];
124         double aZ = Fk[2];
125 
126         xyzEwaldAPlus[0] += delta;
127         updateSources(xyzEwaldAPlus);
128         buildOverlapSource(multipoleTensor);
129         multipoleTensor.generateTensor(xyzEwaldAPlus);
130         double posX = multipoleTensor.multipoleEnergy(mI, mK);
131         xyzEwaldAPlus[0] -= delta2;
132         updateSources(xyzEwaldAPlus);
133         buildOverlapSource(multipoleTensor);
134         multipoleTensor.generateTensor(xyzEwaldAPlus);
135         double negX = multipoleTensor.multipoleEnergy(mI, mK);
136         xyzEwaldAPlus[0] += delta;
137         xyzEwaldAPlus[1] += delta;
138         updateSources(xyzEwaldAPlus);
139         buildOverlapSource(multipoleTensor);
140         multipoleTensor.generateTensor(xyzEwaldAPlus);
141         double posY = multipoleTensor.multipoleEnergy(mI, mK);
142         xyzEwaldAPlus[1] -= delta2;
143         updateSources(xyzEwaldAPlus);
144         buildOverlapSource(multipoleTensor);
145         multipoleTensor.generateTensor(xyzEwaldAPlus);
146         double negY = multipoleTensor.multipoleEnergy(mI, mK);
147         xyzEwaldAPlus[1] += delta;
148         xyzEwaldAPlus[2] += delta;
149         updateSources(xyzEwaldAPlus);
150         buildOverlapSource(multipoleTensor);
151         multipoleTensor.generateTensor(xyzEwaldAPlus);
152         double posZ = multipoleTensor.multipoleEnergy(mI, mK);
153         xyzEwaldAPlus[2] -= delta2;
154         updateSources(xyzEwaldAPlus);
155         buildOverlapSource(multipoleTensor);
156         multipoleTensor.generateTensor(xyzEwaldAPlus);
157         double negZ = multipoleTensor.multipoleEnergy(mI, mK);
158         xyzEwaldAPlus[2] += delta;
159 
160         double expect = aX;
161         double actual = (posX - negX) / delta2;
162         assertEquals(info + " OH Ewald Overlap Force X", expect, actual, fdTolerance);
163         expect = aY;
164         actual = (posY - negY) / delta2;
165         assertEquals(info + " OH Ewald Overlap Force Y", expect, actual, fdTolerance);
166         expect = aZ;
167         actual = (posZ - negZ) / delta2;
168         assertEquals(info + " OH Ewald Overlap Force Z", expect, actual, fdTolerance);
169         Arrays.fill(Fi, 0);
170         Arrays.fill(Fk, 0);
171 
172         // A+ CoreE = Core*(Ewald - Coulomb + Damp)*Mpole -> 2x for i->j j->i energies
173         AmoebaPlusDampTensorGlobal newTensor = new AmoebaPlusDampTensorGlobal(3, alphaI, alphaK, ewaldA);
174         newTensor.setR(xyzEwaldAPlus);
175         newTensor.generateTensor();
176         double eCore = newTensor.coreInteractionAndGradient(mI, mK, Fi, Fk);
177         double answerCore = 1.8003552196793374E-006;
178         assertEquals("Real-Space Core-Mpole Interaction", answerCore, eCore+overlapAns, tolerance);
179 
180         aX = Fk[0];
181         aY = Fk[1];
182         aZ = Fk[2];
183 
184         xyzEwaldAPlus[0] += delta;
185         newTensor.generateTensor(xyzEwaldAPlus);
186         posX = newTensor.coreInteraction(mI, mK);
187         xyzEwaldAPlus[0] -= delta2;
188         newTensor.generateTensor(xyzEwaldAPlus);
189         negX = newTensor.coreInteraction(mI, mK);
190         xyzEwaldAPlus[0] += delta;
191         xyzEwaldAPlus[1] += delta;
192         newTensor.generateTensor(xyzEwaldAPlus);
193         posY = newTensor.coreInteraction(mI, mK);
194         xyzEwaldAPlus[1] -= delta2;
195         newTensor.generateTensor(xyzEwaldAPlus);
196         negY = newTensor.coreInteraction(mI, mK);
197         xyzEwaldAPlus[1] += delta;
198         xyzEwaldAPlus[2] += delta;
199         newTensor.generateTensor(xyzEwaldAPlus);
200         posZ = newTensor.coreInteraction(mI, mK);
201         xyzEwaldAPlus[2] -= delta2;
202         newTensor.generateTensor(xyzEwaldAPlus);
203         negZ = newTensor.coreInteraction(mI, mK);
204         xyzEwaldAPlus[2] += delta;
205 
206         expect = aX;
207         actual = (posX - negX) / delta2;
208         assertEquals(info + " OH Ewald Core Force X", expect, actual, fdTolerance);
209         expect = aY;
210         actual = (posY - negY) / delta2;
211         assertEquals(info + " OH Ewald Core Force Y", expect, actual, fdTolerance);
212         expect = aZ;
213         actual = (posZ - negZ) / delta2;
214         assertEquals(info + " OH Ewald Core Force Z", expect, actual, fdTolerance);
215         Arrays.fill(Fi, 0);
216         Arrays.fill(Fk, 0);
217     }
218 
219     @Test
220     public void polarizationEnergyAndGradientTest() {
221         // Oxygen-Hydrogen A+ Real-Space Direct Pol Energy under Ewald (aplussmall.xyz 1, 536)
222         double[] Fi = new double[3];
223         double[] Ti = new double[3];
224         double[] Tk = new double[3];
225         double[] uI = new double[]{-9.5291976837210871E-002,  -5.1217206248926346E-002,   9.0954068003896951E-002};
226         double[] uK = new double[]{-1.4736767767345298E-002,   1.3165041117687385E-002,  -7.4174987487104807E-005};
227         PolarizableMultipole mI = new PolarizableMultipole(QlAPlusEwal,uI,uI,8.0);
228         PolarizableMultipole mK = new PolarizableMultipole(QmAplusEwal,uK,uK,1.0);
229         CombinedTensorGlobal multipoleTensor = new CombinedTensorGlobal(order);
230         multipoleTensor.setR(xyzEwaldAPlus);
231 
232         // A+ PolE = Dip*(Ewald-TholeDirect)*Mpole
233         multipoleTensor = new CombinedTensorGlobal(order);
234         multipoleTensor.addTermsSeparate(ewaldTerms);
235         multipoleTensor.addTerms(tholeDirect);
236         multipoleTensor.addCoulombMultiplier(-1.0);
237         multipoleTensor.setR(xyzEwaldAPlus);
238         multipoleTensor.generateTensor();
239         double ePolOverlap = multipoleTensor.polarizationEnergyAndGradient(mI, mK, 1.0, 1.0, 0.0, Fi, Ti, Tk);
240         double polAnswer = 4.8757773606071702E-006/2;
241         assertEquals(" OH Real-Space PME PolE interaction", polAnswer, ePolOverlap, tolerance);
242 
243         double[] aplusIndGrad = new double[]{-1.1162913643840848E-005,
244                 6.4235992888231040E-006,   -6.2954494656288314E-006};
245         assertEquals(info + " Polarization GradX", aplusIndGrad[0], Fi[0], tolerance);
246         assertEquals(info + " Polarization GradY", aplusIndGrad[1], Fi[1], tolerance);
247         assertEquals(info + " Polarization GradZ", aplusIndGrad[2], Fi[2], tolerance);
248     }
249 }