ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

remove some more VLAs

Browse Source
This commit is contained in:
Axel Kohlmeyer
2024-09-10 22:12:22 -04:00
parent ac5d2d560d
commit d76e10d2ca
2 changed files with 93 additions and 48 deletions
Download Patch File Download Diff File Expand all files Collapse all files

138
src/QEQ/fix_qeq_ctip.cpp
View File

@ -43,13 +43,13 @@ using namespace FixConst;
FixQEqCTIP::FixQEqCTIP(LAMMPS *lmp, int narg, char **arg) :
FixQEq(lmp, narg, arg), reff(nullptr), reffsq(nullptr), reff4(nullptr), reff7(nullptr),
s2d_self(nullptr)
s2d_self(nullptr), shield(nullptr), shieldcu(nullptr), reffc(nullptr), reffcsq(nullptr),
reffc4(nullptr), reffc7(nullptr), s2d_shift(nullptr), f_shift(nullptr), e_shift(nullptr)
{
cdamp = 0.30;
maxrepeat = 10;
// optional arg
// optional args
int iarg = 8;
while (iarg < narg) {
if (strcmp(arg[iarg], "cdamp") == 0) {
@ -78,6 +78,15 @@ FixQEqCTIP::~FixQEqCTIP()
delete[] reff4;
delete[] reff7;
delete[] s2d_self;
memory->destroy(shield);
memory->destroy(shieldcu);
memory->destroy(reffc);
memory->destroy(reffcsq);
memory->destroy(reffc4);
memory->destroy(reffc7);
memory->destroy(s2d_shift);
memory->destroy(f_shift);
memory->destroy(e_shift);
}
/* ---------------------------------------------------------------------- */
@ -89,7 +98,7 @@ void FixQEqCTIP::init()
neighbor->add_request(this, NeighConst::REQ_FULL);
int ntypes = atom->ntypes;
memory->create(shld,ntypes+1,ntypes+1,"qeq:shielding");
memory->create(shld, ntypes + 1, ntypes + 1, "qeq:shielding");
delete[] reff;
delete[] reffsq;
@ -104,19 +113,65 @@ void FixQEqCTIP::init()
s2d_self = new double[ntypes];
double r = cutoff;
double rsq = r*r;
double r6 = rsq*rsq*rsq;
double rsq = r * r;
double r6 = rsq * rsq * rsq;
double erfcd_cut = exp(-cdamp * cdamp * rsq);
double t_cut = 1.0 / (1.0 + EWALD_P * cdamp * r);
double erfcc_cut = (t_cut * (A1 + t_cut * (A2 + t_cut * (A3 + t_cut * (A4 + t_cut * A5)))) * erfcd_cut) / r;
double erfcc_cut =
(t_cut * (A1 + t_cut * (A2 + t_cut * (A3 + t_cut * (A4 + t_cut * A5)))) * erfcd_cut) / r;
for (int elt1 = 0; elt1 < ntypes; elt1++) {
reff[elt1] = std::cbrt(rsq*r + 1.0/(gamma[elt1+1]*gamma[elt1+1]*gamma[elt1+1]));
reffsq[elt1] = reff[elt1]*reff[elt1];
reff4[elt1] = reffsq[elt1]*reffsq[elt1];
reff7[elt1] = reff4[elt1]*reffsq[elt1]*reff[elt1];
s2d_self[elt1] = 2.0*force->qqr2e*(1.5*erfcc_cut + 2.0*cdamp/MY_PIS*erfcd_cut + cdamp*cdamp*cdamp/MY_PIS*rsq*erfcd_cut + 0.5/reff[elt1] - 1.5/r + r6/reff7[elt1] + cdamp/MY_PIS);
reff[elt1] = std::cbrt(rsq * r + 1.0 / (gamma[elt1 + 1] * gamma[elt1 + 1] * gamma[elt1 + 1]));
reffsq[elt1] = reff[elt1] * reff[elt1];
reff4[elt1] = reffsq[elt1] * reffsq[elt1];
reff7[elt1] = reff4[elt1] * reffsq[elt1] * reff[elt1];
s2d_self[elt1] = 2.0 * force->qqr2e *
(1.5 * erfcc_cut + 2.0 * cdamp / MY_PIS * erfcd_cut +
cdamp * cdamp * cdamp / MY_PIS * rsq * erfcd_cut + 0.5 / reff[elt1] - 1.5 / r +
r6 / reff7[elt1] + cdamp / MY_PIS);
}
memory->destroy(shield);
memory->destroy(shieldcu);
memory->destroy(reffc);
memory->destroy(reffcsq);
memory->destroy(reffc4);
memory->destroy(reffc7);
memory->destroy(s2d_shift);
memory->destroy(f_shift);
memory->destroy(e_shift);
memory->create(shield, ntypes, ntypes, "qeq:shield");
memory->create(shieldcu, ntypes, ntypes, "qeq:shieldcu");
memory->create(reffc, ntypes, ntypes, "qeq:reffc");
memory->create(reffcsq, ntypes, ntypes, "qeq:reffcsq");
memory->create(reffc4, ntypes, ntypes, "qeq:reffc4");
memory->create(reffc7, ntypes, ntypes, "qeq:reffc7");
memory->create(s2d_shift, ntypes, ntypes, "qeq:s2d_shift");
memory->create(f_shift, ntypes, ntypes, "qeq:f_shift");
memory->create(e_shift, ntypes, ntypes, "qeq:e_shift");
double cutoffsq = cutoff * cutoff;
double cutoffcu = cutoffsq * cutoff;
double cutoff4 = cutoffsq * cutoffsq;
double cdampcu = cdamp * cdamp * cdamp;
for (int elt1 = 0; elt1 < atom->ntypes; elt1++) {
for (int elt2 = 0; elt2 < atom->ntypes; elt2++) {
shield[elt1][elt2] = sqrt(gamma[elt1 + 1] * gamma[elt2 + 1]);
shieldcu[elt1][elt2] = shield[elt1][elt2] * shield[elt1][elt2] * shield[elt1][elt2];
reffc[elt1][elt2] = cbrt(cutoffcu + 1 / shieldcu[elt1][elt2]);
reffcsq[elt1][elt2] = reffc[elt1][elt2] * reffc[elt1][elt2];
reffc4[elt1][elt2] = reffcsq[elt1][elt2] * reffcsq[elt1][elt2];
reffc7[elt1][elt2] = reffc4[elt1][elt2] * reffcsq[elt1][elt2] * reffc[elt1][elt2];
s2d_shift[elt1][elt2] = 2.0 * erfcc_cut / cutoffcu +
4.0 * cdamp / MY_PIS * erfcd_cut / cutoffsq + 4.0 * cdampcu / MY_PIS * erfcd_cut -
2 / cutoffcu + 4 * cutoff4 / reffc7[elt1][elt2] - 2 * cutoff / reffc4[elt1][elt2];
f_shift[elt1][elt2] = erfcc_cut / cutoffsq + 2.0 * cdamp / MY_PIS * erfcd_cut / cutoff -
1 / cutoffsq + cutoffsq / reffc4[elt1][elt2];
e_shift[elt1][elt2] = erfcc_cut / cutoff + 1 / reffc[elt1][elt2] - 1 / cutoff;
}
}
}
@ -182,6 +237,7 @@ void FixQEqCTIP::init_matvec()
int inum, ii, i;
int *ilist;
double *q = atom->q, qi;
int *type = atom->type;
double r = cutoff;
double rsq = r*r;
@ -199,15 +255,15 @@ void FixQEqCTIP::init_matvec()
if (atom->mask[i] & groupbit) {
qi=q[i];
if (qi < qmin[atom->type[i]]) {
Hdia_inv[i] = 1. / (eta[atom->type[i]]+2*omega[atom->type[i]]-s2d_self[atom->type[i]-1]);
b_s[i] = -((chi[atom->type[i]]-2*qmin[atom->type[i]]*omega[atom->type[i]]) + chizj[i]);
} else if (qi < qmax[atom->type[i]]) {
Hdia_inv[i] = 1. / (eta[atom->type[i]]-s2d_self[atom->type[i]-1]);
b_s[i] = -(chi[atom->type[i]] + chizj[i]);
if (qi < qmin[type[i]]) {
Hdia_inv[i] = 1. / (eta[type[i]]+2*omega[type[i]]-s2d_self[type[i]-1]);
b_s[i] = -((chi[type[i]]-2*qmin[type[i]]*omega[type[i]]) + chizj[i]);
} else if (qi < qmax[type[i]]) {
Hdia_inv[i] = 1. / (eta[type[i]]-s2d_self[type[i]-1]);
b_s[i] = -(chi[type[i]] + chizj[i]);
} else {
Hdia_inv[i] = 1. / (eta[atom->type[i]]+2*omega[atom->type[i]]-s2d_self[atom->type[i]-1]);
b_s[i] = -((chi[atom->type[i]]-2*qmax[atom->type[i]]*omega[atom->type[i]]) + chizj[i]);
Hdia_inv[i] = 1. / (eta[type[i]]+2*omega[type[i]]-s2d_self[type[i]-1]);
b_s[i] = -((chi[type[i]]-2*qmax[type[i]]*omega[type[i]]) + chizj[i]);
}
b_t[i] = -1.0;
@ -228,16 +284,14 @@ void FixQEqCTIP::compute_H()
{
int inum, jnum, *ilist, *jlist, *numneigh, **firstneigh;
int i, j, ii, jj;
double **x;
double dx, dy, dz, r_sqr, r, reff;
double cutoffsq, cutoffcu, cutoff4, cdampcu, erfcc_cut, erfcd_cut, t_cut;
double erfcc, erfcd, t;
int elt1, elt2;
double shield[atom->ntypes][atom->ntypes], shieldcu[atom->ntypes][atom->ntypes], reffc[atom->ntypes][atom->ntypes], reffcsq[atom->ntypes][atom->ntypes], reffc4[atom->ntypes][atom->ntypes], reffc7[atom->ntypes][atom->ntypes];
double s2d_shift[atom->ntypes][atom->ntypes], f_shift[atom->ntypes][atom->ntypes], e_shift[atom->ntypes][atom->ntypes];
x = atom->x;
double **x = atom->x;
int *mask = atom->mask;
int *type = atom->type;
inum = list->inum;
ilist = list->ilist;
@ -253,19 +307,6 @@ void FixQEqCTIP::compute_H()
t_cut = 1.0 / (1.0 + EWALD_P * cdamp * cutoff);
erfcc_cut = t_cut * (A1 + t_cut * (A2 + t_cut * (A3 + t_cut * (A4 + t_cut * A5)))) * erfcd_cut;
for (elt1 = 0; elt1 < atom->ntypes; elt1++) {
for (elt2 = 0; elt2 < atom->ntypes; elt2++) {
shield[elt1][elt2] = sqrt(gamma[elt1+1] * gamma[elt2+1]);
shieldcu[elt1][elt2] = shield[elt1][elt2] * shield[elt1][elt2] * shield[elt1][elt2];
reffc[elt1][elt2] = cbrt(cutoffcu + 1 / shieldcu[elt1][elt2]);
reffcsq[elt1][elt2] = reffc[elt1][elt2] * reffc[elt1][elt2];
reffc4[elt1][elt2] = reffcsq[elt1][elt2] * reffcsq[elt1][elt2];
reffc7[elt1][elt2] = reffc4[elt1][elt2] * reffcsq[elt1][elt2] * reffc[elt1][elt2];
s2d_shift[elt1][elt2] = 2.0 * erfcc_cut / cutoffcu + 4.0 * cdamp / MY_PIS * erfcd_cut / cutoffsq + 4.0 * cdampcu / MY_PIS * erfcd_cut - 2 / cutoffcu + 4 * cutoff4 / reffc7[elt1][elt2] - 2 * cutoff / reffc4[elt1][elt2];
f_shift[elt1][elt2] = erfcc_cut / cutoffsq + 2.0 * cdamp / MY_PIS * erfcd_cut / cutoff - 1 / cutoffsq + cutoffsq / reffc4[elt1][elt2];
e_shift[elt1][elt2] = erfcc_cut / cutoff + 1 / reffc[elt1][elt2] - 1 / cutoff;
}
}
// fill in the H matrix
m_fill = 0;
@ -288,11 +329,11 @@ void FixQEqCTIP::compute_H()
if (r_sqr <= cutoff_sq) {
H.jlist[m_fill] = j;
r = sqrt(r_sqr);
reff = cbrt(r_sqr * r + 1 / shieldcu[atom->type[i]-1][atom->type[j]-1]);
reff = cbrt(r_sqr * r + 1 / shieldcu[type[i]-1][type[j]-1]);
erfcd = exp(-cdamp * cdamp * r_sqr);
t = 1.0 / (1.0 + EWALD_P * cdamp * r);
erfcc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * erfcd;
H.val[m_fill] = 0.5*force->qqr2e*(erfcc/r+1/reff-1/r-e_shift[atom->type[i]-1][atom->type[j]-1]+f_shift[atom->type[i]-1][atom->type[j]-1]*(r-cutoff)-s2d_shift[atom->type[i]-1][atom->type[j]-1]*0.5*(r-cutoff)*(r-cutoff));
H.val[m_fill] = 0.5*force->qqr2e*(erfcc/r+1/reff-1/r-e_shift[type[i]-1][type[j]-1]+f_shift[type[i]-1][type[j]-1]*(r-cutoff)-s2d_shift[type[i]-1][type[j]-1]*0.5*(r-cutoff)*(r-cutoff));
m_fill++;
}
}
@ -301,8 +342,7 @@ void FixQEqCTIP::compute_H()
}
if (m_fill >= H.m)
error->all(FLERR,"Fix qeq/ctip has insufficient H matrix "
"size: m_fill={} H.m={}\n",m_fill, H.m);
error->all(FLERR,"Fix qeq/ctip has insufficient H matrix size: m_fill={} H.m={}\n",m_fill, H.m);
}
/* ---------------------------------------------------------------------- */
@ -311,6 +351,7 @@ void FixQEqCTIP::sparse_matvec(sparse_matrix *A, double *x, double *b)
{
int i, j, itr_j;
double *q=atom->q, qi;
int *type = atom->type;
nlocal = atom->nlocal;
nall = atom->nlocal + atom->nghost;
@ -318,12 +359,12 @@ void FixQEqCTIP::sparse_matvec(sparse_matrix *A, double *x, double *b)
for (i = 0; i < nlocal; ++i) {
if (atom->mask[i] & groupbit) {
qi=q[i];
if (qi < qmin[atom->type[i]]) {
b[i] = (eta[atom->type[i]]+2*omega[atom->type[i]]-s2d_self[atom->type[i]-1])*x[i];
} else if (qi < qmax[atom->type[i]]) {
b[i] = (eta[atom->type[i]]-s2d_self[atom->type[i]-1]) * x[i];
if (qi < qmin[type[i]]) {
b[i] = (eta[type[i]]+2*omega[type[i]]-s2d_self[type[i]-1])*x[i];
} else if (qi < qmax[type[i]]) {
b[i] = (eta[type[i]]-s2d_self[type[i]-1]) * x[i];
} else {
b[i] = (eta[atom->type[i]]+2*omega[atom->type[i]]-s2d_self[atom->type[i]-1])*x[i];
b[i] = (eta[type[i]]+2*omega[type[i]]-s2d_self[type[i]-1])*x[i];
}
}
}
@ -353,6 +394,7 @@ int FixQEqCTIP::calculate_check_Q()
int *ilist;
double u, s_sum, t_sum;
double *q = atom->q;
int *type = atom->type;
double qi_old,qi_new;
double qi_check1,qi_check2;
double qi_check3;
@ -380,8 +422,8 @@ int FixQEqCTIP::calculate_check_Q()
t_hist[i][0] = t[i];
qi_new = q[i];
qi_check1=(qi_new-qmin[atom->type[i]])*(qi_old-qmin[atom->type[i]]);
qi_check2=(qi_new-qmax[atom->type[i]])*(qi_old-qmax[atom->type[i]]);
qi_check1=(qi_new-qmin[type[i]])*(qi_old-qmin[type[i]]);
qi_check2=(qi_new-qmax[type[i]])*(qi_old-qmax[type[i]]);
if ( qi_check1 < 0.0 || qi_check2 < 0.0 ) {
qi_check3=abs(qi_new-qi_old);
if (qi_check3 > tolerance) n++;

3
src/QEQ/fix_qeq_ctip.h
View File

@ -39,6 +39,9 @@ class FixQEqCTIP : public FixQEq {
void extract_ctip();
int calculate_check_Q();
double *reff, *reffsq, *reff4, *reff7, *s2d_self;
double **shield, **shieldcu, **reffc, **reffcsq, **reffc4, **reffc7;
double **s2d_shift, **f_shift, **e_shift;
double cdamp;
int maxrepeat;
int nout;