View Javadoc
1   // ******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2025.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
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.ParallelTeam;
44  import ffx.numerics.atomic.AtomicDoubleArray3D;
45  import ffx.potential.bonded.Atom;
46  import ffx.potential.parameters.ForceField;
47  import ffx.potential.parameters.MultipoleType.MultipoleFrameDefinition;
48  import ffx.potential.utils.EnergyException;
49  
50  import java.util.Arrays;
51  import java.util.logging.Level;
52  import java.util.logging.Logger;
53  
54  import static java.lang.Double.isInfinite;
55  import static java.lang.Double.isNaN;
56  import static java.lang.String.format;
57  
58  /**
59   * Parallel conversion of torques into forces, and then reduce them.
60   *
61   * @author Michael J. Schnieders
62   * @since 1.0
63   */
64  public class ReduceRegion extends ParallelRegion {
65  
66    private static final Logger logger = Logger.getLogger(ReduceRegion.class.getName());
67  
68    /**
69     * If set to false, multipoles are fixed in their local frame and torques are zero, which is
70     * useful for narrowing down discrepancies between analytic and finite-difference
71     * derivatives(default is true).
72     */
73    private final boolean rotateMultipoles;
74  
75    private final TorqueLoop[] torqueLoop;
76    private final ReduceLoop[] reduceLoop;
77    /**
78     * If lambdaTerm is true, some ligand atom interactions with the environment are being turned
79     * on/off.
80     */
81    private boolean lambdaTerm;
82    /**
83     * If true, compute coordinate gradient.
84     */
85    private boolean gradient;
86    /**
87     * An ordered array of atoms in the system.
88     */
89    private Atom[] atoms;
90    /**
91     * Dimensions of [nsymm][xyz][nAtoms].
92     */
93    private double[][][] coordinates;
94    /**
95     * Multipole frame definition.
96     */
97    private MultipoleFrameDefinition[] frame;
98    /**
99     * Multipole frame defining atoms.
100    */
101   private int[][] axisAtom;
102   /**
103    * Atomic Gradient array.
104    */
105   private AtomicDoubleArray3D grad;
106   /**
107    * Atomic Torque array.
108    */
109   private AtomicDoubleArray3D torque;
110   /**
111    * Partial derivative of the gradient with respect to Lambda.
112    */
113   private AtomicDoubleArray3D lambdaGrad;
114   /**
115    * Partial derivative of the torque with respect to Lambda.
116    */
117   private AtomicDoubleArray3D lambdaTorque;
118 
119   public ReduceRegion(int threadCount, ForceField forceField) {
120     torqueLoop = new TorqueLoop[threadCount];
121     reduceLoop = new ReduceLoop[threadCount];
122     rotateMultipoles = forceField.getBoolean("ROTATE_MULTIPOLES", true);
123   }
124 
125   /**
126    * Execute the ReduceRegion with the passed ParallelTeam.
127    *
128    * @param parallelTeam The ParallelTeam instance to execute with.
129    */
130   public void executeWith(ParallelTeam parallelTeam) {
131     try {
132       parallelTeam.execute(this);
133     } catch (Exception e) {
134       String message = "Exception calculating torques.";
135       logger.log(Level.SEVERE, message, e);
136     }
137   }
138 
139   public void init(
140       boolean lambdaTerm,
141       boolean gradient,
142       Atom[] atoms,
143       double[][][] coordinates,
144       MultipoleFrameDefinition[] frame,
145       int[][] axisAtom,
146       AtomicDoubleArray3D grad,
147       AtomicDoubleArray3D torque,
148       AtomicDoubleArray3D lambdaGrad,
149       AtomicDoubleArray3D lambdaTorque) {
150     this.lambdaTerm = lambdaTerm;
151     this.gradient = gradient;
152     this.atoms = atoms;
153     this.coordinates = coordinates;
154     this.frame = frame;
155     this.axisAtom = axisAtom;
156     this.grad = grad;
157     this.torque = torque;
158     this.lambdaGrad = lambdaGrad;
159     this.lambdaTorque = lambdaTorque;
160   }
161 
162   @Override
163   public void run() throws EnergyException {
164     int nAtoms = atoms.length;
165     try {
166       int threadIndex = getThreadIndex();
167       if (torqueLoop[threadIndex] == null) {
168         torqueLoop[threadIndex] = new TorqueLoop();
169         reduceLoop[threadIndex] = new ReduceLoop();
170       }
171       if (rotateMultipoles) {
172         execute(0, nAtoms - 1, torqueLoop[threadIndex]);
173       }
174       execute(0, nAtoms - 1, reduceLoop[threadIndex]);
175     } catch (Exception e) {
176       if (e instanceof EnergyException) {
177         throw (EnergyException) e;
178       }
179       String message = "Fatal exception computing torque in thread " + getThreadIndex() + "\n";
180       logger.log(Level.SEVERE, message, e);
181     }
182   }
183 
184   private class TorqueLoop extends IntegerForLoop {
185 
186     private Torque torques;
187     private int threadID;
188     double[] trq = new double[3];
189     double[][] g = new double[4][3];
190     int[] frameIndex = new int[4];
191 
192     TorqueLoop() {
193       torques = new Torque();
194     }
195 
196     @Override
197     public void run(int lb, int ub) throws EnergyException {
198       if (gradient) {
199         torque.reduce(lb, ub);
200         for (int i = lb; i <= ub; i++) {
201           // Gradients from torques will exist if the frameIndex is not -1.
202           Arrays.fill(frameIndex, -1);
203           trq[0] = torque.getX(i);
204           trq[1] = torque.getY(i);
205           trq[2] = torque.getZ(i);
206 
207           // Check for undefined torques.
208           if (isNaN(trq[0]) || isInfinite(trq[0])
209               || isNaN(trq[1]) || isInfinite(trq[1])
210               || isNaN(trq[2]) || isInfinite(trq[2])) {
211             Atom a = atoms[i];
212             throw new EnergyException(
213                 format(" Undefined torque (%8.3f,%8.3f,%8.3f) for atom %s.", trq[0], trq[1], trq[2], a));
214           }
215 
216           torques.torque(i, 0, trq, frameIndex, g);
217           for (int j = 0; j < 4; j++) {
218             int index = frameIndex[j];
219             if (index >= 0) {
220               double[] gj = g[j];
221 
222               // Check for undefined torques.
223               if (isNaN(gj[0]) || isInfinite(gj[0])
224                   || isNaN(gj[1]) || isInfinite(gj[1])
225                   || isNaN(gj[2]) || isInfinite(gj[2])) {
226                 Atom ai = atoms[i];
227                 Atom aj = atoms[index];
228                 throw new EnergyException(
229                     format(" Undefined gradient (%8.3f,%8.3f,%8.3f)\n For atom: %s\n From torque of atom %s",
230                         gj[0], gj[1], gj[2], aj, ai));
231               }
232 
233 
234               grad.add(threadID, index, gj[0], gj[1], gj[2]);
235             }
236           }
237         }
238       }
239       if (lambdaTerm) {
240         lambdaTorque.reduce(lb, ub);
241         for (int i = lb; i <= ub; i++) {
242           // Gradients from torques will exist if the frameIndex is not -1.
243           Arrays.fill(frameIndex, -1);
244           trq[0] = lambdaTorque.getX(i);
245           trq[1] = lambdaTorque.getY(i);
246           trq[2] = lambdaTorque.getZ(i);
247           torques.torque(i, 0, trq, frameIndex, g);
248           for (int j = 0; j < 4; j++) {
249             int index = frameIndex[j];
250             if (index >= 0) {
251               double[] gj = g[j];
252               lambdaGrad.add(threadID, index, gj[0], gj[1], gj[2]);
253             }
254           }
255         }
256       }
257     }
258 
259     @Override
260     public IntegerSchedule schedule() {
261       return IntegerSchedule.fixed();
262     }
263 
264     @Override
265     public void start() {
266       threadID = getThreadIndex();
267       torques.init(axisAtom, frame, coordinates);
268     }
269   }
270 
271   private class ReduceLoop extends IntegerForLoop {
272 
273     @Override
274     public void run(int lb, int ub) throws EnergyException {
275       if (gradient) {
276         grad.reduce(lb, ub);
277         for (int i = lb; i <= ub; i++) {
278           Atom ai = atoms[i];
279           ai.addToXYZGradient(grad.getX(i), grad.getY(i), grad.getZ(i));
280         }
281       }
282       if (lambdaTerm) {
283         lambdaGrad.reduce(lb, ub);
284       }
285     }
286 
287     @Override
288     public IntegerSchedule schedule() {
289       return IntegerSchedule.fixed();
290     }
291   }
292 }