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.crystal; 39 40 import java.util.List; 41 import java.util.Objects; 42 import java.util.logging.Level; 43 import java.util.logging.Logger; 44 45 import static java.lang.String.format; 46 47 /** 48 * The ReplicatesCrystal class extends Crystal to generate additional symmetry operators needed to 49 * describe a "replicated" super-cell. 50 * <p> 51 * The replicated crystal cell edges are of length {l*a, m*b, n*c} where l, m and n are integers and 52 * a, b and c are the original unit cell edge lengths. 53 * <br> 54 * The replicates integers l, m and n are chosen large enough for the ReplicatesCrystal to allow 55 * consistent application of the minimum image convention. This is ensured by increasing l, m and/or 56 * n until a sphere of necessary radius fits entirely inside the ReplicatedCrystal. 57 * <br> 58 * 59 * @author Michael J. Schnieders 60 * @see Crystal 61 * @since 1.0 62 */ 63 public class ReplicatesCrystal extends Crystal { 64 65 /** 66 * The logger. 67 */ 68 private static final Logger logger = Logger.getLogger(ReplicatesCrystal.class.getName()); 69 /** 70 * The base unit cell for the system being simulated. 71 */ 72 private final Crystal unitCell; 73 /** 74 * The cut-off distance in Angstroms. 75 */ 76 private final double cutOff; 77 /** 78 * The number of replicates along the a-axis. 79 */ 80 private int l; 81 /** 82 * The number of replicates along the b-axis. 83 */ 84 private int m; 85 /** 86 * The number of replicates along the c-axis. 87 */ 88 private int n; 89 90 /** 91 * Constructor for a ReplicatesCrystal. 92 * 93 * @param unitCell The base unit cell. 94 * @param l Number of replicates along the a-axis. 95 * @param m Number of replicates along the b-axis. 96 * @param n Number of replicates along the c-axis. 97 * @param cutOff2 Twice the cut-off distance. 98 * @since 1.0 99 */ 100 public ReplicatesCrystal(Crystal unitCell, int l, int m, int n, double cutOff2) { 101 super(unitCell.a * l, unitCell.b * m, unitCell.c * n, unitCell.alpha, unitCell.beta, unitCell.gamma, unitCell.spaceGroup.shortName); 102 this.unitCell = unitCell; 103 104 assert (l >= 1); 105 assert (m >= 1); 106 assert (n >= 1); 107 this.l = l; 108 this.m = m; 109 this.n = n; 110 this.cutOff = cutOff2 / 2.0; 111 112 /* 113 At this point, the ReplicatesCrystal references a SpaceGroup instance 114 whose symmetry operators are inconsistent. This is corrected by 115 generating symmetry operators to fill up the ReplicatesCrystal based 116 on the asymmetric unit. 117 */ 118 updateReplicateOperators(); 119 } 120 121 /** 122 * The number of replicates along the a-axis. 123 */ 124 public int getL() { 125 return l; 126 } 127 128 /** 129 * The number of replicates along the b-axis. 130 */ 131 public int getM() { 132 return m; 133 } 134 135 /** 136 * The number of replicates along the c-axis. 137 */ 138 public int getN() { 139 return n; 140 } 141 142 /** 143 * Returns a ReplicatesCrystal large enough to satisfy the minimum image convention for the 144 * specified unit cell and cutoff criteria. If the unit cell is already sufficiently large, then it 145 * is returned. 146 * 147 * @param unitCell The unit cell of the crystal. 148 * @param cutOff2 Two times the cutoff distance. 149 * @return A Crystal or ReplicatesCrystal large enough to satisfy the minimum image convention. 150 */ 151 public static Crystal replicatesCrystalFactory(Crystal unitCell, double cutOff2) { 152 153 return replicatesCrystalFactory(unitCell, cutOff2, new int[3]); 154 } 155 156 /** 157 * Returns a ReplicatesCrystal large enough to satisfy the minimum image convention for the 158 * specified unit cell and cutoff criteria. If the unit cell is already sufficiently large, then it 159 * is returned. 160 * 161 * @param unitCell The unit cell of the crystal. 162 * @param cutOff2 Two times the cutoff distance. 163 * @return A Crystal or ReplicatesCrystal large enough to satisfy the minimum image convention. 164 */ 165 public static Crystal replicatesCrystalFactory(Crystal unitCell, double cutOff2, int[] replicatesVector) { 166 167 if (unitCell == null || unitCell.aperiodic()) { 168 replicatesVector[0] = 0; 169 replicatesVector[1] = 0; 170 replicatesVector[2] = 0; 171 return unitCell; 172 } 173 174 int l = 1; 175 int m = 1; 176 int n = 1; 177 178 double cutOff = cutOff2 / 2.0; 179 180 while (unitCell.interfacialRadiusA * l < cutOff) { 181 l++; 182 } 183 while (unitCell.interfacialRadiusB * m < cutOff) { 184 m++; 185 } 186 while (unitCell.interfacialRadiusC * n < cutOff) { 187 n++; 188 } 189 190 replicatesVector[0] = l; 191 replicatesVector[1] = m; 192 replicatesVector[2] = n; 193 return new ReplicatesCrystal(unitCell, l, m, n, cutOff2); 194 195 } 196 197 /** 198 * Change the cell vectors for the base unit cell, which is followed by an update of the 199 * ReplicateCrystal parameters and possibly the number of replicated cells. 200 * 201 * @param cellVectors 3x3 matrix of cell vectors. 202 * @return True if the perturbation of cell vectors succeeds. 203 */ 204 public boolean setCellVectors(double[][] cellVectors) { 205 206 // First, update the parameters of the unit cell. 207 if (unitCell.setCellVectors(cellVectors)) { 208 209 // Then, update the parameters of the ReplicatesCrystal and possibly the number of replicates. 210 return updateReplicatesDimensions(); 211 } 212 return false; 213 } 214 215 /** 216 * Change the cell vectors and volume for the base unit cell, which is followed by an update of the 217 * ReplicateCrystal parameters and possibly the number of replicated cells. 218 * 219 * @param cellVectors 3x3 matrix of cell vectors. 220 * @param targetAUVolume the target volume for the new cell Vectors. 221 * @return True if the perturbation of cell vectors succeeds. 222 */ 223 public boolean setCellVectorsAndVolume(double[][] cellVectors, double targetAUVolume) { 224 // First, update the parameters of the unit cell. 225 if (unitCell.setCellVectorsAndVolume(cellVectors, targetAUVolume)) { 226 227 // Then, update the parameters of the ReplicatesCrystal and possibly the number of replicates. 228 return updateReplicatesDimensions(); 229 } 230 return false; 231 } 232 233 /** 234 * Change the cell parameters for the base unit cell, which is followed by an update of the 235 * ReplicateCrystal parameters and possibly the number of replicated cells. 236 * 237 * @param a The length of the a-axis for the base unit cell (in Angstroms). 238 * @param b The length of the b-axis for the base unit cell (in Angstroms). 239 * @param c The length of the c-axis for the base unit cell (in Angstroms). 240 * @param alpha The angle between the b-axis and c-axis (in Degrees). 241 * @param beta The angle between the a-axis and c-axis (in Degrees). 242 * @param gamma The angle between the a-axis and b-axis (in Degrees). 243 * @return True is returned if the unit cell and replicates cell are updated successfully. 244 */ 245 @Override 246 public boolean changeUnitCellParameters(double a, double b, double c, double alpha, double beta, double gamma) { 247 248 // First, update the parameters of the unit cell. 249 if (unitCell.changeUnitCellParameters(a, b, c, alpha, beta, gamma)) { 250 251 // Then, update the parameters of the ReplicatesCrystal and possibly the number of replicates. 252 return updateReplicatesDimensions(); 253 } 254 return false; 255 } 256 257 /** 258 * Change the cell parameters for the base unit cell, which is followed by an update of the 259 * ReplicateCrystal parameters and possibly the number of replicated cells. 260 * 261 * @param a The length of the a-axis for the base unit cell (in Angstroms). 262 * @param b The length of the b-axis for the base unit cell (in Angstroms). 263 * @param c The length of the c-axis for the base unit cell (in Angstroms). 264 * @param alpha The angle between the b-axis and c-axis (in Degrees). 265 * @param beta The angle between the a-axis and c-axis (in Degrees). 266 * @param gamma The angle between the a-axis and b-axis (in Degrees). 267 * @param targetAUVolume the target volume for the new cell Vectors. 268 * @return True is returned if the unit cell and replicates cell are updated successfully. 269 */ 270 @Override 271 public boolean changeUnitCellParametersAndVolume(double a, double b, double c, double alpha, double beta, double gamma, double targetAUVolume) { 272 273 // First, update the parameters of the unit cell. 274 if (unitCell.changeUnitCellParametersAndVolume(a, b, c, alpha, beta, gamma, targetAUVolume)) { 275 276 // Then, update the parameters of the ReplicatesCrystal and possibly the number of replicates. 277 return updateReplicatesDimensions(); 278 } 279 return false; 280 } 281 282 283 /** 284 * Two crystals are equal only if all unit cell parameters are exactly the same. 285 */ 286 @Override 287 public boolean equals(Object o) { 288 if (this == o) { 289 return true; 290 } 291 if (o == null || getClass() != o.getClass()) { 292 return false; 293 } 294 ReplicatesCrystal replicatesCrystal = (ReplicatesCrystal) o; 295 return (Objects.equals(unitCell, replicatesCrystal.unitCell) && a == replicatesCrystal.a && b == replicatesCrystal.b && c == replicatesCrystal.c && alpha == replicatesCrystal.alpha && beta == replicatesCrystal.beta && gamma == replicatesCrystal.gamma && spaceGroup.number == replicatesCrystal.spaceGroup.number && spaceGroup.symOps.size() == replicatesCrystal.spaceGroup.symOps.size()); 296 } 297 298 /** 299 * {@inheritDoc} 300 */ 301 @Override 302 public boolean getCheckRestrictions() { 303 return unitCell.getCheckRestrictions(); 304 } 305 306 /** 307 * {@inheritDoc} 308 */ 309 @Override 310 public void setCheckRestrictions(boolean checkRestrictions) { 311 this.checkRestrictions = checkRestrictions; 312 unitCell.setCheckRestrictions(checkRestrictions); 313 } 314 315 /** 316 * Return the density of the ReplicatesCrystal. 317 * 318 * @param mass The mass of the ReplicatesCrystal. 319 * @return The density. 320 */ 321 public double getDensity(double mass) { 322 return unitCell.getDensity(mass); 323 } 324 325 /** 326 * {@inheritDoc} 327 * 328 * <p>Returns the unit cell for this ReplicatesCrystal. This is useful for the reciprocal space 329 * portion of PME that operates on the unit cell even though the real space cutoff requires a 330 * ReplicatesCrystal. 331 */ 332 @Override 333 public Crystal getUnitCell() { 334 return unitCell; 335 } 336 337 /** 338 * Update the ReplicatesCrystal using random parameters with the target density. 339 * 340 * @param density Target density. 341 * @param mass Mass of the ReplicatesCrystal. 342 */ 343 public boolean randomParameters(double density, double mass) { 344 boolean succeed = unitCell.randomParameters(density, mass); 345 if (succeed) { 346 updateReplicatesDimensions(); 347 } 348 return succeed; 349 } 350 351 /** 352 * Update the ReplicatesCrystal dimensions to the target density. 353 * 354 * @param density Target density. 355 * @param mass Mass of the ReplicatesCrystal. 356 */ 357 public void setDensity(double density, double mass) { 358 unitCell.setDensity(density, mass); 359 updateReplicatesDimensions(); 360 } 361 362 /** 363 * {@inheritDoc} 364 * 365 * <p>Include information about the base unit cell and replicates cell. 366 */ 367 @Override 368 public String toString() { 369 StringBuilder sb = new StringBuilder(unitCell.toString()); 370 // Only log the replicates cell if there is more than one replicate. 371 if (l * m * n > 1) { 372 sb.append("\n\n Replicates Cell\n"); 373 sb.append(format(" Dimension: (%3d x%3d x%3d)\n", l, m, n)); 374 sb.append(format(" A-axis: %8.3f\n", a)); 375 sb.append(format(" B-axis: %8.3f\n", b)); 376 sb.append(format(" C-axis: %8.3f\n", c)); 377 sb.append(format(" Alpha: %8.3f\n", alpha)); 378 sb.append(format(" Beta: %8.3f\n", beta)); 379 sb.append(format(" Gamma: %8.3f\n", gamma)); 380 sb.append(format(" Total Symmetry Operators: %8d", spaceGroup.getNumberOfSymOps())); 381 } 382 return sb.toString(); 383 } 384 385 /** 386 * A String containing the replicated unit cell parameters. 387 * 388 * @return A string with the unit cell parameters. 389 */ 390 public String toShortString() { 391 return format("%6.2f %6.2f %6.2f %6.2f %6.2f %6.2f (%3d x%3d x%3d) ", a, b, c, alpha, beta, gamma, l, m, n); 392 } 393 394 private boolean updateReplicatesDimensions() { 395 // Then, update the parameters of the ReplicatesCrystal and possibly the number of replicates. 396 int ll = 1; 397 int mm = 1; 398 int nn = 1; 399 400 while (unitCell.interfacialRadiusA * ll < cutOff) { 401 ll++; 402 } 403 while (unitCell.interfacialRadiusB * mm < cutOff) { 404 mm++; 405 } 406 while (unitCell.interfacialRadiusC * nn < cutOff) { 407 nn++; 408 } 409 if (super.changeUnitCellParameters(unitCell.a * ll, unitCell.b * mm, unitCell.c * nn, unitCell.alpha, unitCell.beta, unitCell.gamma)) { 410 l = ll; 411 m = mm; 412 n = nn; 413 updateReplicateOperators(); 414 return true; 415 } 416 417 return false; 418 } 419 420 /** 421 * Update the list of symmetry operators used to generate the replicates super-cell from the 422 * asymmetric unit. 423 */ 424 private void updateReplicateOperators() { 425 List<SymOp> symOps = spaceGroup.symOps; 426 427 // First, we remove the existing symmetry operators. 428 symOps.clear(); 429 430 /* 431 Now create symmetry operators for each replicate. Note that the first 432 symOp is still equivalent to the asymmetric unit and the first set of 433 symOps are still equivalent to the unit cell. 434 */ 435 double dX = 1.0 / (double) l; 436 double dY = 1.0 / (double) m; 437 double dZ = 1.0 / (double) n; 438 int symOpCount = 0; 439 List<SymOp> ucSymOps = unitCell.spaceGroup.symOps; 440 for (int i = 0; i < l; i++) { 441 for (int j = 0; j < m; j++) { 442 for (int k = 0; k < n; k++) { 443 int ii = 0; 444 for (SymOp symOp : ucSymOps) { 445 double[] repTrans = new double[3]; 446 repTrans[0] = (symOp.tr[0] + i) * dX; 447 repTrans[1] = (symOp.tr[1] + j) * dY; 448 repTrans[2] = (symOp.tr[2] + k) * dZ; 449 SymOp repSymOp = new SymOp(symOp.rot, repTrans, new int[]{i, j, k}); 450 symOps.add(repSymOp); 451 if (logger.isLoggable(Level.FINEST)) { 452 logger.finest(format("\n SymOp %d (%2d,%2d,%2d): %d", symOpCount, i, j, k, ii)); 453 logger.finest(repSymOp.toString()); 454 } 455 ii++; 456 symOpCount++; 457 } 458 } 459 } 460 } 461 } 462 }