From d76e10d2ca613c76dcab46889fba2a79e3315c3e Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 10 Sep 2024 22:12:22 -0400 Subject: [PATCH] remove some more VLAs --- src/QEQ/fix_qeq_ctip.cpp | 138 +++++++++++++++++++++++++-------------- src/QEQ/fix_qeq_ctip.h | 3 + 2 files changed, 93 insertions(+), 48 deletions(-) diff --git a/src/QEQ/fix_qeq_ctip.cpp b/src/QEQ/fix_qeq_ctip.cpp index 56b06dbc81..0649c47412 100644 --- a/src/QEQ/fix_qeq_ctip.cpp +++ b/src/QEQ/fix_qeq_ctip.cpp @@ -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++; diff --git a/src/QEQ/fix_qeq_ctip.h b/src/QEQ/fix_qeq_ctip.h index 4e2af746cf..2c2e201123 100644 --- a/src/QEQ/fix_qeq_ctip.h +++ b/src/QEQ/fix_qeq_ctip.h @@ -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;