View Javadoc
1   // ******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2024.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
36  //
37  // ******************************************************************************
38  package ffx.algorithms.dynamics.thermostats;
39  
40  import static java.lang.String.format;
41  import static org.apache.commons.math3.util.FastMath.exp;
42  import static org.apache.commons.math3.util.FastMath.sqrt;
43  
44  import ffx.potential.SystemState;
45  import ffx.numerics.Constraint;
46  import ffx.numerics.Potential.VARIABLE_TYPE;
47  import java.util.Collections;
48  import java.util.List;
49  import java.util.Random;
50  
51  /**
52   * Thermostat a molecular dynamics trajectory to an external bath using the Bussi, Donadio, and
53   * Parrinello method. This method is similar to Berendsen thermostat, but generates a canonical
54   * distribution.
55   *
56   * @author Michael J. Schnieders
57   * @see <a href="http://dx.doi.org/10.1016/j.cpc.2008.01.006">G. Bussi and M. Parrinello,
58   *     "Stochastic Thermostats: Comparison of Local and Global Schemes", Computer Physics
59   *     Communications, 179, 26-29 (2008)</a>
60   * @since 1.0
61   */
62  public class Bussi extends Thermostat {
63  
64    /** The random number generator used to perturb velocities. */
65    private final Random bussiRandom;
66    /** Bussi thermostat time constant (psec). */
67    private double tau;
68  
69    /**
70     * Constructor for Bussi.
71     *
72     * @param state The MDState to operate on.
73     * @param type the VARIABLE_TYPE of each variable.
74     * @param targetTemperature The target temperature.
75     * @param tau Bussi thermostat time constant (psec).
76     */
77    public Bussi(SystemState state, VARIABLE_TYPE[] type, double targetTemperature, double tau) {
78      this(state, type, targetTemperature, tau, Collections.emptyList());
79    }
80  
81    public Bussi(SystemState state, VARIABLE_TYPE[] type, double targetTemperature, double tau,
82        List<Constraint> constraints) {
83      super(state, type, targetTemperature, constraints);
84      this.name = ThermostatEnum.BUSSI;
85      this.tau = tau;
86      this.bussiRandom = new Random();
87    }
88  
89    /**
90     * Constructor for Bussi.
91     *
92     * @param state The MDState to operate on.
93     * @param type the VARIABLE_TYPE of each variable.
94     * @param targetTemperature a double.
95     */
96    public Bussi(SystemState state, VARIABLE_TYPE[] type, double targetTemperature) {
97      this(state, type, targetTemperature, 0.2e0);
98    }
99  
100   /**
101    * {@inheritDoc}
102    *
103    * <p>Full step velocity modification.
104    */
105   @Override
106   public void fullStep(double dt) {
107     double expTau = exp(-dt / tau);
108     double tempRatio = targetTemperature / state.getTemperature();
109     double rate = (1.0 - expTau) * tempRatio / degreesOfFreedom;
110     double r = bussiRandom.nextGaussian();
111     double s = 0.0;
112     for (int i = 0; i < degreesOfFreedom - 1; i++) {
113       double si = bussiRandom.nextGaussian();
114       s += si * si;
115     }
116     double scale = expTau + (s + r * r) * rate + 2.0 * r * sqrt(expTau * rate);
117     scale = sqrt(scale);
118     if (r + sqrt(expTau / rate) < 0.0) {
119       scale = -scale;
120     }
121     double[] v = state.v();
122     double[] mass = state.getMass();
123     for (int i = 0; i < state.getNumberOfVariables(); i++) {
124       if (mass[i] > 0.0) {
125         v[i] *= scale;
126       }
127     }
128   }
129 
130   /**
131    * Getter for the field <code>tau</code>.
132    *
133    * @return a double.
134    */
135   public double getTau() {
136     return tau;
137   }
138 
139   /**
140    * Setter for the field <code>tau</code>.
141    *
142    * @param tau a double.
143    */
144   public void setTau(double tau) {
145     this.tau = tau;
146   }
147 
148   /**
149    * {@inheritDoc}
150    *
151    * <p>No velocity modifications are made by the Bussi method at the half-step.
152    */
153   @Override
154   public void halfStep(double dt) {
155   }
156 
157   /**
158    * {@inheritDoc}
159    *
160    * <p>Initialize the Random number generator used to apply random forces to the particles.
161    */
162   public void setRandomSeed(long seed) {
163     bussiRandom.setSeed(seed);
164   }
165 
166   /**
167    * Add Thermostat details to the kinetic energy and temperature details.
168    *
169    * @return Description of the thermostat, kinetic energy and temperature.
170    */
171   public String toThermostatString() {
172     return format("\n Bussi Thermostat (tau = %8.3f psec)\n%s", tau, super.toString());
173   }
174 
175   /** {@inheritDoc} */
176   @Override
177   public String toString() {
178     return "Bussi";
179   }
180 }