MolLib/MD/MDGraphite.c

   1 /*
   2  *  MDGraphite.c
   3  *  Molby
   4  *
   5  *  Created by Toshi Nagata on 07/10/07.
   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 "MDGraphite.h"
  19
  20 #include <stdio.h>
  21 #include <stdlib.h>
  22 #include <math.h>
  23 #include <string.h>
  24
  25 /*  Table entry  */
  26 /*  The table should be created for each type of atom of different vdW parameters.
  27         So we need to keep the vdW parameters and the table address.  */
  28 typedef struct MDGraphiteTable {
  29         Double  A, B;    /*  A = 4*pow(r0, 12)*eps, B = 4*pow(r0, 6)*eps  */
  30         Double  offset;  /*  Offset for cutoff truncation */
  31         Double  cutoff2; /*  Square of the cutoff distance  */
  32         Double *cache;  /*  Double[Sx*Sy*Sz*4]  */
  33 } MDGraphiteTable;
  34
  35 struct MDGraphiteArena {
  36         Int natoms_uniq;   /*  Must be the same as arena->natoms_uniq; otherwise, all internal information is rebuilt. */
  37
  38         Int needs_update;   /* If non-zero, then the internal information is rebuilt at the next force calculation  */
  39         Int *table_indices; /* Index to tables[] for each symmetry-unique atom  */
  40         Double *energies;    /* Cache energies of symmetry-unique atoms */
  41         Double last_energy;
  42         Vector *last_forces;
  43         Int ntables;
  44         MDGraphiteTable *tables;
  45         Vector origin;               /*  Origin of the graphite plane in cartesian coordinates */
  46         Vector xaxis, yaxis, zaxis;  /*  The axes of the graphite plane in cartesian coordinates */
  47         Mat33 rcell;                 /*  Convert cartesian to graphite internal coordinates  */
  48
  49         Double R;  /* C-C bond length. Default: 1.23  */
  50         Double C_eps;  /* LJ epsilon for graphite C atoms. Default: 0.0860 * KCAL2INTERNAL  */
  51         Double C_sig;  /* LJ sigma for graphite C atoms. Default: 3.3997  */
  52
  53 /*  The energies and forces at selected points in the asymmetric unit is
  54         cached in graphite->table[]. Sx, Sy, Sz and Sz0 specifies how these points
  55         are selected.
  56         The asymmetric unit is a rectangle (0, 0)-(0.75R, sqrt(3)*0.5R). This rectangle
  57         area (including borders) is split into an (Sx * Sy) lattice, and the lattice
  58         points (i, j) are selected as sample points. For each point, the z-coordinate
  59         in the range [Z_min, Z_lim] is sampled, so that there are Sz0 points span
  60         [Z_min, Z_switch) and (Sz - Sz0) points span [Z_switch, Z_lim]. Thus, the
  61         cartesian coordinates of a sample point (i, j, k) are:
  62                 x = 0.75*R*i/(Sx-1),
  63                 y = sqrt(3)*0.5R*j/(Sy-1),
  64                 z = (k >= Sz0 ? Z_switch+(Z_lim-Z_switch)*(k-Sz0)/(Sz-Sz0-1)
  65                                           : Z_min+(Z_switch-Z_min)*k/Sz0)
  66         The cached value can be accessed by:
  67                 table[((i*Sy + j)*Sz + k)*4 + n]
  68         where n is 0, 1, 2, 3 for the energy, force.x, force.y, force.z, respectively.
  69 */
  70         Int use_table;
  71         Int Sx;  /*  Default: 16  */
  72         Int Sy;  /*  Default: 16  */
  73         Int Sz;  /*  Default: 64  */
  74         Int Sz0; /*  Default: 52  */
  75         Double Z_min;     /*  Default: 1.0  */
  76         Double Z_switch;  /*  Default: 5.0  */
  77         Double Z_lim;     /*  Default: 10.0  */
  78
  79         /*  If periodic_average > 0 and periodic boundary conditions are enabled, then
  80                 the graphite potential will be averaged for periodic_average number of image
  81                 atoms for each periodic cell axes. */
  82         Int periodic_average;
  83
  84 #if DEBUG
  85         FILE *debug_fp;
  86 #endif
  87
  88 };
  89
  90 /*  A = 4*pow(r0, 12)*eps, B = 4*pow(r0, 6)*eps  */
  91 static Double
  92 s_graphite_lj(Vector r, const MDGraphiteTable *tr, Vector *force)
  93 {
  94         Double r2 = r.x*r.x + r.y*r.y + r.z*r.z;
  95         Double w2, w6, w12, k0, k1;
  96         if (r2 >= tr->cutoff2) {
  97                 force->x = force->y = force->z = 0.0;
  98                 return 0.0;
  99         }
 100         w2 = 1.0 / r2;
 101         w6 = w2 * w2 * w2;
 102         w12 = w6 * w6;
 103         k0 = tr->A * w12 - tr->B * w6 + tr->offset;
 104         k1 = 6 * w2 * (2.0 * tr->A * w12 - tr->B * w6);
 105         if (force != NULL) {
 106                 force->x = r.x * k1;
 107                 force->y = r.y * k1;
 108                 force->z = r.z * k1;
 109         }
 110         return k0;
 111 }
 112
 113 /*  Calculate the total potential and force between an atom and the graphite sheet.
 114         The graphite sheet is assumed as a perfect hexagonal array of carbon atoms, and
 115         only van der Waals force is considered.
 116         The origin is at the center of a hexagon, and all carbon atoms lie on the xy plane.
 117         The x-axis points from the origin to one carbon atom. The carbon atoms can be
 118         indexed as:
 119                 (0.5R*(3*i+1), sqrt(3)*0.5R*j), where i+j is odd
 120         This function calculates the vdW potential and force from the graphite
 121         atoms in the range (i,j)=(-10,-15) to (+10,+15).
 122         It is expected that the target atom is close to the z-axis. In the present code,
 123         this function is only invoked for the position in the asymmetric unit, (0,0)..
 124         (1.5R, sqrt(3)*0.5R).
 125 */
 126 static Double
 127 s_graphite_asymunit(MDGraphiteArena *graphite, Vector r, const MDGraphiteTable *tr, Vector *force)
 128 {
 129         Int i, j;
 130         Double e0, en;
 131         Vector r0, f0;
 132         force->x = force->y = force->z = 0.0;
 133         en = 0.0;
 134         for (i = -10; i <= 10; i++) {
 135                 for (j = -15; j <= 15; j++) {
 136                         if ((i + j) % 2 == 0)
 137                                 continue;
 138                         r0.x = r.x - 0.5 * graphite->R * (3 * i + 1);
 139                         r0.y = r.y - graphite->R * 0.866025403784439 * j;
 140                         r0.z = r.z;
 141                         e0 = s_graphite_lj(r0, tr, &f0);
 142                         VecInc(*force, f0);
 143                         en += e0;
 144                         r0.x = r.x - 0.5 * graphite->R * (3 * i - 1);
 145                         e0 = s_graphite_lj(r0, tr, &f0);
 146                         VecInc(*force, f0);
 147                         en += e0;
 148                 }
 149         }
 150         return en;
 151 }
 152
 153 /*  Nearest graphite atoms; x, y coordinates should be multiplied by
 154 R/2 and (sqrt(3)/2)R, respectively  */
 155 /*
 156 static Int sNearestTable[] = {
 157         -1, -3,  1, -3,
 158         -4, -2, -2, -2,  2, -2,  4, -2,
 159         -5, -1, -1, -1,  1, -1,  5, -1,
 160         -4,  0, -2,  0,  2,  0,  4,  0,
 161         -5,  1, -1,  1,  1,  1,  5,  1,
 162         -4,  2, -2,  2,  2,  2,  4,  2,
 163         -1,  3,  1,  3
 164 }; */
 165 static Int sNearestTable[] = {
 166         -2, -4,  2, -4,
 167         -5, -3, -1, -3,  1, -3,  5, -3,
 168         -4, -2, -2, -2,  2, -2,  4, -2,
 169         -7, -1, -5, -1, -1, -1,  1, -1,  5, -1,  7, -1,
 170         -8,  0, -4,  0, -2,  0,  2,  0,  4,  0,  8,  0,
 171         -7,  1, -5,  1, -1,  1,  1,  1,  5,  1,  7,  1,
 172         -4,  2, -2,  2,  2,  2,  4,  2,
 173         -5,  3, -1,  3,  1,  3,  5,  3,
 174         -2,  4,  2, -4
 175 };
 176 #define kSizeOfNearestTable (sizeof(sNearestTable)/sizeof(sNearestTable[0])/2)
 177
 178 /*  The nearest 24 graphite atoms are always explicitly calculated  */
 179 Double
 180 s_graphite_nearest_asymunit(MDGraphiteArena *graphite, Vector r, const MDGraphiteTable *tr, Vector *force)
 181 {
 182         Int i;
 183         Double en;
 184         Vector q, dr, f;
 185         en = force->x = force->y = force->z = 0.0;
 186         for (i = 0; i < kSizeOfNearestTable; i++) {
 187                 q.x = sNearestTable[i*2] * graphite->R * 0.5;
 188                 q.y = sNearestTable[i*2+1] * 0.866025403784439 * graphite->R;
 189                 dr.x = r.x - q.x;
 190                 dr.y = r.y - q.y;
 191                 dr.z = r.z;
 192                 en += s_graphite_lj(dr, tr, &f);
 193                 force->x += f.x;
 194                 force->y += f.y;
 195                 force->z += f.z;
 196         }
 197         return en;
 198 }
 199
 200 static void
 201 s_graphite_make_table(MDGraphiteArena *graphite, MDGraphiteTable *tr)
 202 {
 203         Int i, j, k;
 204         Vector v, f, f0;
 205         Double en, en0, *fp;
 206         k = graphite->Sx * graphite->Sy * graphite->Sz * sizeof(Double) * 4;
 207         if (tr->cache == NULL)
 208                 tr->cache = (Double *)malloc(k);
 209         else
 210                 tr->cache = (Double *)realloc(tr->cache, k);
 211         memset(tr->cache, 0, k);
 212         for (k = 0; k < graphite->Sz; k++) {
 213                 if (k >= graphite->Sz0)
 214                         v.z = graphite->Z_switch + (graphite->Z_lim - graphite->Z_switch) * (k - graphite->Sz0) / (graphite->Sz - graphite->Sz0 - 1);
 215                 else
 216                         v.z = graphite->Z_min + (graphite->Z_switch - graphite->Z_min) * k / graphite->Sz0;  /*  Not (Sz0 - 1)  */
 217                 for (i = 0; i < graphite->Sx; i++) {
 218                         for (j = 0; j < graphite->Sy; j++) {
 219                                 v.x = 0.75 * graphite->R * i / (graphite->Sx - 1.0);
 220                                 v.y = 0.866025403784439 * graphite->R * j / (graphite->Sy - 1.0);
 221                                 en = s_graphite_asymunit(graphite, v, tr, &f);
 222                                 en0 = s_graphite_nearest_asymunit(graphite, v, tr, &f0);
 223                                 fp = tr->cache + ((i * graphite->Sy + j) * graphite->Sz + k) * 4;
 224                                 *fp++ = en - en0;
 225                                 *fp++ = f.x - f0.x;
 226                                 *fp++ = f.y - f0.y;
 227                                 *fp++ = f.z - f0.z;
 228                         }
 229                 }
 230         }
 231 }
 232
 233 static inline Double
 234 s_linear_interpolate(MDGraphiteArena *graphite, Double *base, Double dx, Double dy, Double dz)
 235 {
 236         Int Sy = graphite->Sy;
 237         Int Sz = graphite->Sz;
 238         Double a0 = base[0];
 239         Double a1 = base[Sy * Sz * 4] - a0;
 240         Double a2 = base[Sz * 4] - a0;
 241         Double a3 = base[4] - a0;
 242         Double a4 = base[(Sy + 1) * Sz * 4] - a0 - a1 - a2;
 243         Double a5 = base[(Sz + 1) * 4] - a0 - a2 - a3;
 244         Double a6 = base[(Sy * Sz + 1) * 4] - a0 - a1 - a3;
 245         Double a7 = base[((Sy + 1) * Sz + 1) * 4] - a0 - a1 - a2 - a3 - a4 - a5 - a6;
 246         return a0 + a1 * dx + a2 * dy + a3 * dz + a4 * dx * dy + a5 * dy * dz + a6 * dz * dx + a7 * dx * dy * dz;
 247 }
 248
 249 static Double
 250 s_graphite(MDGraphiteArena *graphite, Vector v, MDGraphiteTable *tr, Vector *f)
 251 {
 252         Int i, j, k;
 253         Double cella, cellb;
 254         Double en, dx, dy, dz, *fp;
 255         Int ix, iy, negx, negy, negz, symop;
 256         Vector v0, f0;
 257
 258         if (v.z < 0) {
 259                 negz = 1;
 260                 v.z = -v.z;
 261         } else negz = 0;
 262         if (v.z > graphite->Z_lim) {
 263                 f->x = f->y = f->z = 0.0;
 264                 return 0.0;
 265         }
 266
 267         cella = 0.75 * graphite->R;
 268         cellb = 0.866025403784439 * graphite->R;
 269         dx = v.x / (cella * 4);
 270         ix = floor(dx + 0.5);
 271         dx -= ix;
 272         dy = v.y / (cellb * 4);
 273         iy = floor(dy + 0.5);
 274         dy -= iy;
 275
 276         /*  Move vector to the asymmetric unit (0,0)-(0.25,0.25) */
 277         if (dx < 0) {
 278                 dx = -dx;
 279                 negx = 1;
 280         } else negx = 0;
 281         if (dy < 0) {
 282                 dy = -dy;
 283                 negy = 1;
 284         } else negy = 0;
 285         if (dy > 0.25) {
 286                 if (dx > 0.25) {
 287                         dx = 0.5 - dx;
 288                         dy -= 0.25;
 289                         symop = 3;
 290                 } else {
 291                         dy = 0.5 - dy;
 292                         symop = 2;
 293                 }
 294         } else if (dx > 0.25) {
 295                 dx = 0.5 - dx;
 296                 dy = 0.25 - dy;
 297                 symop = 1;
 298         } else symop = 0;
 299         v0.x = dx * cella * 4;
 300         v0.y = dy * cellb * 4;
 301         v0.z = v.z;
 302
 303         if (tr->cache != NULL) {
 304                 /* 0 <= dx <= 0.25, 0 <= dy <= 0.25, 1.0 <= v.z <= 10.0 */
 305                 /* 0 <= i <= Sx-1, 0 <= j <= Sy-1, 0 <= k <= Sz-1 */
 306                 dx = dx * (graphite->Sx - 1.0) / 0.25;
 307                 dy = dy * (graphite->Sy - 1.0) / 0.25;
 308                 if (v.z >= graphite->Z_switch)
 309                         dz = (v.z - graphite->Z_switch) * (graphite->Sz - graphite->Sz0 - 1) / (graphite->Z_lim - graphite->Z_switch) + graphite->Sz0;
 310                 else
 311                         dz = (v.z - graphite->Z_min) * graphite->Sz0 / (graphite->Z_switch - graphite->Z_min);
 312                 i = floor(dx);
 313                 dx -= i;
 314                 j = floor(dy);
 315                 dy -= j;
 316                 k = floor(dz);
 317                 dz -= k;
 318                 if (i >= graphite->Sx - 1) {
 319                         i = graphite->Sx - 2;
 320                         dx = 1.0;
 321                 }
 322                 if (j >= graphite->Sy - 1) {
 323                         j = graphite->Sy - 2;
 324                         dy = 1.0;
 325                 }
 326                 if (k < 0) {
 327                         k = 0;
 328                         dz = 0.0;
 329                 } else if (k >= graphite->Sz - 1) {
 330                         k = graphite->Sz - 2;
 331                         dz = 1.0;
 332                 }
 333                 fp = tr->cache + ((i*graphite->Sy + j)*graphite->Sz + k)*4;
 334                 en = s_linear_interpolate(graphite, fp, dx, dy, dz);
 335                 f0.x = s_linear_interpolate(graphite, fp + 1, dx, dy, dz);
 336                 f0.y = s_linear_interpolate(graphite, fp + 2, dx, dy, dz);
 337                 f0.z = s_linear_interpolate(graphite, fp + 3, dx, dy, dz);
 338                 en += s_graphite_nearest_asymunit(graphite, v0, tr, f);
 339                 f->x += f0.x;
 340                 f->y += f0.y;
 341                 f->z += f0.z;
 342         } else {
 343                 en = s_graphite_asymunit(graphite, v0, tr, f);
 344         }
 345
 346         /*  Back transform to original axis system  */
 347         switch (symop) {
 348                 case 3: negx = !negx; break;
 349                 case 2: negy = !negy; break;
 350                 case 1: negx = !negx; negy = !negy; break;
 351                 default: break;
 352         }
 353         if (negx)
 354                 f->x = -f->x;
 355         if (negy)
 356                 f->y = -f->y;
 357         if (negz)
 358                 f->z = -f->z;
 359         return en;
 360 }
 361
 362 MDGraphiteArena *
 363 graphite_new(void)
 364 {
 365         static Vector x = {1, 0, 0}, y = {0, 1, 0}, z = {0, 0, 1};
 366         static Mat33 u = {1, 0, 0, 0, 1, 0, 0, 0, 1};
 367         MDGraphiteArena *graphite = (MDGraphiteArena *)calloc(sizeof(MDGraphiteArena), 1);
 368         graphite->use_table = 1;
 369         graphite->Sx = 16;
 370         graphite->Sy = 16;
 371         graphite->Sz = 64;
 372         graphite->Sz0 = 52;
 373         graphite->Z_min = 1.0;
 374         graphite->Z_switch = 5.0;
 375         graphite->Z_lim = 10.0;
 376         graphite->R = 1.23;
 377         graphite->C_eps = 0.0860 * KCAL2INTERNAL;
 378         graphite->C_sig = 3.3997;
 379         graphite->xaxis = x;
 380         graphite->yaxis = y;
 381         graphite->zaxis = z;
 382         memmove(graphite->rcell, u, sizeof(Mat33));
 383         return graphite;
 384 }
 385
 386 static void
 387 s_graphite_release_tables(MDGraphiteArena *graphite)
 388 {
 389         int i;
 390         if (graphite == NULL || graphite->ntables == 0 || graphite->tables == NULL)
 391                 return;
 392         for (i = 0; i < graphite->ntables; i++) {
 393                 if (graphite->tables[i].cache != NULL)
 394                         free(graphite->tables[i].cache);
 395         }
 396         free(graphite->tables);
 397         graphite->tables = NULL;
 398         graphite->ntables = 0;
 399 }
 400
 401 void
 402 graphite_release(MDGraphiteArena *graphite)
 403 {
 404         if (graphite != NULL) {
 405                 if (graphite->table_indices != NULL)
 406                         free(graphite->table_indices);
 407                 if (graphite->tables != NULL)
 408                         s_graphite_release_tables(graphite);
 409         }
 410         free(graphite);
 411 }
 412
 413 void
 414 graphite_set_origin(MDGraphiteArena *graphite, const Vector *vp)
 415 {
 416         if (graphite == NULL)
 417                 return;
 418         graphite->origin = *vp;
 419 }
 420
 421 /*  (Re)initialize the MDGraphiteArena  */
 422 static void
 423 s_graphite_init(MDGraphiteArena *graphite, MDArena *arena)
 424 {
 425         MDGraphiteTable *tp;
 426         Atom *ap;
 427         int nuniq = arena->natoms_uniq;
 428         int i, j;
 429
 430         /*  Reallocate the internal storage  */
 431         if (graphite->table_indices == NULL)
 432                 graphite->table_indices = (Int *)malloc(sizeof(Int) * nuniq);
 433         else
 434                 graphite->table_indices = (Int *)realloc(graphite->table_indices, sizeof(Int) * nuniq);
 435         if (graphite->energies == NULL)
 436                 graphite->energies = (Double *)malloc(sizeof(Double) * nuniq);
 437         else
 438                 graphite->energies = (Double *)realloc(graphite->energies, sizeof(Double) * nuniq);
 439         if (graphite->last_forces == NULL)
 440                 graphite->last_forces = (Vector *)malloc(sizeof(Vector) * nuniq);
 441         else
 442                 graphite->last_forces = (Vector *)realloc(graphite->last_forces, sizeof(Vector) * nuniq);
 443         if (graphite->tables != NULL)
 444                 s_graphite_release_tables(graphite);
 445
 446         /*  Calculate the tables  */
 447         md_log(arena, "Building the graphite table...\n");
 448
 449         for (i = 0, ap = arena->mol->atoms; i < nuniq; i++, ap++) {
 450                 Int vdw_idx;
 451                 VdwPar *vp;
 452                 Double sig, eps, sigc, epsc, A, B;
 453                 vdw_idx = arena->vdw_par_i[i];
 454                 if (vdw_idx < 0 || vdw_idx >= arena->par->nvdwPars) {
 455                         /*  Out of range  */
 456                         graphite->table_indices[i] = -1;
 457                         continue;
 458                 }
 459                 vp = &(arena->par->vdwPars[vdw_idx]);
 460                 if (fabs(vp->A) < 1e-15) {
 461                         sig = 1.0;
 462                         eps = 0.0;
 463                 } else {
 464                         eps = 0.25 * vp->B * vp->B / vp->A;
 465                         sig = pow(vp->B * 0.25 / eps, 1.0/6.0);
 466                 }
 467                 sigc = (graphite->C_sig + sig) * 0.5;
 468                 epsc = sqrt(graphite->C_eps * eps);
 469                 A =     4.0 * pow(sigc, 12.0) * epsc;
 470                 B = 4.0 * pow(sigc, 6.0) * epsc;
 471                 /*  Are there entries with similar A/B values?  */
 472                 for (j = 0; j < graphite->ntables; j++) {
 473                         tp = graphite->tables + j;
 474                         if (fabs(tp->A - A) < 1e-10 && fabs(tp->B - B) < 1e-10)
 475                                 break;
 476                 }
 477                 if (j >= graphite->ntables) {
 478                         /*  Create a new entry with a table  */
 479                         Vector v;
 480                         j = graphite->ntables;
 481                         tp = (MDGraphiteTable *)AssignArray(&graphite->tables, &graphite->ntables, sizeof(MDGraphiteTable), j, NULL);
 482                         tp->A = A;
 483                         tp->B = B;
 484                         /*  Calculate the offset: set 0 to offset temporarily, and call s_graphite_lj */
 485                         tp->cutoff2 = arena->cutoff * arena->cutoff + 1.0;  /* +1.0 to avoid cutoff */
 486                         tp->offset = 0.0;
 487                         v.x = arena->cutoff;
 488                         v.y = v.z = 0.0;
 489                         tp->offset = -s_graphite_lj(v, tp, &v);
 490                         tp->cutoff2 = arena->cutoff * arena->cutoff; /* This is the real value */
 491                         if (graphite->use_table)
 492                                 s_graphite_make_table(graphite, tp);
 493                         else tp->cache = NULL;
 494                 }
 495                 graphite->table_indices[i] = j;
 496         }
 497         graphite->natoms_uniq = arena->natoms_uniq;
 498         md_log(arena, "%d table(s) were created %s cached data.\n", graphite->ntables, (graphite->use_table ? "with" : "without"));
 499         graphite->needs_update = 0;
 500
 501 #if DEBUG
 502         if (arena->debug_result && arena->debug_output_level > 1) {
 503                 graphite->debug_fp = arena->debug_result;
 504         } else graphite->debug_fp = NULL;
 505 #endif
 506
 507 }
 508
 509 void
 510 graphite_force(MDGraphiteArena *graphite, MDArena *arena, Double *energy, Vector *forces)
 511 {
 512         Molecule *mol = arena->mol;
 513         int i;
 514         int ix, iy, iz, px, py, pz;
 515         Double pxpypz_1;
 516         Vector *avp, *bvp, *cvp;
 517         if (graphite == NULL || mol == NULL || mol->natoms == 0 || arena->natoms_uniq == 0)
 518                 return;
 519         if (graphite->natoms_uniq != arena->natoms_uniq || graphite->needs_update)
 520                 s_graphite_init(graphite, arena);
 521         graphite->last_energy = 0;
 522         memset(graphite->last_forces, 0, sizeof(Vector) * arena->natoms_uniq);
 523         if (graphite->periodic_average && arena->mol->cell != NULL) {
 524                 px = (arena->periodic_a ? graphite->periodic_average : 1);
 525                 py = (arena->periodic_b ? graphite->periodic_average : 1);
 526                 pz = (arena->periodic_c ? graphite->periodic_average : 1);
 527         } else {
 528                 px = py = pz = 1;
 529         }
 530         pxpypz_1 = 1.0 / (px * py * pz);
 531         if (mol->cell != NULL) {
 532                 avp = &(mol->cell->axes[0]);
 533                 bvp = &(mol->cell->axes[1]);
 534                 cvp = &(mol->cell->axes[2]);
 535         }
 536         for (i = 0; i < arena->natoms_uniq; i++) {
 537                 Vector f;
 538                 Double en;
 539                 Vector v = mol->atoms[i].r;
 540                 if (mol->atoms[i].fix_force < 0)
 541                         continue;
 542                 if (mol->atoms[i].occupancy == 0.0)
 543                         continue;
 544                 en = f.x = f.y = f.z = 0.0;
 545                 for (ix = 0; ix < px; ix++) {
 546                         for (iy = 0; iy < py; iy++) {
 547                                 for (iz = 0; iz < pz; iz++) {
 548                                         Vector f0, v0;
 549                                         Double en0;
 550                                         if (mol->cell != NULL) {
 551                                                 v0.x = v.x + avp->x * ix + bvp->x * iy + cvp->x * iz;
 552                                                 v0.y = v.y + avp->y * ix + bvp->y * iy + cvp->y * iz;
 553                                                 v0.z = v.z + avp->z * ix + bvp->z * iy + cvp->z * iz;
 554                                         } else v0 = v;
 555                                         VecDec(v0, graphite->origin);
 556                                         MatrixVec(&v0, graphite->rcell, &v0);
 557                                         en0 = s_graphite(graphite, v0, graphite->tables + graphite->table_indices[i], &f0);
 558                                         en += en0;
 559                                         VecInc(f, f0);
 560                                 }
 561                         }
 562                 }
 563                 en *= pxpypz_1;
 564                 VecScaleSelf(f, pxpypz_1);
 565                 graphite->energies[i] = en;
 566                 graphite->last_energy += en;
 567                 graphite->last_forces[i] = f;
 568         }
 569         for (i = arena->natoms_uniq; i < mol->natoms; i++) {
 570                 Symop symop = mol->atoms[i].symop;
 571                 if (mol->atoms[i].fix_force < 0)
 572                         continue;
 573                 if (symop.alive) {
 574                         int n = mol->atoms[i].symbase;
 575                         if (n >= 0 && n < arena->natoms_uniq)
 576                                 graphite->last_energy += graphite->energies[n];
 577                 }
 578         }
 579         if (energy != NULL)
 580                 *energy += graphite->last_energy;
 581         if (forces != NULL) {
 582                 for (i = 0; i < arena->natoms_uniq; i++)
 583                         VecInc(forces[i], graphite->last_forces[i]);
 584         }
 585 }
 586
 587 void
 588 graphite_get_axes(MDGraphiteArena *graphite, Vector *op, Vector *xp, Vector *yp, Vector *zp, Double *rp)
 589 {
 590         if (graphite != NULL) {
 591                 if (op != NULL)
 592                         *op = graphite->origin;
 593                 if (xp != NULL)
 594                         *xp = graphite->xaxis;
 595                 if (yp != NULL)
 596                         *yp = graphite->yaxis;
 597                 if (zp != NULL)
 598                         *zp = graphite->zaxis;
 599                 if (rp != NULL)
 600                         *rp = graphite->R;
 601         }
 602 }
 603