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.nonbonded.GeneralizedKirkwood;
47  
48  import java.util.logging.Level;
49  import java.util.logging.Logger;
50  
51  import static ffx.potential.parameters.MultipoleType.t001;
52  import static ffx.potential.parameters.MultipoleType.t010;
53  import static ffx.potential.parameters.MultipoleType.t100;
54  
55  /**
56   * Parallel summation and reduction of components of the induced dipole field at each atom.
57   *
58   * @author Michael J. Schnieders
59   * @since 1.0
60   */
61  public class InducedDipoleFieldReduceRegion extends ParallelRegion {
62  
63    private static final Logger logger = Logger.getLogger(DirectRegion.class.getName());
64    private final InducedDipoleFieldReduceLoop[] inducedDipoleFieldReduceLoop;
65    /**
66     * Dimensions of [nsymm][nAtoms][3]
67     */
68    public double[][][] inducedDipole;
69  
70    public double[][][] inducedDipoleCR;
71    /**
72     * An ordered array of atoms in the system.
73     */
74    private Atom[] atoms;
75  
76    private double[][] cartesianDipolePhi;
77    private double[][] cartesianDipolePhiCR;
78    /**
79     * Field array.
80     */
81    private AtomicDoubleArray3D field;
82    /**
83     * Chain rule field array.
84     */
85    private AtomicDoubleArray3D fieldCR;
86    /**
87     * Flag to indicate use of generalized Kirkwood.
88     */
89    private boolean generalizedKirkwoodTerm;
90  
91    private GeneralizedKirkwood generalizedKirkwood;
92    private double aewald;
93    private double aewald3;
94  
95    private double soluteDielectric;
96  
97    public InducedDipoleFieldReduceRegion(int nt) {
98      inducedDipoleFieldReduceLoop = new InducedDipoleFieldReduceLoop[nt];
99    }
100 
101   /**
102    * Execute the InducedDipoleFieldReduceRegion with the passed ParallelTeam.
103    *
104    * @param parallelTeam The ParallelTeam instance to execute with.
105    */
106   public void executeWith(ParallelTeam parallelTeam) {
107     try {
108       parallelTeam.execute(this);
109     } catch (Exception e) {
110       String message = " Exception computing induced dipole field.\n";
111       logger.log(Level.WARNING, message, e);
112     }
113   }
114 
115   public void init(
116       Atom[] atoms,
117       double[][][] inducedDipole,
118       double[][][] inducedDipoleCR,
119       boolean generalizedKirkwoodTerm,
120       GeneralizedKirkwood generalizedKirkwood,
121       EwaldParameters ewaldParameters,
122       double soluteDielectric,
123       double[][] cartesianDipolePhi,
124       double[][] cartesianDipolePhiCR,
125       AtomicDoubleArray3D field,
126       AtomicDoubleArray3D fieldCR) {
127     // Input
128     this.atoms = atoms;
129     this.inducedDipole = inducedDipole;
130     this.inducedDipoleCR = inducedDipoleCR;
131     this.generalizedKirkwoodTerm = generalizedKirkwoodTerm;
132     this.generalizedKirkwood = generalizedKirkwood;
133     this.aewald = ewaldParameters.aewald;
134     this.aewald3 = ewaldParameters.aewald3;
135     this.soluteDielectric = soluteDielectric;
136     this.cartesianDipolePhi = cartesianDipolePhi;
137     this.cartesianDipolePhiCR = cartesianDipolePhiCR;
138     // Output
139     this.field = field;
140     this.fieldCR = fieldCR;
141   }
142 
143   @Override
144   public void run() throws Exception {
145     try {
146       int threadID = getThreadIndex();
147       if (inducedDipoleFieldReduceLoop[threadID] == null) {
148         inducedDipoleFieldReduceLoop[threadID] = new InducedDipoleFieldReduceLoop();
149       }
150       int nAtoms = atoms.length;
151       execute(0, nAtoms - 1, inducedDipoleFieldReduceLoop[threadID]);
152     } catch (Exception e) {
153       String message = " Fatal exception computing the mutual induced dipoles in thread "
154           + getThreadIndex() + "\n";
155       logger.log(Level.SEVERE, message, e);
156     }
157   }
158 
159   private class InducedDipoleFieldReduceLoop extends IntegerForLoop {
160 
161     @Override
162     public void run(int lb, int ub) throws Exception {
163       int threadID = getThreadIndex();
164 
165       final double[][] induced0 = inducedDipole[0];
166       final double[][] inducedCR0 = inducedDipoleCR[0];
167       // Add the PME self and reciprocal space fields to the real space field.
168       if (aewald > 0.0) {
169         for (int i = lb; i <= ub; i++) {
170           double[] dipolei = induced0[i];
171           double[] dipoleCRi = inducedCR0[i];
172           final double[] phii = cartesianDipolePhi[i];
173           final double[] phiCRi = cartesianDipolePhiCR[i];
174           double fx = aewald3 * dipolei[0] - phii[t100];
175           double fy = aewald3 * dipolei[1] - phii[t010];
176           double fz = aewald3 * dipolei[2] - phii[t001];
177           double fxCR = aewald3 * dipoleCRi[0] - phiCRi[t100];
178           double fyCR = aewald3 * dipoleCRi[1] - phiCRi[t010];
179           double fzCR = aewald3 * dipoleCRi[2] - phiCRi[t001];
180           field.add(threadID, i, fx, fy, fz);
181           fieldCR.add(threadID, i, fxCR, fyCR, fzCR);
182         }
183       }
184 
185       // Reduce the total direct field.
186       field.reduce(lb, ub);
187       fieldCR.reduce(lb, ub);
188 
189       // Scale the total induced field by the inverse dielectric.
190       if (soluteDielectric > 1.0) {
191         double inverseDielectric = 1.0 / soluteDielectric;
192         for (int i = lb; i <= ub; i++) {
193           field.scale(0, i, inverseDielectric);
194           fieldCR.scale(0, i, inverseDielectric);
195         }
196       }
197 
198       // Add the GK reaction field.
199       if (generalizedKirkwoodTerm) {
200         AtomicDoubleArray3D fieldGK = generalizedKirkwood.getFieldGK();
201         AtomicDoubleArray3D fieldGKCR = generalizedKirkwood.getFieldGKCR();
202         for (int i = lb; i <= ub; i++) {
203           field.add(0, i, fieldGK.getX(i), fieldGK.getY(i), fieldGK.getZ(i));
204           fieldCR.add(0, i, fieldGKCR.getX(i), fieldGKCR.getY(i), fieldGKCR.getZ(i));
205         }
206       }
207     }
208 
209     @Override
210     public IntegerSchedule schedule() {
211       return IntegerSchedule.fixed();
212     }
213   }
214 }