struct gmx_ewald_tab_t
{
int nx, ny, nz, kmax;
- cvec **eir;
+ cvec** eir;
t_complex *tab_xy, *tab_qxyz;
};
-void init_ewald_tab(struct gmx_ewald_tab_t **et, const t_inputrec *ir, FILE *fp)
+void init_ewald_tab(struct gmx_ewald_tab_t** et, const t_inputrec* ir, FILE* fp)
{
snew(*et, 1);
if (fp)
fprintf(fp, "Will do ordinary reciprocal space Ewald sum.\n");
}
- (*et)->nx = ir->nkx+1;
- (*et)->ny = ir->nky+1;
- (*et)->nz = ir->nkz+1;
+ (*et)->nx = ir->nkx + 1;
+ (*et)->ny = ir->nky + 1;
+ (*et)->nz = ir->nkz + 1;
(*et)->kmax = std::max((*et)->nx, std::max((*et)->ny, (*et)->nz));
(*et)->eir = nullptr;
(*et)->tab_xy = nullptr;
//! Calculates wave vectors.
static void calc_lll(const rvec box, rvec lll)
{
- lll[XX] = 2.0*M_PI/box[XX];
- lll[YY] = 2.0*M_PI/box[YY];
- lll[ZZ] = 2.0*M_PI/box[ZZ];
+ lll[XX] = 2.0 * M_PI / box[XX];
+ lll[YY] = 2.0 * M_PI / box[YY];
+ lll[ZZ] = 2.0 * M_PI / box[ZZ];
}
//! Make tables for the structure factor parts
-static void tabulateStructureFactors(int natom, const rvec x[], int kmax, cvec **eir, const rvec lll)
+static void tabulateStructureFactors(int natom, const rvec x[], int kmax, cvec** eir, const rvec lll)
{
- int i, j, m;
+ int i, j, m;
if (kmax < 1)
{
for (m = 0; (m < 3); m++)
{
- eir[1][i][m].re = std::cos(x[i][m]*lll[m]);
- eir[1][i][m].im = std::sin(x[i][m]*lll[m]);
+ eir[1][i][m].re = std::cos(x[i][m] * lll[m]);
+ eir[1][i][m].im = std::sin(x[i][m] * lll[m]);
}
for (j = 2; (j < kmax); j++)
{
for (m = 0; (m < 3); m++)
{
- eir[j][i][m] = cmul(eir[j-1][i][m], eir[1][i][m]);
+ eir[j][i][m] = cmul(eir[j - 1][i][m], eir[1][i][m]);
}
}
}
}
-real do_ewald(const t_inputrec *ir,
+real do_ewald(const t_inputrec* ir,
const rvec x[],
rvec f[],
const real chargeA[],
const real chargeB[],
const matrix box,
- const t_commrec *cr,
+ const t_commrec* cr,
int natoms,
matrix lrvir,
real ewaldcoeff,
real lambda,
- real *dvdlambda,
- gmx_ewald_tab_t *et)
+ real* dvdlambda,
+ gmx_ewald_tab_t* et)
{
- real factor = -1.0/(4*ewaldcoeff*ewaldcoeff);
- const real *charge;
+ real factor = -1.0 / (4 * ewaldcoeff * ewaldcoeff);
+ const real* charge;
real energy_AB[2], energy;
rvec lll;
int lowiy, lowiz, ix, iy, iz, n, q;
}
/* 1/(Vol*e0) */
- real scaleRecip = 4.0*M_PI/(boxDiag[XX]*boxDiag[YY]*boxDiag[ZZ])*ONE_4PI_EPS0/ir->epsilon_r;
+ real scaleRecip = 4.0 * M_PI / (boxDiag[XX] * boxDiag[YY] * boxDiag[ZZ]) * ONE_4PI_EPS0 / ir->epsilon_r;
if (!et->eir) /* allocate if we need to */
{
energy_AB[q] = 0;
for (ix = 0; ix < et->nx; ix++)
{
- mx = ix*lll[XX];
+ mx = ix * lll[XX];
for (iy = lowiy; iy < et->ny; iy++)
{
- my = iy*lll[YY];
+ my = iy * lll[YY];
if (iy >= 0)
{
for (n = 0; n < natoms; n++)
}
for (iz = lowiz; iz < et->nz; iz++)
{
- mz = iz*lll[ZZ];
- m2 = mx*mx+my*my+mz*mz;
- ak = std::exp(m2*factor)/m2;
- akv = 2.0*ak*(1.0/m2-factor);
+ mz = iz * lll[ZZ];
+ m2 = mx * mx + my * my + mz * mz;
+ ak = std::exp(m2 * factor) / m2;
+ akv = 2.0 * ak * (1.0 / m2 - factor);
if (iz >= 0)
{
for (n = 0; n < natoms; n++)
{
- et->tab_qxyz[n] = rcmul(charge[n], cmul(et->tab_xy[n],
- et->eir[iz][n][ZZ]));
+ et->tab_qxyz[n] = rcmul(charge[n], cmul(et->tab_xy[n], et->eir[iz][n][ZZ]));
}
}
else
{
for (n = 0; n < natoms; n++)
{
- et->tab_qxyz[n] = rcmul(charge[n], cmul(et->tab_xy[n],
- conjugate(et->eir[-iz][n][ZZ])));
+ et->tab_qxyz[n] = rcmul(
+ charge[n], cmul(et->tab_xy[n], conjugate(et->eir[-iz][n][ZZ])));
}
}
cs += et->tab_qxyz[n].re;
ss += et->tab_qxyz[n].im;
}
- energy_AB[q] += ak*(cs*cs+ss*ss);
- tmp = scale*akv*(cs*cs+ss*ss);
- lrvir[XX][XX] -= tmp*mx*mx;
- lrvir[XX][YY] -= tmp*mx*my;
- lrvir[XX][ZZ] -= tmp*mx*mz;
- lrvir[YY][YY] -= tmp*my*my;
- lrvir[YY][ZZ] -= tmp*my*mz;
- lrvir[ZZ][ZZ] -= tmp*mz*mz;
+ energy_AB[q] += ak * (cs * cs + ss * ss);
+ tmp = scale * akv * (cs * cs + ss * ss);
+ lrvir[XX][XX] -= tmp * mx * mx;
+ lrvir[XX][YY] -= tmp * mx * my;
+ lrvir[XX][ZZ] -= tmp * mx * mz;
+ lrvir[YY][YY] -= tmp * my * my;
+ lrvir[YY][ZZ] -= tmp * my * mz;
+ lrvir[ZZ][ZZ] -= tmp * mz * mz;
for (n = 0; n < natoms; n++)
{
/*tmp=scale*ak*(cs*tab_qxyz[n].im-ss*tab_qxyz[n].re);*/
- tmp = scale*ak*(cs*et->tab_qxyz[n].im-ss*et->tab_qxyz[n].re);
- f[n][XX] += tmp*mx*2*scaleRecip;
- f[n][YY] += tmp*my*2*scaleRecip;
- f[n][ZZ] += tmp*mz*2*scaleRecip;
+ tmp = scale * ak * (cs * et->tab_qxyz[n].im - ss * et->tab_qxyz[n].re);
+ f[n][XX] += tmp * mx * 2 * scaleRecip;
+ f[n][YY] += tmp * my * 2 * scaleRecip;
+ f[n][ZZ] += tmp * mz * 2 * scaleRecip;
#if 0
f[n][XX] += tmp*mx;
f[n][YY] += tmp*my;
f[n][ZZ] += tmp*mz;
#endif
}
- lowiz = 1-et->nz;
+ lowiz = 1 - et->nz;
}
- lowiy = 1-et->ny;
+ lowiy = 1 - et->ny;
}
}
}
}
else
{
- energy = (1.0 - lambda)*energy_AB[0] + lambda*energy_AB[1];
- *dvdlambda += scaleRecip*(energy_AB[1] - energy_AB[0]);
+ energy = (1.0 - lambda) * energy_AB[0] + lambda * energy_AB[1];
+ *dvdlambda += scaleRecip * (energy_AB[1] - energy_AB[0]);
}
- lrvir[XX][XX] = -0.5*scaleRecip*(lrvir[XX][XX]+energy);
- lrvir[XX][YY] = -0.5*scaleRecip*(lrvir[XX][YY]);
- lrvir[XX][ZZ] = -0.5*scaleRecip*(lrvir[XX][ZZ]);
- lrvir[YY][YY] = -0.5*scaleRecip*(lrvir[YY][YY]+energy);
- lrvir[YY][ZZ] = -0.5*scaleRecip*(lrvir[YY][ZZ]);
- lrvir[ZZ][ZZ] = -0.5*scaleRecip*(lrvir[ZZ][ZZ]+energy);
+ lrvir[XX][XX] = -0.5 * scaleRecip * (lrvir[XX][XX] + energy);
+ lrvir[XX][YY] = -0.5 * scaleRecip * (lrvir[XX][YY]);
+ lrvir[XX][ZZ] = -0.5 * scaleRecip * (lrvir[XX][ZZ]);
+ lrvir[YY][YY] = -0.5 * scaleRecip * (lrvir[YY][YY] + energy);
+ lrvir[YY][ZZ] = -0.5 * scaleRecip * (lrvir[YY][ZZ]);
+ lrvir[ZZ][ZZ] = -0.5 * scaleRecip * (lrvir[ZZ][ZZ] + energy);
lrvir[YY][XX] = lrvir[XX][YY];
lrvir[ZZ][XX] = lrvir[XX][ZZ];
return energy;
}
-real ewald_charge_correction(const t_commrec *cr,
- const t_forcerec *fr,
+real ewald_charge_correction(const t_commrec* cr,
+ const t_forcerec* fr,
const real lambda,
const matrix box,
- real *dvdlambda,
+ real* dvdlambda,
tensor vir)
{
if (MASTER(cr))
{
/* Apply charge correction */
- vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
+ vol = box[XX][XX] * box[YY][YY] * box[ZZ][ZZ];
- fac = M_PI*ONE_4PI_EPS0/(fr->ic->epsilon_r*2.0*vol*vol*gmx::square(fr->ic->ewaldcoeff_q));
+ fac = M_PI * ONE_4PI_EPS0
+ / (fr->ic->epsilon_r * 2.0 * vol * vol * gmx::square(fr->ic->ewaldcoeff_q));
- qs2A = fr->qsum[0]*fr->qsum[0];
- qs2B = fr->qsum[1]*fr->qsum[1];
+ qs2A = fr->qsum[0] * fr->qsum[0];
+ qs2B = fr->qsum[1] * fr->qsum[1];
- vc = (qs2A*(1 - lambda) + qs2B*lambda)*fac;
+ vc = (qs2A * (1 - lambda) + qs2B * lambda) * fac;
- enercorr = -vol*vc;
+ enercorr = -vol * vc;
- *dvdlambda += -vol*(qs2B - qs2A)*fac;
+ *dvdlambda += -vol * (qs2B - qs2A) * fac;
for (d = 0; d < DIM; d++)
{
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff