View Javadoc
1   //******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2025.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
36  //
37  //******************************************************************************
38  package ffx.algorithms.commands;
39  
40  import ffx.algorithms.cli.AlgorithmsCommand;
41  import ffx.algorithms.cli.DynamicsOptions;
42  import ffx.algorithms.cli.LambdaParticleOptions;
43  import ffx.algorithms.cli.MultiDynamicsOptions;
44  import ffx.algorithms.cli.OSTOptions;
45  import ffx.algorithms.cli.ThermodynamicsOptions;
46  import ffx.algorithms.thermodynamics.OrthogonalSpaceTempering;
47  import ffx.numerics.Potential;
48  import ffx.potential.ForceFieldEnergy;
49  import ffx.utilities.FFXBinding;
50  import org.apache.commons.io.FilenameUtils;
51  import picocli.CommandLine.Command;
52  import picocli.CommandLine.Option;
53  import picocli.CommandLine.Parameters;
54  
55  import java.io.File;
56  import java.io.FileWriter;
57  import java.io.IOException;
58  import java.util.Collections;
59  import java.util.List;
60  
61  /**
62   * The Histogram script prints out a Orthogonal Space histogram from a *.his file.
63   * <br>
64   * Usage:
65   * <br>
66   * ffxc Histogram [options] &lt;filename&gt;
67   */
68  @Command(description = " Evaluate the Orthogonal Space Histogram.", name = "Histogram")
69  public class Histogram extends AlgorithmsCommand {
70  
71    /**
72     * -s or --save Save the histogram, PMF and 2D bias to files.
73     */
74    @Option(names = {"-s", "--save"}, paramLabel = "false",
75        description = "Save the bias histogram to histogram.txt, the total PMF to pmf.txt, and 2D PMF to pmf.2D.txt")
76    private boolean save = false;
77  
78    /**
79     * -b or --bias By default, the PMF is saved. This flag flips the sign to give the OST bias.
80     */
81    @Option(names = {"-b", "--bias"}, paramLabel = "false",
82        description = "By default, the PMF is saved. This flag flips the sign to give the OST bias.")
83    private boolean bias = false;
84  
85    /**
86     * An XYZ or PDB input file.
87     */
88    @Parameters(arity = "1", paramLabel = "file",
89        description = "XYZ or PDB input file.")
90    private String filename;
91  
92    private OrthogonalSpaceTempering orthogonalSpaceTempering;
93    private File saveDir = null;
94  
95    /**
96     * Histogram Constructor.
97     */
98    public Histogram() {
99      super();
100   }
101 
102   /**
103    * Histogram Constructor.
104    *
105    * @param binding The Binding to use.
106    */
107   public Histogram(FFXBinding binding) {
108     super(binding);
109   }
110 
111   /**
112    * Histogram constructor that sets the command line arguments.
113    *
114    * @param args Command line arguments.
115    */
116   public Histogram(String[] args) {
117     super(args);
118   }
119 
120   /**
121    * {@inheritDoc}
122    */
123   @Override
124   public Histogram run() {
125 
126     // Init the context and bind variables.
127     if (!init()) {
128       return this;
129     }
130 
131     // Load the MolecularAssembly.
132     activeAssembly = getActiveAssembly(filename);
133     if (activeAssembly == null) {
134       logger.info(helpString());
135       return this;
136     }
137 
138     // Set the filename.
139     filename = activeAssembly.getFile().getAbsolutePath();
140 
141     logger.info("\n Evaluating Histogram for " + filename);
142 
143     File structureFile = new File(FilenameUtils.normalize(filename));
144     structureFile = new File(structureFile.getAbsolutePath());
145     String baseFilename = FilenameUtils.removeExtension(structureFile.getName());
146     File histogramRestart = new File(baseFilename + ".his");
147     File lambdaRestart = null;
148 
149     // Get a reference to the active system's ForceFieldEnergy and atom array.
150     ForceFieldEnergy energy = activeAssembly.getPotentialEnergy();
151 
152     // Print the current energy
153     energy.energy(true, true);
154     logger.info("");
155 
156     if (saveDir == null || !saveDir.exists() || !saveDir.isDirectory() || !saveDir.canWrite()) {
157       saveDir = new File(FilenameUtils.getFullPath(filename));
158     }
159 
160     // Construct some options with defaults.
161     DynamicsOptions dynamicsOptions = new DynamicsOptions();
162     LambdaParticleOptions lambdaParticleOptions = new LambdaParticleOptions();
163     MultiDynamicsOptions multiDynamicsOptions = new MultiDynamicsOptions();
164     ThermodynamicsOptions thermodynamicsOptions = new ThermodynamicsOptions();
165     OSTOptions ostOptions = new OSTOptions();
166 
167     // Construct the Thermodynamics instance.
168     orthogonalSpaceTempering = ostOptions.constructOST(energy,
169         lambdaRestart, histogramRestart, activeAssembly, null,
170         dynamicsOptions, thermodynamicsOptions, lambdaParticleOptions,
171         algorithmListener, !multiDynamicsOptions.isSynchronous());
172 
173     if (save) {
174       orthogonalSpaceTempering.setMolecularAssembly(activeAssembly);
175       OrthogonalSpaceTempering.Histogram histogram = orthogonalSpaceTempering.getHistogram();
176       histogram.updateFreeEnergyDifference(false, true);
177       StringBuffer sb = histogram.evaluateTotalOSTBias(bias);
178 
179       String dirName = saveDir.toString() + File.separator;
180       String file = dirName + "pmf.txt";
181       logger.info(" Writing " + file);
182       try (FileWriter fileWriter = new FileWriter(file)) {
183         fileWriter.write(sb.toString());
184       } catch (IOException e) {
185         logger.severe("Error writing file " + file + ": " + e.getMessage());
186       }
187 
188       sb = histogram.evaluate2DOSTBias(bias);
189       file = dirName + "pmf.2D.txt";
190       logger.info(" Writing " + file);
191       try (FileWriter fileWriter = new FileWriter(file)) {
192         fileWriter.write(sb.toString());
193       } catch (IOException e) {
194         logger.severe("Error writing file " + file + ": " + e.getMessage());
195       }
196     }
197     return this;
198   }
199 
200   /**
201    * {@inheritDoc}
202    */
203   @Override
204   public List<Potential> getPotentials() {
205     List<Potential> potentials;
206     if (orthogonalSpaceTempering == null) {
207       potentials = Collections.emptyList();
208     } else {
209       potentials = Collections.singletonList(orthogonalSpaceTempering);
210     }
211     return potentials;
212   }
213 
214 }