An error in PME energy calculation is corrected. MDArean#ewald_order is implemented.

git-svn-id: svn+ssh://svn.sourceforge.jp/svnroot/molby/trunk@339 a2be9bc6-48de-4e38-9406-05402d4bc13c

diff --git a/MolLib/MD/MDCore.c b/MolLib/MD/MDCore.c
index 10e656e..d06921b 100644 (file)
--- a/MolLib/MD/MDCore.c
+++ b/MolLib/MD/MDCore.c
@@ -3505,6 +3505,7 @@ md_set_default(MDArena *arena)
        arena->ewald_grid_y = 16;
        arena->ewald_grid_z = 16;
        arena->ewald_freq = 2;
+       arena->ewald_order = 8;
        
        /*      arena->pressure_coupling = 0.4;
        arena->pressure_trial_width = 0.01;

diff --git a/MolLib/MD/MDCore.h b/MolLib/MD/MDCore.h
index 12108f1..eae78dc 100644 (file)
--- a/MolLib/MD/MDCore.h
+++ b/MolLib/MD/MDCore.h
@@ -279,8 +279,9 @@ typedef struct MDArena {
        Int    use_ewald;       /*  1 if Ewald technique is to be used  */
        Double ewald_beta;      /*  beta for Ewald sum  */
        Int    ewald_grid_x, ewald_grid_y, ewald_grid_z;  /*  Number of grids for each direction; if all zero, then direct Ewald is used  */
-       Int    ewald_freq;
-       
+       Int    ewald_freq;      /*  How frequent should Ewald sum be calculated? (default: 2)  */
+       Int    ewald_order;     /*  Order of cardinal B-spline for approximating structure factor (default: 8)  */
+
        /*  Runtime fields  */
        
        /*  Initialize flag  */

