MolLib/MD/MDPressure.c

   1 /*
   2  *  MDPressure.c
   3  *  Molby
   4  *
   5  *  Created by Toshi Nagata on 07/10/12.
   6  *  Copyright 2007 Toshi Nagata. All rights reserved.
   7  *
   8  This program is free software; you can redistribute it and/or modify
   9  it under the terms of the GNU General Public License as published by
  10  the Free Software Foundation version 2 of the License.
  11
  12  This program is distributed in the hope that it will be useful,
  13  but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  GNU General Public License for more details.
  16  */
  17
  18 #include "MDPressure.h"
  19 #include "MDForce.h"
  20
  21 #include <stdlib.h>
  22 #include <math.h>
  23 #include <string.h>
  24
  25 /*  Calculate the lattice deformation stress  */
  26 static Double
  27 s_lattice_deformation_stress(const Mat33 apply, const Transform celltr, const Transform old_celltr)
  28 {
  29         Mat33 grad, grad_tr, grad_inv, grad_inv_tr;
  30         Mat33 strain, tension;
  31         Double eng, volume;
  32         Int i;
  33
  34         /*  Some of the calculations duplicate with md_scale_cell(), but it will not cause
  35                 a serious slowdown  */
  36         volume = fabs(MatrixDeterminant(celltr));
  37
  38         /*  Deformation gradient (temp) = celltr * old_rcelltr  */
  39         MatrixInvert(grad, old_celltr);
  40         MatrixMul(grad, celltr, grad);
  41         MatrixTranspose(grad_tr, grad);
  42         MatrixInvert(grad_inv, grad);
  43         MatrixTranspose(grad_inv_tr, grad_inv);
  44         MatrixMul(strain, grad_tr, grad);
  45         strain[0] -= 1.0;
  46         strain[4] -= 1.0;
  47         strain[8] -= 1.0;
  48         MatrixScale(strain, strain, 0.5);
  49         MatrixMul(tension, apply, grad_inv_tr);
  50         MatrixMul(tension, grad_inv, tension);
  51         MatrixScale(tension, tension, MatrixDeterminant(grad));
  52         eng = 0.0;
  53         for (i = 0; i < 9; i++)
  54                 eng += strain[i] * tension[i];
  55         eng *= volume * BAR2INTERNALPRESSURE;
  56         return eng;
  57 }
  58
  59 MDPressureArena *
  60 pressure_new(void)
  61 {
  62         MDPressureArena *pressure = (MDPressureArena *)calloc(sizeof(MDPressureArena), 1);
  63         if (pressure == NULL)
  64                 return NULL;
  65         pressure->freq = 20;
  66         pressure->coupling = 0.3;
  67         pressure->trial_width = 0.001;
  68         pressure->control_algorithm = 0;
  69         pressure->fluctuate_cell_origin = 0.01;
  70         pressure->fluctuate_cell_orientation = 0.01;
  71         pressure->apply[0] = 1;
  72         pressure->apply[1] = 0;
  73         pressure->apply[2] = 0;
  74         pressure->apply[3] = 0;
  75         pressure->apply[4] = 1;
  76         pressure->apply[5] = 0;
  77         pressure->apply[6] = 0;
  78         pressure->apply[7] = 0;
  79         pressure->apply[8] = 1;
  80         pressure->cell_flexibility[0] = -1;
  81         pressure->cell_flexibility[1] = -1;
  82         pressure->cell_flexibility[2] = -1;
  83         pressure->cell_flexibility[3] = 0;
  84         pressure->cell_flexibility[4] = 0;
  85         pressure->cell_flexibility[5] = 0;
  86         pressure->cell_flexibility[6] = 0;
  87         pressure->cell_flexibility[7] = 0;
  88         return pressure;
  89 }
  90
  91 void
  92 pressure_release(MDPressureArena *pressure)
  93 {
  94         if (pressure == NULL)
  95                 return;
  96         free(pressure);
  97 }
  98
  99 void
 100 pressure_prepare(MDArena *arena)
 101 {
 102         MDPressureArena *pressure = arena->pressure;
 103         if (pressure == NULL)
 104                 return;
 105         pressure->mc_accept = pressure->mc_reject = 0;
 106         if (pressure->temporary_energies != NULL)
 107                 free(pressure->temporary_energies);
 108         pressure->temporary_energies = (Double *)calloc(sizeof(Double), arena->natoms_uniq);
 109         if (pressure->temporary_energies == NULL)
 110                 md_panic(arena, "Low memory");
 111         if (pressure->temporary_velocities != NULL)
 112                 free(pressure->temporary_velocities);
 113         pressure->temporary_velocities = (Vector *)calloc(sizeof(Vector), arena->natoms_uniq);
 114         if (pressure->temporary_velocities == NULL)
 115                 md_panic(arena, "Low memory");
 116 }
 117
 118 /*  Pressure control  */
 119 void
 120 pressure_control(MDArena *arena)
 121 {
 122         Vector cella_save, cellb_save, cellc_save, cello_save;
 123         Double de, total_energy_save, energies_save[kEndIndex];
 124         MDPressureArena *pressure = arena->pressure;
 125         Atom *ap;
 126         Transform tf;
 127         Vector v;
 128         Vector cello_new;
 129         Transform celltr_save;
 130         Double w, w0;
 131         int i;
 132         int needs_cell_recalculate = 0;
 133         XtalCell *cell = arena->mol->cell;
 134
 135         if (pressure == NULL || (!arena->periodic_a && !arena->periodic_b && !arena->periodic_c))
 136                 return;
 137         if (pressure->freq == 0 || arena->step % pressure->freq != 0)
 138                 return;
 139
 140         while (md_rand() < pressure->coupling) {
 141
 142                 Double flex;
 143                 Transform tf0;
 144
 145                 /*  Trial move of the periodic cell  */
 146                 memset(tf, 0, sizeof(tf));
 147                 tf[0] = tf[4] = tf[8] = 1;
 148                 memmove(tf0, tf, sizeof(tf));
 149                 w0 = (md_rand() - 0.5) * pressure->trial_width;
 150                 for (i = 0; i < 3; i++) {
 151                         static char idx[] = {4, 8, 7, 5, 0, 8, 6, 2, 0, 4, 3, 1};
 152                         flex = pressure->cell_flexibility[i + 3];
 153                         if (flex > 0) {
 154                                 w = w0 = (md_rand() + md_rand() - 1.0) * pressure->trial_width * 0.5;
 155                         } else if (flex < 0) {
 156                                 w = w0;
 157                         } else w = 0;
 158                         w *= fabs(flex);
 159                         tf[idx[i * 4]] = tf[idx[i * 4 + 1]] = sqrt(1 - w * w);
 160                         tf[idx[i * 4 + 2]] = tf[idx[i * 4 + 3]] = w;
 161                         TransformMul(tf0, tf0, tf);
 162                         tf[idx[i * 4]] = tf[idx[i * 4 + 1]] = 1;
 163                         tf[idx[i * 4 + 2]] = tf[idx[i * 4 + 3]] = 0;
 164                 }
 165                 memmove(tf, tf0, sizeof(tf));
 166                 w0 = (md_rand() - 0.5) * pressure->trial_width;
 167                 for (i = 0; i < 3; i++) {
 168                         int j;
 169                         flex = pressure->cell_flexibility[i];
 170                         if (flex > 0) {
 171                                 w = w0 = (md_rand() - 0.5) * pressure->trial_width;
 172                         } else if (flex < 0) {
 173                                 w = w0;
 174                         } else w = 0;
 175                         for (j = 0; j < 3; j++) {
 176                                 tf[i * 3 + j] *= 1 + w * fabs(flex);
 177                         }
 178                 }
 179
 180                 /*  The new cell vector (bij) = (tik)(akj), so B = A*T  */
 181                 /*  The cell transform matrix R = A*T*(A^-1)  */
 182                 MatrixMul(tf, cell->tr, tf);
 183                 MatrixMul(tf, tf, cell->rtr);
 184         /*      MatrixInvert(mat, arena->celltr);
 185                 MatrixMul(tf, tf, mat); */
 186
 187                 cella_save = cell->axes[0];
 188                 cellb_save = cell->axes[1];
 189                 cellc_save = cell->axes[2];
 190                 cello_save = cell->origin;
 191                 memmove(celltr_save, cell->tr, sizeof(celltr_save));
 192
 193                 /*  Save a snapshot  */
 194                 md_snapshot(arena, 0);
 195                 total_energy_save = arena->total_energy;
 196                 memmove(energies_save, arena->energies, sizeof(Double) * kEndIndex);
 197
 198                 /*  Transform the cell  */
 199                 md_scale_cell(arena, tf, (pressure->use_atomic_move ? 1 : 2));
 200
 201                 /*  Fluctuate cell origin and orientation if requested  */
 202                 if (pressure->cell_flexibility[6] != 0) {
 203                         cello_new = cell->origin;
 204                         cello_new.x += (md_rand() - 0.5) * pressure->fluctuate_cell_origin;
 205                         cello_new.y += (md_rand() - 0.5) * pressure->fluctuate_cell_origin;
 206                         cello_new.z += (md_rand() - 0.5) * pressure->fluctuate_cell_origin;
 207                         cell->origin = cello_new;
 208                         needs_cell_recalculate = 1;
 209                 }
 210                 if (pressure->cell_flexibility[7] != 0) {
 211                         Vector axis;
 212                         Mat33 mat2;
 213                         switch ((int)pressure->cell_flexibility[7]) {
 214                                 case 1: axis.x = 1.0; axis.y = axis.z = 0.0; break;
 215                                 case 2: axis.y = 1.0; axis.x = axis.z = 0.0; break;
 216                                 case 3: axis.z = 1.0; axis.x = axis.y = 0.0; break;
 217                                 case 4:
 218                                         axis.x = md_rand() - 0.5;
 219                                         axis.y = md_rand() - 0.5;
 220                                         axis.z = md_rand() - 0.5;
 221                                         w = 1.0 / VecLength(axis);
 222                                         VecScaleSelf(axis, w);
 223                                         break;
 224                                 default:
 225                                         axis.x = axis.y = axis.z = 0.0; /* Just to keep compiler happy */
 226                                         break;
 227                         }
 228                         MatrixRotation(mat2, &axis, (md_rand() - 0.5) * pressure->fluctuate_cell_orientation);
 229                         MatrixVec(&cell->axes[0], mat2, &cell->axes[0]);
 230                         MatrixVec(&cell->axes[1], mat2, &cell->axes[1]);
 231                         MatrixVec(&cell->axes[2], mat2, &cell->axes[2]);
 232                         needs_cell_recalculate = 1;
 233                 }
 234                 if (needs_cell_recalculate)
 235                         md_update_cell(arena);
 236
 237                 /*  Recalc the energies  */
 238                 md_amend_by_symmetry(arena);
 239                 calc_force(arena);
 240                 md_calc_kinetic_energy(arena);
 241                 pressure->mc_delta_potential = arena->total_energy - total_energy_save;
 242                 pressure->mc_delta_kinetic = arena->energies[kKineticIndex] - energies_save[kKineticIndex];
 243                 pressure->mc_stress = -s_lattice_deformation_stress(pressure->apply, cell->tr, celltr_save);
 244                 de = pressure->mc_delta_potential + pressure->mc_delta_kinetic + pressure->mc_stress;
 245                 pressure->mc_delta_volume = fabs(MatrixDeterminant(cell->tr)) - fabs(MatrixDeterminant(celltr_save));
 246
 247                 /*  Metropolis scheme  */
 248                 if (de > 0 && md_rand() > exp(-de / (BOLTZMANN * arena->temperature))) {
 249                         /*  Reject  */
 250                         cell->axes[0] = cella_save;
 251                         cell->axes[1] = cellb_save;
 252                         cell->axes[2] = cellc_save;
 253                         cell->origin = cello_save;
 254                         md_update_cell(arena);
 255                         if (pressure->control_algorithm != 1) {
 256                                 md_restore(arena, 0);
 257                                 arena->total_energy = total_energy_save;
 258                                 memmove(arena->energies, energies_save, sizeof(Double) * kEndIndex);
 259                         }
 260                         pressure->mc_reject++;
 261                 } else {
 262                         pressure->mc_accept++;
 263                         if (pressure->control_algorithm == 1) {
 264                                 Vector dr;
 265                                 int j, k;
 266                                 /*  Update positions and velocities  */
 267                                 for (i = 0, ap = arena->mol->atoms; i < arena->mol->natoms; i++, ap++) {
 268                                         if (ap->periodic_exclude)
 269                                                 continue;
 270                                         if (i < arena->natoms_uniq) {
 271                                                 k = arena->fragment_indices[i];
 272                                                 v = pressure->temporary_velocities[i];
 273                                         } else if (ap->symop.alive && (j = ap->symbase) >= 0 && j < arena->natoms_uniq) {
 274                                                 k = arena->fragment_indices[j];
 275                                                 v = pressure->temporary_velocities[j];
 276                                         } else continue;
 277                                         dr = arena->fragment_info[k].pos;
 278                                         if (ap->symop.alive) {
 279                                                 md_transform_vec_by_symmetry(arena, &dr, &dr, ap->symop, 1);
 280                                                 md_transform_vec_by_symmetry(arena, &v, &v, ap->symop, 1);
 281                                         }
 282                                         VecInc(ap->r, dr);
 283                                         VecInc(arena->verlets_dr[i], dr);
 284                                         ap->v = v;
 285                                 }
 286                         }
 287                 }
 288         /*      if (pressure->mc_callback != NULL) {
 289                         if (Tcl_EvalObjEx((Tcl_Interp *)(arena->tcl_interp), pressure->mc_callback, 0) != TCL_OK) {
 290                                 arena->request_abort = 1;
 291                                 return;
 292                         }
 293                 } */
 294         } /* end while (mrand() < pressure->coupling) */
 295 }
 296