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 }