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-2024.
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.numerics.multipole;
39  
40  import static ffx.numerics.multipole.GKSource.GK_MULTIPOLE_ORDER.DIPOLE;
41  import static ffx.numerics.multipole.GKSource.GK_MULTIPOLE_ORDER.MONOPOLE;
42  import static ffx.numerics.multipole.GKSource.GK_MULTIPOLE_ORDER.QUADRUPOLE;
43  import static ffx.numerics.multipole.GKSource.GK_TENSOR_MODE.BORN;
44  import static ffx.numerics.multipole.GKSource.GK_TENSOR_MODE.POTENTIAL;
45  import static java.util.Arrays.fill;
46  
47  public class GKEnergyGlobal {
48  
49    private final GKSource gkSource;
50    private final GKTensorGlobal gkMonopole;
51    private final GKTensorGlobal gkDipole;
52    private final GKTensorGlobal gkQuadrupole;
53  
54    /**
55     * Constructor for GKEnergyGlobal.
56     *
57     * @param gkc The GK generalizing function constant.
58     * @param epsilon The solvent dielectric.
59     * @param gradient If true, compute the gradient and torque.
60     */
61    public GKEnergyGlobal(double gkc, double epsilon, boolean gradient) {
62      int monopoleOrder = 2;
63      int dipoleOrder = 3;
64      int quadrupoleOrder = 4;
65      if (gradient) {
66        monopoleOrder = 3;
67        dipoleOrder = 4;
68        quadrupoleOrder = 5;
69      }
70      gkSource = new GKSource(quadrupoleOrder, gkc);
71      gkMonopole = new GKTensorGlobal(MONOPOLE, monopoleOrder, gkSource, 1.0, epsilon);
72      gkDipole = new GKTensorGlobal(DIPOLE, dipoleOrder, gkSource, 1.0, epsilon);
73      gkQuadrupole = new GKTensorGlobal(QUADRUPOLE, quadrupoleOrder, gkSource, 1.0, epsilon);
74    }
75  
76    public void initPotential(double[] r, double r2, double rbi, double rbk) {
77      gkSource.generateSource(POTENTIAL, QUADRUPOLE, r2, rbi, rbk);
78      gkMonopole.setR(r);
79      gkDipole.setR(r);
80      gkQuadrupole.setR(r);
81      gkMonopole.generateTensor();
82      gkDipole.generateTensor();
83      gkQuadrupole.generateTensor();
84    }
85  
86    public void initBorn(double[] r, double r2, double rbi, double rbk) {
87      gkSource.generateSource(BORN, QUADRUPOLE, r2, rbi, rbk);
88      gkMonopole.setR(r);
89      gkDipole.setR(r);
90      gkQuadrupole.setR(r);
91      gkMonopole.generateTensor();
92      gkDipole.generateTensor();
93      gkQuadrupole.generateTensor();
94    }
95  
96    public double multipoleEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
97      double em = gkMonopole.multipoleEnergy(mI, mK);
98      double ed = gkDipole.multipoleEnergy(mI, mK);
99      double eq = gkQuadrupole.multipoleEnergy(mI, mK);
100     return em + ed + eq;
101   }
102 
103   public double polarizationEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
104     double emp = gkMonopole.polarizationEnergy(mI, mK);
105     double edp = gkDipole.polarizationEnergy(mI, mK);
106     double eqp = gkQuadrupole.polarizationEnergy(mI, mK);
107     return emp + edp + eqp;
108   }
109 
110   public double multipoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
111       double[] gradI, double[] torqueI, double[] torqueK) {
112     double[] gI = new double[3];
113     double[] gK = new double[3];
114     double[] tI = new double[3];
115     double[] tK = new double[3];
116     double em = gkMonopole.multipoleEnergyAndGradient(mI, mK, gI, gK, tI, tK);
117     for (int j = 0; j < 3; j++) {
118       gradI[j] = gI[j];
119       torqueI[j] = tI[j];
120       torqueK[j] = tK[j];
121     }
122     fill(gI, 0.0);
123     fill(gK, 0.0);
124     fill(tI, 0.0);
125     fill(tK, 0.0);
126     double ed = gkDipole.multipoleEnergyAndGradient(mI, mK, gI, gK, tI, tK);
127     for (int j = 0; j < 3; j++) {
128       gradI[j] += gI[j];
129       torqueI[j] += tI[j];
130       torqueK[j] += tK[j];
131     }
132     fill(gI, 0.0);
133     fill(gK, 0.0);
134     fill(tI, 0.0);
135     fill(tK, 0.0);
136     double eq = gkQuadrupole.multipoleEnergyAndGradient(mI, mK, gI, gK, tI, tK);
137     for (int j = 0; j < 3; j++) {
138       gradI[j] += gI[j];
139       torqueI[j] += tI[j];
140       torqueK[j] += tK[j];
141     }
142 
143     return em + ed + eq;
144   }
145 
146   public double polarizationEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
147       double mutualMask,
148       double[] gradI, double[] torqueI, double[] torqueK) {
149     double[] gI = new double[3];
150     double[] tI = new double[3];
151     double[] tK = new double[3];
152     double emp = gkMonopole.polarizationEnergyAndGradient(mI, mK, 1.0, 1.0, mutualMask, gI, tI, tK);
153     for (int j = 0; j < 3; j++) {
154       gradI[j] = gI[j];
155       torqueI[j] = tI[j];
156       torqueK[j] = tK[j];
157     }
158     fill(gI, 0.0);
159     fill(tI, 0.0);
160     fill(tK, 0.0);
161     double edp = gkDipole.polarizationEnergyAndGradient(mI, mK, 1.0, 1.0, mutualMask, gI, tI, tK);
162     for (int j = 0; j < 3; j++) {
163       gradI[j] += gI[j];
164       torqueI[j] += tI[j];
165       torqueK[j] += tK[j];
166     }
167     fill(gI, 0.0);
168     fill(tI, 0.0);
169     fill(tK, 0.0);
170     double eqp = gkQuadrupole.polarizationEnergyAndGradient(mI, mK, 1.0, 1.0, mutualMask, gI, tI,
171         tK);
172     for (int j = 0; j < 3; j++) {
173       gradI[j] += gI[j];
174       torqueI[j] += tI[j];
175       torqueK[j] += tK[j];
176     }
177 
178     // Sum the GK polarization interaction energy.
179     return emp + edp + eqp;
180   }
181 
182   public double multipoleEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
183     double db = gkMonopole.multipoleEnergyBornGrad(mI, mK);
184     db += gkDipole.multipoleEnergyBornGrad(mI, mK);
185     db += gkQuadrupole.multipoleEnergyBornGrad(mI, mK);
186     return db;
187   }
188 
189   public double polarizationEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK,
190       boolean mutual) {
191     // Compute the GK polarization Born chain-rule term.
192     double db = gkMonopole.polarizationEnergyBornGrad(mI, mK);
193     db += gkDipole.polarizationEnergyBornGrad(mI, mK);
194     db += gkQuadrupole.polarizationEnergyBornGrad(mI, mK);
195     // Add the mutual polarization contribution to Born chain-rule term.
196     if (mutual) {
197       db += gkDipole.mutualPolarizationEnergyBornGrad(mI, mK);
198     }
199     return db;
200   }
201 
202 }