1734 lines
48 KiB
C++
1734 lines
48 KiB
C++
// clang-format off
|
|
/* ----------------------------------------------------------------------
|
|
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
|
https://www.lammps.org/ Sandia National Laboratories
|
|
LAMMPS development team: developers@lammps.org
|
|
|
|
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
|
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
|
certain rights in this software. This software is distributed under
|
|
the GNU General Public License.
|
|
|
|
See the README file in the top-level LAMMPS directory.
|
|
------------------------------------------------------------------------- */
|
|
|
|
#include "pair_amoeba.h"
|
|
|
|
#include "amoeba_convolution.h"
|
|
#include "atom.h"
|
|
#include "comm.h"
|
|
#include "domain.h"
|
|
#include "error.h"
|
|
#include "fix_store_atom.h"
|
|
#include "math_const.h"
|
|
#include "math_special.h"
|
|
#include "my_page.h"
|
|
#include "neigh_list.h"
|
|
#include "timer.h"
|
|
|
|
#include <cmath>
|
|
|
|
using namespace LAMMPS_NS;
|
|
using namespace MathConst;
|
|
|
|
using MathSpecial::cube;
|
|
|
|
enum{INDUCE,RSD,SETUP_AMOEBA,SETUP_HIPPO,KMPOLE,AMGROUP}; // forward comm
|
|
enum{FIELD,ZRSD,TORQUE,UFLD}; // reverse comm
|
|
enum{VDWL,REPULSE,QFER,DISP,MPOLE,POLAR,USOLV,DISP_LONG,MPOLE_LONG,POLAR_LONG};
|
|
enum{MUTUAL,OPT,TCG,DIRECT};
|
|
enum{GEAR,ASPC,LSQR};
|
|
enum{BUILD,APPLY};
|
|
enum{GORDON1,GORDON2};
|
|
|
|
static constexpr double DEBYE = 4.80321; // conversion factor from q-Angs (real units) to Debye
|
|
|
|
/* ----------------------------------------------------------------------
|
|
induce = induced dipole moments via pre-conditioned CG solver
|
|
adapted from Tinker induce0a() routine
|
|
------------------------------------------------------------------------- */
|
|
|
|
void PairAmoeba::induce()
|
|
{
|
|
bool done;
|
|
int i,j,m,itype,iter;
|
|
double polmin;
|
|
double eps,epsold;
|
|
double epsd,epsp;
|
|
double udsum,upsum;
|
|
double a,ap,b,bp;
|
|
double sum,sump,term;
|
|
double reduce[4],allreduce[4];
|
|
|
|
// set cutoffs, taper coeffs, and PME params
|
|
|
|
if (use_ewald) choose(POLAR_LONG);
|
|
else choose(POLAR);
|
|
|
|
// owned atoms
|
|
|
|
int nlocal = atom->nlocal;
|
|
|
|
// zero out the induced dipoles at each site
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
for (j = 0; j < 3; j++) {
|
|
uind[i][j] = 0.0;
|
|
uinp[i][j] = 0.0;
|
|
}
|
|
}
|
|
|
|
// get the electrostatic field due to permanent multipoles
|
|
|
|
dfield0c(field,fieldp);
|
|
|
|
// reverse comm to sum field,fieldp from ghost atoms to owned atoms
|
|
|
|
crstyle = FIELD;
|
|
comm->reverse_comm(this);
|
|
|
|
// set induced dipoles to polarizability times direct field
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
itype = amtype[i];
|
|
for (j = 0; j < 3; j++) {
|
|
udir[i][j] = polarity[itype] * field[i][j];
|
|
udirp[i][j] = polarity[itype] * fieldp[i][j];
|
|
if (pcgguess) {
|
|
uind[i][j] = udir[i][j];
|
|
uinp[i][j] = udirp[i][j];
|
|
}
|
|
}
|
|
}
|
|
|
|
// get induced dipoles via the OPT extrapolation method
|
|
// NOTE: could rewrite these loops to avoid allocating
|
|
// uopt,uoptp with a optorder+1 dimension, just optorder
|
|
// since no need to store optorder+1 values after these loops
|
|
|
|
if (poltyp == OPT) {
|
|
for (i = 0; i < nlocal; i++) {
|
|
for (j = 0; j < 3; j++) {
|
|
uopt[i][0][j] = udir[i][j];
|
|
uoptp[i][0][j] = udirp[i][j];
|
|
}
|
|
}
|
|
|
|
for (m = 1; m <= optorder; m++) {
|
|
optlevel = m - 1; // used in umutual1() for fopt,foptp
|
|
|
|
cfstyle = INDUCE;
|
|
comm->forward_comm(this);
|
|
|
|
ufield0c(field,fieldp);
|
|
|
|
crstyle = FIELD;
|
|
comm->reverse_comm(this);
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
itype = amtype[i];
|
|
for (j = 0; j < 3; j++) {
|
|
uopt[i][m][j] = polarity[itype] * field[i][j];
|
|
uoptp[i][m][j] = polarity[itype] * fieldp[i][j];
|
|
uind[i][j] = uopt[i][m][j];
|
|
uinp[i][j] = uoptp[i][m][j];
|
|
}
|
|
}
|
|
}
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
for (j = 0; j < 3; j++) {
|
|
uind[i][j] = 0.0;
|
|
uinp[i][j] = 0.0;
|
|
usum[i][j] = 0.0;
|
|
usump[i][j] = 0.0;
|
|
for (m = 0; m <= optorder; m++) {
|
|
usum[i][j] += uopt[i][m][j];
|
|
usump[i][j] += uoptp[i][m][j];
|
|
uind[i][j] += copt[m]*usum[i][j];
|
|
uinp[i][j] += copt[m]*usump[i][j];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// set tolerances for computation of mutual induced dipoles
|
|
|
|
if (poltyp == MUTUAL) {
|
|
done = false;
|
|
iter = 0;
|
|
polmin = 0.00000001;
|
|
eps = 100.0;
|
|
|
|
// estimate induced dipoles using a polynomial predictor
|
|
|
|
if (use_pred && nualt == maxualt) {
|
|
ulspred();
|
|
|
|
double ***udalt = fixudalt->tstore;
|
|
double ***upalt = fixupalt->tstore;
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
for (j = 0; j < 3; j++) {
|
|
udsum = 0.0;
|
|
upsum = 0.0;
|
|
for (m = 0; m < nualt; m++) {
|
|
udsum += bpred[m]*udalt[i][m][j];
|
|
upsum += bpredp[m]*upalt[i][m][j];
|
|
}
|
|
uind[i][j] = udsum;
|
|
uinp[i][j] = upsum;
|
|
}
|
|
}
|
|
}
|
|
|
|
// estimate induced dipoles via inertial extended Lagrangian
|
|
// not supported for now
|
|
// requires uaux,upaux to persist with each atom
|
|
// also requires a velocity vector(s) to persist
|
|
// also requires updating uaux,upaux in the Verlet integration
|
|
|
|
/*
|
|
if (use_ielscf) {
|
|
for (i = 0; i < nlocal; i++) {
|
|
for (j = 0; j < 3; j++) {
|
|
uind[i][j] = uaux[i][j];
|
|
uinp[i][j] = upaux[i][j];
|
|
}
|
|
}
|
|
}
|
|
*/
|
|
|
|
// get the electrostatic field due to induced dipoles
|
|
|
|
cfstyle = INDUCE;
|
|
comm->forward_comm(this);
|
|
|
|
ufield0c(field,fieldp);
|
|
|
|
crstyle = FIELD;
|
|
comm->reverse_comm(this);
|
|
|
|
// set initial conjugate gradient residual and conjugate vector
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
itype = amtype[i];
|
|
|
|
poli[i] = MAX(polmin,polarity[itype]);
|
|
for (j = 0; j < 3; j++) {
|
|
if (pcgguess) {
|
|
rsd[i][j] = (udir[i][j]-uind[i][j])/poli[i] + field[i][j];
|
|
rsdp[i][j] = (udirp[i][j]-uinp[i][j])/poli[i] + fieldp[i][j];
|
|
} else {
|
|
rsd[i][j] = udir[i][j] / poli[i];
|
|
rsdp[i][j] = udirp[i][j] / poli[i];
|
|
}
|
|
zrsd[i][j] = rsd[i][j];
|
|
zrsdp[i][j] = rsdp[i][j];
|
|
}
|
|
}
|
|
|
|
if (pcgprec) {
|
|
cfstyle = RSD;
|
|
comm->forward_comm(this);
|
|
uscale0b(BUILD,rsd,rsdp,zrsd,zrsdp);
|
|
uscale0b(APPLY,rsd,rsdp,zrsd,zrsdp);
|
|
crstyle = ZRSD;
|
|
comm->reverse_comm(this);
|
|
}
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
for (j = 0; j < 3; j++) {
|
|
conj[i][j] = zrsd[i][j];
|
|
conjp[i][j] = zrsdp[i][j];
|
|
}
|
|
}
|
|
|
|
// conjugate gradient iteration of the mutual induced dipoles
|
|
|
|
while (!done) {
|
|
iter++;
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
for (j = 0; j < 3; j++) {
|
|
vec[i][j] = uind[i][j];
|
|
vecp[i][j] = uinp[i][j];
|
|
uind[i][j] = conj[i][j];
|
|
uinp[i][j] = conjp[i][j];
|
|
}
|
|
}
|
|
|
|
cfstyle = INDUCE;
|
|
comm->forward_comm(this);
|
|
|
|
ufield0c(field,fieldp);
|
|
|
|
crstyle = FIELD;
|
|
comm->reverse_comm(this);
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
for (j = 0; j < 3; j++) {
|
|
uind[i][j] = vec[i][j];
|
|
uinp[i][j] = vecp[i][j];
|
|
vec[i][j] = conj[i][j]/poli[i] - field[i][j];
|
|
vecp[i][j] = conjp[i][j]/poli[i] - fieldp[i][j];
|
|
}
|
|
}
|
|
|
|
a = 0.0;
|
|
ap = 0.0;
|
|
sum = 0.0;
|
|
sump = 0.0;
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
for (j = 0; j < 3; j++) {
|
|
a += conj[i][j]*vec[i][j];
|
|
ap += conjp[i][j]*vecp[i][j];
|
|
sum += rsd[i][j]*zrsd[i][j];
|
|
sump += rsdp[i][j]*zrsdp[i][j];
|
|
}
|
|
}
|
|
|
|
reduce[0] = a;
|
|
reduce[1] = ap;
|
|
reduce[2] = sum;
|
|
reduce[3] = sump;
|
|
MPI_Allreduce(reduce,allreduce,4,MPI_DOUBLE,MPI_SUM,world);
|
|
a = allreduce[0];
|
|
ap = allreduce[1];
|
|
sum = allreduce[2];
|
|
sump = allreduce[3];
|
|
|
|
if (a != 0.0) a = sum / a;
|
|
if (ap != 0.0) ap = sump / ap;
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
for (j = 0; j < 3; j++) {
|
|
uind[i][j] = uind[i][j] + a*conj[i][j];
|
|
uinp[i][j] = uinp[i][j] + ap*conjp[i][j];
|
|
rsd[i][j] = rsd[i][j] - a*vec[i][j];
|
|
rsdp[i][j] = rsdp[i][j] - ap*vecp[i][j];
|
|
zrsd[i][j] = rsd[i][j];
|
|
zrsdp[i][j] = rsdp[i][j];
|
|
}
|
|
}
|
|
|
|
if (pcgprec) {
|
|
cfstyle = RSD;
|
|
comm->forward_comm(this);
|
|
uscale0b(APPLY,rsd,rsdp,zrsd,zrsdp);
|
|
crstyle = ZRSD;
|
|
comm->reverse_comm(this);
|
|
}
|
|
|
|
b = 0.0;
|
|
bp = 0.0;
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
for (j = 0; j < 3; j++) {
|
|
b += rsd[i][j]*zrsd[i][j];
|
|
bp += rsdp[i][j]*zrsdp[i][j];
|
|
}
|
|
}
|
|
|
|
reduce[0] = b;
|
|
reduce[1] = bp;
|
|
MPI_Allreduce(reduce,allreduce,4,MPI_DOUBLE,MPI_SUM,world);
|
|
b = allreduce[0];
|
|
bp = allreduce[1];
|
|
|
|
if (sum != 0.0) b /= sum;
|
|
if (sump != 0.0) bp /= sump;
|
|
|
|
epsd = 0.0;
|
|
epsp = 0.0;
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
for (j = 0; j < 3; j++) {
|
|
conj[i][j] = zrsd[i][j] + b*conj[i][j];
|
|
conjp[i][j] = zrsdp[i][j] + bp*conjp[i][j];
|
|
epsd += rsd[i][j]*rsd[i][j];
|
|
epsp += rsdp[i][j]*rsdp[i][j];
|
|
}
|
|
}
|
|
|
|
reduce[0] = epsd;
|
|
reduce[1] = epsp;
|
|
MPI_Allreduce(reduce,allreduce,4,MPI_DOUBLE,MPI_SUM,world);
|
|
epsd = allreduce[0];
|
|
epsp = allreduce[1];
|
|
|
|
// check the convergence of the mutual induced dipoles
|
|
|
|
epsold = eps;
|
|
eps = MAX(epsd,epsp);
|
|
eps = DEBYE * sqrt(eps/atom->natoms);
|
|
|
|
if (eps < poleps) done = true;
|
|
// also commented out in induce.f of Tinker
|
|
// if (eps > epsold) done = true;
|
|
if (iter >= politer) done = true;
|
|
|
|
// apply a "peek" iteration to the mutual induced dipoles
|
|
|
|
if (done) {
|
|
for (i = 0; i < nlocal; i++) {
|
|
term = pcgpeek * poli[i];
|
|
for (j = 0; j < 3; j++) {
|
|
uind[i][j] += term*rsd[i][j];
|
|
uinp[i][j] += term*rsdp[i][j];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// terminate the calculation if dipoles failed to converge
|
|
// NOTE: could make this an error
|
|
|
|
if (iter >= politer || eps > epsold)
|
|
if (comm->me == 0)
|
|
error->warning(FLERR,"AMOEBA induced dipoles did not converge");
|
|
}
|
|
|
|
// update the lists of previous induced dipole values
|
|
// shift previous m values up to m+1, add new values at m = 0
|
|
// only when preconditioner is used
|
|
|
|
if (use_pred) {
|
|
double ***udalt = fixudalt->tstore;
|
|
double ***upalt = fixupalt->tstore;
|
|
|
|
nualt = MIN(nualt+1,maxualt);
|
|
for (i = 0; i < nlocal; i++) {
|
|
for (j = 0; j < 3; j++) {
|
|
for (m = nualt-1; m > 0; m--) {
|
|
udalt[i][m][j] = udalt[i][m-1][j];
|
|
upalt[i][m][j] = upalt[i][m-1][j];
|
|
}
|
|
udalt[i][0][j] = uind[i][j];
|
|
upalt[i][0][j] = uinp[i][j];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
ulspred = induced dipole prediction coeffs
|
|
|
|
ulspred uses standard extrapolation or a least squares fit
|
|
to set coefficients of an induced dipole predictor polynomial
|
|
literature references:
|
|
|
|
J. Kolafa, "Time-Reversible Always Stable Predictor-Corrector
|
|
Method for Molecular Dynamics of Polarizable Molecules", Journal
|
|
of Computational Chemistry, 25, 335-342 (2004)
|
|
|
|
W. Wang and R. D. Skeel, "Fast Evaluation of Polarizable Forces",
|
|
Journal of Chemical Physics, 123, 164107 (2005)
|
|
------------------------------------------------------------------------- */
|
|
|
|
void PairAmoeba::ulspred()
|
|
{
|
|
int i,j,k,m;
|
|
double coeff,udk,upk;
|
|
double amax,apmax;
|
|
|
|
// set the Gear predictor binomial coefficients
|
|
|
|
if (polpred == GEAR) {
|
|
for (i = 0; i < nualt; i++) {
|
|
coeff = gear[i];
|
|
bpred[i] = coeff;
|
|
bpredp[i] = coeff;
|
|
bpreds[i] = coeff;
|
|
bpredps[i] = coeff;
|
|
}
|
|
|
|
// set always stable predictor-corrector (ASPC) coefficients
|
|
|
|
} else if (polpred == ASPC) {
|
|
for (i = 0; i < nualt; i++) {
|
|
coeff = aspc[i];
|
|
bpred[i] = coeff;
|
|
bpredp[i] = coeff;
|
|
bpreds[i] = coeff;
|
|
bpredps[i] = coeff;
|
|
}
|
|
|
|
// derive normal equations corresponding to least squares fit
|
|
|
|
} else if (polpred == LSQR) {
|
|
double ***udalt = fixudalt->tstore;
|
|
double ***upalt = fixupalt->tstore;
|
|
|
|
for (k = 0; k < nualt; k++) {
|
|
b_ualt[k] = 0.0;
|
|
bp_ualt[k] = 0.0;
|
|
for (m = k; m < nualt; m++) {
|
|
c_ualt[k][m] = 0.0;
|
|
cp_ualt[k][m] = 0.0;
|
|
}
|
|
}
|
|
|
|
int nlocal = atom->nlocal;
|
|
for (i = 0; i < nlocal; i++) {
|
|
for (j = 0; j < 3; j++) {
|
|
for (k = 0; k < nualt; k++) {
|
|
udk = udalt[i][k][j];
|
|
upk = upalt[i][k][j];
|
|
for (m = k; m < nualt; m++) {
|
|
c_ualt[k][m] += udk*udalt[i][m][j];
|
|
cp_ualt[k][m] += upk*upalt[i][m][j];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
i = 0;
|
|
for (k = 1; k < nualt; k++) {
|
|
b_ualt[k-1] = c_ualt[0][k];
|
|
bp_ualt[k-1] = cp_ualt[0][k];
|
|
for (m = k; m < nualt; m++) {
|
|
a_ualt[i] = c_ualt[k][m];
|
|
ap_ualt[i] = cp_ualt[k][m];
|
|
i++;
|
|
}
|
|
}
|
|
|
|
// check for nonzero coefficients and solve normal equations
|
|
|
|
k = nualt - 1;
|
|
amax = 0.0;
|
|
apmax = 0.0;
|
|
for (i = 0; i < k*(k+1)/2; i++) {
|
|
amax = MAX(amax,a_ualt[i]);
|
|
apmax = MAX(apmax,ap_ualt[i]);
|
|
}
|
|
if (amax != 0.0) cholesky(nualt-1,a_ualt,b_ualt);
|
|
if (apmax != 0.0) cholesky(nualt-1,ap_ualt,bp_ualt);
|
|
|
|
// transfer the final solution to the coefficient vector
|
|
|
|
for (k = 0; k < nualt-1; k++) {
|
|
bpred[k] = b_ualt[k];
|
|
bpredp[k] = bp_ualt[k];
|
|
bpreds[k] = b_ualt[k];
|
|
bpredps[k] = bp_ualt[k];
|
|
}
|
|
bpred[nualt-1] = 0.0;
|
|
bpredp[nualt-1] = 0.0;
|
|
bpreds[nualt-1] = 0.0;
|
|
bpredps[nualt-1] = 0.0;
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
ufield0c = mutual induction via Ewald sum
|
|
ufield0c computes the mutual electrostatic field due to
|
|
induced dipole moments via Ewald summation
|
|
------------------------------------------------------------------------- */
|
|
|
|
void PairAmoeba::ufield0c(double **field, double **fieldp)
|
|
{
|
|
int i,j;
|
|
double term;
|
|
|
|
// zero field,fieldp for owned and ghost atoms
|
|
|
|
int nlocal = atom->nlocal;
|
|
int nall = nlocal + atom->nghost;
|
|
|
|
for (i = 0; i < nall; i++) {
|
|
for (j = 0; j < 3; j++) {
|
|
field[i][j] = 0.0;
|
|
fieldp[i][j] = 0.0;
|
|
}
|
|
}
|
|
|
|
double time0, time1, time2;
|
|
if (timer->has_sync()) MPI_Barrier(world);
|
|
time0 = platform::walltime();
|
|
|
|
// get the real space portion of the mutual field
|
|
|
|
if (polar_rspace_flag) umutual2b(field,fieldp);
|
|
time1 = platform::walltime();
|
|
|
|
// get the reciprocal space part of the mutual field
|
|
|
|
if (polar_kspace_flag) umutual1(field,fieldp);
|
|
time2 = platform::walltime();
|
|
|
|
// add the self-energy portion of the mutual field
|
|
|
|
term = (4.0/3.0) * aewald*aewald*aewald / MY_PIS;
|
|
for (i = 0; i < nlocal; i++) {
|
|
for (j = 0; j < 3; j++) {
|
|
field[i][j] += term*uind[i][j];
|
|
fieldp[i][j] += term*uinp[i][j];
|
|
}
|
|
}
|
|
|
|
// accumulate timing information
|
|
|
|
time_mutual_rspace += time1 - time0;
|
|
time_mutual_kspace += time2 - time1;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
uscale0b = dipole preconditioner via neigh list
|
|
uscale0b builds and applies a preconditioner for the conjugate
|
|
gradient induced dipole solver using a neighbor pair list
|
|
------------------------------------------------------------------------- */
|
|
|
|
void PairAmoeba::uscale0b(int mode, double **rsd, double **rsdp,
|
|
double **zrsd, double **zrsdp)
|
|
{
|
|
int i,j,itype,jtype,iclass,jclass,igroup,jgroup;
|
|
int ii,jj;
|
|
double xi,yi,zi;
|
|
double xr,yr,zr;
|
|
double r,r2,rr3,rr5;
|
|
double pdi,pti;
|
|
double polmin;
|
|
double poli,polik;
|
|
double alphai,alphak;
|
|
double damp,expdamp;
|
|
double pgamma;
|
|
double scale3,scale5;
|
|
double m1,m2,m3;
|
|
double m4,m5,m6;
|
|
double factor_uscale,factor_wscale;
|
|
double dmpik[5];
|
|
|
|
// owned atoms
|
|
|
|
double **x = atom->x;
|
|
int nlocal = atom->nlocal;
|
|
int nall = nlocal + atom->nghost;
|
|
|
|
// neighbor list info
|
|
|
|
int inum,jnum;
|
|
int *ilist,*jlist;
|
|
double *pclist;
|
|
|
|
inum = list->inum;
|
|
ilist = list->ilist;
|
|
|
|
// ------------------------------------------------
|
|
// apply the preconditioning matrix to the current residual
|
|
// ------------------------------------------------
|
|
|
|
if (mode == APPLY) {
|
|
|
|
// use diagonal preconditioner elements as first approximation
|
|
|
|
polmin = 0.00000001;
|
|
for (i = 0; i < nlocal; i++) {
|
|
itype = amtype[i];
|
|
poli = udiag * MAX(polmin,polarity[itype]);
|
|
for (j = 0; j < 3; j++) {
|
|
zrsd[i][j] = poli * rsd[i][j];
|
|
zrsdp[i][j] = poli * rsdp[i][j];
|
|
}
|
|
}
|
|
|
|
// zero zrsd,zrsdp for ghost atoms only
|
|
|
|
for (i = nlocal; i < nall; i++) {
|
|
for (j = 0; j < 3; j++) {
|
|
zrsd[i][j] = 0.0;
|
|
zrsdp[i][j] = 0.0;
|
|
}
|
|
}
|
|
|
|
// use the off-diagonal preconditioner elements in second phase
|
|
|
|
for (ii = 0; ii < inum; ii++) {
|
|
i = ilist[ii];
|
|
jlist = firstneigh_precond[i];
|
|
jnum = numneigh_precond[i];
|
|
pclist = firstneigh_pcpc[i];
|
|
|
|
for (jj = 0; jj < jnum; jj++) {
|
|
j = jlist[jj];
|
|
j &= NEIGHMASK15;
|
|
|
|
m1 = pclist[0];
|
|
m2 = pclist[1];
|
|
m3 = pclist[2];
|
|
m4 = pclist[3];
|
|
m5 = pclist[4];
|
|
m6 = pclist[5];
|
|
|
|
zrsd[i][0] += m1*rsd[j][0] + m2*rsd[j][1] + m3*rsd[j][2];
|
|
zrsd[i][1] += m2*rsd[j][0] + m4*rsd[j][1] + m5*rsd[j][2];
|
|
zrsd[i][2] += m3*rsd[j][0] + m5*rsd[j][1] + m6*rsd[j][2];
|
|
zrsd[j][0] += m1*rsd[i][0] + m2*rsd[i][1] + m3*rsd[i][2];
|
|
zrsd[j][1] += m2*rsd[i][0] + m4*rsd[i][1] + m5*rsd[i][2];
|
|
zrsd[j][2] += m3*rsd[i][0] + m5*rsd[i][1] + m6*rsd[i][2];
|
|
zrsdp[i][0] += m1*rsdp[j][0] + m2*rsdp[j][1] + m3*rsdp[j][2];
|
|
zrsdp[i][1] += m2*rsdp[j][0] + m4*rsdp[j][1] + m5*rsdp[j][2];
|
|
zrsdp[i][2] += m3*rsdp[j][0] + m5*rsdp[j][1] + m6*rsdp[j][2];
|
|
zrsdp[j][0] += m1*rsdp[i][0] + m2*rsdp[i][1] + m3*rsdp[i][2];
|
|
zrsdp[j][1] += m2*rsdp[i][0] + m4*rsdp[i][1] + m5*rsdp[i][2];
|
|
zrsdp[j][2] += m3*rsdp[i][0] + m5*rsdp[i][1] + m6*rsdp[i][2];
|
|
|
|
pclist += 6;
|
|
}
|
|
}
|
|
|
|
return;
|
|
}
|
|
|
|
// ------------------------------------------------
|
|
// build the off-diagonal elements of preconditioning matrix
|
|
// ------------------------------------------------
|
|
|
|
dpage_pcpc->reset();
|
|
|
|
// determine the off-diagonal elements of the preconditioner
|
|
|
|
for (ii = 0; ii < inum; ii++) {
|
|
i = ilist[ii];
|
|
itype = amtype[i];
|
|
iclass = amtype2class[itype];
|
|
igroup = amgroup[i];
|
|
|
|
jlist = firstneigh_precond[i];
|
|
jnum = numneigh_precond[i];
|
|
|
|
xi = x[i][0];
|
|
yi = x[i][1];
|
|
zi = x[i][2];
|
|
|
|
poli = polarity[itype];
|
|
if (amoeba) {
|
|
pdi = pdamp[itype];
|
|
pti = thole[itype];
|
|
} else {
|
|
alphai = palpha[iclass];
|
|
}
|
|
|
|
// evaluate all sites in induce neigh list, no cutoff
|
|
// store results in plist for re-use in APPLY
|
|
|
|
pclist = dpage_pcpc->get(6*jnum);
|
|
firstneigh_pcpc[i] = pclist;
|
|
|
|
for (jj = 0; jj < jnum; jj++) {
|
|
j = jlist[jj];
|
|
factor_wscale = special_polar_wscale[sbmask15(j)];
|
|
j &= NEIGHMASK15;
|
|
|
|
xr = x[j][0] - xi;
|
|
yr = x[j][1] - yi;
|
|
zr = x[j][2] - zi;
|
|
r2 = xr*xr + yr* yr + zr*zr;
|
|
r = sqrt(r2);
|
|
|
|
jtype = amtype[j];
|
|
jclass = amtype2class[jtype];
|
|
jgroup = amgroup[j];
|
|
|
|
if (igroup == jgroup) factor_uscale = polar_uscale;
|
|
else factor_uscale = 1.0;
|
|
|
|
if (amoeba) {
|
|
scale3 = factor_uscale;
|
|
scale5 = factor_uscale;
|
|
damp = pdi * pdamp[jtype];
|
|
if (damp != 0.0) {
|
|
pgamma = MIN(pti,thole[jtype]);
|
|
damp = -pgamma * cube(r/damp);
|
|
if (damp > -50.0) {
|
|
expdamp = exp(damp);
|
|
scale3 *= 1.0 - expdamp;
|
|
scale5 *= 1.0 - expdamp*(1.0-damp);
|
|
}
|
|
}
|
|
} else {
|
|
alphak = palpha[jclass];
|
|
dampmut(r,alphai,alphak,dmpik);
|
|
scale3 = factor_wscale * dmpik[2];
|
|
scale5 = factor_wscale * dmpik[4];
|
|
}
|
|
|
|
polik = poli * polarity[jtype];
|
|
rr3 = scale3 * polik / (r*r2);
|
|
rr5 = 3.0 * scale5 * polik / (r*r2*r2);
|
|
|
|
pclist[0] = rr5*xr*xr - rr3;
|
|
pclist[1] = rr5*xr*yr;
|
|
pclist[2] = rr5*xr*zr;
|
|
pclist[3] = rr5*yr*yr - rr3;
|
|
pclist[4] = rr5*yr*zr;
|
|
pclist[5] = rr5*zr*zr - rr3;
|
|
pclist += 6;
|
|
}
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
dfield0c = direct induction via Ewald sum
|
|
dfield0c computes the mutual electrostatic field due to
|
|
permanent multipole moments via Ewald summation
|
|
------------------------------------------------------------------------- */
|
|
|
|
void PairAmoeba::dfield0c(double **field, double **fieldp)
|
|
{
|
|
int i,j;
|
|
double term;
|
|
|
|
// zero out field,fieldp for owned and ghost atoms
|
|
|
|
int nlocal = atom->nlocal;
|
|
int nall = nlocal + atom->nghost;
|
|
|
|
for (i = 0; i < nall; i++) {
|
|
for (j = 0; j < 3; j++) {
|
|
field[i][j] = 0.0;
|
|
fieldp[i][j] = 0.0;
|
|
}
|
|
}
|
|
|
|
// get the reciprocal space part of the permanent field
|
|
|
|
double time0, time1, time2;
|
|
if (timer->has_sync()) MPI_Barrier(world);
|
|
time0 = platform::walltime();
|
|
|
|
if (polar_kspace_flag) udirect1(field);
|
|
time1 = platform::walltime();
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
for (j = 0; j < 3; j++) {
|
|
fieldp[i][j] = field[i][j];
|
|
}
|
|
}
|
|
|
|
// get the real space portion of the permanent field
|
|
|
|
if (polar_rspace_flag) udirect2b(field,fieldp);
|
|
time2 = platform::walltime();
|
|
|
|
// get the self-energy portion of the permanent field
|
|
|
|
term = (4.0/3.0) * aewald*aewald*aewald / MY_PIS;
|
|
for (i = 0; i < nlocal; i++) {
|
|
for (j = 0; j < 3; j++) {
|
|
field[i][j] += term*rpole[i][j+1];
|
|
fieldp[i][j] += term*rpole[i][j+1];
|
|
}
|
|
}
|
|
|
|
// accumulate timing information
|
|
|
|
time_direct_kspace += time1 - time0;
|
|
time_direct_rspace += time2 - time1;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
umutual1 = Ewald recip mutual induced field
|
|
umutual1 computes the reciprocal space contribution of the
|
|
induced atomic dipole moments to the field
|
|
------------------------------------------------------------------------- */
|
|
|
|
void PairAmoeba::umutual1(double **field, double **fieldp)
|
|
{
|
|
int i,j,k,m,n;
|
|
int nxlo,nxhi,nylo,nyhi,nzlo,nzhi;
|
|
double term;
|
|
double a[3][3]; // indices not flipped vs Fortran
|
|
|
|
// return if the Ewald coefficient is zero
|
|
|
|
if (aewald < 1.0e-6) return;
|
|
|
|
// convert Cartesian dipoles to fractional coordinates
|
|
|
|
for (j = 0; j < 3; j++) {
|
|
a[0][j] = nfft1 * recip[0][j];
|
|
a[1][j] = nfft2 * recip[1][j];
|
|
a[2][j] = nfft3 * recip[2][j];
|
|
}
|
|
|
|
int nlocal = atom->nlocal;
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
for (j = 0; j < 3; j++) {
|
|
fuind[i][j] = a[j][0]*uind[i][0] + a[j][1]*uind[i][1] + a[j][2]*uind[i][2];
|
|
fuinp[i][j] = a[j][0]*uinp[i][0] + a[j][1]*uinp[i][1] + a[j][2]*uinp[i][2];
|
|
}
|
|
}
|
|
|
|
double time0, time1;
|
|
|
|
// gridpre = my portion of 4d grid in brick decomp w/ ghost values
|
|
|
|
FFT_SCALAR ****gridpre = (FFT_SCALAR ****) ic_kspace->zero();
|
|
|
|
// map 2 values to grid
|
|
|
|
if (timer->has_sync()) MPI_Barrier(world);
|
|
time0 = platform::walltime();
|
|
|
|
grid_uind(fuind,fuinp,gridpre);
|
|
|
|
time1 = platform::walltime();
|
|
time_grid_uind += (time1 - time0);
|
|
|
|
// pre-convolution operations including forward FFT
|
|
// gridfft = my portion of complex 3d grid in FFT decomposition
|
|
|
|
FFT_SCALAR *gridfft = ic_kspace->pre_convolution();
|
|
|
|
// ---------------------
|
|
// convolution operation
|
|
// ---------------------
|
|
|
|
nxlo = ic_kspace->nxlo_fft;
|
|
nxhi = ic_kspace->nxhi_fft;
|
|
nylo = ic_kspace->nylo_fft;
|
|
nyhi = ic_kspace->nyhi_fft;
|
|
nzlo = ic_kspace->nzlo_fft;
|
|
nzhi = ic_kspace->nzhi_fft;
|
|
|
|
// use qfac values stored in udirect1()
|
|
|
|
m = n = 0;
|
|
for (k = nzlo; k <= nzhi; k++) {
|
|
for (j = nylo; j <= nyhi; j++) {
|
|
for (i = nxlo; i <= nxhi; i++) {
|
|
term = qfac[m++];
|
|
gridfft[n] *= term;
|
|
gridfft[n+1] *= term;
|
|
n += 2;
|
|
}
|
|
}
|
|
}
|
|
|
|
// post-convolution operations including backward FFT
|
|
// gridppost = my portion of 4d grid in brick decomp w/ ghost values
|
|
|
|
FFT_SCALAR ****gridpost = (FFT_SCALAR ****) ic_kspace->post_convolution();
|
|
|
|
// get potential
|
|
|
|
if (timer->has_sync()) MPI_Barrier(world);
|
|
time0 = platform::walltime();
|
|
|
|
fphi_uind(gridpost,fdip_phi1,fdip_phi2,fdip_sum_phi);
|
|
|
|
time1 = platform::walltime();
|
|
time_fphi_uind += (time1 - time0);
|
|
|
|
// store fractional reciprocal potentials for OPT method
|
|
|
|
if (poltyp == OPT) {
|
|
for (i = 0; i < nlocal; i++) {
|
|
for (j = 0; j < 10; j++) {
|
|
fopt[i][optlevel][j] = fdip_phi1[i][j];
|
|
foptp[i][optlevel][j] = fdip_phi2[i][j];
|
|
}
|
|
}
|
|
}
|
|
|
|
// convert the dipole fields from fractional to Cartesian
|
|
|
|
for (i = 0; i < 3; i++) {
|
|
a[0][i] = nfft1 * recip[0][i];
|
|
a[1][i] = nfft2 * recip[1][i];
|
|
a[2][i] = nfft3 * recip[2][i];
|
|
}
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
for (j = 0; j < 3; j++) {
|
|
dipfield1[i][j] = a[j][0]*fdip_phi1[i][1] +
|
|
a[j][1]*fdip_phi1[i][2] + a[j][2]*fdip_phi1[i][3];
|
|
dipfield2[i][j] = a[j][0]*fdip_phi2[i][1] +
|
|
a[j][1]*fdip_phi2[i][2] + a[j][2]*fdip_phi2[i][3];
|
|
}
|
|
}
|
|
|
|
// increment the field at each multipole site
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
for (j = 0; j < 3; j++) {
|
|
field[i][j] -= dipfield1[i][j];
|
|
fieldp[i][j] -= dipfield2[i][j];
|
|
}
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
umutual2b = Ewald real mutual field via list
|
|
umutual2b computes the real space contribution of the induced
|
|
atomic dipole moments to the field via a neighbor list
|
|
------------------------------------------------------------------------- */
|
|
|
|
void PairAmoeba::umutual2b(double **field, double **fieldp)
|
|
{
|
|
int i,j,m,ii,jj,jnum;
|
|
double fid[3],fkd[3];
|
|
double fip[3],fkp[3];
|
|
double *uindi,*uindj,*uinpi,*uinpj;
|
|
|
|
// neigh list
|
|
|
|
int inum = list->inum;
|
|
int *ilist = list->ilist;
|
|
int *jlist;
|
|
double *tdipdip;
|
|
|
|
// loop over owned atoms and neighs
|
|
// compute field terms for each pairwise interaction
|
|
// using tdipdip values stored by udirect2b()
|
|
|
|
for (ii = 0; ii < inum; ii++) {
|
|
i = ilist[ii];
|
|
uindi = uind[i];
|
|
uinpi = uinp[i];
|
|
jlist = firstneigh_dipole[i];
|
|
tdipdip = firstneigh_dipdip[i];
|
|
jnum = numneigh_dipole[i];
|
|
|
|
for (jj = 0; jj < jnum; jj++) {
|
|
j = jlist[jj];
|
|
uindj = uind[j];
|
|
uinpj = uinp[j];
|
|
|
|
fid[0] = tdipdip[0]*uindj[0] + tdipdip[1]*uindj[1] + tdipdip[2]*uindj[2];
|
|
fid[1] = tdipdip[1]*uindj[0] + tdipdip[3]*uindj[1] + tdipdip[4]*uindj[2];
|
|
fid[2] = tdipdip[2]*uindj[0] + tdipdip[4]*uindj[1] + tdipdip[5]*uindj[2];
|
|
|
|
fkd[0] = tdipdip[0]*uindi[0] + tdipdip[1]*uindi[1] + tdipdip[2]*uindi[2];
|
|
fkd[1] = tdipdip[1]*uindi[0] + tdipdip[3]*uindi[1] + tdipdip[4]*uindi[2];
|
|
fkd[2] = tdipdip[2]*uindi[0] + tdipdip[4]*uindi[1] + tdipdip[5]*uindi[2];
|
|
|
|
fip[0] = tdipdip[0]*uinpj[0] + tdipdip[1]*uinpj[1] + tdipdip[2]*uinpj[2];
|
|
fip[1] = tdipdip[1]*uinpj[0] + tdipdip[3]*uinpj[1] + tdipdip[4]*uinpj[2];
|
|
fip[2] = tdipdip[2]*uinpj[0] + tdipdip[4]*uinpj[1] + tdipdip[5]*uinpj[2];
|
|
|
|
fkp[0] = tdipdip[0]*uinpi[0] + tdipdip[1]*uinpi[1] + tdipdip[2]*uinpi[2];
|
|
fkp[1] = tdipdip[1]*uinpi[0] + tdipdip[3]*uinpi[1] + tdipdip[4]*uinpi[2];
|
|
fkp[2] = tdipdip[2]*uinpi[0] + tdipdip[4]*uinpi[1] + tdipdip[5]*uinpi[2];
|
|
|
|
tdipdip += 6;
|
|
|
|
// increment the field at each site due to this interaction
|
|
|
|
for (m = 0; m < 3; m++) {
|
|
field[i][m] += fid[m];
|
|
field[j][m] += fkd[m];
|
|
fieldp[i][m] += fip[m];
|
|
fieldp[j][m] += fkp[m];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
udirect1 = Ewald recip direct induced field
|
|
udirect1 computes the reciprocal space contribution of the
|
|
permanent atomic multipole moments to the field
|
|
since corresponding values in empole and epolar are different
|
|
------------------------------------------------------------------------- */
|
|
|
|
void PairAmoeba::udirect1(double **field)
|
|
{
|
|
int i,j,k,m,n;
|
|
int nhalf1,nhalf2,nhalf3;
|
|
int nxlo,nxhi,nylo,nyhi,nzlo,nzhi;
|
|
double r1,r2,r3;
|
|
double h1,h2,h3;
|
|
double volterm,denom;
|
|
double hsq,expterm;
|
|
double term,pterm;
|
|
|
|
// return if the Ewald coefficient is zero
|
|
|
|
if (aewald < 1.0e-6) return;
|
|
|
|
pterm = (MY_PI/aewald) * (MY_PI/aewald);
|
|
double volbox = domain->prd[0] * domain->prd[1] * domain->prd[2];
|
|
volterm = MY_PI * volbox;
|
|
|
|
// FFT moduli pre-computations
|
|
// set igrid for each atom and its B-spline coeffs
|
|
|
|
nfft1 = i_kspace->nx;
|
|
nfft2 = i_kspace->ny;
|
|
nfft3 = i_kspace->nz;
|
|
bsorder = i_kspace->order;
|
|
|
|
moduli();
|
|
bspline_fill();
|
|
|
|
// copy the multipole moments into local storage areas
|
|
|
|
int nlocal = atom->nlocal;
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
cmp[i][0] = rpole[i][0];
|
|
cmp[i][1] = rpole[i][1];
|
|
cmp[i][2] = rpole[i][2];
|
|
cmp[i][3] = rpole[i][3];
|
|
cmp[i][4] = rpole[i][4];
|
|
cmp[i][5] = rpole[i][8];
|
|
cmp[i][6] = rpole[i][12];
|
|
cmp[i][7] = 2.0 * rpole[i][5];
|
|
cmp[i][8] = 2.0 * rpole[i][6];
|
|
cmp[i][9] = 2.0 * rpole[i][9];
|
|
}
|
|
|
|
// convert Cartesian multipoles to fractional coordinates
|
|
|
|
cmp_to_fmp(cmp,fmp);
|
|
|
|
// gridpre = my portion of 3d grid in brick decomp w/ ghost values
|
|
// zeroed by setup()
|
|
|
|
FFT_SCALAR ***gridpre = (FFT_SCALAR ***) i_kspace->zero();
|
|
|
|
// map multipole moments to grid
|
|
|
|
grid_mpole(fmp,gridpre);
|
|
|
|
// pre-convolution operations including forward FFT
|
|
// gridfft = my 1d portion of complex 3d grid in FFT decomp
|
|
|
|
FFT_SCALAR *gridfft = i_kspace->pre_convolution();
|
|
|
|
// ---------------------
|
|
// convolution operation
|
|
// ---------------------
|
|
|
|
nhalf1 = (nfft1+1) / 2;
|
|
nhalf2 = (nfft2+1) / 2;
|
|
nhalf3 = (nfft3+1) / 2;
|
|
|
|
nxlo = i_kspace->nxlo_fft;
|
|
nxhi = i_kspace->nxhi_fft;
|
|
nylo = i_kspace->nylo_fft;
|
|
nyhi = i_kspace->nyhi_fft;
|
|
nzlo = i_kspace->nzlo_fft;
|
|
nzhi = i_kspace->nzhi_fft;
|
|
|
|
m = n = 0;
|
|
for (k = nzlo; k <= nzhi; k++) {
|
|
for (j = nylo; j <= nyhi; j++) {
|
|
for (i = nxlo; i <= nxhi; i++) {
|
|
r1 = (i >= nhalf1) ? i-nfft1 : i;
|
|
r2 = (j >= nhalf2) ? j-nfft2 : j;
|
|
r3 = (k >= nhalf3) ? k-nfft3 : k;
|
|
h1 = recip[0][0]*r1 + recip[0][1]*r2 + recip[0][2]*r3; // matvec
|
|
h2 = recip[1][0]*r1 + recip[1][1]*r2 + recip[1][2]*r3;
|
|
h3 = recip[2][0]*r1 + recip[2][1]*r2 + recip[2][2]*r3;
|
|
hsq = h1*h1 + h2*h2 + h3*h3;
|
|
term = -pterm * hsq;
|
|
expterm = 0.0;
|
|
if (term > -50.0 && hsq != 0.0) {
|
|
denom = volterm*hsq*bsmod1[i]*bsmod2[j]*bsmod3[k];
|
|
expterm = exp(term) / denom;
|
|
}
|
|
qfac[m++] = expterm;
|
|
gridfft[n] *= expterm;
|
|
gridfft[n+1] *= expterm;
|
|
n += 2;
|
|
}
|
|
}
|
|
}
|
|
|
|
// post-convolution operations including backward FFT
|
|
// gridppost = my portion of 3d grid in brick decomp w/ ghost values
|
|
|
|
FFT_SCALAR ***gridpost = (FFT_SCALAR ***) i_kspace->post_convolution();
|
|
|
|
// get potential
|
|
|
|
fphi_mpole(gridpost,fphi);
|
|
|
|
// convert the field from fractional to Cartesian
|
|
|
|
fphi_to_cphi(fphi,cphi);
|
|
|
|
// increment the field at each multipole site
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
field[i][0] -= cphi[i][1];
|
|
field[i][1] -= cphi[i][2];
|
|
field[i][2] -= cphi[i][3];
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
udirect2b = Ewald real direct field via list
|
|
udirect2b computes the real space contribution of the permanent
|
|
atomic multipole moments to the field via a neighbor list
|
|
------------------------------------------------------------------------- */
|
|
|
|
void PairAmoeba::udirect2b(double **field, double **fieldp)
|
|
{
|
|
int i,j,m,n,ii,jj,jextra,ndip,itype,jtype,iclass,jclass,igroup,jgroup;
|
|
double xr,yr,zr,r,r2;
|
|
double rr1,rr2,rr3;
|
|
double rr5,rr7;
|
|
double rr3i,rr5i,rr7i;
|
|
double rr3k,rr5k,rr7k;
|
|
double rr3ik,rr5ik;
|
|
double bfac,exp2a;
|
|
double ci,dix,diy,diz;
|
|
double qixx,qiyy,qizz;
|
|
double qixy,qixz,qiyz;
|
|
double ck,dkx,dky,dkz;
|
|
double qkxx,qkyy,qkzz;
|
|
double qkxy,qkxz,qkyz;
|
|
double dir,dkr;
|
|
double qix,qiy,qiz,qir;
|
|
double qkx,qky,qkz,qkr;
|
|
double corei,corek;
|
|
double vali,valk;
|
|
double alphai,alphak;
|
|
double ralpha,aefac;
|
|
double aesq2,aesq2n;
|
|
double pdi,pti,ddi;
|
|
double pgamma;
|
|
double damp,expdamp;
|
|
double scale3,scale5;
|
|
double scale7,scalek;
|
|
double bn[4],bcn[3];
|
|
double fid[3],fkd[3];
|
|
double fip[3],fkp[3];
|
|
double dmpi[7],dmpk[7];
|
|
double dmpik[5];
|
|
double factor_dscale,factor_pscale,factor_uscale,factor_wscale;
|
|
|
|
int inum,jnum;
|
|
int *ilist,*jlist,*numneigh,**firstneigh;
|
|
|
|
// owned atoms
|
|
|
|
double **x = atom->x;
|
|
double *pval = atom->dvector[index_pval];
|
|
|
|
// neigh list
|
|
|
|
inum = list->inum;
|
|
ilist = list->ilist;
|
|
numneigh = list->numneigh;
|
|
firstneigh = list->firstneigh;
|
|
|
|
aesq2 = 2.0 * aewald * aewald;
|
|
aesq2n = 0.0;
|
|
if (aewald > 0.0) aesq2n = 1.0 / (MY_PIS*aewald);
|
|
|
|
// rebuild dipole-dipole pair list and store pairwise dipole matrices
|
|
// done one atom at a time in real-space double loop over atoms & neighs
|
|
|
|
int *neighptr;
|
|
double *tdipdip;
|
|
|
|
// compute the real space portion of the Ewald summation
|
|
|
|
ipage_dipole->reset();
|
|
dpage_dipdip->reset();
|
|
|
|
for (ii = 0; ii < inum; ii++) {
|
|
i = ilist[ii];
|
|
itype = amtype[i];
|
|
iclass = amtype2class[itype];
|
|
igroup = amgroup[i];
|
|
jlist = firstneigh[i];
|
|
jnum = numneigh[i];
|
|
|
|
n = ndip = 0;
|
|
neighptr = ipage_dipole->vget();
|
|
tdipdip = dpage_dipdip->vget();
|
|
|
|
ci = rpole[i][0];
|
|
dix = rpole[i][1];
|
|
diy = rpole[i][2];
|
|
diz = rpole[i][3];
|
|
qixx = rpole[i][4];
|
|
qixy = rpole[i][5];
|
|
qixz = rpole[i][6];
|
|
qiyy = rpole[i][8];
|
|
qiyz = rpole[i][9];
|
|
qizz = rpole[i][12];
|
|
if (amoeba) {
|
|
pdi = pdamp[itype];
|
|
pti = thole[itype];
|
|
ddi = dirdamp[itype];
|
|
} else {
|
|
corei = pcore[iclass];
|
|
alphai = palpha[iclass];
|
|
vali = pval[i];
|
|
}
|
|
|
|
// evaluate all sites within the cutoff distance
|
|
|
|
for (jj = 0; jj < jnum; jj++) {
|
|
jextra = jlist[jj];
|
|
j = jextra & NEIGHMASK15;
|
|
|
|
xr = x[j][0] - x[i][0];
|
|
yr = x[j][1] - x[i][1];
|
|
zr = x[j][2] - x[i][2];
|
|
r2 = xr*xr + yr* yr + zr*zr;
|
|
if (r2 > off2) continue;
|
|
|
|
jtype = amtype[j];
|
|
jclass = amtype2class[jtype];
|
|
jgroup = amgroup[j];
|
|
|
|
if (amoeba) {
|
|
factor_wscale = special_polar_wscale[sbmask15(jextra)];
|
|
if (igroup == jgroup) {
|
|
factor_pscale = special_polar_piscale[sbmask15(jextra)];
|
|
factor_dscale = polar_dscale;
|
|
factor_uscale = polar_uscale;
|
|
} else {
|
|
factor_pscale = special_polar_pscale[sbmask15(jextra)];
|
|
factor_dscale = factor_uscale = 1.0;
|
|
}
|
|
|
|
} else {
|
|
factor_wscale = special_polar_wscale[sbmask15(jextra)];
|
|
if (igroup == jgroup) {
|
|
factor_dscale = factor_pscale = special_polar_piscale[sbmask15(jextra)];
|
|
factor_uscale = polar_uscale;
|
|
} else {
|
|
factor_dscale = factor_pscale = special_polar_pscale[sbmask15(jextra)];
|
|
factor_uscale = 1.0;
|
|
}
|
|
}
|
|
|
|
r = sqrt(r2);
|
|
rr1 = 1.0 / r;
|
|
rr2 = rr1 * rr1;
|
|
rr3 = rr2 * rr1;
|
|
rr5 = 3.0 * rr2 * rr3;
|
|
rr7 = 5.0 * rr2 * rr5;
|
|
ck = rpole[j][0];
|
|
dkx = rpole[j][1];
|
|
dky = rpole[j][2];
|
|
dkz = rpole[j][3];
|
|
qkxx = rpole[j][4];
|
|
qkxy = rpole[j][5];
|
|
qkxz = rpole[j][6];
|
|
qkyy = rpole[j][8];
|
|
qkyz = rpole[j][9];
|
|
qkzz = rpole[j][12];
|
|
|
|
// intermediates involving moments and separation distance
|
|
|
|
dir = dix*xr + diy*yr + diz*zr;
|
|
qix = qixx*xr + qixy*yr + qixz*zr;
|
|
qiy = qixy*xr + qiyy*yr + qiyz*zr;
|
|
qiz = qixz*xr + qiyz*yr + qizz*zr;
|
|
qir = qix*xr + qiy*yr + qiz*zr;
|
|
dkr = dkx*xr + dky*yr + dkz*zr;
|
|
qkx = qkxx*xr + qkxy*yr + qkxz*zr;
|
|
qky = qkxy*xr + qkyy*yr + qkyz*zr;
|
|
qkz = qkxz*xr + qkyz*yr + qkzz*zr;
|
|
qkr = qkx*xr + qky*yr + qkz*zr;
|
|
|
|
// calculate the real space Ewald error function terms
|
|
|
|
ralpha = aewald * r;
|
|
bn[0] = erfc(ralpha) * rr1;
|
|
exp2a = exp(-ralpha*ralpha);
|
|
aefac = aesq2n;
|
|
for (m = 1; m <= 3; m++) {
|
|
bfac = m+m-1;
|
|
aefac = aesq2 * aefac;
|
|
bn[m] = (bfac*bn[m-1]+aefac*exp2a) * rr2;
|
|
}
|
|
|
|
// find the field components for Thole polarization damping
|
|
|
|
if (amoeba) {
|
|
scale3 = 1.0;
|
|
scale5 = 1.0;
|
|
scale7 = 1.0;
|
|
damp = pdi * pdamp[jtype];
|
|
if (damp != 0.0) {
|
|
pgamma = MIN(ddi,dirdamp[jtype]);
|
|
if (pgamma != 0.0) {
|
|
damp = pgamma * pow(r/damp,1.5);
|
|
if (damp < 50.0) {
|
|
expdamp = exp(-damp) ;
|
|
scale3 = 1.0 - expdamp ;
|
|
scale5 = 1.0 - expdamp*(1.0+0.5*damp);
|
|
scale7 = 1.0 - expdamp*(1.0+0.65*damp + 0.15*damp*damp);
|
|
}
|
|
} else {
|
|
pgamma = MIN(pti,thole[jtype]);
|
|
damp = pgamma * cube(r/damp);
|
|
if (damp < 50.0) {
|
|
expdamp = exp(-damp);
|
|
scale3 = 1.0 - expdamp;
|
|
scale5 = 1.0 - expdamp*(1.0+damp);
|
|
scale7 = 1.0 - expdamp*(1.0+damp + 0.6*damp*damp);
|
|
}
|
|
}
|
|
}
|
|
|
|
scalek = factor_dscale;
|
|
bcn[0] = bn[1] - (1.0-scalek*scale3)*rr3;
|
|
bcn[1] = bn[2] - (1.0-scalek*scale5)*rr5;
|
|
bcn[2] = bn[3] - (1.0-scalek*scale7)*rr7;
|
|
fid[0] = -xr*(bcn[0]*ck-bcn[1]*dkr+bcn[2]*qkr) -
|
|
bcn[0]*dkx + 2.0*bcn[1]*qkx;
|
|
fid[1] = -yr*(bcn[0]*ck-bcn[1]*dkr+bcn[2]*qkr) -
|
|
bcn[0]*dky + 2.0*bcn[1]*qky;
|
|
fid[2] = -zr*(bcn[0]*ck-bcn[1]*dkr+bcn[2]*qkr) -
|
|
bcn[0]*dkz + 2.0*bcn[1]*qkz;
|
|
fkd[0] = xr*(bcn[0]*ci+bcn[1]*dir+bcn[2]*qir) -
|
|
bcn[0]*dix - 2.0*bcn[1]*qix;
|
|
fkd[1] = yr*(bcn[0]*ci+bcn[1]*dir+bcn[2]*qir) -
|
|
bcn[0]*diy - 2.0*bcn[1]*qiy;
|
|
fkd[2] = zr*(bcn[0]*ci+bcn[1]*dir+bcn[2]*qir) -
|
|
bcn[0]*diz - 2.0*bcn[1]*qiz;
|
|
|
|
scalek = factor_pscale;
|
|
bcn[0] = bn[1] - (1.0-scalek*scale3)*rr3;
|
|
bcn[1] = bn[2] - (1.0-scalek*scale5)*rr5;
|
|
bcn[2] = bn[3] - (1.0-scalek*scale7)*rr7;
|
|
fip[0] = -xr*(bcn[0]*ck-bcn[1]*dkr+bcn[2]*qkr) -
|
|
bcn[0]*dkx + 2.0*bcn[1]*qkx;
|
|
fip[1] = -yr*(bcn[0]*ck-bcn[1]*dkr+bcn[2]*qkr) -
|
|
bcn[0]*dky + 2.0*bcn[1]*qky;
|
|
fip[2] = -zr*(bcn[0]*ck-bcn[1]*dkr+bcn[2]*qkr) -
|
|
bcn[0]*dkz + 2.0*bcn[1]*qkz;
|
|
fkp[0] = xr*(bcn[0]*ci+bcn[1]*dir+bcn[2]*qir) -
|
|
bcn[0]*dix - 2.0*bcn[1]*qix;
|
|
fkp[1] = yr*(bcn[0]*ci+bcn[1]*dir+bcn[2]*qir) -
|
|
bcn[0]*diy - 2.0*bcn[1]*qiy;
|
|
fkp[2] = zr*(bcn[0]*ci+bcn[1]*dir+bcn[2]*qir) -
|
|
bcn[0]*diz - 2.0*bcn[1]*qiz;
|
|
|
|
// find terms needed later to compute mutual polarization
|
|
|
|
if (poltyp != DIRECT) {
|
|
scale3 = 1.0;
|
|
scale5 = 1.0;
|
|
damp = pdi * pdamp[jtype];
|
|
if (damp != 0.0) {
|
|
pgamma = MIN(pti,thole[jtype]);
|
|
damp = pgamma * cube(r/damp);
|
|
if (damp < 50.0) {
|
|
expdamp = exp(-damp);
|
|
scale3 = 1.0 - expdamp;
|
|
scale5 = 1.0 - expdamp*(1.0+damp);
|
|
}
|
|
}
|
|
scalek = factor_uscale;
|
|
bcn[0] = bn[1] - (1.0-scalek*scale3)*rr3;
|
|
bcn[1] = bn[2] - (1.0-scalek*scale5)*rr5;
|
|
|
|
neighptr[n++] = j;
|
|
tdipdip[ndip++] = -bcn[0] + bcn[1]*xr*xr;
|
|
tdipdip[ndip++] = bcn[1]*xr*yr;
|
|
tdipdip[ndip++] = bcn[1]*xr*zr;
|
|
tdipdip[ndip++] = -bcn[0] + bcn[1]*yr*yr;
|
|
tdipdip[ndip++] = bcn[1]*yr*zr;
|
|
tdipdip[ndip++] = -bcn[0] + bcn[1]*zr*zr;
|
|
}
|
|
|
|
// find the field components for charge penetration damping
|
|
|
|
} else {
|
|
corek = pcore[jclass];
|
|
alphak = palpha[jclass];
|
|
valk = pval[j];
|
|
dampdir(r,alphai,alphak,dmpi,dmpk);
|
|
|
|
scalek = factor_dscale;
|
|
rr3i = bn[1] - (1.0-scalek*dmpi[2])*rr3;
|
|
rr5i = bn[2] - (1.0-scalek*dmpi[4])*rr5;
|
|
rr7i = bn[3] - (1.0-scalek*dmpi[6])*rr7;
|
|
rr3k = bn[1] - (1.0-scalek*dmpk[2])*rr3;
|
|
rr5k = bn[2] - (1.0-scalek*dmpk[4])*rr5;
|
|
rr7k = bn[3] - (1.0-scalek*dmpk[6])*rr7;
|
|
rr3 = bn[1] - (1.0-scalek)*rr3;
|
|
fid[0] = -xr*(rr3*corek + rr3k*valk - rr5k*dkr + rr7k*qkr) -
|
|
rr3k*dkx + 2.0*rr5k*qkx;
|
|
fid[1] = -yr*(rr3*corek + rr3k*valk - rr5k*dkr + rr7k*qkr) -
|
|
rr3k*dky + 2.0*rr5k*qky;
|
|
fid[2] = -zr*(rr3*corek + rr3k*valk - rr5k*dkr + rr7k*qkr) -
|
|
rr3k*dkz + 2.0*rr5k*qkz;
|
|
fkd[0] = xr*(rr3*corei + rr3i*vali + rr5i*dir + rr7i*qir) -
|
|
rr3i*dix - 2.0*rr5i*qix;
|
|
fkd[1] = yr*(rr3*corei + rr3i*vali + rr5i*dir + rr7i*qir) -
|
|
rr3i*diy - 2.0*rr5i*qiy;
|
|
fkd[2] = zr*(rr3*corei + rr3i*vali + rr5i*dir + rr7i*qir) -
|
|
rr3i*diz - 2.0*rr5i*qiz;
|
|
|
|
scalek = factor_pscale;
|
|
rr3 = rr2 * rr1;
|
|
rr3i = bn[1] - (1.0-scalek*dmpi[2])*rr3;
|
|
rr5i = bn[2] - (1.0-scalek*dmpi[4])*rr5;
|
|
rr7i = bn[3] - (1.0-scalek*dmpi[6])*rr7;
|
|
rr3k = bn[1] - (1.0-scalek*dmpk[2])*rr3;
|
|
rr5k = bn[2] - (1.0-scalek*dmpk[4])*rr5;
|
|
rr7k = bn[3] - (1.0-scalek*dmpk[6])*rr7;
|
|
rr3 = bn[1] - (1.0-scalek)*rr3;
|
|
fip[0] = -xr*(rr3*corek + rr3k*valk - rr5k*dkr + rr7k*qkr) -
|
|
rr3k*dkx + 2.0*rr5k*qkx;
|
|
fip[1] = -yr*(rr3*corek + rr3k*valk - rr5k*dkr + rr7k*qkr) -
|
|
rr3k*dky + 2.0*rr5k*qky;
|
|
fip[2] = -zr*(rr3*corek + rr3k*valk - rr5k*dkr + rr7k*qkr) -
|
|
rr3k*dkz + 2.0*rr5k*qkz;
|
|
fkp[0] = xr*(rr3*corei + rr3i*vali + rr5i*dir + rr7i*qir) -
|
|
rr3i*dix - 2.0*rr5i*qix;
|
|
fkp[1] = yr*(rr3*corei + rr3i*vali + rr5i*dir + rr7i*qir) -
|
|
rr3i*diy - 2.0*rr5i*qiy;
|
|
fkp[2] = zr*(rr3*corei + rr3i*vali + rr5i*dir + rr7i*qir) -
|
|
rr3i*diz - 2.0*rr5i*qiz;
|
|
|
|
// find terms needed later to compute mutual polarization
|
|
|
|
if (poltyp != DIRECT) {
|
|
dampmut(r,alphai,alphak,dmpik);
|
|
scalek = factor_wscale;
|
|
rr3 = rr2 * rr1;
|
|
rr3ik = bn[1] - (1.0-scalek*dmpik[2])*rr3;
|
|
rr5ik = bn[2] - (1.0-scalek*dmpik[4])*rr5;
|
|
|
|
neighptr[n++] = j;
|
|
tdipdip[ndip++] = -rr3ik + rr5ik*xr*xr;
|
|
tdipdip[ndip++] = rr5ik*xr*yr;
|
|
tdipdip[ndip++] = rr5ik*xr*zr;
|
|
tdipdip[ndip++] = -rr3ik + rr5ik*yr*yr;
|
|
tdipdip[ndip++] = rr5ik*yr*zr;
|
|
tdipdip[ndip++] = -rr3ik + rr5ik*zr*zr;
|
|
}
|
|
}
|
|
|
|
// increment the field at each site due to this interaction
|
|
|
|
for (m = 0; m < 3; m++) {
|
|
field[i][m] += fid[m];
|
|
field[j][m] += fkd[m];
|
|
fieldp[i][m] += fip[m];
|
|
fieldp[j][m] += fkp[m];
|
|
}
|
|
}
|
|
|
|
firstneigh_dipole[i] = neighptr;
|
|
firstneigh_dipdip[i] = tdipdip;
|
|
numneigh_dipole[i] = n;
|
|
ipage_dipole->vgot(n);
|
|
dpage_dipdip->vgot(ndip);
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
dampmut = mutual field damping coefficents
|
|
dampmut generates coefficients for the mutual field damping
|
|
function for powers of the interatomic distance
|
|
------------------------------------------------------------------------- */
|
|
|
|
void PairAmoeba::dampmut(double r, double alphai, double alphak, double *dmpik)
|
|
{
|
|
double termi,termk;
|
|
double termi2,termk2;
|
|
double alphai2,alphak2;
|
|
double eps,diff;
|
|
double expi,expk;
|
|
double dampi,dampk;
|
|
double dampi2,dampi3;
|
|
double dampi4,dampi5;
|
|
double dampk2,dampk3;
|
|
|
|
// compute tolerance and exponential damping factors
|
|
|
|
eps = 0.001;
|
|
diff = fabs(alphai-alphak);
|
|
dampi = alphai * r;
|
|
dampk = alphak * r;
|
|
expi = exp(-dampi);
|
|
expk = exp(-dampk);
|
|
|
|
// valence-valence charge penetration damping for Gordon f1 (HIPPO)
|
|
|
|
dampi2 = dampi * dampi;
|
|
dampi3 = dampi * dampi2;
|
|
if (diff < eps) {
|
|
dampi4 = dampi2 * dampi2;
|
|
dampi5 = dampi2 * dampi3;
|
|
dmpik[2] = 1.0 - (1.0 + dampi + 0.5*dampi2 +
|
|
7.0*dampi3/48.0 + dampi4/48.0)*expi;
|
|
dmpik[4] = 1.0 - (1.0 + dampi + 0.5*dampi2 + dampi3/6.0 +
|
|
dampi4/24.0 + dampi5/144.0)*expi;
|
|
} else {
|
|
dampk2 = dampk * dampk;
|
|
dampk3 = dampk * dampk2;
|
|
alphai2 = alphai * alphai;
|
|
alphak2 = alphak * alphak;
|
|
termi = alphak2 / (alphak2-alphai2);
|
|
termk = alphai2 / (alphai2-alphak2);
|
|
termi2 = termi * termi;
|
|
termk2 = termk * termk;
|
|
dmpik[2] = 1.0 - termi2*(1.0+dampi+0.5*dampi2)*expi -
|
|
termk2*(1.0+dampk+0.5*dampk2)*expk -
|
|
2.0*termi2*termk*(1.0+dampi)*expi - 2.0*termk2*termi*(1.0+dampk)*expk;
|
|
dmpik[4] = 1.0 - termi2*(1.0+dampi+0.5*dampi2 + dampi3/6.0)*expi -
|
|
termk2*(1.0+dampk+0.5*dampk2 + dampk3/6.00)*expk -
|
|
2.0*termi2*termk *(1.0+dampi+dampi2/3.0)*expi -
|
|
2.0*termk2*termi *(1.0+dampk+dampk2/3.0)*expk;
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
dampdir = direct field damping coefficents
|
|
dampdir generates coefficients for the direct field damping
|
|
function for powers of the interatomic distance
|
|
------------------------------------------------------------------------- */
|
|
|
|
void PairAmoeba::dampdir(double r, double alphai, double alphak,
|
|
double *dmpi, double *dmpk)
|
|
{
|
|
double eps,diff;
|
|
double expi,expk;
|
|
double dampi,dampk;
|
|
double dampi2,dampk2;
|
|
double dampi3,dampk3;
|
|
double dampi4,dampk4;
|
|
|
|
// compute tolerance and exponential damping factors
|
|
|
|
eps = 0.001;
|
|
diff = fabs(alphai-alphak);
|
|
dampi = alphai * r;
|
|
dampk = alphak * r;
|
|
expi = exp(-dampi);
|
|
expk = exp(-dampk);
|
|
|
|
// core-valence charge penetration damping for Gordon f1 (HIPPO)
|
|
|
|
dampi2 = dampi * dampi;
|
|
dampi3 = dampi * dampi2;
|
|
dampi4 = dampi2 * dampi2;
|
|
dmpi[2] = 1.0 - (1.0 + dampi + 0.5*dampi2)*expi;
|
|
dmpi[4] = 1.0 - (1.0 + dampi + 0.5*dampi2 + dampi3/6.0)*expi;
|
|
dmpi[6] = 1.0 - (1.0 + dampi + 0.5*dampi2 + dampi3/6.0 + dampi4/30.0)*expi;
|
|
if (diff < eps) {
|
|
dmpk[2] = dmpi[2];
|
|
dmpk[4] = dmpi[4];
|
|
dmpk[6] = dmpi[6];
|
|
} else {
|
|
dampk2 = dampk * dampk;
|
|
dampk3 = dampk * dampk2;
|
|
dampk4 = dampk2 * dampk2;
|
|
dmpk[2] = 1.0 - (1.0 + dampk + 0.5*dampk2)*expk;
|
|
dmpk[4] = 1.0 - (1.0 + dampk + 0.5*dampk2 + dampk3/6.0)*expk;
|
|
dmpk[6] = 1.0 - (1.0 + dampk + 0.5*dampk2 + dampk3/6.0 + dampk4/30.0)*expk;
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
cholesky = modified Cholesky linear solver
|
|
cholesky uses a modified Cholesky method to solve the linear
|
|
system Ax = b, returning "x" in "b"; "A" is a real symmetric
|
|
positive definite matrix with its upper triangle (including the
|
|
diagonal) stored by rows
|
|
literature reference:
|
|
R. S. Martin, G. Peters and J. H. Wilkinson, "Symmetric
|
|
Decomposition of a Positive Definite Matrix", Numerische
|
|
Mathematik, 7, 362-383 (1965)
|
|
------------------------------------------------------------------------- */
|
|
|
|
void PairAmoeba::cholesky(int nvar, double *a, double *b)
|
|
{
|
|
int i,j,k;
|
|
int ii,ij,ik,ki,kk;
|
|
int im,jk,jm;
|
|
double r,s,t;
|
|
|
|
// all code in this method is exact translation from Fortran version
|
|
// decrement pointers so can access vectors like Fortran does from 1:N
|
|
|
|
a--;
|
|
b--;
|
|
|
|
// Cholesky factorization to reduce A to (L)(D)(L transpose)
|
|
// L has a unit diagonal; store 1.0/D on the diagonal of A
|
|
|
|
ii = 1;
|
|
for (i = 1; i <= nvar; i++) {
|
|
im = i - 1;
|
|
if (i != 1) {
|
|
ij = i;
|
|
for (j = 1; j <= im; j++) {
|
|
r = a[ij];
|
|
if (j != 1) {
|
|
ik = i;
|
|
jk = j;
|
|
jm = j - 1;
|
|
for (k = 1; k <= jm; k++) {
|
|
r = r - a[ik]*a[jk];
|
|
ik = nvar - k + ik;
|
|
jk = nvar - k + jk;
|
|
}
|
|
}
|
|
a[ij] = r;
|
|
ij = nvar - j + ij;
|
|
}
|
|
}
|
|
|
|
r = a[ii];
|
|
if (i != 1) {
|
|
kk = 1;
|
|
ik = i;
|
|
for (k = 1; k <= im; k++) {
|
|
s = a[ik];
|
|
t = s * a[kk];
|
|
a[ik] = t;
|
|
r = r - s*t;
|
|
ik = nvar - k + ik;
|
|
kk = nvar - k + 1 + kk;
|
|
}
|
|
}
|
|
|
|
a[ii] = 1.0 / r;
|
|
ii = nvar - i + 1 + ii;
|
|
}
|
|
|
|
// solve linear equations; first solve Ly = b for y
|
|
|
|
for (i = 1; i <= nvar; i++) {
|
|
if (i != 1) {
|
|
ik = i;
|
|
im = i - 1;
|
|
r = b[i];
|
|
for (k = 1; k <= im; k++) {
|
|
r = r - b[k]*a[ik];
|
|
ik = nvar - k + ik;
|
|
}
|
|
b[i] = r;
|
|
}
|
|
}
|
|
|
|
// finally, solve (D)(L transpose)(x) = y for x
|
|
|
|
ii = nvar*(nvar+1)/2;
|
|
for (j = 1; j <= nvar; j++) {
|
|
i = nvar + 1 - j;
|
|
r = b[i] * a[ii];
|
|
if (j != 1) {
|
|
im = i + 1;
|
|
ki = ii + 1;
|
|
for (k = im; k <= nvar; k++) {
|
|
r = r - a[ki]*b[k];
|
|
ki = ki + 1;
|
|
}
|
|
}
|
|
b[i] = r;
|
|
ii = ii - j - 1;
|
|
}
|
|
}
|