From c2ed63d25da2b5430b5c3270b625692260752b73 Mon Sep 17 00:00:00 2001 From: Jacob M Miller Date: Wed, 11 Jun 2025 16:41:00 -0500 Subject: [PATCH 1/6] Fixed BAR snapshot selection. --- .../src/main/groovy/ffx/algorithms/groovy/BAR.groovy | 4 ++-- .../src/main/java/ffx/potential/parsers/BARFilter.java | 9 +++++---- 2 files changed, 7 insertions(+), 6 deletions(-) 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/potential/src/main/java/ffx/potential/parsers/BARFilter.java b/modules/potential/src/main/java/ffx/potential/parsers/BARFilter.java index 2fea799c6b..d71857f002 100644 --- a/modules/potential/src/main/java/ffx/potential/parsers/BARFilter.java +++ b/modules/potential/src/main/java/ffx/potential/parsers/BARFilter.java @@ -150,7 +150,6 @@ public boolean readFile() { try (BufferedReader br = new BufferedReader(new FileReader(barFile))) { String data; while ((data = br.readLine()) != null) { - count++; String[] tokens = data.trim().split(" +"); int numTokens = tokens.length; if (data.contains(".xyz") || data.contains(".pdb") || numTokens < 3) { @@ -163,6 +162,7 @@ public boolean readFile() { temperature2 = parseDouble(tokens[1]); } } else if (endingSnap != 0) { + count++; snapshots = (endingSnap - startingSnap) + 1; if (count >= 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])); } From de9b47ab7b006ecc3359e8e0db55e67817b203d4 Mon Sep 17 00:00:00 2001 From: "Michael J. Schnieders" Date: Thu, 12 Jun 2025 11:57:44 -0500 Subject: [PATCH 2/6] Add the pj.vt property to turn on the use of virtual threads instead of platform threads. --- .../src/main/java/edu/rit/pj/KillRegion.java | 2 +- .../main/java/edu/rit/pj/PJProperties.java | 20 ++++++++++ .../java/edu/rit/pj/ParallelConstruct.java | 15 +++----- .../java/edu/rit/pj/ParallelTeamThread.java | 38 ++++++++++++++++--- 4 files changed, 59 insertions(+), 16 deletions(-) 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..05bbf56752 100644 --- a/modules/pj/src/main/java/edu/rit/pj/KillRegion.java +++ b/modules/pj/src/main/java/edu/rit/pj/KillRegion.java @@ -58,6 +58,6 @@ public class KillRegion extends ParallelRegion { /** {@inheritDoc} */ @Override public void run() throws Exception { - // Does precisely nothing save exist. + // Empty run method. } } 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/ParallelTeamThread.java b/modules/pj/src/main/java/edu/rit/pj/ParallelTeamThread.java index 5895158449..9909e22add 100644 --- a/modules/pj/src/main/java/edu/rit/pj/ParallelTeamThread.java +++ b/modules/pj/src/main/java/edu/rit/pj/ParallelTeamThread.java @@ -54,9 +54,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,26 +86,52 @@ class ParallelTeamThread // Exception thrown while setting up a parallel construct, or null if none. volatile Throwable myConstructException; + // The thread instance + private final Thread thread; + // 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(); + + boolean isVirtual = PJProperties.getPjVt(); + if (isVirtual) { + // Create a virtual thread + thread = Thread.ofVirtual() + .name("ParallelTeamThread-" + theIndex) + .unstarted(this); + } else { + // Create a new thread for this parallel team thread. + thread = new Thread(this, "ParallelTeamThread-" + theIndex); + } + + thread.setDaemon(true); + thread.start(); } // Exported operations. + + /** + * Get the underlying thread. + * + * @return The Thread instance + */ + public Thread getThread() { + return thread; + } + /** * Run this parallel team thread. */ From 1967e8d29555097694ec524ec711faf688736c68 Mon Sep 17 00:00:00 2001 From: "Michael J. Schnieders" Date: Sat, 14 Jun 2025 15:04:38 -0500 Subject: [PATCH 3/6] Shut down ParallelTeam instances in Complex3DParallelTest and Real3DParallelTest; Shut down the ParallelTeam that creates the vacuum N^2 neighbor-list in ParticleMeshEwald --- .../numerics/fft/Complex3DParallelTest.java | 6 +++++ .../ffx/numerics/fft/Real3DParallelTest.java | 6 +++++ .../src/main/java/edu/rit/pj/KillRegion.java | 5 ++++ .../main/java/edu/rit/pj/ParallelTeam.java | 2 +- .../java/edu/rit/pj/ParallelTeamThread.java | 24 ++++++++++++++----- .../groovy/PrepareSpaceGroups.groovy | 4 ++-- .../nonbonded/ParticleMeshEwald.java | 7 ++++-- pom.xml | 1 - 8 files changed, 43 insertions(+), 12 deletions(-) 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 05bbf56752..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 { // 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/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 9909e22add..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; /** @@ -89,6 +90,11 @@ class ParallelTeamThread // 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; @@ -106,18 +112,21 @@ public ParallelTeamThread(ParallelTeam theTeam, myTeam = theTeam; myIndex = theIndex; + String name = "ParallelTeamThread-" + theIndex; + boolean isVirtual = PJProperties.getPjVt(); if (isVirtual) { - // Create a virtual thread - thread = Thread.ofVirtual() - .name("ParallelTeamThread-" + theIndex) - .unstarted(this); + // 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 = new Thread(this, "ParallelTeamThread-" + theIndex); + thread = Thread.ofPlatform().name(name).unstarted(this); + thread.setDaemon(true); } - thread.setDaemon(true); + long numThreads = totalThreads.incrementAndGet(); + // System.out.printf(" Creating team %s thread %s out of %d.\n", myTeam, name, numThreads); + thread.start(); } @@ -135,6 +144,7 @@ public Thread getThread() { /** * Run this parallel team thread. */ + @Override public void run() { for (;;) { // Wait until released by the main thread. @@ -146,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/nonbonded/ParticleMeshEwald.java b/modules/potential/src/main/java/ffx/potential/nonbonded/ParticleMeshEwald.java index d0f4ee6a5b..c73dd3536c 100755 --- a/modules/potential/src/main/java/ffx/potential/nonbonded/ParticleMeshEwald.java +++ b/modules/potential/src/main/java/ffx/potential/nonbonded/ParticleMeshEwald.java @@ -1484,8 +1484,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; diff --git a/pom.xml b/pom.xml index 238bdd50e1..fa8453ddd7 100644 --- a/pom.xml +++ b/pom.xml @@ -922,7 +922,6 @@ -5 - https://ffx.biochem.uiowa.edu/download.html From db86e41842854c567ff6383553e8238b7682f77f Mon Sep 17 00:00:00 2001 From: "Michael J. Schnieders" Date: Mon, 16 Jun 2025 14:53:37 -0500 Subject: [PATCH 4/6] Added the uniform background charge correction term to PME --- ipynb-kotlin/ffx.json | 1 + .../ffx/numerics/fft/Complex3DParallel.java | 34 +++++++++++------ .../nonbonded/ParticleMeshEwald.java | 28 +++++++++----- .../nonbonded/pme/RealSpaceEnergyRegion.java | 1 + .../nonbonded/pme/ReciprocalEnergyRegion.java | 37 +++++++++++++++++-- .../java/ffx/potential/groovy/EnergyTest.java | 10 ++--- pom.xml | 2 +- 7 files changed, 84 insertions(+), 29 deletions(-) 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/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/potential/src/main/java/ffx/potential/nonbonded/ParticleMeshEwald.java b/modules/potential/src/main/java/ffx/potential/nonbonded/ParticleMeshEwald.java index c73dd3536c..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; } @@ -1963,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) { @@ -1976,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; @@ -2134,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; @@ -2145,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..9afd2669b7 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,30 @@ 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) { + 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); + } + chargeCorrectionEnergy = chargeCorrectionEnergy * permanentScale; + } } public double getInducedDipoleReciprocalEnergy() { @@ -237,6 +261,10 @@ public double getPermanentSelfEnergy() { return permanentSelfEnergy; } + public double getPermanentChargeCorrectionEnergy() { + return chargeCorrectionEnergy; + } + public void init( Atom[] atoms, Crystal crystal, @@ -343,6 +371,7 @@ public void start() { private class PermanentReciprocalEnergyLoop extends IntegerForLoop { int threadID; + double totalCharge; double eSelf; double eRecip; @@ -353,12 +382,13 @@ public void finish() { } @Override - public void run(int lb, int ub) throws Exception { + public void run(int lb, int ub) { // Permanent multipole self energy and gradient. for (int i = lb; i <= ub; i++) { if (use[i]) { double[] in = globalMultipole[0][i]; + totalCharge += in[t000]; double cii = in[t000] * in[t000]; double dii = in[t100] * in[t100] + in[t010] * in[t010] + in[t001] * in[t001]; double qii = in[t200] * in[t200] + in[t020] * in[t020] + in[t002] * in[t002] @@ -495,6 +525,7 @@ public IntegerSchedule schedule() { @Override public void start() { + totalCharge = 0.0; eSelf = 0.0; eRecip = 0.0; threadID = getThreadIndex(); 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 fa8453ddd7..2e41f14284 100644 --- a/pom.xml +++ b/pom.xml @@ -1081,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 From d693f32ee8674d572e5bef7e9c17f4dc3d713a69 Mon Sep 17 00:00:00 2001 From: acthiel212 Date: Wed, 18 Jun 2025 14:51:06 -0500 Subject: [PATCH 5/6] Added derviaties of the volume Ewald correction for titrating systems. --- .../nonbonded/pme/ReciprocalEnergyRegion.java | 10 ++++++++++ 1 file changed, 10 insertions(+) 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 9afd2669b7..b76cdb61ac 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 @@ -241,6 +241,16 @@ public void finish() { shareddEdLambda.addAndGet(chargeCorrectionEnergy * dlPowPerm * dEdLSign); sharedd2EdLambda2.addAndGet(chargeCorrectionEnergy * d2lPowPerm * dEdLSign); } + if (esvTerm){ + for(int i=0; i Date: Thu, 19 Jun 2025 14:56:15 -0500 Subject: [PATCH 6/6] Added conditional evaluation to total charge correction to check if it should calculate the correction when an ExtendedSystem is found. The default behavior is not to calculate it when the ExtendedSystem is present. --- .../src/main/java/ffx/potential/extended/ExtendedSystem.java | 2 ++ .../ffx/potential/nonbonded/pme/ReciprocalEnergyRegion.java | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) 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/pme/ReciprocalEnergyRegion.java b/modules/potential/src/main/java/ffx/potential/nonbonded/pme/ReciprocalEnergyRegion.java index b76cdb61ac..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 @@ -226,7 +226,7 @@ public void finish() { permanentReciprocalEnergy += permanentReciprocalEnergyLoop[i].eRecip; } - if (totalCharge == 0.0) { + if (totalCharge == 0.0 || (esvTerm && !extendedSystem.useTotalChargeCorrection)) { chargeCorrectionEnergy = 0.0; } else { Crystal unitCell = crystal.getUnitCell();