diff --git a/MolLib/MD/MDEwald.c b/MolLib/MD/MDEwald.c
index b2ffef8..fe90378 100755 (executable)
--- a/MolLib/MD/MDEwald.c
+++ b/MolLib/MD/MDEwald.c
@@ -197,17 +197,23 @@ s_prepare_direct_ewald(MDArena *arena)
                if (pass == 1)
                        ip = arena->pme->kxyz;
                for (n = ninit; n <= 50; n += nstep) {
-                       Int nx, ny, nz;
+                       Int nx, ny, nz, nx0, ny0, nz0;
                        Vector vx, vy, vm;
                        Double len;
                        dc = 0;
-                       for (nx = -n + 1; nx < n; nx++) {
+                       nx0 = (arena->periodic_a ? n : 1);
+                       ny0 = (arena->periodic_b ? n : 1);
+                       nz0 = (arena->periodic_c ? n : 1);
+                       for (nx = -nx0 + 1; nx < nx0; nx++) {
                                VecScale(vx, raxes[0], nx);
-                               for (ny = -n + 1; ny < n; ny++) {
+                               for (ny = -ny0 + 1; ny < ny0; ny++) {
                                        VecScale(vy, raxes[1], ny);
-                                       for (nz = -n + 1; nz < n; nz++) {
-                                               if (n > ninit && nx > -n + nstep && nx < n - nstep && ny > -n + nstep && ny < n - nstep && nz > -n + nstep && nz < n - nstep)
-                                                       continue;
+                                       for (nz = -nz0 + 1; nz < nz0; nz++) {
+                                               if (n > ninit)
+                                                       if (arena->periodic_a == 0 || (nx > -n + nstep && nx < n - nstep))
+                                                               if (arena->periodic_b == 0 || (ny > -n + nstep && ny < n - nstep))
+                                                                       if (arena->periodic_c == 0 || (nz > -n + nstep && nz < n - nstep))
+                                                                               continue;
                                                if (nx == 0 && ny == 0 && nz == 0)
                                                        continue;
                                                VecScale(vm, raxes[2], nz);
@@ -247,29 +253,33 @@ calc_direct_electrostatic(MDArena *arena)
        Vector *forces, *f;
        Double energy;
        Double coul;
-       Int i, j, n, nx, ny, nz, zero;
+       Int i, j, n, nx, ny, nz, nx0, ny0, nz0, zero, ninit, nstep;
        Atom *ap, *apj;
        
-       forces = (Vector *)calloc(sizeof(Vector), mol->natoms);
+       forces = &arena->forces[kPMEIndex * arena->mol->natoms];
        energy = 0.0;
        coul = COULOMBIC / arena->dielectric;
-       for (n = 5; n < 50; n += 5) {
+       ninit = 3;
+       nstep = 1;
+       for (n = ninit; n <= 50; n += nstep) {
                Vector vx, vy, vz;
                Double de;
                de = 0.0;
-               for (nx = -n + 1; nx < n; nx++) {
-                       if (nx > -n + 5 && nx < n - 5)
-                               continue;
+               nx0 = (arena->periodic_a ? n : 1);
+               ny0 = (arena->periodic_b ? n : 1);
+               nz0 = (arena->periodic_c ? n : 1);
+               for (nx = -nx0 + 1; nx < nx0; nx++) {
                        VecScale(vx, cell->axes[0], nx);
-                       for (ny = -n + 1; ny < n; ny++) {
-                               if (ny > -n + 5 && ny < n - 5)
-                                       continue;
+                       for (ny = -ny0 + 1; ny < ny0; ny++) {
                                VecScale(vy, cell->axes[1], ny);
-                               for (nz = -n + 1; nz < n; nz++) {
-                                       if (nz > -n + 5 && nz < n - 5)
-                                               continue;
-                                       VecScale(vz, cell->axes[2], nz);
+                               for (nz = -nz0 + 1; nz < nz0; nz++) {
+                                       if (n > ninit)
+                                               if (arena->periodic_a == 0 || (nx > -n + nstep && nx < n - nstep))
+                                                       if (arena->periodic_b == 0 || (ny > -n + nstep && ny < n - nstep))
+                                                               if (arena->periodic_c == 0 || (nz > -n + nstep && nz < n - nstep))
+                                                                       continue;
                                        zero = (nx == 0 && ny == 0 && nz == 0);
+                                       VecScale(vz, cell->axes[2], nz);
                                        for (i = 0, ap = mol->atoms, f = forces; i < mol->natoms; i++, ap = ATOM_NEXT(ap), f++) {
                                                Vector ri, rj, rij;
                                                ri = ap->r;
@@ -297,16 +307,22 @@ calc_direct_electrostatic(MDArena *arena)
                }
                de *= 0.5 * coul;
                energy += de;
-               if (arena->debug_result != NULL) {
+               if (arena->debug_result != NULL && arena->debug_output_level > 0) {
                        fprintf(arena->debug_result, "Direct electrostatic energy\n");
                        fprintf(arena->debug_result, "N = %d energy = %.9g %.9g\n", n, energy * INTERNAL2KCAL, de * INTERNAL2KCAL);
-                       for (i = 0; i < mol->natoms; i++) {
-                               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);
+                       if (arena->debug_output_level > 1) {
+                               for (i = 0; i < mol->natoms; i++) {
+                                       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);
+                               }
                        }
                }
                if (fabs(de) < 1e-7 || fabs(de / energy) < 1e-6)
                        break;
        }
+       if (arena->debug_result != NULL) {
+               fprintf(arena->debug_result, "Direct electrostatic energy = %g\n", energy);
+       }
+       arena->energies[kPMEIndex] += energy;
 }
 
 /*  Direct Ewald Sum  */
@@ -379,9 +395,12 @@ calc_direct_ewald_force(MDArena *arena)
                VecScaleSelf(forces[i], coul_v);
        }
        arena->energies[kPMEIndex] += energy_rec;
-       if (arena->debug_result != NULL && arena->debug_output_level > 0) {
-               for (i = 0; i < mol->natoms; i++) {
-                       fprintf(arena->debug_result, "  reciprocal force %d = {%g,%g,%g}\n", i, forces[i].x * INTERNAL2KCAL, forces[i].y * INTERNAL2KCAL, forces[i].z * INTERNAL2KCAL);
+       if (arena->debug_result != NULL) {
+               fprintf(arena->debug_result, "Reciprocal energy = %g\n", energy_rec);
+               if (arena->debug_output_level > 0) {
+                       for (i = 0; i < mol->natoms; i++) {
+                               fprintf(arena->debug_result, "  reciprocal force %d = {%g,%g,%g}\n", i, forces[i].x * INTERNAL2KCAL, forces[i].y * INTERNAL2KCAL, forces[i].z * INTERNAL2KCAL);
+                       }
                }
        }
 }
@@ -431,39 +450,54 @@ calc_pme_force(MDArena *arena)
                r.z *= grid_z;
                mp = pme->marray + (grid_x + grid_y + grid_z) * 2 * i;
                for (k = 0; k < grid_x; k++) {
-                       n0 = floor((r.x - order - k) / grid_x) + 1;
-                       n1 = floor((r.x - k) / grid_x);
-                       m0 = m1 = 0.0;
-                       for (n = n0; n <= n1; n++) {
-                               t = r.x - k - n * grid_x;
-                               m0 += s_calc_spline(order, t);
-                               m1 += (s_calc_spline(order - 1, t) - s_calc_spline(order - 1, t - 1)) * grid_x;
+                       if (grid_x == 1) {
+                               m0 = 1.0;
+                               m1 = 1.0;
+                       } else {
+                               n0 = floor((r.x - order - k) / grid_x) + 1;
+                               n1 = floor((r.x - k) / grid_x);
+                               m0 = m1 = 0.0;
+                               for (n = n0; n <= n1; n++) {
+                                       t = r.x - k - n * grid_x;
+                                       m0 += s_calc_spline(order, t);
+                                       m1 += (s_calc_spline(order - 1, t) - s_calc_spline(order - 1, t - 1)) * grid_x;
+                               }
                        }
                        mp[k * 2] = m0;
                        mp[k * 2 + 1] = m1;
                        if (debug_level > 1) fprintf(arena->debug_result, "mx[%d] (%g, %g)\n", k, m0, m1);
                }
                for (k = 0; k < grid_y; k++) {
-                       n0 = floor((r.y - order - k) / grid_y) + 1;
-                       n1 = floor((r.y - k) / grid_y);
-                       m0 = m1 = 0.0;
-                       for (n = n0; n <= n1; n++) {
-                               t = r.y - k - n * grid_y;
-                               m0 += s_calc_spline(order, t);
-                               m1 += (s_calc_spline(order - 1, t) - s_calc_spline(order - 1, t - 1)) * grid_y;
+                       if (grid_y == 1) {
+                               m0 = 1.0;
+                               m1 = 1.0;
+                       } else {
+                               n0 = floor((r.y - order - k) / grid_y) + 1;
+                               n1 = floor((r.y - k) / grid_y);
+                               m0 = m1 = 0.0;
+                               for (n = n0; n <= n1; n++) {
+                                       t = r.y - k - n * grid_y;
+                                       m0 += s_calc_spline(order, t);
+                                       m1 += (s_calc_spline(order - 1, t) - s_calc_spline(order - 1, t - 1)) * grid_y;
+                               }
                        }
                        mp[(k + grid_x) * 2] = m0;
                        mp[(k + grid_x) * 2 + 1] = m1;
                        if (debug_level > 1) fprintf(arena->debug_result, "my[%d] (%g, %g)\n", k, m0, m1);
                }
                for (k = 0; k < grid_z; k++) {
-                       n0 = floor((r.z - order - k) / grid_z) + 1;
-                       n1 = floor((r.z - k) / grid_z);
-                       m0 = m1 = 0.0;
-                       for (n = n0; n <= n1; n++) {
-                               t = r.z - k - n * grid_z;
-                               m0 += s_calc_spline(order, t);
-                               m1 += (s_calc_spline(order - 1, t) - s_calc_spline(order - 1, t - 1)) * grid_z;
+                       if (grid_z == 1) {
+                               m0 = 1.0;
+                               m1 = 1.0;
+                       } else {
+                               n0 = floor((r.z - order - k) / grid_z) + 1;
+                               n1 = floor((r.z - k) / grid_z);
+                               m0 = m1 = 0.0;
+                               for (n = n0; n <= n1; n++) {
+                                       t = r.z - k - n * grid_z;
+                                       m0 += s_calc_spline(order, t);
+                                       m1 += (s_calc_spline(order - 1, t) - s_calc_spline(order - 1, t - 1)) * grid_z;
+                               }
                        }
                        mp[(k + grid_x + grid_y) * 2] = m0;
                        mp[(k + grid_x + grid_y) * 2 + 1] = m1;
@@ -598,7 +632,7 @@ calc_pme_force(MDArena *arena)
                        }
                }
        }
-       *energy *= dielec_r / 0.5;
+       *energy *= dielec_r * 0.5;
        
        /*  qfarray is transformed by FFT into qqarray  */
        fftw_execute(pme->plan2);
@@ -677,6 +711,8 @@ calc_ewald_force(MDArena *arena)
                        calc_pme_force(arena);
                else if (arena->use_ewald == 2)
                        calc_direct_ewald_force(arena);
+               else if (arena->use_ewald == -1)   /*  Direct sum of electrostatic (very slow)  */
+                       calc_direct_electrostatic(arena);
                else return;
                arena->energies[kPMEIndex] += arena->pme->self_correction;
                if (arena->debug_result != NULL) {
@@ -735,11 +771,11 @@ s_pme_spline_test(MDArena *arena)
 #define DIM 32
 #define ORDER 8
        Int kk = DIM;
-       fftw_complex a1[DIM], a2[DIM], a3[DIM], a4[DIM], b[DIM], c1, c2, c3;
+       fftw_complex a1[DIM], a2[DIM], a3[DIM], a4[DIM], b[DIM], c1, c2;
        double b2[DIM], c[DIM];
        double d, w, x;
        double len = 10.0;  /*  Unit cell dimension  */
-       Int i, j, k, m, n0, n1;
+       Int k, m, n0, n1;
        x = 3.7;
        for (m = 0; m < DIM; m++) {
                k = (m <= DIM / 2 ? m : m - DIM);
@@ -907,7 +943,13 @@ pme_init(MDArena *arena)
                                pme->plan1 = fftw_plan_dft_3d(pme->grid_x, pme->grid_y, pme->grid_z, pme->qarray, pme->qfarray, FFTW_BACKWARD, FFTW_ESTIMATE);
                                pme->plan2 = fftw_plan_dft_3d(pme->grid_x, pme->grid_y, pme->grid_z, pme->qfarray, pme->qqarray, FFTW_FORWARD, FFTW_ESTIMATE);
                        }
-                       pme->order = 8;
+                       pme->order = arena->ewald_order;
+                       /*  Should be even numbers in [4..MAX_DIM_SPLINE]  */
+                       if (pme->order < 4)
+                               pme->order = 4;
+                       if (pme->order > MAX_DIM_SPLINE)
+                               pme->order = MAX_DIM_SPLINE;
+                       pme->order -= pme->order % 2;
                        /*  Initialize B array  */
                        for (i = 0; i < 3; i++) {
                                int idx_ofs, grid_i;
@@ -919,8 +961,9 @@ pme_init(MDArena *arena)
                                for (j = 0; j < grid_i; j++) {
                                        fftw_complex b1, b2;
                                        double t;
-                                       if (0 && j == grid_i / 2) {
-                                               b1[0] = b1[1] = 0.0;
+                                       if (grid_i == 1) {
+                                               b1[0] = 1.0;
+                                               b1[1] = 0.0;
                                        } else {
                                                t = 2 * PI * (pme->order - 1) * j / grid_i;
                                                b1[0] = cos(t);

diff --git a/MolLib/Ruby_bind/ruby_bind.c b/MolLib/Ruby_bind/ruby_bind.c
index c99066e..54bc819 100644 (file)
--- a/MolLib/Ruby_bind/ruby_bind.c
+++ b/MolLib/Ruby_bind/ruby_bind.c
@@ -4430,7 +4430,7 @@ static int
 s_Molecule_AtomIndexFromValue(Molecule *mol, VALUE val)
 {
        int n;
-       char *p;
+       char *p = "";
        if (FIXNUM_P(val)) {
                n = FIX2INT(val);
                if (n >= 0 && n < mol->natoms)
@@ -4442,12 +4442,13 @@ s_Molecule_AtomIndexFromValue(Molecule *mol, VALUE val)
        }
        if (n >= 0 && n < mol->natoms)
                return n;
+       p = StringValuePtr(val);
        if (n == -1)
-               rb_raise(rb_eMolbyError, "no such atom: %s", StringValuePtr(val));
+               rb_raise(rb_eMolbyError, "no such atom: %s", p);
        else if (n == -2)
-               rb_raise(rb_eMolbyError, "bad format of atom specification: %s", StringValuePtr(val));
+               rb_raise(rb_eMolbyError, "bad format of atom specification: %s", p);
        else
-               rb_raise(rb_eMolbyError, "error in converting value to atom index: %s", StringValuePtr(val));
+               rb_raise(rb_eMolbyError, "error in converting value to atom index: %s", p);
        return 0; /* Not reached */
 }
 

diff --git a/MolLib/Ruby_bind/ruby_md.c b/MolLib/Ruby_bind/ruby_md.c
index bc12290..967a1b2 100644 (file)
--- a/MolLib/Ruby_bind/ruby_md.c
+++ b/MolLib/Ruby_bind/ruby_md.c
@@ -334,7 +334,7 @@ s_DielectricSym, s_GradientConvergenceSym, s_CoordinateConvergenceSym, s_UseXplo
 s_Scale14VdwSym, s_Scale14ElectSym, s_RelocateCenterSym, 
 s_SurfaceProbeRadiusSym, s_SurfaceTensionSym, s_SurfacePotentialFreqSym, s_UseGraphiteSym,
 s_AlchemicalLambdaSym, s_AlchemicalDeltaLambdaSym, s_AlchemicalEnergySym, s_MinimizeCellSym,
-s_UseEwaldSym, s_EwaldBetaSym, s_EwaldGridSym, s_EwaldFreqSym;
+s_UseEwaldSym, s_EwaldBetaSym, s_EwaldGridSym, s_EwaldFreqSym, s_EwaldOrderSym;
 
 struct s_MDArenaAttrDef {
        char *name;
@@ -384,8 +384,9 @@ static struct s_MDArenaAttrDef s_MDArenaAttrDefTable[] = {
        {"minimize_cell",     &s_MinimizeCellSym,     0, 0, 'b', offsetof(MDArena, minimize_cell)},
        {"use_ewald",         &s_UseEwaldSym,         0, 0, 'i', offsetof(MDArena, use_ewald)},
        {"ewald_beta",        &s_EwaldBetaSym,        0, 0, 'f', offsetof(MDArena, ewald_beta)},
-       {"ewald_grid",        &s_EwaldGridSym,        0, 0, 0, offsetof(MDArena, ewald_grid_x)},
+       {"ewald_grid",        &s_EwaldGridSym,        0, 0, 0,   offsetof(MDArena, ewald_grid_x)},
        {"ewald_freq",        &s_EwaldFreqSym,        0, 0, 'i', offsetof(MDArena, ewald_freq)},
+       {"ewald_order",       &s_EwaldOrderSym,       0, 0, 'i', offsetof(MDArena, ewald_order)},
        {NULL} /* Sentinel */
 };