MolLib/MD/MDEwald.c

   1 /*
   2  *  MDEwald.c
   3  *
   4  *  Created by Toshi Nagata on 2012/11/03.
   5  *  Copyright 2012 Toshi Nagata. All rights reserved.
   6  *
   7  This program is free software; you can redistribute it and/or modify
   8  it under the terms of the GNU General Public License as published by
   9  the Free Software Foundation version 2 of the License.
  10
  11  This program is distributed in the hope that it will be useful,
  12  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  GNU General Public License for more details.
  15  */
  16
  17 #include "MDCore.h"
  18 #include "MDEwald.h"
  19
  20 #include <stdlib.h>
  21 #include <string.h>
  22 #include <math.h>
  23
  24 #include <fftw3.h>
  25
  26 /*  Temporary structure for Particle Mesh Ewald (Some fields are also used for Direct Ewald)  */
  27 struct MDPME {
  28         Int grid_x, grid_y, grid_z;  /*  Grid size for each cell direction (should be even numbers)  */
  29         Int natoms;                  /*  Number of atoms (necessary to indicate the size of allocated memory)  */
  30         fftw_complex *qarray;        /*  Q array; size grid_x * grid_y * grid_z  */
  31         fftw_complex *qfarray;       /*  Finv(Q) (Finv is inverse FFT); size grid_x * grid_y * grid_z  */
  32         fftw_complex *qqarray;       /*  F(B*C*Finv(Q)) (F is FFT); size grid_x * grid_y * grid_z  */
  33         fftw_complex *b0array;       /*  Complex B array (for debug); size grid_x + grid_y + grid_z  */
  34         double       *marray;        /*  M array; size (grid_x + grid_y + grid_z) * 2 * natoms  */
  35         double       *barray;        /*  B array; size grid_x + grid_y + grid_z  */
  36         double       *bcarray;       /*  B*C array; size grid_x * grid_y * grid_z  */
  37         double       *carray;        /*  C array (for debug); size grid_x * grid_y * grid_z  */
  38         fftw_plan plan1;             /*  Inverse FFT of Q */
  39         fftw_plan plan2;             /*  FFT of (B*C*Finv(Q))  */
  40         Int       order;             /*  Order of cardinal B spline (must be even and >=4 and <= MAX_DIM_SPLINE )  */
  41         Double direct_limit;         /*  Maximum (k/beta) for direct Ewald summation (default = 1.29; this is where erfc(lim*PI) = 1e-8  */
  42         Int *kxyz;                   /*  Reciprocal lattice indices for direct Ewald summation  */
  43         Int nkxyz;
  44         Double self_correction;      /*  Self correction term, -beta/sqrt(PI)*sigma(charge**2)  */
  45 };
  46
  47 static Double *s_spline_coeffs;
  48
  49 #define SPLINE_COEFF(_n, _m, _k) s_spline_coeffs[(_k) + MAX_DIM_SPLINE * ((_m) + MAX_DIM_SPLINE * ((_n) - 2))]
  50
  51 /*  Calculate the coefficients of the cardinal B-spline functions  */
  52 static void
  53 s_initialize_spline_coeffs(void)
  54 {
  55         Int n, m, k;
  56         if (s_spline_coeffs != NULL && SPLINE_COEFF(2, 0, 1) == 1.0 && SPLINE_COEFF(2, 1, 1) == -1.0)
  57                 return;  /*  Already initialized  */
  58         s_spline_coeffs = (Double *)calloc(sizeof(Double), (MAX_DIM_SPLINE - 1) * MAX_DIM_SPLINE * MAX_DIM_SPLINE);
  59         SPLINE_COEFF(2, 0, 0) = 0.0;
  60         SPLINE_COEFF(2, 0, 1) = 1.0;
  61         SPLINE_COEFF(2, 1, 0) = 2.0;
  62         SPLINE_COEFF(2, 1, 1) = -1.0;
  63         for (n = 3; n <= MAX_DIM_SPLINE; n++) {
  64                 for (m = 0; m < n; m++) {
  65                         Double a, b, b_prev;
  66                         b_prev = 0.0;
  67                         for (k = 0; k < n; k++) {
  68                                 Int sgn, j, bin;
  69                                 b = 0.0;
  70                                 bin = 1;
  71                                 sgn = 1;
  72                                 if (m > 0) {
  73                                         for (j = k; j < n - 1; j++) {
  74                                                 b += SPLINE_COEFF(n - 1, m - 1, j) * bin * sgn;
  75                                                 sgn *= -1;
  76                                                 bin = bin * (j + 1) / (j + 1 - k); /* Always dividable; this is a binomial coefficient (j+1 k)  */
  77                                         }
  78                                 }
  79                                 if (k > 0 && m < n - 1)
  80                                         a = SPLINE_COEFF(n - 1, m, k - 1);
  81                                 else a = 0.0;
  82                                 a = (a + n * b - b_prev) / (n - 1);
  83                                 b_prev = b;
  84                                 SPLINE_COEFF(n, m, k) = a;
  85                         }
  86                 }
  87         }
  88 }
  89
  90 static Double
  91 s_calc_spline(Int order, Double x)
  92 {
  93         Int m, i;
  94         Double y;
  95         if (x <= 0 || x >= order)
  96                 return 0.0;
  97         m = floor(x);
  98         y = SPLINE_COEFF(order, m, order - 1);
  99         for (i = order - 2; i >= 0; i--)
 100                 y = y * x + SPLINE_COEFF(order, m, i);
 101         return y;
 102 }
 103
 104 inline static void
 105 s_complex_add(fftw_complex dst, fftw_complex src1, fftw_complex src2)
 106 {
 107         dst[0] = src1[0] + src2[0];
 108         dst[1] = src1[1] + src2[1];
 109 }
 110
 111 inline static void
 112 s_complex_scale_add(fftw_complex dst, fftw_complex src1, fftw_complex src2, double scale)
 113 {
 114         dst[0] = src1[0] + src2[0] * scale;
 115         dst[1] = src1[1] + src2[1] * scale;
 116 }
 117
 118 inline static void
 119 s_complex_mul(fftw_complex dst, fftw_complex src1, fftw_complex src2)
 120 {
 121         double real;
 122         real = src1[0] * src2[0] - src1[1] * src2[1];
 123         dst[1] = src1[0] * src2[1] + src1[1] * src2[0];
 124         dst[0] = real;
 125 }
 126
 127 /*  For debug  */
 128 static void
 129 s_output_array_as_cube(MDArena *arena, fftw_complex *ary, const char *fname)
 130 {
 131         /*  For debug  */
 132         /*  Output array for visualize  */
 133         int i, j, k, n;
 134         double max, min, d;
 135         FILE *fp = fopen(fname, "w");
 136         Molecule *mol = arena->xmol;
 137         XtalCell *cp = mol->cell;
 138         MDPME *pme = arena->pme;
 139
 140         fprintf(fp, "PME test\n");
 141         fprintf(fp, "Output of the array\n");
 142         fprintf(fp, "%5d %11.6f %11.6f %11.6f\n", -(mol->natoms), cp->origin.x, cp->origin.y, cp->origin.z);
 143         fprintf(fp, "%5d %11.6f %11.6f %11.6f\n", pme->grid_x, cp->axes[0].x / pme->grid_x, cp->axes[0].y / pme->grid_x, cp->axes[0].z / pme->grid_x);
 144         fprintf(fp, "%5d %11.6f %11.6f %11.6f\n", pme->grid_y, cp->axes[1].x / pme->grid_y, cp->axes[1].y / pme->grid_y, cp->axes[1].z / pme->grid_y);
 145         fprintf(fp, "%5d %11.6f %11.6f %11.6f\n", pme->grid_z, cp->axes[2].x / pme->grid_z, cp->axes[2].y / pme->grid_z, cp->axes[2].z / pme->grid_z);
 146
 147         /*  Atomic information  */
 148         for (i = 0; i < mol->natoms; i++) {
 149                 Atom *ap = ATOM_AT_INDEX(mol->atoms, i);
 150                 fprintf(fp, "%5d %11.6f %11.6f %11.6f %11.6f\n", ap->atomicNumber, (double)ap->charge, ap->r.x, ap->r.y, ap->r.z);
 151         }
 152         fprintf(fp, "%5d%5d\n", 1, 1);
 153         max = -1e10;
 154         min = 1e10;
 155         for (i = 0; i < pme->grid_x * pme->grid_y * pme->grid_z; i++) {
 156                 d = sqrt(ary[i][0] * ary[i][0] + ary[i][1] * ary[i][1]);
 157                 if (d > max)
 158                         max = d;
 159                 if (d < min)
 160                         min = d;
 161         }
 162         max = (max > -min ? max : -min);
 163         for (i = n = 0; i < pme->grid_x; i++) {
 164                 for (j = 0; j < pme->grid_y; j++) {
 165                         for (k = 0; k < pme->grid_z; k++) {
 166                                 fprintf(fp, " %12.5e", sqrt(ary[n][0] * ary[n][0] + ary[n][1] * ary[n][1]) / max);
 167                                 n++;
 168                                 if (k == pme->grid_z - 1 || k % 6 == 5)
 169                                         fprintf(fp, "\n");
 170                         }
 171                 }
 172         }
 173         fclose(fp);
 174 }
 175
 176 /*  Prepare kxyz for direct Ewald sum  */
 177 void
 178 s_prepare_direct_ewald(MDArena *arena)
 179 {
 180         Vector raxes[3];
 181         XtalCell *cell = arena->mol->cell;
 182         Int pass, count, dc, n, ninit, nstep;
 183         raxes[0].x = cell->rtr[0];
 184         raxes[0].y = cell->rtr[3];
 185         raxes[0].z = cell->rtr[6];
 186         raxes[1].x = cell->rtr[1];
 187         raxes[1].y = cell->rtr[4];
 188         raxes[1].z = cell->rtr[7];
 189         raxes[2].x = cell->rtr[2];
 190         raxes[2].y = cell->rtr[5];
 191         raxes[2].z = cell->rtr[8];
 192         ninit = 3;
 193         nstep = 1;
 194         count = 0;
 195         for (pass = 0; pass < 2; pass++) {
 196                 Int *ip;
 197                 if (pass == 1)
 198                         ip = arena->pme->kxyz;
 199                 for (n = ninit; n <= 50; n += nstep) {
 200                         Int nx, ny, nz, nx0, ny0, nz0;
 201                         Vector vx, vy, vm;
 202                         Double len;
 203                         dc = 0;
 204                         nx0 = (arena->periodic_a ? n : 1);
 205                         ny0 = (arena->periodic_b ? n : 1);
 206                         nz0 = (arena->periodic_c ? n : 1);
 207                         for (nx = -nx0 + 1; nx < nx0; nx++) {
 208                                 VecScale(vx, raxes[0], nx);
 209                                 for (ny = -ny0 + 1; ny < ny0; ny++) {
 210                                         VecScale(vy, raxes[1], ny);
 211                                         for (nz = -nz0 + 1; nz < nz0; nz++) {
 212                                                 if (n > ninit)
 213                                                         if (arena->periodic_a == 0 || (nx > -n + nstep && nx < n - nstep))
 214                                                                 if (arena->periodic_b == 0 || (ny > -n + nstep && ny < n - nstep))
 215                                                                         if (arena->periodic_c == 0 || (nz > -n + nstep && nz < n - nstep))
 216                                                                                 continue;
 217                                                 if (nx == 0 && ny == 0 && nz == 0)
 218                                                         continue;
 219                                                 VecScale(vm, raxes[2], nz);
 220                                                 VecInc(vm, vx);
 221                                                 VecInc(vm, vy);
 222                                                 len = VecLength(vm);
 223                                                 if (len / arena->ewald_beta <= arena->pme->direct_limit) {
 224                                                         if (pass == 0) {
 225                                                                 dc++;
 226                                                                 count++;
 227                                                         } else {
 228                                                                 *ip++ = nx;
 229                                                                 *ip++ = ny;
 230                                                                 *ip++ = nz;
 231                                                                 count--;
 232                                                         }
 233                                                 }
 234                                         }
 235                                 }
 236                         }
 237                         if ((pass == 0 && dc == 0) || (pass == 1 && count <= 0))
 238                                 break;
 239                 }
 240                 if (pass == 0) {
 241                         arena->pme->kxyz = (Int *)realloc(arena->pme->kxyz, sizeof(Int) * 3 * count);
 242                         arena->pme->nkxyz = count;
 243                 }
 244         }
 245 }
 246
 247 /*  For debug  */
 248 void
 249 calc_direct_electrostatic(MDArena *arena)
 250 {
 251         Molecule *mol = arena->mol;
 252         XtalCell *cell = arena->mol->cell;
 253         Vector *forces, *f;
 254         Double energy;
 255         Double coul;
 256         Int i, j, n, nx, ny, nz, nx0, ny0, nz0, zero, ninit, nstep;
 257         Atom *ap, *apj;
 258
 259         forces = &arena->forces[kPMEIndex * arena->mol->natoms];
 260         energy = 0.0;
 261         coul = COULOMBIC / arena->dielectric;
 262         ninit = 3;
 263         nstep = 1;
 264         for (n = ninit; n <= 50; n += nstep) {
 265                 Vector vx, vy, vz;
 266                 Double de;
 267                 de = 0.0;
 268                 nx0 = (arena->periodic_a ? n : 1);
 269                 ny0 = (arena->periodic_b ? n : 1);
 270                 nz0 = (arena->periodic_c ? n : 1);
 271                 for (nx = -nx0 + 1; nx < nx0; nx++) {
 272                         VecScale(vx, cell->axes[0], nx);
 273                         for (ny = -ny0 + 1; ny < ny0; ny++) {
 274                                 VecScale(vy, cell->axes[1], ny);
 275                                 for (nz = -nz0 + 1; nz < nz0; nz++) {
 276                                         if (n > ninit)
 277                                                 if (arena->periodic_a == 0 || (nx > -n + nstep && nx < n - nstep))
 278                                                         if (arena->periodic_b == 0 || (ny > -n + nstep && ny < n - nstep))
 279                                                                 if (arena->periodic_c == 0 || (nz > -n + nstep && nz < n - nstep))
 280                                                                         continue;
 281                                         zero = (nx == 0 && ny == 0 && nz == 0);
 282                                         VecScale(vz, cell->axes[2], nz);
 283                                         for (i = 0, ap = mol->atoms, f = forces; i < mol->natoms; i++, ap = ATOM_NEXT(ap), f++) {
 284                                                 Vector ri, rj, rij;
 285                                                 ri = ap->r;
 286                                                 for (j = 0, apj = mol->atoms; j < mol->natoms; j++, apj = ATOM_NEXT(apj)) {
 287                                                         Double w1, w2, w3, qiqj;
 288                                                         if (i == j && zero)
 289                                                                 continue;
 290                                                         qiqj = ap->charge * apj->charge;
 291                                                         rj = apj->r;
 292                                                         VecSub(rij, ri, rj);
 293                                                         VecInc(rij, vx);
 294                                                         VecInc(rij, vy);
 295                                                         VecInc(rij, vz);
 296                                                         w1 = 1.0 / VecLength2(rij);
 297                                                         w2 = sqrt(w1);
 298                                                         de += w2 * qiqj;
 299                                                         w3 = w1 * w2 * qiqj;
 300                                                         f->x += rij.x * w3;
 301                                                         f->y += rij.y * w3;
 302                                                         f->z += rij.z * w3;
 303                                                 }
 304                                         }
 305                                 }
 306                         }
 307                 }
 308                 de *= 0.5 * coul;
 309                 energy += de;
 310                 if (arena->debug_result != NULL && arena->debug_output_level > 0) {
 311                         fprintf(arena->debug_result, "Direct electrostatic energy\n");
 312                         fprintf(arena->debug_result, "N = %d energy = %.9g %.9g\n", n, energy * INTERNAL2KCAL, de * INTERNAL2KCAL);
 313                         if (arena->debug_output_level > 1) {
 314                                 for (i = 0; i < mol->natoms; i++) {
 315                                         fprintf(arena->debug_result, "  %d: %.9g %.9g %.9g\n", i, forces[i].x * coul * INTERNAL2KCAL, forces[i].y * coul * INTERNAL2KCAL, forces[i].z * coul * INTERNAL2KCAL);
 316                                 }
 317                         }
 318                 }
 319                 if (fabs(de) < 1e-7 || fabs(de / energy) < 1e-6)
 320                         break;
 321         }
 322         if (arena->debug_result != NULL) {
 323                 fprintf(arena->debug_result, "Direct electrostatic energy = %g\n", energy);
 324         }
 325         arena->energies[kPMEIndex] += energy;
 326 }
 327
 328 /*  Direct Ewald Sum  */
 329 void
 330 calc_direct_ewald_force(MDArena *arena)
 331 {
 332         Molecule *mol = arena->mol;
 333         XtalCell *cell = arena->mol->cell;
 334         Vector raxes[3], *forces;
 335         Double energy_rec, coul_v;
 336         Int i, n, nx, ny, nz, *ip;
 337         Atom *ap;
 338
 339         /*  Reciprocal cell axes  */
 340         raxes[0].x = cell->rtr[0];
 341         raxes[0].y = cell->rtr[3];
 342         raxes[0].z = cell->rtr[6];
 343         raxes[1].x = cell->rtr[1];
 344         raxes[1].y = cell->rtr[4];
 345         raxes[1].z = cell->rtr[7];
 346         raxes[2].x = cell->rtr[2];
 347         raxes[2].y = cell->rtr[5];
 348         raxes[2].z = cell->rtr[8];
 349         energy_rec = 0.0;
 350         coul_v = COULOMBIC / arena->dielectric / (2 * PI * TransformDeterminant(cell->tr));
 351         if (arena->debug_result != NULL)
 352                 fprintf(arena->debug_result, "Direct Ewald sum\n");
 353         forces = &arena->forces[kPMEIndex * arena->mol->natoms];
 354         for (n = 0, ip = arena->pme->kxyz; n < arena->pme->nkxyz; n++, ip += 3) {
 355                 Vector vm, f;
 356                 Double de, w, ww, m2;
 357                 Double s_real, s_imag, *mp;
 358                 nx = ip[0];
 359                 ny = ip[1];
 360                 nz = ip[2];
 361                 de = 0.0;
 362                 VecScale(vm, raxes[0], nx);
 363                 VecScaleInc(vm, raxes[1], ny);
 364                 VecScaleInc(vm, raxes[2], nz);
 365                 m2 = VecLength2(vm);
 366                 s_real = s_imag = 0;
 367                 mp = arena->pme->marray;
 368                 for (i = 0, ap = mol->atoms; i < mol->natoms; i++, ap = ATOM_NEXT(ap)) {
 369                         Double mr = VecDot(vm, ap->r) * 2 * PI;
 370                         mp[0] = ap->charge * cos(mr);
 371                         mp[1] = ap->charge * sin(mr);
 372                         s_real += mp[0];
 373                         s_imag += mp[1];
 374                         mp += 2;
 375                 }
 376                 w = PI * PI * m2 / (arena->ewald_beta * arena->ewald_beta);
 377                 ww = exp(-w) / m2;
 378                 de += ww * (s_real * s_real + s_imag * s_imag);
 379                 mp = arena->pme->marray;
 380                 for (i = 0; i < mol->natoms; i++) {
 381                         Double w2 = 4 * PI * ww * (s_imag * mp[0] - s_real * mp[1]);
 382                         f.x = w2 * (nx * raxes[0].x + ny * raxes[1].x + nz * raxes[2].x);
 383                         f.y = w2 * (nx * raxes[0].y + ny * raxes[1].y + nz * raxes[2].y);
 384                         f.z = w2 * (nx * raxes[0].z + ny * raxes[1].z + nz * raxes[2].z);
 385                         VecInc(forces[i], f);
 386                         mp += 2;
 387                 }
 388                 if (arena->debug_result != NULL && arena->debug_output_level > 0) {
 389                         fprintf(arena->debug_result, "(%d %d %d) %.9g %.9g %.6g\n", nx, ny, nz, (ww * (s_real * s_real + s_imag * s_imag)) *coul_v * INTERNAL2KCAL, de * coul_v * INTERNAL2KCAL, sqrt(w) / PI);
 390                 }
 391                 de *= coul_v;
 392                 energy_rec += de;
 393         }
 394         for (i = 0; i < mol->natoms; i++) {
 395                 VecScaleSelf(forces[i], coul_v);
 396         }
 397         arena->energies[kPMEIndex] += energy_rec;
 398         if (arena->debug_result != NULL) {
 399                 fprintf(arena->debug_result, "Reciprocal energy = %g\n", energy_rec);
 400                 if (arena->debug_output_level > 0) {
 401                         for (i = 0; i < mol->natoms; i++) {
 402                                 fprintf(arena->debug_result, "  reciprocal force %d = {%g,%g,%g}\n", i, forces[i].x * INTERNAL2KCAL, forces[i].y * INTERNAL2KCAL, forces[i].z * INTERNAL2KCAL);
 403                         }
 404                 }
 405         }
 406 }
 407
 408 /*  Particle Mesh Ewald Sum  */
 409 void
 410 calc_pme_force(MDArena *arena)
 411 {
 412         Molecule *mol;
 413         MDPME *pme;
 414         Atom *ap;
 415         Vector r;
 416         Int grid_x, grid_y, grid_z, order;
 417         Int i, j, k, kx, ky, kz, n, n0, n1, n2;
 418         Int debug_level = -1;
 419         double m0, m1, t, *mp, *bp, *bcp;
 420         fftw_complex *cp;
 421         double volume;
 422         double beta2;
 423         double dielec_r;
 424         Double *energy = &arena->energies[kPMEIndex];
 425         Vector *forces = &arena->forces[kPMEIndex * arena->mol->natoms];
 426
 427         if (arena == NULL || (mol = arena->mol) == NULL || (pme = arena->pme) == NULL)
 428                 return;
 429
 430         if (arena->debug_result != NULL) {
 431                 debug_level = arena->debug_output_level;
 432         }
 433
 434         forces = &arena->forces[kPMEIndex * arena->mol->natoms];
 435         dielec_r = COULOMBIC / arena->dielectric;
 436
 437         /*  Fill M and Q  */
 438         grid_x = pme->grid_x;
 439         grid_y = pme->grid_y;
 440         grid_z = pme->grid_z;
 441         order = pme->order;
 442         cp = pme->qarray;
 443         beta2 = arena->ewald_beta * arena->ewald_beta;
 444         memset(cp, 0, sizeof(fftw_complex) * grid_x * grid_y * grid_z);
 445         for (i = 0, ap = mol->atoms; i < mol->natoms; i++, ap = ATOM_NEXT(ap)) {
 446                 /*  Calculate fractional coordinate  */
 447                 TransformVec(&r, mol->cell->rtr, &ap->r);
 448                 r.x *= grid_x;
 449                 r.y *= grid_y;
 450                 r.z *= grid_z;
 451                 mp = pme->marray + (grid_x + grid_y + grid_z) * 2 * i;
 452                 for (k = 0; k < grid_x; k++) {
 453                         if (grid_x == 1) {
 454                                 m0 = 1.0;
 455                                 m1 = 1.0;
 456                         } else {
 457                                 n0 = floor((r.x - order - k) / grid_x) + 1;
 458                                 n1 = floor((r.x - k) / grid_x);
 459                                 m0 = m1 = 0.0;
 460                                 for (n = n0; n <= n1; n++) {
 461                                         t = r.x - k - n * grid_x;
 462                                         m0 += s_calc_spline(order, t);
 463                                         m1 += (s_calc_spline(order - 1, t) - s_calc_spline(order - 1, t - 1)) * grid_x;
 464                                 }
 465                         }
 466                         mp[k * 2] = m0;
 467                         mp[k * 2 + 1] = m1;
 468                         if (debug_level > 1) fprintf(arena->debug_result, "mx[%d] (%g, %g)\n", k, m0, m1);
 469                 }
 470                 for (k = 0; k < grid_y; k++) {
 471                         if (grid_y == 1) {
 472                                 m0 = 1.0;
 473                                 m1 = 1.0;
 474                         } else {
 475                                 n0 = floor((r.y - order - k) / grid_y) + 1;
 476                                 n1 = floor((r.y - k) / grid_y);
 477                                 m0 = m1 = 0.0;
 478                                 for (n = n0; n <= n1; n++) {
 479                                         t = r.y - k - n * grid_y;
 480                                         m0 += s_calc_spline(order, t);
 481                                         m1 += (s_calc_spline(order - 1, t) - s_calc_spline(order - 1, t - 1)) * grid_y;
 482                                 }
 483                         }
 484                         mp[(k + grid_x) * 2] = m0;
 485                         mp[(k + grid_x) * 2 + 1] = m1;
 486                         if (debug_level > 1) fprintf(arena->debug_result, "my[%d] (%g, %g)\n", k, m0, m1);
 487                 }
 488                 for (k = 0; k < grid_z; k++) {
 489                         if (grid_z == 1) {
 490                                 m0 = 1.0;
 491                                 m1 = 1.0;
 492                         } else {
 493                                 n0 = floor((r.z - order - k) / grid_z) + 1;
 494                                 n1 = floor((r.z - k) / grid_z);
 495                                 m0 = m1 = 0.0;
 496                                 for (n = n0; n <= n1; n++) {
 497                                         t = r.z - k - n * grid_z;
 498                                         m0 += s_calc_spline(order, t);
 499                                         m1 += (s_calc_spline(order - 1, t) - s_calc_spline(order - 1, t - 1)) * grid_z;
 500                                 }
 501                         }
 502                         mp[(k + grid_x + grid_y) * 2] = m0;
 503                         mp[(k + grid_x + grid_y) * 2 + 1] = m1;
 504                         if (debug_level > 1) fprintf(arena->debug_result, "mz[%d] (%g, %g)\n", k, m0, m1);
 505                 }
 506                 for (kx = 0; kx < grid_x; kx++) {
 507                         for (ky = 0; ky < grid_y; ky++) {
 508                                 for (kz = 0; kz < grid_z; kz++) {
 509                                         cp[kz + grid_z * (ky + grid_y * kx)][0] += ap->charge * mp[kx * 2] * mp[(ky + grid_x) * 2] * mp[(kz + grid_x + grid_y) * 2];
 510                                 }
 511                         }
 512                 }
 513         }
 514         if (debug_level > 1) {
 515                 for (kx = 0; kx < grid_x; kx++) {
 516                         for (ky = 0; ky < grid_y; ky++) {
 517                                 for (kz = 0; kz < grid_z; kz++) {
 518                                         j = kz + grid_z * (ky + grid_y * kx);
 519                                         fprintf(arena->debug_result, "Q[%d,%d,%d] = (%g, %g)\n", kx, ky, kz, pme->qarray[j][0], pme->qarray[j][1]);
 520                                 }
 521                         }
 522                 }
 523         }
 524
 525         /*  Q is transformed by inverse FFT into qfarray */
 526         fftw_execute(pme->plan1);
 527
 528         /*  Dump Finv(Q) and structure factor */
 529         if (debug_level > 1) {
 530                 fprintf(arena->debug_result, "kx ky kz FB_re FB_im S_re S_im Sa_re Sa_im\n");
 531                 for (kx = 0; kx < grid_x; kx++) {
 532                         for (ky = 0; ky < grid_y; ky++) {
 533                                 for (kz = 0; kz < grid_z; kz++) {
 534                                         fftw_complex bf;
 535                                         j = kz + grid_z * (ky + grid_y * kx);
 536                                         fprintf(arena->debug_result, "%2d %2d %2d ", kx, ky, kz);
 537                                         bf[0] = pme->qfarray[j][0];
 538                                         bf[1] = pme->qfarray[j][1];
 539                                         s_complex_mul(bf, bf, pme->b0array[kx]);
 540                                         s_complex_mul(bf, bf, pme->b0array[ky + grid_x]);
 541                                         s_complex_mul(bf, bf, pme->b0array[kz + grid_x + grid_y]);
 542                                         fprintf(arena->debug_result, "%g %g ", bf[0], bf[1]);
 543                                         m0 = m1 = 0.0;
 544                                         bf[0] = bf[1] = 0;
 545                                         for (i = 0, ap = mol->atoms; i < mol->natoms; i++, ap = ATOM_NEXT(ap)) {
 546                                                 /*  Calc S(m)  */
 547                                                 fftw_complex cw1, cw2;
 548                                                 Int kx0, ky0, kz0;
 549                                                 TransformVec(&r, mol->cell->rtr, &ap->r);
 550                                                 kx0 = (kx <= grid_x / 2 ? kx : kx - grid_x);
 551                                                 ky0 = (ky <= grid_y / 2 ? ky : ky - grid_y);
 552                                                 kz0 = (kz <= grid_z / 2 ? kz : kz - grid_z);
 553                                                 t = 2 * PI * (kx0 * r.x + ky0 * r.y + kz0 * r.z);
 554                                                 m0 += ap->charge * cos(t);
 555                                                 m1 += ap->charge * sin(t);
 556                                                 /*  Calc interpolated S(m)  */
 557                                                 n0 = floor(r.x * grid_x - order) + 1;
 558                                                 n1 = floor(r.x * grid_x);
 559                                                 cw1[0] = cw1[1] = 0.0;
 560                                                 for (j = n0; j <= n1; j++) {
 561                                                         double w = 2 * PI * kx * j / grid_x;
 562                                                         t = s_calc_spline(order, r.x * grid_x - j);
 563                                                         cw1[0] += t * cos(w);
 564                                                         cw1[1] += t * sin(w);
 565                                                 }
 566                                                 s_complex_mul(cw2, cw1, pme->b0array[kx]);
 567                                                 n0 = floor(r.y * grid_y - order) + 1;
 568                                                 n1 = floor(r.y * grid_y);
 569                                                 cw1[0] = cw1[1] = 0.0;
 570                                                 for (j = n0; j <= n1; j++) {
 571                                                         double w = 2 * PI * ky * j / grid_y;
 572                                                         t = s_calc_spline(order, r.y * grid_y - j);
 573                                                         cw1[0] += t * cos(w);
 574                                                         cw1[1] += t * sin(w);
 575                                                 }
 576                                                 s_complex_mul(cw1, cw1, pme->b0array[ky + grid_x]);
 577                                                 s_complex_mul(cw2, cw2, cw1);
 578                                                 n0 = floor(r.z * grid_z - order) + 1;
 579                                                 n1 = floor(r.z * grid_z);
 580                                                 cw1[0] = cw1[1] = 0.0;
 581                                                 for (j = n0; j <= n1; j++) {
 582                                                         double w = 2 * PI * kz * j / grid_z;
 583                                                         t = s_calc_spline(order, r.z * grid_z - j);
 584                                                         cw1[0] += t * cos(w);
 585                                                         cw1[1] += t * sin(w);
 586                                                 }
 587                                                 s_complex_mul(cw1, cw1, pme->b0array[kz + grid_x + grid_y]);
 588                                                 s_complex_mul(cw2, cw2, cw1);
 589                                                 bf[0] += ap->charge * cw2[0];
 590                                                 bf[1] += ap->charge * cw2[1];
 591                                         }
 592                                         fprintf(arena->debug_result, "%g %g %g %g\n", m0, m1, bf[0], bf[1]);
 593                                 }
 594                         }
 595                 }
 596         }
 597
 598         /*  C is built and B*C is multiplied into qfarray  */
 599         volume = TransformDeterminant(mol->cell->tr);
 600         bp = pme->barray;
 601         bcp = pme->bcarray;
 602         *energy = 0.0;
 603         for (kx = 0; kx < grid_x; kx++) {
 604                 n0 = (kx <= grid_x / 2 ? kx : kx - grid_x);
 605                 for (ky = 0; ky < grid_y; ky++) {
 606                         n1 = (ky <= grid_y / 2 ? ky : ky - grid_y);
 607                         for (kz = 0; kz < grid_z; kz++) {
 608                                 Transform *rtr = &(mol->cell->rtr);
 609                                 Vector mv;
 610                                 double d;
 611                                 n2 = (kz <= grid_z / 2 ? kz : kz - grid_z);
 612                                 /*  rtr[0,3,6], rtr[1,4,7], rtr[2,5,8] are the reciprocal unit cell vectors  */
 613                                 if (n0 == 0 && n1 == 0 && n2 == 0)
 614                                         d = 0.0;
 615                                 else {
 616                                         mv.x = (*rtr)[0] * n0 + (*rtr)[1] * n1 + (*rtr)[2] * n2;
 617                                         mv.y = (*rtr)[3] * n0 + (*rtr)[4] * n1 + (*rtr)[5] * n2;
 618                                         mv.z = (*rtr)[6] * n0 + (*rtr)[7] * n1 + (*rtr)[8] * n2;
 619                                         d = VecLength2(mv);
 620                                         d = exp(-PI * PI * d / beta2) / (PI * volume * d);
 621                                 }
 622                                 k = kz + grid_z * (ky + grid_y * kx);
 623                                 pme->carray[k] = d;
 624                                 bcp[k] = bp[kx] * bp[grid_x + ky] * bp[grid_x + grid_y + kz] * d;
 625                                 m0 = pme->qfarray[k][0];
 626                                 m1 = pme->qfarray[k][1];
 627                                 pme->qarray[k][0] = m0; /* for debug */
 628                                 pme->qarray[k][1] = m1; /* for debug */
 629                                 *energy += (m0 * m0 + m1 * m1) * bcp[k];
 630                                 pme->qfarray[k][0] = m0 * bcp[k];
 631                                 pme->qfarray[k][1] = m1 * bcp[k];
 632                         }
 633                 }
 634         }
 635         *energy *= dielec_r * 0.5;
 636
 637         /*  qfarray is transformed by FFT into qqarray  */
 638         fftw_execute(pme->plan2);
 639
 640         if (debug_level > 1) {
 641                 fprintf(arena->debug_result, "kx ky kz FB_re FB_im F2B2C b c bc QF_re QF_im\n");
 642                 for (kx = 0; kx < grid_x; kx++) {
 643                         for (ky = 0; ky < grid_y; ky++) {
 644                                 for (kz = 0; kz < grid_z; kz++) {
 645                                         fftw_complex bf;
 646                                         j = kz + grid_z * (ky + grid_y * kx);
 647                                         fprintf(arena->debug_result, "%2d %2d %2d ", kx, ky, kz);
 648                                         bf[0] = pme->qarray[j][0];
 649                                         bf[1] = pme->qarray[j][1];
 650                                         s_complex_mul(bf, bf, pme->b0array[kx]);
 651                                         s_complex_mul(bf, bf, pme->b0array[ky + grid_x]);
 652                                         s_complex_mul(bf, bf, pme->b0array[kz + grid_x + grid_y]);
 653                                         fprintf(arena->debug_result, "%g %g ", bf[0], bf[1]);
 654                                         fprintf(arena->debug_result, "%g ", (bf[0] * bf[0] + bf[1] * bf[1]) * pme->carray[j]);
 655                                 //      fprintf(arena->debug_result, "%g %g ", pme->qarray[j][0] * pme->bcarray[j], pme->qarray[j][1] * pme->bcarray[j]);
 656                                         fprintf(arena->debug_result, "%g %g %g ", pme->barray[kx] * pme->barray[ky + grid_x] * pme->barray[kz + grid_x + grid_y], pme->carray[j], pme->bcarray[j]);
 657                                         fprintf(arena->debug_result, "%g %g\n", pme->qqarray[j][0], pme->qqarray[j][1]);
 658                                 }
 659                         }
 660                 }
 661         }
 662
 663         /*  Calculate force  */
 664         memset(forces, 0, sizeof(Vector) * mol->natoms);
 665         for (kx = 0; kx < grid_x; kx++) {
 666                 for (ky = 0; ky < grid_y; ky++) {
 667                         for (kz = 0; kz < grid_z; kz++) {
 668                                 k = kz + grid_z * (ky + grid_y * kx);
 669                                 mp = pme->marray;
 670                                 for (i = 0; i < mol->natoms; i++) {
 671                                         t = ATOM_AT_INDEX(mol->atoms, i)->charge * pme->qqarray[k][0];
 672                                         forces[i].x += t * mp[kx * 2 + 1] * mp[(ky + grid_x) * 2] * mp[(kz + grid_x + grid_y) * 2];
 673                                         forces[i].y += t * mp[kx * 2] * mp[(ky + grid_x) * 2 + 1] * mp[(kz + grid_x + grid_y) * 2];
 674                                         forces[i].z += t * mp[kx * 2] * mp[(ky + grid_x) * 2] * mp[(kz + grid_x + grid_y) * 2 + 1];
 675                                         mp += (grid_x + grid_y + grid_z) * 2;
 676                                 }
 677                         }
 678                 }
 679         }
 680         {
 681                 /*  Scale the force vectors and transform into cartesian  */
 682                 Transform tf;
 683                 memmove(&tf, &mol->cell->rtr, sizeof(Transform));
 684                 tf[9] = tf[10] = tf[11] = 0.0;
 685                 for (i = 0; i < mol->natoms; i++) {
 686                         Vector v;
 687                         VecScale(v, forces[i], dielec_r);
 688                         TransformVec(&forces[i], tf, &v);
 689                 }
 690         }
 691
 692         if (debug_level > 0) {
 693                 fprintf(arena->debug_result, "PME reciprocal energy = %.9g\n", *energy * INTERNAL2KCAL);
 694                 fprintf(arena->debug_result, "PME reciprocal forces:\n");
 695                 for (i = 0; i < mol->natoms; i++)
 696                         fprintf(arena->debug_result, "%d %.9g %.9g %.9g\n", i, forces[i].x * INTERNAL2KCAL, forces[i].y * INTERNAL2KCAL, forces[i].z * INTERNAL2KCAL);
 697                 if (debug_level > 1) {
 698                         fprintf(arena->debug_result, "Grid values:\n");
 699                         for (i = 0; i < (grid_x + grid_y + grid_z); i++) {
 700                                 fprintf(arena->debug_result, "%d %g %g\n", i, pme->marray[i * 2], pme->marray[i * 2 + 1]);
 701                         }
 702                 }
 703         }
 704 }
 705
 706 void
 707 calc_ewald_force(MDArena *arena)
 708 {
 709         if (arena->use_ewald != 0) {
 710                 if (arena->use_ewald == 1)
 711                         calc_pme_force(arena);
 712                 else if (arena->use_ewald == 2)
 713                         calc_direct_ewald_force(arena);
 714                 else if (arena->use_ewald == -1)   /*  Direct sum of electrostatic (very slow)  */
 715                         calc_direct_electrostatic(arena);
 716                 else return;
 717                 arena->energies[kPMEIndex] += arena->pme->self_correction;
 718                 if (arena->debug_result != NULL) {
 719                         fprintf(arena->debug_result, "%s Ewald reciprocal energy = %.9g\n", (arena->use_ewald == 1 ? "Particle Mesh" : "Direct"), (arena->energies[kPMEIndex] - arena->pme->self_correction) * INTERNAL2KCAL);
 720                         fprintf(arena->debug_result, "Self correction energy = %.9g\n", arena->pme->self_correction * INTERNAL2KCAL);
 721                         fprintf(arena->debug_result, "Net Ewald indirect energy = %.9g\n", arena->energies[kPMEIndex] * INTERNAL2KCAL);
 722                 }
 723         }
 724 }
 725
 726 /*  For debug  */
 727 int
 728 pme_test(MDArena *arena)
 729 {
 730         s_initialize_spline_coeffs();
 731         calc_pme_force(arena);
 732         return 0;
 733 }
 734
 735 void
 736 pme_release(MDArena *arena)
 737 {
 738         if (arena == NULL || arena->pme == NULL)
 739                 return;
 740         if (arena->pme->plan1 != 0)
 741                 fftw_destroy_plan(arena->pme->plan1);
 742         if (arena->pme->plan2 != 0)
 743                 fftw_destroy_plan(arena->pme->plan2);
 744         if (arena->pme->qarray != NULL)
 745                 fftw_free(arena->pme->qarray);
 746         if (arena->pme->qqarray != NULL)
 747                 fftw_free(arena->pme->qqarray);
 748         if (arena->pme->qfarray != NULL)
 749                 fftw_free(arena->pme->qfarray);
 750         if (arena->pme->marray != NULL)
 751                 free(arena->pme->marray);
 752         if (arena->pme->bcarray != NULL)
 753                 free(arena->pme->bcarray);
 754         if (arena->pme->barray != NULL)
 755                 free(arena->pme->barray);
 756         if (arena->pme->kxyz)
 757                 free(arena->pme->kxyz);
 758         free(arena->pme);
 759         arena->pme = NULL;
 760 }
 761
 762 static void
 763 s_pme_spline_test(MDArena *arena)
 764 {
 765         /*  Compare:
 766             exp(2*PI*i*m*x/K)                            ; target (S(m))
 767                 b(m)*Sum(k, spline(x - k)*exp(2*PI*i*m*k/K)) ; interpolation
 768             b(m)*Fourier(k, m, Sum(n, spline(x - k - n*K))) ; Fourier
 769       b(m) = exp(2*PI*i*(order-1)*m/K)/Sum(k=0..order-2, spline(k+1)*exp(2*PI*i*m*k/K))
 770           Fourier(k, m, f(k)) = Sum(k, f(k)*exp(2*PI*i*m*k/K))  */
 771 #define DIM 32
 772 #define ORDER 8
 773         Int kk = DIM;
 774         fftw_complex a1[DIM], a2[DIM], a3[DIM], a4[DIM], b[DIM], c1, c2;
 775         double b2[DIM], c[DIM];
 776         double d, w, x;
 777         double len = 10.0;  /*  Unit cell dimension  */
 778         Int k, m, n0, n1;
 779         x = 3.7;
 780         for (m = 0; m < DIM; m++) {
 781                 k = (m <= DIM / 2 ? m : m - DIM);
 782                 d = 2 * PI * k * x / kk;
 783                 a1[m][0] = cos(d);
 784                 a1[m][1] = sin(d);
 785         }
 786         for (m = 0; m < DIM; m++) {
 787                 c1[0] = c1[1] = 0.0;
 788                 for (k = 0; k < ORDER - 2; k++) {
 789                         d = 2 * PI * m * k / kk;
 790                         w = s_calc_spline(ORDER, k + 1.0);
 791                         c1[0] += w * cos(d);
 792                         c1[1] += w * sin(d);
 793                 }
 794                 c1[1] *= -1;
 795                 d = c1[0] * c1[0] + c1[1] * c1[1];
 796                 c1[0] /= d;
 797                 c1[1] /= d;
 798                 d = 2 * PI * (ORDER - 1) * m / kk;
 799                 c2[0] = cos(d);
 800                 c2[1] = sin(d);
 801                 s_complex_mul(b[m], c2, c1);
 802         }
 803         for (m = 0; m < DIM; m++) {
 804                 n0 = floor(x - ORDER) + 1;
 805                 n1 = floor(x);
 806                 c1[0] = c1[1] = 0.0;
 807                 for (k = n0; k <= n1; k++) {
 808                         d = 2 * PI * m * k / kk;
 809                         w = s_calc_spline(ORDER, x - k);
 810                         c1[0] += w * cos(d);
 811                         c1[1] += w * sin(d);
 812                 }
 813                 s_complex_mul(a2[m], b[m], c1);
 814         }
 815 //      b(m)*Fourier(k, m, Sum(n, spline(x - k - n*K))) ; Fourier
 816         for (m = 0; m < DIM; m++) {
 817                 n0 = floor((x - m - ORDER) / kk) + 1;
 818                 n1 = floor((x - m) / kk);
 819                 w = 0.0;
 820                 for (k = n0; k <= n1; k++) {
 821                         w += s_calc_spline(ORDER, x - m - k * kk);
 822                 }
 823                 a3[m][0] = w;
 824                 a3[m][1] = 0.0;
 825         }
 826         //  Fourier(k, m, f(k)) = Sum(k, f(k)*exp(2*PI*i*m*k/K))
 827         for (m = 0; m < DIM; m++) {
 828                 c1[0] = c1[1] = 0.0;
 829                 for (k = 0; k < DIM; k++) {
 830                         d = -2 * PI * m * k / kk;  /*  Inverse Fourier  */
 831                         c1[0] += a3[k][0] * cos(d);
 832                         c1[1] += a3[k][0] * sin(d);
 833                 }
 834                 b[m][1] *= -1;
 835                 s_complex_mul(a4[m], b[m], c1);  /*  Finv(f(k))*b_conjugate should be S_conjugate  */
 836                 b[m][1] *= -1;
 837         }
 838         //  b2 = |b|^2
 839         //  c = exp(-PI^2*|m|^2/(beta^2))/|m|^2)
 840         for (m = 0; m < DIM; m++) {
 841                 b2[m] = b[m][0] * b[m][0] + b[m][1] * b[m][1];
 842                 w = (m <= DIM / 2 ? m : m - DIM) / len;
 843                 c[m] = exp(-PI * PI * w * w / (0.5 * 0.5)) / (w * w);
 844         }
 845         fprintf(stderr, "a1_re a1_im a2_re a2_im a3 a4_re a4_im b2c\n");
 846         for (m = 0; m < DIM; m++) {
 847                 fprintf(stderr, "%.6g %.6g %.6g %.6g %.6g %.6g %.6g %.6g\n", a1[m][0], a1[m][1], a2[m][0], a2[m][1], a3[m][0], a4[m][0], a4[m][1], b2[m] * c[m]);
 848         }
 849         d = 0.0;
 850         w = 0.0;
 851         for (m = 1; m < DIM; m++) {
 852                 d += b2[m] * c[m] * (a1[m][0] * a1[m][0] + a1[m][1] * a1[m][1]);
 853                 w += b2[m] * c[m] * (a4[m][0] * a4[m][0] + a4[m][1] * a4[m][1]);
 854         }
 855         fprintf(stderr, "\nEnergy = %g %g\n", d, w);
 856 }
 857
 858 void
 859 pme_init(MDArena *arena)
 860 {
 861         MDPME *pme;
 862         Int xyz, x_y_z, i, j, k;
 863         if (arena == NULL || arena->mol == NULL || arena->mol->cell == NULL) {
 864                 if (arena->pme != NULL)
 865                         pme_release(arena);
 866                 return;
 867         }
 868         if (arena->use_ewald == 2) {
 869                 /*  Direct Ewald  */
 870                 if (arena->pme != NULL)
 871                         pme_release(arena);
 872                 arena->pme = (MDPME *)calloc(sizeof(MDPME), 1);
 873                 if (arena->pme == NULL)
 874                         return;
 875                 /*  Only marray and kxyz is allocated (for force calculation)  */
 876                 arena->pme->marray = (double *)malloc(sizeof(double) * 2 * arena->mol->natoms);
 877                 arena->pme->direct_limit = 1.29;
 878                 /*  kxyz is initialized  */
 879                 s_prepare_direct_ewald(arena);
 880         } else {
 881                 s_initialize_spline_coeffs();
 882 #if 0
 883                 {
 884                         /*  spline test  */
 885                         int n;
 886                         double x;
 887                         for (n = 3; n <= MAX_DIM_SPLINE; n++) {
 888                                 for (x = 0.0; x <= n; x += 0.1) {
 889                                         double y = s_calc_spline(n, x);
 890                                         double yy = (x * s_calc_spline(n - 1, x) + (n - x) * s_calc_spline(n - 1, x - 1.0)) / (n - 1);
 891                                         if (fabs(y - yy) > 1e-8) {
 892                                                 fprintf(stderr, "spline test failed n = %d x = %g y = %g yy = %g\n", n, x, y, yy);
 893                                                 break;
 894                                         }
 895                                 }
 896                         }
 897                 }
 898                 s_pme_spline_test(arena);
 899 #endif
 900
 901                 xyz = arena->ewald_grid_x * arena->ewald_grid_y * arena->ewald_grid_z;
 902                 x_y_z = arena->ewald_grid_x + arena->ewald_grid_y + arena->ewald_grid_z;
 903                 if (xyz == 0) {
 904                         if (arena->pme != NULL)
 905                                 pme_release(arena);
 906                         return;
 907                 }
 908                 if (arena->pme == NULL) {
 909                         arena->pme = (MDPME *)calloc(sizeof(MDPME), 1);
 910                         if (arena->pme == NULL)
 911                                 return;
 912                 }
 913                 pme = arena->pme;
 914                 if (pme->grid_x != arena->ewald_grid_x || pme->grid_y != arena->ewald_grid_y || pme->grid_z != arena->ewald_grid_z || pme->natoms != arena->mol->natoms) {
 915                         /*  Need to rebuild internal table  */
 916                         if (pme->grid_x * pme->grid_y * pme->grid_z < xyz) {
 917                                 if (pme->qarray != NULL)
 918                                         fftw_free(pme->qarray);
 919                                 pme->qarray = fftw_alloc_complex(xyz);
 920                                 if (pme->qfarray != NULL)
 921                                         fftw_free(pme->qfarray);
 922                                 pme->qfarray = fftw_alloc_complex(xyz);
 923                                 if (pme->qqarray != NULL)
 924                                         fftw_free(pme->qqarray);
 925                                 pme->qqarray = fftw_alloc_complex(xyz);
 926                                 pme->bcarray = (double *)realloc(pme->bcarray, sizeof(double) * xyz);
 927                                 pme->carray = (double *)realloc(pme->carray, sizeof(double) * xyz);
 928                         }
 929                         if (pme->b0array != NULL)
 930                                 fftw_free(pme->b0array);
 931                         pme->b0array = fftw_alloc_complex(x_y_z);
 932                         pme->marray = (double *)realloc(pme->marray, sizeof(double) * x_y_z * 2 * arena->mol->natoms);
 933                         pme->barray = (double *)realloc(pme->barray, sizeof(double) * x_y_z);
 934                         pme->grid_x = arena->ewald_grid_x;
 935                         pme->grid_y = arena->ewald_grid_y;
 936                         pme->grid_z = arena->ewald_grid_z;
 937                         pme->natoms = arena->mol->natoms;
 938                         if (xyz != 0 && x_y_z != 0) {
 939                                 if (pme->plan1 != NULL)
 940                                         fftw_destroy_plan(pme->plan1);
 941                                 if (pme->plan2 != NULL)
 942                                         fftw_destroy_plan(pme->plan2);
 943                                 pme->plan1 = fftw_plan_dft_3d(pme->grid_x, pme->grid_y, pme->grid_z, pme->qarray, pme->qfarray, FFTW_BACKWARD, FFTW_ESTIMATE);
 944                                 pme->plan2 = fftw_plan_dft_3d(pme->grid_x, pme->grid_y, pme->grid_z, pme->qfarray, pme->qqarray, FFTW_FORWARD, FFTW_ESTIMATE);
 945                         }
 946                         pme->order = arena->ewald_order;
 947                         /*  Should be even numbers in [4..MAX_DIM_SPLINE]  */
 948                         if (pme->order < 4)
 949                                 pme->order = 4;
 950                         if (pme->order > MAX_DIM_SPLINE)
 951                                 pme->order = MAX_DIM_SPLINE;
 952                         pme->order -= pme->order % 2;
 953                         /*  Initialize B array  */
 954                         for (i = 0; i < 3; i++) {
 955                                 int idx_ofs, grid_i;
 956                                 switch (i) {
 957                                         case 0: grid_i = pme->grid_x; idx_ofs = 0; break;
 958                                         case 1: grid_i = pme->grid_y; idx_ofs = pme->grid_x; break;
 959                                         case 2: grid_i = pme->grid_z; idx_ofs = pme->grid_x + pme->grid_y; break;
 960                                 }
 961                                 for (j = 0; j < grid_i; j++) {
 962                                         fftw_complex b1, b2;
 963                                         double t;
 964                                         if (grid_i == 1) {
 965                                                 b1[0] = 1.0;
 966                                                 b1[1] = 0.0;
 967                                         } else {
 968                                                 t = 2 * PI * (pme->order - 1) * j / grid_i;
 969                                                 b1[0] = cos(t);
 970                                                 b1[1] = sin(t);
 971                                                 b2[0] = b2[1] = 0.0;
 972                                                 for (k = 0; k < pme->order - 1; k++) {
 973                                                         double s;
 974                                                         t = 2 * PI * j * k / grid_i;
 975                                                         s = s_calc_spline(pme->order, k + 1);
 976                                                         b2[0] += cos(t) * s;
 977                                                         b2[1] += sin(t) * s;
 978                                                 }
 979                                                 /*  Invert b2  */
 980                                                 t = 1.0 / (b2[0] * b2[0] + b2[1] * b2[1]);
 981                                                 b2[0] *= t;
 982                                                 b2[1] *= -t;
 983                                                 s_complex_mul(b1, b1, b2);
 984                                         }
 985                                         pme->b0array[idx_ofs + j][0] = b1[0];
 986                                         pme->b0array[idx_ofs + j][1] = b1[1];
 987                                         pme->barray[idx_ofs + j] = b1[0] * b1[0] + b1[1] * b1[1];
 988                                 }
 989                         }
 990                 }
 991         }
 992         /*  Calculate self-correction term  */
 993         {
 994                 Double en;
 995                 Atom *ap;
 996                 en = 0.0;
 997                 for (i = 0, ap = arena->mol->atoms; i < arena->mol->natoms; i++, ap = ATOM_NEXT(ap)) {
 998                         en += ap->charge * ap->charge;
 999                 }
1000                 arena->pme->self_correction = -(COULOMBIC/arena->dielectric) * arena->ewald_beta * PI2R * en;
1001         }
1002 }