diff --git a/modules/potential/src/main/groovy/ffx/potential/groovy/BFactors.groovy b/modules/potential/src/main/groovy/ffx/potential/groovy/BFactors.groovy new file mode 100644 index 0000000000..a3e83e554a --- /dev/null +++ b/modules/potential/src/main/groovy/ffx/potential/groovy/BFactors.groovy @@ -0,0 +1,715 @@ +//****************************************************************************** +// +// Title: Force Field X. +// Description: Force Field X - Software for Molecular Biophysics. +// Copyright: Copyright (c) Michael J. Schnieders 2001-2025. +// +// This file is part of Force Field X. +// +// Force Field X is free software; you can redistribute it and/or modify it +// under the terms of the GNU General Public License version 3 as published by +// the Free Software Foundation. +// +// Force Field X is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +// details. +// +// You should have received a copy of the GNU General Public License along with +// Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple +// Place, Suite 330, Boston, MA 02111-1307 USA +// +// Linking this library statically or dynamically with other modules is making a +// combined work based on this library. Thus, the terms and conditions of the +// GNU General Public License cover the whole combination. +// +// As a special exception, the copyright holders of this library give you +// permission to link this library with independent modules to produce an +// executable, regardless of the license terms of these independent modules, and +// to copy and distribute the resulting executable under terms of your choice, +// provided that you also meet, for each linked independent module, the terms +// and conditions of the license of that module. An independent module is a +// module which is not derived from or based on this library. If you modify this +// library, you may extend this exception to your version of the library, but +// you are not obligated to do so. If you do not wish to do so, delete this +// exception statement from your version. +// +//****************************************************************************** +package ffx.potential.groovy + +import ffx.potential.bonded.Atom +import ffx.potential.bonded.Residue +import ffx.potential.cli.PotentialScript +import ffx.potential.utils.GetProteinFeatures +import picocli.CommandLine +import picocli.CommandLine.Command +import picocli.CommandLine.Parameters +import picocli.CommandLine.Option +import org.apache.commons.io.FilenameUtils + +@Command(description = " Extract B-factors/confidence scores and structural features from PDB files", name = "BFactors") +class BFactors extends PotentialScript { + + @Option(names = ["-c", "--csv"], paramLabel = "output.csv", + description = "Optional CSV file to append B-factors to. If not provided, a new CSV will be created.") + private String csvFilename = null + + @Option(names = ["-d", "--delimiter"], paramLabel = ",", + description = "Delimiter for CSV file (default is comma)") + private String delimiter = "," + + @Option(names = ["-b", "--batch"], paramLabel = "false", + description = "Process a directory of PDB files and output features to PDB_Features.csv") + private boolean batchMode = false + + @Option(names = ["-m", "--maestro"], paramLabel = "false", + description = "Calculate Instability Heat Score (IHS) using MAESTRO predictions") + private boolean calculateIHS = false + + @Option(names = ["-s", "--stabilityDir"], paramLabel = "stabilityPredictions", + description = "Directory containing MAESTRO prediction files (.maestro.txt)") + private String stabilityDir = "stabilityPredictions" + + @Parameters(arity = "1", paramLabel = "file", + description = 'The atomic coordinate file in PDB format or directory containing PDB files') + private String filename = null + + BFactors() { this(new Binding()) } + BFactors(Binding binding) { super(binding) } + + @Override + BFactors run() { + if (!init()) return null + + if (batchMode) { + processDirectory() + } else { + processSingleFile() + } + + return this + } + + private void processSingleFile() { + activeAssembly = getActiveAssembly(filename) + if (activeAssembly == null) { + logger.info(helpString()) + return + } + + List residues = activeAssembly.getResidueList() + List> bFactorData = collectBFactors(residues) + writeBFactorsToCSV(bFactorData) + + // Calculate and print median B-factor + printMedianBFactor(bFactorData) + + // Add IHS calculation for single file + if (calculateIHS) { + String pdbId = FilenameUtils.getBaseName(filename) + double ihs = calculateIHS(residues, pdbId) + logger.info("\nFinal IHS for ${pdbId}: ${String.format("%.2f", ihs)}\n") + + // You could also write this to a CSV if needed + Map result = [ + pdbId: pdbId, + medianBfactor: String.format("%.2f", calculateMedian(bFactorData)), + totalSASA: "0.0", // Would need to calculate SASA as done in processDirectory() + numResidues: residues.size(), + IHS: String.format("%.2f", ihs) + ] + writeExtendedResultsToCSV([result]) + } + } + + private void processDirectory() { + File directory = new File(filename) + if (!directory.isDirectory()) { + logger.severe("The provided path is not a directory: ${filename}") + return + } + + List> fileResults = [] + List> geneHeatScores = [] + File[] pdbFiles = directory.listFiles({ File f -> f.name.toLowerCase().endsWith('.pdb') } as FileFilter) + + if (pdbFiles == null || pdbFiles.length == 0) { + logger.severe("No PDB files found in directory: ${filename}") + return + } + + // Create PDB to gene mapping + Map pdbToGene = createPdbToGeneMap() + + logger.info("\nProcessing ${pdbFiles.length} PDB files in directory: ${filename}") + int processedCount = 0 + int errorCount = 0 + + // Check stability directory exists if IHS calculation is requested + File stabilityDirectory = null + if (calculateIHS) { + stabilityDirectory = new File(stabilityDir) + if (!stabilityDirectory.exists() || !stabilityDirectory.isDirectory()) { + logger.severe("Stability predictions directory not found: ${stabilityDir}") + return + } + } + + pdbFiles.each { pdbFile -> + try { + logger.info("Processing file: ${pdbFile.name}") + + // Create a new binding with required properties + Binding fileBinding = new Binding() + fileBinding.setVariable("args", [pdbFile.absolutePath] as String[]) + fileBinding.setVariable("csvFilename", csvFilename) + fileBinding.setVariable("delimiter", delimiter) + fileBinding.setVariable("batchMode", batchMode) + fileBinding.setVariable("calculateIHS", calculateIHS) + fileBinding.setVariable("stabilityDir", stabilityDir) + + // Create new script instance with initialized binding + BFactors fileProcessor = new BFactors(fileBinding) + fileProcessor.run() + + // Check if assembly was created successfully + if (fileProcessor.activeAssembly == null) { + logger.warning("Could not create assembly for file: ${pdbFile.name}") + errorCount++ + return // continues to next file + } + + // Get structural features + List residues = fileProcessor.activeAssembly.residueList + List> bFactorData = fileProcessor.collectBFactors(residues) + + if (bFactorData.isEmpty()) { + logger.warning("No B-factors collected from file: ${pdbFile.name}") + errorCount++ + return + } + + // Initialize GetProteinFeatures and calculate SASA + GetProteinFeatures proteinFeatures = new GetProteinFeatures() + residues.each { residue -> + double residueSA = GetProteinFeatures.standardSurfaceArea.getOrDefault( + residue.getAminoAcid3(), 0.0) + proteinFeatures.saveFeatures(residue, residueSA, false, false) + } + double totalSASA = proteinFeatures.getTotalSurfaceArea() + + // Calculate median B-factor + double median = calculateMedian(bFactorData) + String pdbId = FilenameUtils.getBaseName(pdbFile.name) + + // Calculate IHS if requested + double ihs = 0.0 + if (calculateIHS) { + ihs = calculateIHS(residues, pdbId) + logger.info("\nFinal IHS for ${pdbId}: ${String.format("%.2f", ihs)}\n") + } + + // Collect all results + Map result = [ + pdbId: pdbId, + medianBfactor: String.format("%.2f", median), + totalSASA: String.format("%.2f", totalSASA), + numResidues: residues.size() + ] + + if (calculateIHS) { + result["IHS"] = String.format("%.2f", ihs) + } + + // After calculating IHS, add to geneHeatScores if available + if (calculateIHS) { + String cleanPdbId = FilenameUtils.getBaseName(pdbFile.name).toUpperCase() + String geneName = pdbToGene.get(cleanPdbId, "Unknown") + + geneHeatScores << [ + pdbId: cleanPdbId, + geneName: geneName, + heatScore: String.format("%.2f", ihs) + ] + } + + fileResults << result + + logger.info(String.format(" Processed %s - Median B-factor: %.2f, Total SASA: %.2f, Residues: %d", + pdbId, median, totalSASA, residues.size())) + processedCount++ + + } catch (Exception e) { + logger.warning("Error processing file ${pdbFile.name}: ${e.message}") + errorCount++ + } + } + + writeExtendedResultsToCSV(fileResults) + // Write gene heat scores to separate CSV if IHS was calculated + if (calculateIHS && !geneHeatScores.isEmpty()) { + writeGeneHeatScoresToCSV(geneHeatScores) + } + + logger.info("\nProcessing complete:") + logger.info(String.format(" Successfully processed: %d files", processedCount)) + logger.info(String.format(" Files with errors: %d", errorCount)) + logger.info(" Results written to PDB_Features.csv") + if (calculateIHS) { + logger.info(" Gene heat scores written to Gene_Heat_Scores.csv") + } + } + + private void writeGeneHeatScoresToCSV(List> geneHeatScores) { + String outputFilename = "Gene_Heat_Scores.csv" + + try { + File outputFile = new File(outputFilename) + outputFile.withWriter { writer -> + // Write header + writer.write("PDB_ID${delimiter}Gene_Name${delimiter}Heat_Score\n") + + // Write data sorted by heat score (descending) + geneHeatScores.sort { -it.heatScore.toDouble() }.each { result -> + writer.write("${result.pdbId}${delimiter}${result.geneName}${delimiter}${result.heatScore}\n") + } + } + logger.info("\nSuccessfully wrote gene heat scores to ${outputFilename}") + } catch (Exception e) { + logger.severe("Error writing gene heat scores to CSV: ${e.message}") + } + } + private Map createPdbToGeneMap() { + return [ + "1WKU": "ACTN3", + "2R0O": "ACTN4", + "3MY0": "ACVRL1", + "3PH9": "AGR3", + "2BP1": "AKR7A2", + "2X18": "AKT3", + "1CPU": "AMY2A", + "3L81": "AP4M1", + "1E3G": "AR", + "5NP0": "ATM", + "2HLQ": "BMPR2", + "1K2P": "BTK", + "5C9J": "CD1C", + "4AAA": "CDKL2", + "1KWM": "CPB1", + "2JDG": "CRYGB", + "1YB5": "CRYZ", + "1FH0": "CTSV", + "2B5L": "CUL4B", + "2I4I": "DDX3X", + "3ZFM": "EPHB2", + "5HYN": "EZH2", + "3ZGQ": "FIT5", + "2C6Y": "FOXK2", + "2A2C": "GALK2", + "1F5N": "GBP1", + "1V4S": "GCK", + "1A22": "GH1", + "2QKH": "GIPR", + "5ER7": "GJB2", + "2H8R": "HNF1B", + "1T0P": "ICAM3", + "2M1W": "IL18RAP", + "4UXT": "KIF5A", + "3KXZ": "LCK", + "1N7D": "LDLR", + "2P4E": "LRP8", + "5TJA": "MCOLN1", + "3U84": "MEN1", + "4XI6": "MIB1", + "1JAP": "MMP8", + "3THX": "MSH3", + "4PA0": "MYH7", + "1H4R": "NF2", + "1IKN": "NFKBIA", + "5UHV": "NRAS", + "4FZV": "NSUN2", + "3Q91": "NUDT14", + "4NHX": "OGFOD1", + "2DEW": "PADI4", + "1DMW": "PAH", + "1PDV": "PARK7", + "3ELO": "PLA3G1B", + "3IKM": "POLG", + "4OR9": "PPP3CB", + "5BZZ": "PTEN", + "2SHP": "PTPN11", + "4PJU": "RAD21", + "1XAP": "RARB", + "2F8X": "RBPJ", + "4MTH": "REG3A", + "6GXZ": "RPAP3", + "4D9T": "RPS6KA3", + "4UYB": "SEC14L3", + "1ANT": "SERPINC1", + "4MHX": "SGSH", + "4PYP": "SLC2A1", + "3VFD": "SPAST", + "1YVL": "STAT1", + "3BFX": "SULT1C2", + "4LG9": "TBL1XR1", + "3KAS": "TFRC", + "2VFJ": "TNFAIP3", + "4YFI": "TNNI3K", + "1QBK": "TNPO1", + "1A36": "TOP1", + "1UOL": "TP53", + "1FLK": "TRAF3", + "4OLI": "TYK2", + "1UOU": "TYMP", + "1Y8Q": "UBA2", + "3HU3": "VCP", + "4NS5": "ZMYND11" + ] + } + + /** + * Map PDB ID to UniProt ID using UniProt's REST API + */ + private List mapPdbToUniprot(String pdbId) { + logger.info("Attempting to map PDB ${pdbId} to UniProt IDs using REST API...") + List uniprotIds = [] + + try { + // Step 1: Submit ID mapping job + def submitUrl = "https://rest.uniprot.org/idmapping/run" + def submitCmd = ["curl", "--request", "POST", submitUrl, + "--form", "ids=${pdbId}", + "--form", "from=PDB", + "--form", "to=UniProtKB"] + + logger.info("Executing: ${submitCmd.join(' ')}") + def submitProcess = submitCmd.execute() + submitProcess.waitFor() + + if (submitProcess.exitValue() != 0) { + logger.warning("Submit command failed with exit code ${submitProcess.exitValue()}") + logger.warning("Error output: ${submitProcess.err.text}") + return uniprotIds + } + + def submitResponse = submitProcess.text + logger.info("Submit response: ${submitResponse}") + + // Parse job ID from response + def jobId = (submitResponse =~ /"jobId":"([^"]+)"/)?.with { it[0][1] } + if (!jobId) { + logger.warning("Could not extract job ID from response") + return uniprotIds + } + logger.info("Job ID: ${jobId}") + + // Step 2: Poll for job completion + def statusUrl = "https://rest.uniprot.org/idmapping/status/${jobId}" + int attempts = 0 + int maxAttempts = 10 + boolean jobComplete = false + + while (attempts < maxAttempts && !jobComplete) { + sleep(3000) // Wait 3 seconds between checks + attempts++ + + def statusCmd = ["curl", "-s", statusUrl] + def statusProcess = statusCmd.execute() + statusProcess.waitFor() + + if (statusProcess.exitValue() != 0) { + logger.warning("Status check failed for attempt ${attempts}") + continue + } + + def statusResponse = statusProcess.text + if (statusResponse.contains('"jobStatus":"FINISHED"')) { + jobComplete = true + } + } + + if (!jobComplete) { + logger.warning("Job did not complete within ${maxAttempts} attempts") + return uniprotIds + } + + // Step 3: Get results + def resultsUrl = "https://rest.uniprot.org/idmapping/uniprotkb/results/stream/${jobId}?format=tsv" + def resultsCmd = ["curl", "-s", resultsUrl] + logger.info("Executing: ${resultsCmd.join(' ')}") + + def resultsProcess = resultsCmd.execute() + resultsProcess.waitFor() + + if (resultsProcess.exitValue() != 0) { + logger.warning("Results command failed with exit code ${resultsProcess.exitValue()}") + logger.warning("Error output: ${resultsProcess.err.text}") + return uniprotIds + } + + def results = resultsProcess.text + logger.info("Results:\n${results}") + + // Parse TSV results - collect all UniProt IDs + def lines = results.readLines() + if (lines.size() > 1) { + lines[1..-1].each { line -> + def parts = line.split("\t") + if (parts.size() >= 2) { + String uniprotId = parts[1].trim() + if (uniprotId && !uniprotIds.contains(uniprotId)) { + uniprotIds.add(uniprotId) + } + } + } + } + + if (uniprotIds) { + logger.info("Successfully mapped to UniProt IDs: ${uniprotIds.join(', ')}") + } else { + logger.warning("No UniProt IDs found in results") + } + + } catch (Exception e) { + logger.warning("Error mapping PDB to UniProt: ${e.message}") + e.printStackTrace() + } + + return uniprotIds + } + + /** + * Find MAESTRO file using UniProt ID pattern with detailed logging + */ + private File findMaestroFile(String pdbId, File maestroDir) { + logger.info("\nSearching for MAESTRO file matching PDB: ${pdbId}") + + // First try direct PDB ID match + logger.info("Attempting direct PDB ID match...") + File directMatch = maestroDir.listFiles().find { + it.name.startsWith(pdbId) && it.name.endsWith(".maestro.txt") + } + if (directMatch) { + logger.info("Found direct match: ${directMatch.name}") + return directMatch + } + + // If no direct match, try UniProt mapping + logger.info("No direct match found, attempting UniProt mapping...") + List uniprotIds = mapPdbToUniprot(pdbId) + if (!uniprotIds) { + logger.warning("Could not map PDB ID ${pdbId} to any UniProt IDs") + return null + } + logger.info("Mapped ${pdbId} to UniProt IDs: ${uniprotIds.join(', ')}") + + // Try all possible patterns for each UniProt ID + for (String uniprotId : uniprotIds) { + logger.info("Searching for files matching UniProt ID: ${uniprotId}") + + // Define patterns for this specific uniprotId + def patterns = [ + ~/(?i)AF-${uniprotId}-F1.*\.maestro\.txt/, + ~/(?i)${uniprotId}\.maestro\.txt/, + ~/(?i)${uniprotId}-.*\.maestro\.txt/ + ] + + // Try all pattern variations + for (pattern in patterns) { + def potentialFiles = maestroDir.listFiles().findAll { file -> + file.name ==~ pattern + } + + if (potentialFiles) { + logger.info("Found ${potentialFiles.size()} matching files with pattern ${pattern}:") + potentialFiles.each { f -> logger.info(" ${f.name}") } + return potentialFiles[0] // Return first match + } + } + } + + logger.warning("No MAESTRO files found matching any UniProt patterns") + return null + } + + /** + * Calculate Instability Heat Score (IHS) for a protein structure using MAESTRO predictions + * + * @param residues List of residues in the protein + * @param pdbId The PDB ID of the protein + * @return Calculated IHS value + */ + private double calculateIHS(List residues, String pdbId) { + logger.info("\nStarting IHS calculation for PDB: ${pdbId}") + + // Clean PDB ID (remove .pdb extension if present) + String cleanPdbId = pdbId.replaceAll(/(?i)\.pdb$/, "") + logger.info("Cleaned PDB ID: ${cleanPdbId}") + + File maestroDir = new File(stabilityDir) + if (!maestroDir.exists()) { + logger.warning("Stability directory not found: ${maestroDir.absolutePath}") + return 0.0 + } + + File maestroFile = findMaestroFile(cleanPdbId, maestroDir) + if (!maestroFile) { + logger.warning("No matching MAESTRO file found for ${cleanPdbId}") + logger.warning("Searched in: ${maestroDir.absolutePath}") + return 0.0 + } + logger.info("Found MAESTRO file: ${maestroFile.absolutePath}") + + List significantDDGs = [] + int totalMutations = 0 + int linesProcessed = 0 + + // Parse the MAESTRO file + try { + logger.info("Parsing MAESTRO file...") + maestroFile.eachLine { line -> + linesProcessed++ + if (!line.startsWith("#") && !line.trim().isEmpty()) { + String[] parts = line.trim().split("\\s+") + if (parts.size() >= 7) { // Ensure we have all expected columns + String mutation = parts[3] + if (mutation != "wildtype") { + totalMutations++ + try { + double ddG = Math.abs(parts[6].toDouble()) + if (ddG > 0.5) { + significantDDGs.add(ddG) + } + } catch (NumberFormatException e) { + logger.warning("Could not parse ddG value in line ${linesProcessed}: ${line}") + } + } + } + } + } + + logger.info("MAESTRO file parsing complete") + logger.info("Total lines processed: ${linesProcessed}") + logger.info("Total mutations found: ${totalMutations}") + logger.info("Significant mutations (|ΔΔG| > 0.5): ${significantDDGs.size()}") + + // Calculate IHS according to the formula: (1/|S|) * sum(|ΔΔG|) + if (significantDDGs.isEmpty()) { + logger.info("No significant mutations found for IHS calculation") + return 0.0 + } + + double sumDDG = significantDDGs.sum() + double ihs = sumDDG / significantDDGs.size() + + logger.info("IHS calculation results:") + logger.info(" Sum of |ΔΔG| values: ${String.format("%.2f", sumDDG)}") + logger.info(" Number of significant mutations: ${significantDDGs.size()}") + logger.info(" Calculated IHS: ${String.format("%.2f", ihs)}") + + return ihs + + } catch (Exception e) { + logger.warning("Error parsing MAESTRO file ${maestroFile.name}: ${e.message}") + e.printStackTrace() + return 0.0 + } + } + + + private List> collectBFactors(List residues) { + List> bFactorData = [] + + residues.each { residue -> + residue.getAtomList().each { atom -> + bFactorData << [ + resNum: residue.getResidueNumber(), + resName: residue.getName(), + atomName: atom.getName(), + chain: residue.getChainID(), + bFactor: atom.getTempFactor() + ] + } + } + return bFactorData + } + + private void writeBFactorsToCSV(List> bFactorData) { + String baseName = FilenameUtils.getBaseName(filename) + String outputFilename = csvFilename ?: "b_factors_${baseName}.csv" + + try { + File outputFile = new File(outputFilename) + boolean fileExists = outputFile.exists() + + outputFile.withWriter { writer -> + if (!fileExists || !csvFilename) { + writer.write("Residue_Number${delimiter}Residue_Name${delimiter}Atom_Name${delimiter}Chain${delimiter}B-factor\n") + } + + bFactorData.each { data -> + writer.write("${data.resNum}${delimiter}${data.resName}${delimiter}${data.atomName}${delimiter}${data.chain}${delimiter}${data.bFactor}\n") + } + } + + logger.info("\n Successfully wrote B-factors to ${outputFilename}") + } catch (Exception e) { + logger.severe("Error writing B-factors to CSV: ${e.message}") + } + } + + private double calculateMedian(List> bFactorData) { + List bFactors = bFactorData.collect { it.bFactor } + bFactors.sort() + + int size = bFactors.size() + if (size == 0) { + return Double.NaN + } + + if (size % 2 == 0) { + return (bFactors[size/2 - 1] + bFactors[size/2]) / 2.0 + } else { + return bFactors[size/2] + } + } + + private void printMedianBFactor(List> bFactorData) { + List bFactors = bFactorData.collect { it.bFactor } + bFactors.sort() + double median = calculateMedian(bFactorData) + + logger.info("\n B-factor statistics:") + logger.info(String.format(" Total atoms: %d", bFactors.size())) + logger.info(String.format(" Minimum B-factor: %.2f", bFactors.min())) + logger.info(String.format(" Maximum B-factor: %.2f", bFactors.max())) + logger.info(String.format(" Median B-factor: %.2f", median)) + } + + private void writeExtendedResultsToCSV(List> results) { + String outputFilename = "PDB_Features.csv" + + try { + File outputFile = new File(outputFilename) + outputFile.withWriter { writer -> + // Write header + if (calculateIHS) { + writer.write("PDB_ID${delimiter}Median_Bfactor${delimiter}Total_SASA${delimiter}Num_Residues${delimiter}IHS\n") + } else { + writer.write("PDB_ID${delimiter}Median_Bfactor${delimiter}Total_SASA${delimiter}Num_Residues\n") + } + + // Write data + results.each { result -> + if (calculateIHS) { + writer.write("${result.pdbId}${delimiter}${result.medianBfactor}${delimiter}${result.totalSASA}${delimiter}${result.numResidues}${delimiter}${result.IHS}\n") + } else { + writer.write("${result.pdbId}${delimiter}${result.medianBfactor}${delimiter}${result.totalSASA}${delimiter}${result.numResidues}\n") + } + } + } + } catch (Exception e) { + logger.severe("Error writing results to CSV: ${e.message}") + } + } +} \ No newline at end of file diff --git a/modules/potential/src/main/groovy/ffx/potential/groovy/NucleicAcidAnalysis.groovy b/modules/potential/src/main/groovy/ffx/potential/groovy/NucleicAcidAnalysis.groovy new file mode 100644 index 0000000000..14db90f54d --- /dev/null +++ b/modules/potential/src/main/groovy/ffx/potential/groovy/NucleicAcidAnalysis.groovy @@ -0,0 +1,300 @@ +//****************************************************************************** +// +// Title: Force Field X. +// Description: Force Field X - Software for Molecular Biophysics. +// Copyright: Copyright (c) Michael J. Schnieders 2001-2024. +// +// This file is part of Force Field X. +// +// Force Field X is free software; you can redistribute it and/or modify it +// under the terms of the GNU General Public License version 3 as published by +// the Free Software Foundation. +// +// Force Field X is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +// details. +// +// You should have received a copy of the GNU General Public License along with +// Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple +// Place, Suite 330, Boston, MA 02111-1307 USA +// +// Linking this library statically or dynamically with other modules is making a +// combined work based on this library. Thus, the terms and conditions of the +// GNU General Public License cover the whole combination. +// +// As a special exception, the copyright holders of this library give you +// permission to link this library with independent modules to produce an +// executable, regardless of the license terms of these independent modules, and +// to copy and distribute the resulting executable under terms of your choice, +// provided that you also meet, for each linked independent module, the terms +// and conditions of the license of that module. An independent module is a +// module which is not derived from or based on this library. If you modify this +// library, you may extend this exception to your version of the library, but +// you are not obligated to do so. If you do not wish to do so, delete this +// exception statement from your version. +// +//****************************************************************************** + +package ffx.potential.groovy + +import ffx.numerics.math.DoubleMath +import ffx.potential.ForceFieldEnergy +import ffx.potential.bonded.Atom +import ffx.potential.bonded.NucleicAcidUtils +import ffx.potential.bonded.Residue +import ffx.potential.cli.PotentialScript +import ffx.potential.parsers.PDBFilter +import ffx.potential.parsers.SystemFilter +import ffx.potential.parsers.XYZFilter +import ffx.crystal.Crystal +import picocli.CommandLine.Command +import picocli.CommandLine.Parameters + +import static java.lang.String.format +import static org.apache.commons.io.FilenameUtils.getBaseName +import static org.apache.commons.math3.util.FastMath.toDegrees +import static org.apache.commons.math3.util.FastMath.abs +import static org.apache.commons.math3.util.FastMath.atan +import static org.apache.commons.math3.util.FastMath.sin +import static org.apache.commons.math3.util.FastMath.cos + +@Command(description = "Nucleic Acid Analysis", name = "NucleicAcidAnalysis") +class NucleicAcidAnalysis extends PotentialScript { + + @Parameters(arity = "1", paramLabel = "file", + description = 'The atomic coordinate file in PDB or XYZ or ARC format.') + private List filenames = null + + private List residues + + NucleicAcidAnalysis() { + this(new Binding()) + } + + NucleicAcidAnalysis(Binding binding) { + super(binding) + } + + @Override + NucleicAcidAnalysis run() { + if (!init()) { + return null + } + + activeAssembly = getActiveAssembly(filenames[0]) + if (activeAssembly == null) { + logger.info(helpString()) + return null + } + + ForceFieldEnergy forceFieldEnergy = activeAssembly.getPotentialEnergy() + double[] x = new double[forceFieldEnergy.getNumberOfVariables()] + forceFieldEnergy.getCoordinates(x) + forceFieldEnergy.energy(x) + + // Create NucleicAcidUtils instance once + def naUtils = new NucleicAcidUtils() + + // Filter for nucleic acid residues + residues = activeAssembly.getResidueList().findAll { residue -> + def name = residue.name.toUpperCase() + name in ["DA", "DC", "DG", "DT", "DU", "A", "C", "G", "T", "U", + "DAD", "DCY", "DGU", "DTY", "URA"] + } + + SystemFilter systemFilter = potentialFunctions.getFilter() + if (systemFilter instanceof XYZFilter || systemFilter instanceof PDBFilter) { + int numModels = systemFilter.countNumModels() + + if (numModels == 1) { + // Single frame - output all information + printHeader() + analyzeAndPrintResidues() + } else { + // Multiple frames (ARC file) - output only pseudorotation and sugar pucker per frame + int frameNumber = 1 + while (systemFilter.readNext()) { + Crystal crystal = activeAssembly.getCrystal() + forceFieldEnergy.setCrystal(crystal) + forceFieldEnergy.getCoordinates(x) + forceFieldEnergy.energy(x) + + // Filter residues again in case the model changed + residues = activeAssembly.getResidueList().findAll { residue -> + def name = residue.name.toUpperCase() + name in ["DA", "DC", "DG", "DT", "DU", "A", "C", "G", "T", "U", + "DAD", "DCY", "DGU", "DTY", "URA"] + } + + println("\nFrame $frameNumber:") + println("Residue Name P Sugar pucker") + println("-------------------------------------------") + + for (Residue residue : residues) { + def v0 = getDihedral(residue, "C2'", "C1'", "O4'", "C4'") + def v1 = getDihedral(residue, "O4'", "C1'", "C2'", "C3'") + def v2 = getDihedral(residue, "C1'", "C2'", "C3'", "C4'") + def v3 = getDihedral(residue, "O4'", "C4'", "C3'", "C2'") + def v4 = getDihedral(residue, "C1'", "O4'", "C4'", "C3'") + + Double P = (v0 != null && v1 != null && v3 != null && v4 != null && v2 != null) ? + calculateP(v0, v1, v2, v3, v4) : null + String sugarPucker = determineSugarPucker(P) + + println(String.format("%-10s %-8s %-10s %-14s", + residue.getResidueNumber(), + residue.name, + formatValue(P), + sugarPucker + )) + } + frameNumber++ + } + } + } + + return this + } + + private void printHeader() { + println("Residue Name V0 V1 V2 V3 V4 P νmax χ γ α β δ ε ζ Sugar pucker Stage") + println("------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------") + } + + private void analyzeAndPrintResidues() { + for (Residue residue : residues) { + def v0 = getDihedral(residue, "C2'", "C1'", "O4'", "C4'") + def v1 = getDihedral(residue, "O4'", "C1'", "C2'", "C3'") + def v2 = getDihedral(residue, "C1'", "C2'", "C3'", "C4'") + def v3 = getDihedral(residue, "O4'", "C4'", "C3'", "C2'") + def v4 = getDihedral(residue, "C1'", "O4'", "C4'", "C3'") + def chi = getDihedral(residue, "O4'", "C1'", "N9", "C4") ?: getDihedral(residue, "O4'", "C1'", "N1", "C2") + def gamma = getDihedral(residue, "O5'", "C5'", "C4'", "C3'") + def alpha = getDihedral(residue, "O3'", "P", "O5'", "C5'") + def beta = getDihedral(residue, "P", "O5'", "C5'", "C4'") + def delta = getDihedral(residue, "C5'", "C4'", "C3'", "O3'") + def epsilon = getDihedral(residue, "C4'", "C3'", "O3'", "P") + def zeta = getDihedral(residue, "C3'", "O3'", "P", "O5'") + + Double P = (v0 != null && v1 != null && v3 != null && v4 != null && v2 != null) ? calculateP(v0, v1, v2, v3, v4) : null + double nuMax = (v2 != null && P != null) ? Math.abs(v2 / Math.cos(Math.toRadians(P))) : null + String type = determineType(residue, P) + String stage = determineStage(delta, chi, P) + + println(String.format("%-10s %-8s %-10s %-10s %-10s %-10s %-10s %-10s %-10s %-10s %-10s %-10s %-10s %-10s %-10s %-10s %-14s %-18s", + residue.getResidueNumber(), + residue.name, + formatValue(v0), formatValue(v1), formatValue(v2), + formatValue(v3), formatValue(v4), + formatValue(P), formatValue(nuMax), + formatValue(chi), formatValue(gamma), + formatValue(alpha), formatValue(beta), + formatValue(delta), formatValue(epsilon), + formatValue(zeta), + type, + stage + )) + } + } + + private Double getDihedral(Residue residue, String atom1, String atom2, String atom3, String atom4) { + try { + Atom a1 = residue.getAtomByName(atom1, true) + Atom a2 = residue.getAtomByName(atom2, true) + Atom a3 = residue.getAtomByName(atom3, true) + Atom a4 = residue.getAtomByName(atom4, true) + + if (a1 != null && a2 != null && a3 != null && a4 != null) { + return toDegrees(DoubleMath.dihedralAngle(a1.getXYZ(null), a2.getXYZ(null), a3.getXYZ(null), a4.getXYZ(null))) + } + } catch (Exception e) { + logger.warning("Could not calculate dihedral for ${residue} atoms $atom1-$atom2-$atom3-$atom4") + } + return null + } + + private String formatValue(Double value) { + return value != null ? String.format("%.2f", value) : "N/A" + } + + private static double calculateP(double v0, double v1, double v2, double v3, double v4) { + try { + double sin36 = Math.sin(Math.toRadians(36)) + double sin72 = Math.sin(Math.toRadians(72)) + double denominator = 2 * v2 * (sin36 + sin72) + if (Math.abs(denominator) < 1e-10) { + return Double.NaN + } + double p = ((v4 - v0) - (v3 - v1)) / denominator + + double P; + if (v2 < 0) { + P = Math.toDegrees(Math.atan(p)) + 180.0 + } else if (p < 0) { + P = Math.toDegrees(Math.atan(p)) + 360.0 + } else { + P = Math.toDegrees(Math.atan(p)) + } + return P + } catch (Exception e) { + return Double.NaN + } + } + + private String determineType(Residue residue, Double P) { + if (P == null) return "Unknown" + + String base = switch (residue.name) { + case "DAD" -> "Ade" + case "DGU" -> "Gua" + case "DCY" -> "Cyt" + case "DTY" -> "Thy" + case "URA" -> "Ura" + default -> "Unknown" + } + + String sugarPucker = determineSugarPucker(P) + return base + ", " + sugarPucker + } + + private String determineSugarPucker(Double P) { + if (P == null) return "Unknown" + + if (P >= 0 && P < 36) { + return "C3'-endo" + } else if (P >= 36 && P < 72) { + return "C4'-endo" + } else if (P >= 72 && P < 108) { + return "O4'-endo" + } else if (P >= 108 && P < 144) { + return "C1'-exo" + } else if (P >= 144 && P < 180) { + return "C2'-endo" + } else if (P >= 180 && P < 216) { + return "C3'-exo" + } else if (P >= 216 && P < 252) { + return "C4'-exo" + } else if (P >= 252 && P < 288) { + return "O4'-exo" + } else if (P >= 288 && P < 324) { + return "C1'-endo" + } else if (P >= 324 || P < 0) { + return "C2'-endo" + } + return "Unknown" + } + + private String determineStage(Double delta, Double chi, Double P) { + if (delta == null || chi == null || P == null) return "Unknown" + + if (delta > 120.0 && chi > -160) { + return "Stage 1: Slide and Roll Change" + } else if (delta < 120.0 && P >= 0 && P < 180.0) { + return "Stage 2: Backbone and Sugar Pucker Change" + } else if (delta < 120.0 && P >= 180.0 && P < 360.0) { + return "Stage 3: Roll and Inclination to Finish Transition" + } + return "Unknown" + } +} \ No newline at end of file diff --git a/modules/potential/src/main/groovy/ffx/potential/groovy/Solvator.groovy b/modules/potential/src/main/groovy/ffx/potential/groovy/Solvator.groovy index f075988b31..1edde7182f 100644 --- a/modules/potential/src/main/groovy/ffx/potential/groovy/Solvator.groovy +++ b/modules/potential/src/main/groovy/ffx/potential/groovy/Solvator.groovy @@ -63,6 +63,9 @@ import java.util.stream.Collectors import static org.apache.commons.math3.util.FastMath.round import static java.lang.String.format +import ffx.potential.extended.ExtendedSystem +import ffx.potential.utils.StructureMetrics + /** * The Solvator script puts a box of solvent around a solute. *
@@ -77,7 +80,7 @@ class Solvator extends PotentialScript { * --sFi or --solventFile specifies the file to read for the solvent. There is currently no default. */ @Option(names = ['--sFi', '--solventFile'], paramLabel = "water", - description = 'A file containing the solvent box to be used. There is currently no default.') + description = 'A file containing the solvent box to be used. There is currently no default.') private String solventFileName = null /** @@ -85,7 +88,7 @@ class Solvator extends PotentialScript { * with format ReferenceAtomStart ReferenceAtomEnd Concentration [Charge]. See forcefieldx/examples/nacl.ions. */ @Option(names = ['--iFi', '--ionFile'], paramLabel = "ions", - description = "Name of the file containing ions. Must also have a .ions file (e.g. nacl.pdb must also have nacl.ions). Default: no ions.") + description = "Name of the file containing ions. Must also have a .ions file (e.g. nacl.pdb must also have nacl.ions). Default: no ions.") private String ionFileName = null /** @@ -93,21 +96,21 @@ class Solvator extends PotentialScript { * but is not recommended for simulations long enough to see solute rotation. */ @Option(names = ['-r', '--rectangular'], paramLabel = "false", defaultValue = "false", - description = "Use a rectangular prism rather than a cube for solvation.") + description = "Use a rectangular prism rather than a cube for solvation.") private boolean rectangular = false /** * -p or --padding Sets the minimum amount of solvent padding around the solute. */ @Option(names = ['-p', '--padding'], paramLabel = "9.0", defaultValue = "9.0", - description = "Sets the minimum amount of solvent padding around the solute.") + description = "Sets the minimum amount of solvent padding around the solute.") private double padding = 9.0 /** * -b or --boundary Delete solvent molecules that infringe closer than this to the solute. */ @Option(names = ['-b', '--boundary'], paramLabel = "2.5", defaultValue = "2.5", - description = "Delete solvent molecules that infringe closer than this to the solute.") + description = "Delete solvent molecules that infringe closer than this to the solute.") private double boundary = 2.5 /** * -x or --translate Move solute molecules to center of box. @@ -120,18 +123,22 @@ class Solvator extends PotentialScript { * --abc or --boxLengths Specify unit cell box lengths, instead of calculating them. */ @Option(names = ['--abc', '--boxLengths'], paramLabel = "a,b,c", - description = "Specify a comma-separated set of unit cell box lengths, instead of calculating them (a,b,c)") + description = "Specify a comma-separated set of unit cell box lengths, instead of calculating them (a,b,c)") private String manualBox = null @Option(names = ['-s', '--randomSeed'], paramLabel = "auto", - description = "Specify a random seed for ion placement.") + description = "Specify a random seed for ion placement.") private String seedString = null + @Option(names = ['--pH', '--constantPH'], paramLabel = '7.4', + description = 'pH value for the system. If set, titration states will be initialized based on this pH.') + private Double pH = null + /** * The final argument is an atomic coordinate file in PDB or XYZ format. */ @Parameters(arity = "1", paramLabel = "file", - description = 'The atomic coordinate file in PDB or XYZ format.') + description = 'The atomic coordinate file in PDB or XYZ format.') String filename = null private MolecularAssembly solute @@ -193,10 +200,57 @@ class Solvator extends PotentialScript { return this } + // Calculate the monopole moment (total charge) of the system + ForceFieldEnergy forceFieldEnergy = activeAssembly.getPotentialEnergy() + + // If pH is specified, create and attach an ExtendedSystem + if (pH != null) { + logger.info("\n Initializing titration states for pH " + pH) + + ExtendedSystem esvSystem = new ExtendedSystem(activeAssembly, pH, null) + forceFieldEnergy.attachExtendedSystem(esvSystem) + + // Set titration states based on pH + esvSystem.reGuessLambdas() + logger.info(esvSystem.getLambdaList()) + + } + + // Coordinates and energy call to ensure PME and other potentials are fully initialized + int nVars = forceFieldEnergy.getNumberOfVariables() + double[] x = new double[nVars] + forceFieldEnergy.getCoordinates(x) + forceFieldEnergy.energy(x, true) + + // Now safely compute the monopole (total charge) + if (forceFieldEnergy.getPmeNode() == null) { + logger.severe("PME (electrostatics) is not initialized. Check force field settings.") + return this + } + + Atom[] activeAtoms = activeAssembly.getActiveAtomArray() + if (activeAtoms == null || activeAtoms.length == 0) { + logger.severe("No active atoms found in the assembly.") + return this + } + + double[] moments + try { + moments = forceFieldEnergy.getPmeNode().computeMoments(activeAtoms, false) + if (moments == null || moments.length == 0) { + logger.severe("Failed to compute moments (null or empty array returned).") + return this + } + //double totalCharge = moments[0] // First element is net charge + logger.info(format("\n System has a total charge of %8.4g before solvation", moments[0])) + } catch (Exception e) { + logger.severe("Error computing moments: " + e.getMessage()) + return this + } + solute = activeAssembly ForceFieldEnergy soluteEnergy = activeAssembly.getPotentialEnergy() - solvent = potentialFunctions.open(solventFileName) Atom[] baseSolventAtoms = solvent.getActiveAtomArray() @@ -223,7 +277,7 @@ class Solvator extends PotentialScript { if (!soluteCrystal.aperiodic()) { if (!soluteCrystal.spaceGroup.shortName.equalsIgnoreCase("P1")){ - logger.severe(" Solute must be aperiodic or strictly P1 periodic") + logger.severe(" Solute must be aperiodic or strictly P1 periodic") } } @@ -231,7 +285,6 @@ class Solvator extends PotentialScript { logger.severe(" Solvent box must be periodic (and P1)!") } - for (int i = 0; i < nSolute; i++) { Atom ati = soluteAtoms[i] double[] xyzi = new double[3] @@ -292,7 +345,7 @@ class Solvator extends PotentialScript { } } else { logger.severe( - " The manualBox option requires either 1 box length, or 3 comma-separated values.") + " The manualBox option requires either 1 box length, or 3 comma-separated values.") } } else { if (rectangular) { @@ -307,8 +360,8 @@ class Solvator extends PotentialScript { } logger.info( - format(" Molecule will be solvated in a periodic box of size %10.4g, %10.4g, %10.4g", - newBox[0], newBox[1], newBox[2])) + format(" Molecule will be solvated in a periodic box of size %10.4g, %10.4g, %10.4g", + newBox[0], newBox[1], newBox[2])) if(translate){ double[] soluteTranslate = new double[3] @@ -334,8 +387,8 @@ class Solvator extends PotentialScript { } logger.info( - format(" Solute size: %12.4g, %12.4g, %12.4g", soluteBoundingBox[0], soluteBoundingBox[1], - soluteBoundingBox[2])) + format(" Solute size: %12.4g, %12.4g, %12.4g", soluteBoundingBox[0], soluteBoundingBox[1], + soluteBoundingBox[2])) int[] solventReplicas = new int[3] double[] solventBoxVectors = new double[3] @@ -348,7 +401,7 @@ class Solvator extends PotentialScript { } Crystal newCrystal = new Crystal(newBox[0], newBox[1], newBox[2], 90, 90, 90, "P1") - soluteEnergy.setCrystal(newCrystal, true) + forceFieldEnergy.setCrystal(newCrystal, true) List bondedNodes = solvent.getAllBondedEntities() MSNode[] solventEntities = bondedNodes.toArray(new MSNode[bondedNodes.size()]) @@ -381,10 +434,10 @@ class Solvator extends PotentialScript { double[] xyzOffset = new double[3] int currXYZIndex = Arrays.stream(soluteAtoms). - mapToInt({ - it.getXyzIndex() - }) - .max().getAsInt() + mapToInt({ + it.getXyzIndex() + }) + .max().getAsInt() int currResSeq = 1 char[] possibleChains = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X', 'Y', 'Z'] as char[] @@ -442,13 +495,13 @@ class Solvator extends PotentialScript { comi[j] = xyzOffset[j] + solventCOMs[i][j] if (comi[j] < 0) { logger.warning(format( - " Skipping a copy of solvent molecule %d for violating minimum boundary 0,0,0. This should not occur!", - i)) + " Skipping a copy of solvent molecule %d for violating minimum boundary 0,0,0. This should not occur!", + i)) continue MoleculePlace } else if (comi[j] > newBox[j]) { logger.fine(format( - " Skipping a copy of solvent molecule %d for violating maximum boundary %12.4g,%12.4g,%12.4g", - i, newBox[0], newBox[1], newBox[2])) + " Skipping a copy of solvent molecule %d for violating maximum boundary %12.4g,%12.4g,%12.4g", + i, newBox[0], newBox[1], newBox[2])) continue MoleculePlace } } @@ -466,10 +519,10 @@ class Solvator extends PotentialScript { newXYZ[j] += xyzOffset[j] } Atom newAtom = new Atom(++currXYZIndex, parentAtom, newXYZ, currResSeq, solventChain, - Character.toString(solventChain)) + Character.toString(solventChain)) logger.fine( - format(" New atom %s at chain %c on residue %s-%d", newAtom, newAtom.getChainID(), - newAtom.getResidueName(), newAtom.getResidueNumber())) + format(" New atom %s at chain %c on residue %s-%d", newAtom, newAtom.getChainID(), + newAtom.getResidueName(), newAtom.getResidueNumber())) newAtom.setHetero(!moli instanceof Polymer) newAtom.setResName(moli.getName()) if (newAtom.getAltLoc() == null) { @@ -488,7 +541,7 @@ class Solvator extends PotentialScript { if (overlapFound) { logger.fine( - format(" Skipping a copy of molecule %d for overlapping with the solute.", i)) + format(" Skipping a copy of molecule %d for overlapping with the solute.", i)) continue MoleculePlace } @@ -522,7 +575,7 @@ class Solvator extends PotentialScript { BufferedReader br Pattern ionicPattern = - Pattern.compile("^\\s*([0-9]+) +([0-9]+) +([0-9]+(?:\\.[0-9]*)?|NEUT\\S*)") + Pattern.compile("^\\s*([0-9]+) +([0-9]+) +([0-9]+(?:\\.[0-9]*)?|NEUT\\S*)") Pattern concPatt = Pattern.compile("^[0-9]+(?:\\.[0-9]*)?") // Parse .ions file to figure out which ions need to be added. try { @@ -560,7 +613,7 @@ class Solvator extends PotentialScript { neutAnion = ia } else { logger.severe(format(" Specified a neutralizing ion %s with no net charge!", - ia.toString())) + ia.toString())) } } else { byConc.add(ia) @@ -576,7 +629,8 @@ class Solvator extends PotentialScript { int ionResSeq = 0 // Begin swapping waters for ions. - double initialCharge = activeAssembly.getCharge(false) + // Use the monopole moment (total charge) we calculated earlier + double initialCharge = moments[0] // @formatter:off logger.info(" Charge before addition of ions is ${initialCharge}") // @formatter:on @@ -587,7 +641,7 @@ class Solvator extends PotentialScript { if (nIons > newMolecules.size()) { // @formatter:off logger.severe( - " Insufficient solvent molecules remain (${newMolecules.size()}) to add ${nIons} ions!") + " Insufficient solvent molecules remain (${newMolecules.size()}) to add ${nIons} ions!") // @formatter:on } logger.info(format(" Number of ions to place: %d\n", nIons)) @@ -604,7 +658,7 @@ class Solvator extends PotentialScript { if (initialCharge > 0) { if (neutAnion == null) { logger.info( - " No counter-anion specified; system will be cationic at charge ${initialCharge}") + " No counter-anion specified; system will be cationic at charge ${initialCharge}") } else { logger.info(" Neutralizing system with ${neutAnion.toString()}") double charge = neutAnion.charge @@ -612,7 +666,7 @@ class Solvator extends PotentialScript { double netCharge = initialCharge + (nAnions * charge) if (nAnions > newMolecules.size()) { logger.severe( - " Insufficient solvent molecules remain (${newMolecules.size()}) to add ${nAnions} counter-anions!") + " Insufficient solvent molecules remain (${newMolecules.size()}) to add ${nAnions} counter-anions!") } for (int i = 0; i < nAnions; i++) { Atom[] newIon = swapIon(newMolecules, neutAnion, currXYZIndex, ionChain, ionResSeq++) @@ -620,13 +674,13 @@ class Solvator extends PotentialScript { addedIons.add(newIon) } logger.info( - format(" System neutralized to %8.4g charge with %d counter-anions", netCharge, - nAnions)) + format(" System neutralized to %8.4g charge with %d counter-anions", netCharge, + nAnions)) } } else if (initialCharge < 0) { if (neutCation == null) { logger.info( - " No counter-cation specified; system will be anionic at charge ${initialCharge}") + " No counter-cation specified; system will be anionic at charge ${initialCharge}") } else { logger.info(" Neutralizing system with ${neutCation.toString()}") double charge = neutCation.charge @@ -634,7 +688,7 @@ class Solvator extends PotentialScript { double netCharge = initialCharge + (nCations * charge) if (nCations > newMolecules.size()) { logger.severe( - " Insufficient solvent molecules remain (${newMolecules.size()}) to add ${nCations} counter-cations!") + " Insufficient solvent molecules remain (${newMolecules.size()}) to add ${nCations} counter-cations!") // @formatter:on } for (int i = 0; i < nCations; i++) { @@ -643,8 +697,8 @@ class Solvator extends PotentialScript { addedIons.add(newIon) } logger.info( - format(" System neutralized to %8.4g charge with %d counter-cations", netCharge, - nCations)) + format(" System neutralized to %8.4g charge with %d counter-cations", netCharge, + nCations)) } } else { logger.info(" System is neutral; no neutralizing ions needed.") @@ -822,7 +876,7 @@ class Solvator extends PotentialScript { } } logger.fine(format(" Constructing domain %3d-%3d-%3d with %d atoms.", i, j, k, - neighborAtoms.size())) + neighborAtoms.size())) logger.fine(format(" Neighbors along X: %s", neighborsI.toListString())) logger.fine(format(" Neighbors along Y: %s", neighborsJ.toListString())) logger.fine(format(" Neighbors along Z: %s\n", neighborsK.toListString())) @@ -836,8 +890,8 @@ class Solvator extends PotentialScript { double avgAtoms = ((double) soluteAtoms.length) / nBoxes time += System.nanoTime() logger.info(format( - " Decomposed the solute into %d domains of side length %11.4g, averaging %10.3g atoms apiece in %8.3g sec.", - nBoxes, domainLength, avgAtoms, (time * 1E-9))) + " Decomposed the solute into %d domains of side length %11.4g, averaging %10.3g atoms apiece in %8.3g sec.", + nBoxes, domainLength, avgAtoms, (time * 1E-9))) } /** @@ -854,27 +908,27 @@ class Solvator extends PotentialScript { int k = (int) (xyz[2] / domainLength) if (i >= nX) { logger.fine( - format(" Atom %s violates the X boundary at %12.7g Angstroms", a.toString(), xyz[0])) + format(" Atom %s violates the X boundary at %12.7g Angstroms", a.toString(), xyz[0])) i = 0 } if (j >= nY) { logger.fine( - format(" Atom %s violates the Y boundary at %12.7g Angstroms", a.toString(), xyz[1])) + format(" Atom %s violates the Y boundary at %12.7g Angstroms", a.toString(), xyz[1])) j = 0 } if (k >= nZ) { logger.fine( - format(" Atom %s violates the Y boundary at %12.7g Angstroms", a.toString(), xyz[2])) + format(" Atom %s violates the Y boundary at %12.7g Angstroms", a.toString(), xyz[2])) k = 0 } return Arrays.stream(withAdjacent[i][j][k]). - anyMatch({Atom s -> - double[] xyzS = new double[3] - xyzS = s.getXYZ(xyzS) - double dist = crystal.minDistOverSymOps(xyz, xyzS) - return dist < threshold - }) + anyMatch({Atom s -> + double[] xyzS = new double[3] + xyzS = s.getXYZ(xyzS) + double dist = crystal.minDistOverSymOps(xyz, xyzS) + return dist < threshold + }) } } @@ -958,9 +1012,9 @@ class Solvator extends PotentialScript { newXYZ[j] += atomOffsets[i][j] } Atom newAtom = new Atom(++currXYZIndex, fromAtom, newXYZ, resSeq, chain, - Character.toString(chain)) + Character.toString(chain)) logger.fine(format(" New ion atom %s at chain %c on ion molecule %s-%d", newAtom, - newAtom.getChainID(), newAtom.getResidueName(), newAtom.getResidueNumber())) + newAtom.getChainID(), newAtom.getResidueName(), newAtom.getResidueNumber())) newAtom.setHetero(true) newAtom.setResName(fromAtom.getResidueName()) if (newAtom.getAltLoc() == null) { @@ -974,8 +1028,8 @@ class Solvator extends PotentialScript { @Override String toString() { StringBuilder sb = new StringBuilder( - format(" Ion addition with %d atoms per ion, concentration %10.3g mM, " + - "charge %6.2f, used to neutralize %b", atoms.length, conc, charge, toNeutralize)) + format(" Ion addition with %d atoms per ion, concentration %10.3g mM, " + + "charge %6.2f, used to neutralize %b", atoms.length, conc, charge, toNeutralize)) for (Atom atom : atoms) { sb.append("\n").append(" Includes atom ").append(atom.toString()) } @@ -1018,7 +1072,7 @@ class Solvator extends PotentialScript { * @return Newly created set of ion atoms. */ private static Atom[] swapIon(List solvent, IonAddition ia, int currXYZIndex, char chain, - int resSeq) { + int resSeq) { int nSolv = solvent.size() - 1 Atom[] lastSolvent = solvent.get(nSolv) solvent.remove(nSolv) @@ -1026,4 +1080,4 @@ class Solvator extends PotentialScript { return ia.createIon(com, currXYZIndex, chain, resSeq) } -} +} \ No newline at end of file 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..bce59b8bd7 100755 --- a/modules/potential/src/main/java/ffx/potential/nonbonded/ParticleMeshEwald.java +++ b/modules/potential/src/main/java/ffx/potential/nonbonded/ParticleMeshEwald.java @@ -2726,7 +2726,7 @@ protected void assignPolarizationGroups() { * @param forceEnergy Force calculation of the electrostatic energy (rotate multipoles, perform * SCF). */ - public void computeMoments(Atom[] activeAtoms, boolean forceEnergy) { + public double[] computeMoments(Atom[] activeAtoms, boolean forceEnergy) { // Zero out total charge, dipole and quadrupole components. var netchg = 0.0; var netdpl = 0.0; @@ -2870,6 +2870,7 @@ public void computeMoments(Atom[] activeAtoms, boolean forceEnergy) { // Eigenvalues are returned in descending order, but logged below in ascending order. var netqdp = e.getRealEigenvalues(); + // Log the results (unchanged from original) logger.info("\n Electric Moments\n"); logger.info(format(" Total Electric Charge: %13.5f Electrons\n", netchg)); logger.info(format(" Dipole Moment Magnitude: %13.5f Debye\n", netdpl)); @@ -2878,8 +2879,13 @@ public void computeMoments(Atom[] activeAtoms, boolean forceEnergy) { logger.info(format(" (Buckinghams) %13.5f %13.5f %13.5f", yxqdp, yyqdp, yzqdp)); logger.info(format(" %13.5f %13.5f %13.5f\n", zxqdp, zyqdp, zzqdp)); logger.info( - format( - " Principal Axes Quadrupole %13.5f %13.5f %13.5f\n", netqdp[2], netqdp[1], netqdp[0])); + format( + " Principal Axes Quadrupole %13.5f %13.5f %13.5f\n", netqdp[2], netqdp[1], netqdp[0])); + + // Return the moments as a double array: [charge, dipole_magnitude, dipole_x, dipole_y, dipole_z, ...] + return new double[] { + netchg, netdpl, xdpl, ydpl, zdpl, xxqdp, xyqdp, xzqdp, yxqdp, yyqdp, yzqdp, zxqdp, zyqdp, zzqdp + }; } /**