diff --git a/ipynb-kotlin/ffx.json b/ipynb-kotlin/ffx.json
index ab5baef550..776511c3a8 100644
--- a/ipynb-kotlin/ffx.json
+++ b/ipynb-kotlin/ffx.json
@@ -100,6 +100,7 @@
"jakarta.activation-api-2.1.3.jar",
"jakarta.xml.bind-api-4.0.2.jar",
"javaparser-core-3.26.4.jar",
+ "javaparser-core-3.27.0.jar",
"javax.annotation-api-1.3.2.jar",
"jaxb-core-4.0.3.jar",
"jaxb-runtime-4.0.3.jar",
diff --git a/modules/algorithms/src/main/groovy/ffx/algorithms/groovy/BAR.groovy b/modules/algorithms/src/main/groovy/ffx/algorithms/groovy/BAR.groovy
index 640a75397b..27bf5be910 100644
--- a/modules/algorithms/src/main/groovy/ffx/algorithms/groovy/BAR.groovy
+++ b/modules/algorithms/src/main/groovy/ffx/algorithms/groovy/BAR.groovy
@@ -112,11 +112,11 @@ class BAR extends AlgorithmsScript {
private boolean sortedArc = false
@Option(names = ["--ss", "--startSnapshot"], paramLabel = "0",
- description = "Start at this snapshot when reading in Tinker BAR files.")
+ description = "Start at this snapshot when reading in Tinker BAR files (indexed from 0).")
private int startingSnapshot = 0
@Option(names = ["--es", "--endSnapshot"], paramLabel = "0",
- description = "End at this snapshot when reading in Tinker BAR files.")
+ description = "End at this snapshot when reading in Tinker BAR files (indexed from 0).")
private int endingSnapshot = 0
/**
diff --git a/modules/numerics/src/main/java/ffx/numerics/fft/Complex3DParallel.java b/modules/numerics/src/main/java/ffx/numerics/fft/Complex3DParallel.java
index 850055ba38..224256da26 100644
--- a/modules/numerics/src/main/java/ffx/numerics/fft/Complex3DParallel.java
+++ b/modules/numerics/src/main/java/ffx/numerics/fft/Complex3DParallel.java
@@ -331,6 +331,7 @@ public Complex3DParallel(int nX, int nY, int nZ, ParallelTeam parallelTeam,
useSIMD = false;
}
+ // Do not pack FFTs by default for now.
packFFTs = false;
String pack = System.getProperty("fft.packFFTs", Boolean.toString(packFFTs));
try {
@@ -360,11 +361,13 @@ public Complex3DParallel(int nX, int nY, int nZ, ParallelTeam parallelTeam,
fftZ[i] = new Complex(nZ, dataLayout1D, internalImZ, nY);
fftZ[i].setUseSIMD(useSIMD);
}
- if (!localZTranspose) {
- work3D = new double[2 * nX * nY * nZ];
- } else {
+
+ if (localZTranspose) {
work3D = null;
+ } else {
+ work3D = new double[2 * nX * nY * nZ];
}
+
fftRegion = new FFTRegion();
ifftRegion = new IFFTRegion();
convRegion = new ConvolutionRegion();
@@ -488,6 +491,7 @@ public void initTiming() {
/**
* Get the timing string.
+ *
* @return The timing string.
*/
public String timingString() {
@@ -1030,11 +1034,15 @@ private class TransposeLoop extends IntegerForLoop {
public void run(final int lb, final int ub) {
for (int x = lb; x <= ub; x++) {
for (int z = 0; z < nZ; z++) {
+ int inputOffset = x * nextX + z * nextZ;
+ int workOffset = x * trNextX + z * trNextZ;
for (int y = 0; y < nY; y++) {
- double real = input[x * nextX + y * nextY + z * nextZ];
- double imag = input[x * nextX + y * nextY + z * nextZ + im];
- work3D[y * trNextY + z * trNextZ + x * trNextX] = real;
- work3D[y * trNextY + z * trNextZ + x * trNextX + internalImZ] = imag;
+ int inputIndex = inputOffset + y * nextY;
+ double real = input[inputIndex];
+ double imag = input[inputIndex + im];
+ int workIndex = workOffset + y * trNextY;
+ work3D[workIndex] = real;
+ work3D[workIndex + internalImZ] = imag;
}
}
}
@@ -1082,11 +1090,15 @@ private class UnTransposeLoop extends IntegerForLoop {
public void run(final int lb, final int ub) {
for (int x = 0; x < nX; x++) {
for (int y = 0; y < nY; y++) {
+ int workOffset = x * trNextX + y * trNextY;
+ int inputOffset = x * nextX + y * nextY;
for (int z = lb; z <= ub; z++) {
- double real = work3D[y * trNextY + z * trNextZ + x * trNextX];
- double imag = work3D[y * trNextY + z * trNextZ + x * trNextX + internalImZ];
- input[x * nextX + y * nextY + z * nextZ] = real;
- input[x * nextX + y * nextY + z * nextZ + im] = imag;
+ int workIndex = workOffset + z * trNextZ;
+ double real = work3D[workIndex];
+ double imag = work3D[workIndex + internalImZ];
+ int inputIndex = inputOffset + z * nextZ;
+ input[inputIndex] = real;
+ input[inputIndex + im] = imag;
}
}
}
diff --git a/modules/numerics/src/test/java/ffx/numerics/fft/Complex3DParallelTest.java b/modules/numerics/src/test/java/ffx/numerics/fft/Complex3DParallelTest.java
index 725ac93293..0804b3ce47 100755
--- a/modules/numerics/src/test/java/ffx/numerics/fft/Complex3DParallelTest.java
+++ b/modules/numerics/src/test/java/ffx/numerics/fft/Complex3DParallelTest.java
@@ -45,6 +45,7 @@
import java.util.Random;
import ffx.utilities.FFXTest;
+import org.junit.After;
import org.junit.Before;
import org.junit.Test;
import org.junit.runner.RunWith;
@@ -101,6 +102,11 @@ public void setUp() {
}
}
+ @After
+ public void after() throws Exception {
+ parallelTeam.shutdown();
+ }
+
/** Test of convolution method, of class Complex3D. */
@Test
public void testConvolution() {
diff --git a/modules/numerics/src/test/java/ffx/numerics/fft/Real3DParallelTest.java b/modules/numerics/src/test/java/ffx/numerics/fft/Real3DParallelTest.java
index 4dc1d0a533..41f204d718 100755
--- a/modules/numerics/src/test/java/ffx/numerics/fft/Real3DParallelTest.java
+++ b/modules/numerics/src/test/java/ffx/numerics/fft/Real3DParallelTest.java
@@ -45,6 +45,7 @@
import java.util.Random;
import ffx.utilities.FFXTest;
+import org.junit.After;
import org.junit.Before;
import org.junit.Test;
import org.junit.runner.RunWith;
@@ -108,6 +109,11 @@ public void setUp() {
}
}
+ @After
+ public void after() throws Exception {
+ parallelTeam.shutdown();
+ }
+
/** Test of convolution method, of class Real3DParallel. */
@Test
public void testConvolution() {
diff --git a/modules/pj/src/main/java/edu/rit/pj/KillRegion.java b/modules/pj/src/main/java/edu/rit/pj/KillRegion.java
index 3a5fac60fe..0db008478b 100644
--- a/modules/pj/src/main/java/edu/rit/pj/KillRegion.java
+++ b/modules/pj/src/main/java/edu/rit/pj/KillRegion.java
@@ -55,9 +55,14 @@
* @author Jacob Litman
*/
public class KillRegion extends ParallelRegion {
+
/** {@inheritDoc} */
@Override
public void run() throws Exception {
- // Does precisely nothing save exist.
+ // Empty run method.
+ ParallelTeamThread currentThread = getCurrentThread();
+ Thread thread = currentThread.getThread();
+ long totalThreads = ParallelTeamThread.totalThreads.getAndDecrement();
+ // System.out.printf(" Killing team %s thread %s of %d\n", currentThread.myTeam, thread, totalThreads);
}
}
diff --git a/modules/pj/src/main/java/edu/rit/pj/PJProperties.java b/modules/pj/src/main/java/edu/rit/pj/PJProperties.java
index 8bfa33536e..fa746b5fd5 100644
--- a/modules/pj/src/main/java/edu/rit/pj/PJProperties.java
+++ b/modules/pj/src/main/java/edu/rit/pj/PJProperties.java
@@ -233,6 +233,26 @@ public static int getPjNt() {
return k;
}
+ /**
+ * Determine whether a {@linkplain ParallelTeam} should use virtual threads.
+ *
+ * If the "pj.vt" Java system property is true, then virtual threads
+ * will be used. If the property is false or not specified, then the
+ * program will use platform threads.
+ *
+ * The default is false, meaning that the program will use platform threads.
+ *
+ * @return true if the program should use virtual threads, false otherwise.
+ */
+ public static boolean getPjVt() {
+ String pjVtProperty = System.getProperty("pj.vt");
+ if (pjVtProperty != null) {
+ return Boolean.parseBoolean(pjVtProperty);
+ } else {
+ return false; // Default is false
+ }
+ }
+
/**
* Determine the schedule for a parallel loop in an SMP parallel program.
* This is the schedule that will be used if the parallel for loop has a
diff --git a/modules/pj/src/main/java/edu/rit/pj/ParallelConstruct.java b/modules/pj/src/main/java/edu/rit/pj/ParallelConstruct.java
index a9df6fb78c..24fb58297c 100644
--- a/modules/pj/src/main/java/edu/rit/pj/ParallelConstruct.java
+++ b/modules/pj/src/main/java/edu/rit/pj/ParallelConstruct.java
@@ -148,16 +148,13 @@ ParallelTeamThread getCurrentThread() {
if (myTeam == null) {
throw new IllegalStateException("ParallelConstruct.getCurrentThread(): No parallel team executing");
}
- try {
- ParallelTeamThread current = (ParallelTeamThread) Thread.currentThread();
- if (current.myTeam != this.myTeam) {
- throw new IllegalStateException("ParallelConstruct.getCurrentThread(): Current thread is not executing this parallel construct");
+ Thread currentThread = Thread.currentThread();
+ for (ParallelTeamThread teamThread : myTeam.myThread) {
+ if (teamThread.getThread() == currentThread) {
+ return teamThread;
}
- return current;
- } catch (ClassCastException exc) {
- throw new IllegalStateException("ParallelConstruct.getCurrentThread(): Current thread is not a parallel team thread",
- exc);
}
+ throw new IllegalStateException("ParallelConstruct.getCurrentThread(): Current thread is not executing this parallel construct");
}
-}
+}
\ No newline at end of file
diff --git a/modules/pj/src/main/java/edu/rit/pj/ParallelTeam.java b/modules/pj/src/main/java/edu/rit/pj/ParallelTeam.java
index 548e92fc57..53cb3d5939 100644
--- a/modules/pj/src/main/java/edu/rit/pj/ParallelTeam.java
+++ b/modules/pj/src/main/java/edu/rit/pj/ParallelTeam.java
@@ -162,7 +162,7 @@ public final void execute(ParallelRegion theRegion)
// Record parallel region.
myRegion = theRegion;
- myExceptionMap = new ConcurrentHashMap(K, 0.75f, K);
+ myExceptionMap = new ConcurrentHashMap<>(K, 0.75f, K);
theRegion.myTeam = this;
try {
diff --git a/modules/pj/src/main/java/edu/rit/pj/ParallelTeamThread.java b/modules/pj/src/main/java/edu/rit/pj/ParallelTeamThread.java
index 5895158449..92d3ab2af3 100644
--- a/modules/pj/src/main/java/edu/rit/pj/ParallelTeamThread.java
+++ b/modules/pj/src/main/java/edu/rit/pj/ParallelTeamThread.java
@@ -44,6 +44,7 @@
//******************************************************************************
package edu.rit.pj;
+import edu.rit.pj.reduction.SharedLong;
import java.util.concurrent.Semaphore;
/**
@@ -54,9 +55,9 @@
* @version 20-Dec-2007
*/
class ParallelTeamThread
- extends Thread {
+ implements Runnable {
-// Hidden data members.
+ // Hidden data members.
// Reference to the parallel team.
ParallelTeam myTeam;
@@ -86,29 +87,64 @@ class ParallelTeamThread
// Exception thrown while setting up a parallel construct, or null if none.
volatile Throwable myConstructException;
+ // The thread instance
+ private final Thread thread;
+
+ /**
+ * Number of threads created by the ParallelTeamThread constructor.
+ */
+ public static SharedLong totalThreads = new SharedLong(0);
+
// 128 bytes of extra padding to avert cache interference.
private long p0, p1, p2, p3, p4, p5, p6, p7;
private long p8, p9, pa, pb, pc, pd, pe, pf;
// Exported constructors.
+
/**
* Construct a new parallel team thread.
*
- * @param theTeam Parallel team to which this thread belongs.
+ * @param theTeam Parallel team to which this thread belongs.
* @param theIndex Index of this thread within the team.
*/
public ParallelTeamThread(ParallelTeam theTeam,
- int theIndex) {
+ int theIndex) {
myTeam = theTeam;
myIndex = theIndex;
- setDaemon(true);
- start();
+
+ String name = "ParallelTeamThread-" + theIndex;
+
+ boolean isVirtual = PJProperties.getPjVt();
+ if (isVirtual) {
+ // Create a virtual thread (all virtual threads are daemon threads).
+ thread = Thread.ofVirtual().name(name).unstarted(this);
+ } else {
+ // Create a new thread for this parallel team thread.
+ thread = Thread.ofPlatform().name(name).unstarted(this);
+ thread.setDaemon(true);
+ }
+
+ long numThreads = totalThreads.incrementAndGet();
+ // System.out.printf(" Creating team %s thread %s out of %d.\n", myTeam, name, numThreads);
+
+ thread.start();
}
// Exported operations.
+
+ /**
+ * Get the underlying thread.
+ *
+ * @return The Thread instance
+ */
+ public Thread getThread() {
+ return thread;
+ }
+
/**
* Run this parallel team thread.
*/
+ @Override
public void run() {
for (;;) {
// Wait until released by the main thread.
@@ -120,6 +156,8 @@ public void run() {
} catch (Exception ex) {
System.err.printf("Error exiting parallel team thread: %s%n", ex);
}
+
+ // Exit the for loop and run method, which will allow threads to be garbage collected.
break;
}
diff --git a/modules/potential/src/main/groovy/ffx/potential/groovy/PrepareSpaceGroups.groovy b/modules/potential/src/main/groovy/ffx/potential/groovy/PrepareSpaceGroups.groovy
index 14c04927f7..cd7d380b01 100644
--- a/modules/potential/src/main/groovy/ffx/potential/groovy/PrepareSpaceGroups.groovy
+++ b/modules/potential/src/main/groovy/ffx/potential/groovy/PrepareSpaceGroups.groovy
@@ -222,8 +222,8 @@ class PrepareSpaceGroups extends PotentialScript {
continue
}
- logger.info(format("\n Preparing %s (CSD percent: %7.4f, PDB Rank: %d)",
- spacegroup2.shortName, getCCDCPercent(num), getPDBRank(spacegroup)))
+ logger.info(format("\n Preparing %d %s (CSD percent: %7.4f, PDB Rank: %d)",
+ num, spacegroup2.shortName, getCCDCPercent(num), getPDBRank(spacegroup)))
// Create the directory.
String sgDirName = spacegroup2.shortName.replace('/', '_')
diff --git a/modules/potential/src/main/java/ffx/potential/extended/ExtendedSystem.java b/modules/potential/src/main/java/ffx/potential/extended/ExtendedSystem.java
index 0db26eff95..ed3300bee6 100644
--- a/modules/potential/src/main/java/ffx/potential/extended/ExtendedSystem.java
+++ b/modules/potential/src/main/java/ffx/potential/extended/ExtendedSystem.java
@@ -149,6 +149,7 @@ public class ExtendedSystem implements Potential {
*/
private final double LYStitrBiasMag;
private final List constraints;
+ public final boolean useTotalChargeCorrection;
private final boolean doBias;
private final boolean doElectrostatics;
private final boolean doPolarization;
@@ -317,6 +318,7 @@ public ExtendedSystem(MolecularAssembly mola, double pH, final File esvFile) {
doElectrostatics = properties.getBoolean("esv.elec", true);
doBias = properties.getBoolean("esv.bias", true);
doPolarization = properties.getBoolean("esv.polarization", true);
+ useTotalChargeCorrection = properties.getBoolean("esv.total.charge.correction", false);
thetaFriction = properties.getDouble("esv.friction", ExtendedSystem.THETA_FRICTION);
thetaMass = properties.getDouble("esv.mass", ExtendedSystem.THETA_MASS);
diff --git a/modules/potential/src/main/java/ffx/potential/nonbonded/ParticleMeshEwald.java b/modules/potential/src/main/java/ffx/potential/nonbonded/ParticleMeshEwald.java
index d0f4ee6a5b..c2622185e2 100755
--- a/modules/potential/src/main/java/ffx/potential/nonbonded/ParticleMeshEwald.java
+++ b/modules/potential/src/main/java/ffx/potential/nonbonded/ParticleMeshEwald.java
@@ -320,6 +320,7 @@ public class ParticleMeshEwald implements LambdaInterface {
protected double permanentRealSpaceEnergy;
protected double permanentSelfEnergy;
protected double permanentReciprocalEnergy;
+ protected double permanentChargeCorrectionEnergy;
/**
* Polarization energy = inducedRealSpaceEnergy + inducedSelfEnergy + inducedReciprocalEnergy.
*/
@@ -435,27 +436,27 @@ public class ParticleMeshEwald implements LambdaInterface {
*/
private AtomicDoubleArrayImpl atomicDoubleArrayImpl;
/**
- * Field array for each thread. [threadID][X/Y/Z][atomID]
+ * Field array for each thread.
*/
private AtomicDoubleArray3D field;
/**
- * Chain rule field array for each thread. [threadID][X/Y/Z][atomID]
+ * Chain rule field array for each thread.
*/
private AtomicDoubleArray3D fieldCR;
/**
- * Gradient array for each thread. [threadID][X/Y/Z][atomID]
+ * Gradient array for each thread.
*/
private AtomicDoubleArray3D grad;
/**
- * Torque array for each thread. [threadID][X/Y/Z][atomID]
+ * Torque array for each thread.
*/
private AtomicDoubleArray3D torque;
/**
- * Partial derivative of the gradient with respect to Lambda. [threadID][X/Y/Z][atomID]
+ * Partial derivative of the gradient with respect to Lambda.
*/
private AtomicDoubleArray3D lambdaGrad;
/**
- * Partial derivative of the torque with respect to Lambda. [threadID][X/Y/Z][atomID]
+ * Partial derivative of the torque with respect to Lambda.
*/
private AtomicDoubleArray3D lambdaTorque;
@@ -827,6 +828,7 @@ public double energy(boolean gradient, boolean print) {
permanentRealSpaceEnergy = 0.0;
permanentSelfEnergy = 0.0;
permanentReciprocalEnergy = 0.0;
+ permanentChargeCorrectionEnergy = 0.0;
polarizationEnergy = 0.0;
inducedRealSpaceEnergy = 0.0;
inducedSelfEnergy = 0.0;
@@ -1079,6 +1081,10 @@ public double getPermRecipEnergy() {
return permanentReciprocalEnergy;
}
+ public double getPermanentChargeCorrectionEnergy() {
+ return permanentChargeCorrectionEnergy;
+ }
+
public double getPolarEps() {
return poleps;
}
@@ -1484,8 +1490,11 @@ private void initSoftCore() {
alchemicalParameters.vaporPermanentSchedule = vacuumNeighborList.getPairwiseSchedule();
alchemicalParameters.vaporEwaldSchedule = alchemicalParameters.vaporPermanentSchedule;
alchemicalParameters.vacuumRanges = new Range[maxThreads];
- vacuumNeighborList.setDisableUpdates(
- forceField.getBoolean("DISABLE_NEIGHBOR_UPDATES", false));
+ try {
+ vacuumNeighborList.destroy();
+ } catch (Exception ex) {
+ logger.warning(" Exception in destroying vacuumNeighborList");
+ }
} else {
alchemicalParameters.vaporCrystal = null;
alchemicalParameters.vaporLists = null;
@@ -1960,6 +1969,7 @@ private double computeEnergy(boolean print) {
double eself = 0.0;
double erecip = 0.0;
+ double ecorrect = 0.0;
double eselfi = 0.0;
double erecipi = 0.0;
if (reciprocalSpaceTerm && ewaldParameters.aewald > 0.0) {
@@ -1973,6 +1983,7 @@ private double computeEnergy(boolean print) {
reciprocalEnergyRegion.executeWith(parallelTeam);
eself = reciprocalEnergyRegion.getPermanentSelfEnergy();
erecip = reciprocalEnergyRegion.getPermanentReciprocalEnergy();
+ ecorrect = reciprocalEnergyRegion.getPermanentChargeCorrectionEnergy();
eselfi = reciprocalEnergyRegion.getInducedDipoleSelfEnergy();
erecipi = reciprocalEnergyRegion.getInducedDipoleReciprocalEnergy();
interactions += nAtoms;
@@ -2131,10 +2142,11 @@ private double computeEnergy(boolean print) {
permanentRealSpaceEnergy += ereal;
permanentSelfEnergy += eself;
permanentReciprocalEnergy += erecip;
+ permanentChargeCorrectionEnergy += ecorrect;
inducedRealSpaceEnergy += ereali;
inducedSelfEnergy += eselfi;
inducedReciprocalEnergy += erecipi;
- permanentMultipoleEnergy += eself + erecip + ereal;
+ permanentMultipoleEnergy += eself + erecip + ereal + ecorrect;
polarizationEnergy += eselfi + erecipi + ereali;
totalMultipoleEnergy += ereal + eself + erecip + ereali + eselfi + erecipi;
@@ -2142,11 +2154,12 @@ private double computeEnergy(boolean print) {
if (logger.isLoggable(Level.FINE)) {
StringBuilder sb = new StringBuilder();
sb.append(format("\n Global Cartesian PME, lambdaMode=%s\n", lambdaMode.toString()));
- sb.append(format(" Multipole Self-Energy: %16.8f\n", eself));
sb.append(format(" Multipole Reciprocal: %16.8f\n", erecip));
+ sb.append(format(" Multipole Self-Energy: %16.8f\n", eself));
+ sb.append(format(" Multipole Correction: %16.8f\n", ecorrect));
sb.append(format(" Multipole Real Space: %16.8f\n", ereal));
- sb.append(format(" Polarization Self-Energy:%16.8f\n", eselfi));
sb.append(format(" Polarization Reciprocal: %16.8f\n", erecipi));
+ sb.append(format(" Polarization Self-Energy:%16.8f\n", eselfi));
sb.append(format(" Polarization Real Space: %16.8f\n", ereali));
if (generalizedKirkwoodTerm) {
sb.append(format(" Generalized Kirkwood: %16.8f\n", solvationEnergy));
diff --git a/modules/potential/src/main/java/ffx/potential/nonbonded/pme/RealSpaceEnergyRegion.java b/modules/potential/src/main/java/ffx/potential/nonbonded/pme/RealSpaceEnergyRegion.java
index a691f0598a..6882908bec 100644
--- a/modules/potential/src/main/java/ffx/potential/nonbonded/pme/RealSpaceEnergyRegion.java
+++ b/modules/potential/src/main/java/ffx/potential/nonbonded/pme/RealSpaceEnergyRegion.java
@@ -267,6 +267,7 @@ public class RealSpaceEnergyRegion extends ParallelRegion implements MaskingInte
* Set polarizationScale to L for part 1. Set polarizationScale to (1-L) for parts 2 & 3.
*/
private double polarizationScale = 1.0;
+
// *************************************************************************
// Mutable Particle Mesh Ewald constants.
private double aewald;
diff --git a/modules/potential/src/main/java/ffx/potential/nonbonded/pme/ReciprocalEnergyRegion.java b/modules/potential/src/main/java/ffx/potential/nonbonded/pme/ReciprocalEnergyRegion.java
index b0e46152f2..9cf2d3375c 100644
--- a/modules/potential/src/main/java/ffx/potential/nonbonded/pme/ReciprocalEnergyRegion.java
+++ b/modules/potential/src/main/java/ffx/potential/nonbonded/pme/ReciprocalEnergyRegion.java
@@ -57,6 +57,8 @@
import static ffx.potential.parameters.MultipoleType.t201;
import static ffx.potential.parameters.MultipoleType.t210;
import static ffx.potential.parameters.MultipoleType.t300;
+import static java.lang.Math.PI;
+import static java.lang.String.format;
import static org.apache.commons.math3.util.FastMath.sqrt;
import edu.rit.pj.IntegerForLoop;
@@ -83,9 +85,10 @@ public class ReciprocalEnergyRegion extends ParallelRegion {
private static final Logger logger = Logger.getLogger(ReciprocalEnergyRegion.class.getName());
- private static final double SQRT_PI = sqrt(Math.PI);
+ private static final double SQRT_PI = sqrt(PI);
private static final int tensorCount = MultipoleUtilities.tensorCount(3);
private final double electric;
+ private final double aewald;
private final double aewald1;
private final double aewald2;
private final double aewald3;
@@ -110,6 +113,7 @@ public class ReciprocalEnergyRegion extends ParallelRegion {
private double[][] fracMultipolePhi;
private double[][] fracInducedDipolePhi;
private double[][] fracInducedDipolePhiCR;
+ private double chargeCorrectionEnergy;
private double permanentSelfEnergy;
private double permanentReciprocalEnergy;
/** An ordered array of atoms in the system. */
@@ -182,8 +186,8 @@ public ReciprocalEnergyRegion(int nt, double aewald, double electric) {
inducedDipoleSelfEnergy = new SharedDouble();
inducedDipoleRecipEnergy = new SharedDouble();
maxThreads = nt;
-
this.electric = electric;
+ this.aewald = aewald;
aewald1 = -electric * aewald / SQRT_PI;
aewald2 = 2.0 * aewald * aewald;
aewald3 = -2.0 / 3.0 * electric * aewald * aewald * aewald / SQRT_PI;
@@ -215,10 +219,40 @@ public void finish() {
*/
permanentSelfEnergy = 0.0;
permanentReciprocalEnergy = 0.0;
+ double totalCharge = 0.0;
for (int i = 0; i < maxThreads; i++) {
+ totalCharge += permanentReciprocalEnergyLoop[i].totalCharge;
permanentSelfEnergy += permanentReciprocalEnergyLoop[i].eSelf;
permanentReciprocalEnergy += permanentReciprocalEnergyLoop[i].eRecip;
}
+
+ if (totalCharge == 0.0 || (esvTerm && !extendedSystem.useTotalChargeCorrection)) {
+ chargeCorrectionEnergy = 0.0;
+ } else {
+ Crystal unitCell = crystal.getUnitCell();
+ int nSymm = unitCell.getNumSymOps();
+ // Compute the total charge for the unit cell from the asymmetric unit total.
+ totalCharge = totalCharge * nSymm;
+ double denom = unitCell.volume * aewald * aewald;
+ chargeCorrectionEnergy = -0.5 * electric * PI * (totalCharge * totalCharge) / denom;
+ // Normalize the unit cell correction to the asymmetric unit.
+ chargeCorrectionEnergy = chargeCorrectionEnergy / nSymm;
+ if (lambdaTerm) {
+ shareddEdLambda.addAndGet(chargeCorrectionEnergy * dlPowPerm * dEdLSign);
+ sharedd2EdLambda2.addAndGet(chargeCorrectionEnergy * d2lPowPerm * dEdLSign);
+ }
+ if (esvTerm){
+ for(int i=0; i= startingSnap + 1 && count <= endingSnap + 1) {
if (count <= snaps1) {
@@ -174,7 +174,7 @@ public boolean readFile() {
} else {
logger.warning(format(" BAR entry of (%3d) is larger than total entries (%3d).", count, snaps1));
}
- } else if (count >= snaps1 + startingSnap + 2 && count <= snaps1 + endingSnap + 2) {
+ } else if (count >= snaps1 + startingSnap + 1 && count <= snaps1 + endingSnap + 1) {
if (count <= snaps1 + snaps2 + 1) {
if (numTokens == 4) {
vol2.add(parseDouble(tokens[3]));
@@ -186,13 +186,14 @@ public boolean readFile() {
}
}
} else {
- if (count <= snaps1 + 1 && count != 1) {
+ count++;
+ if (count <= snaps1) {
if (numTokens == 4) {
vol1.add(parseDouble(tokens[3]));
}
ens1lam1.add(parseDouble(tokens[1]));
ens1lam2.add(parseDouble(tokens[2]));
- } else if (count > snaps1 + 2) {
+ } else {
if (numTokens == 4) {
vol2.add(parseDouble(tokens[3]));
}
diff --git a/modules/potential/src/test/java/ffx/potential/groovy/EnergyTest.java b/modules/potential/src/test/java/ffx/potential/groovy/EnergyTest.java
index b96038bd26..24043c3346 100644
--- a/modules/potential/src/test/java/ffx/potential/groovy/EnergyTest.java
+++ b/modules/potential/src/test/java/ffx/potential/groovy/EnergyTest.java
@@ -99,7 +99,7 @@ public static Collection data() {
0,
4003.27183547,
741643,
- -12303.93157019,
+ -12305.43400643136,
332114,
-2811.45683671,
332114,
@@ -136,7 +136,7 @@ public static Collection data() {
0,
16013.08734188,
2966572,
- -49215.72628076,
+ -49221.7360257262,
1328456,
-11245.82734685,
1328456,
@@ -172,7 +172,7 @@ public static Collection data() {
0,
31718.84803063,
3480619,
- -78751.36662993,
+ -78752.17402965919,
1463589,
-31790.15161303,
1463589,
@@ -388,7 +388,7 @@ public static Collection data() {
0,
12862.34983956,
1482946,
- -32737.08497141565, // or -32736.92207773 (with no FFT factors of 6 & 7)
+ -32737.101745561464,
623627,
-12934.94219259, // or -12934.93512829 (with no FFT factors of 6 & 7)
623627,
@@ -424,7 +424,7 @@ public static Collection data() {
0,
4860.41072437,
2279196,
- -40036.04804903,
+ -40036.07471652395,
3562034,
0.0,
3562034,
diff --git a/pom.xml b/pom.xml
index 238bdd50e1..2e41f14284 100644
--- a/pom.xml
+++ b/pom.xml
@@ -922,7 +922,6 @@
-5
-
https://ffx.biochem.uiowa.edu/download.html
@@ -1082,7 +1081,7 @@
true
3.4.2
3.1.0
- 3.26.4
+ 3.27.0
1.7.1
3.11.2
4.0.2