From fe502c71d379e18296d4fb2a3cc9b5223042ca94 Mon Sep 17 00:00:00 2001 From: phankl Date: Fri, 8 Apr 2022 16:53:12 +0100 Subject: [PATCH] moved contents of compute to mesolj function for future modularity --- src/MESONT/pair_mesocnt.cpp | 339 +++++++++++++++++++----------------- src/MESONT/pair_mesocnt.h | 13 +- 2 files changed, 184 insertions(+), 168 deletions(-) diff --git a/src/MESONT/pair_mesocnt.cpp b/src/MESONT/pair_mesocnt.cpp index c5d33c656e..20de200cc0 100644 --- a/src/MESONT/pair_mesocnt.cpp +++ b/src/MESONT/pair_mesocnt.cpp @@ -98,6 +98,159 @@ PairMesoCNT::~PairMesoCNT() /* ---------------------------------------------------------------------- */ void PairMesoCNT::compute(int eflag, int vflag) +{ + ev_init(eflag, vflag); + + mesolj(); + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- */ + +void PairMesoCNT::allocate() +{ + allocated = 1; + int ntypes = atom->ntypes; + int np1 = ntypes + 1; + int init_size = 1; + + memory->create(cutsq, np1, np1, "pair:cutsq"); + memory->create(setflag, np1, np1, "pair:setflag"); + for (int i = 1; i <= ntypes; i++) + for (int j = i; j <= ntypes; j++) setflag[i][j] = 0; + + memory->create(end_types, nend_types, "pair:end_types"); + + memory->create(uinf_coeff, uinf_points, 4, "pair:uinf_coeff"); + memory->create(gamma_coeff, gamma_points, 4, "pair:gamma_coeff"); + memory->create(phi_coeff, phi_points, phi_points, 4, 4, "pair:phi_coeff"); + memory->create(usemi_coeff, usemi_points, usemi_points, 4, 4, "pair:usemi_coeff"); + + memory->create(numchainlist, init_size, "pair:numchainlist"); + memory->create(nchainlist, init_size, init_size, "pair:nchainlist"); + memory->create(endlist, init_size, init_size, "pair:endlist"); + memory->create(chainlist, init_size, init_size, init_size, "pair:chainlist"); + + memory->create(w, init_size, "pair:w"); + memory->create(wnode, init_size, "pair:wnode"); + memory->create(dq_w, init_size, 3, "pair:dq_w"); + memory->create(q1_dq_w, init_size, 3, 3, "pair:q1_dq_w"); + memory->create(q2_dq_w, init_size, 3, 3, "pair:q2_dq_w"); + + memory->create(param, 7, "pair:param"); + + memory->create(flocal, 2, 3, "pair:flocal"); + memory->create(fglobal, 4, 3, "pair:fglobal"); + memory->create(basis, 3, 3, "pair:basis"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairMesoCNT::settings(int narg, char ** /* arg */) +{ + if (narg != 0) error->all(FLERR, "Illegal pair_style command"); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairMesoCNT::coeff(int narg, char **arg) +{ + if (narg < 4) error->all(FLERR, "Incorrect args for pair coefficients"); + read_file(arg[2]); + + nend_types = narg - 3; + + if (!allocated) allocate(); + + // end atom types + for (int i = 3; i < narg; i++) end_types[i - 3] = utils::inumeric(FLERR, arg[i], false, lmp); + + // units, eV to energy unit conversion + ang = force->angstrom; + ang_inv = 1.0 / ang; + if (strcmp(update->unit_style, "lj") == 0) + error->all(FLERR, "Pair style mesocnt does not support lj units"); + else if (strcmp(update->unit_style, "real") == 0) + eunit = 23.06054966; + else if (strcmp(update->unit_style, "metal") == 0) + eunit = 1.0; + else if (strcmp(update->unit_style, "si") == 0) + eunit = 1.6021765e-19; + else if (strcmp(update->unit_style, "cgs") == 0) + eunit = 1.6021765e-12; + else if (strcmp(update->unit_style, "electron") == 0) + eunit = 3.674932248e-2; + else if (strcmp(update->unit_style, "micro") == 0) + eunit = 1.6021765e-4; + else if (strcmp(update->unit_style, "nano") == 0) + eunit = 1.6021765e2; + else + error->all(FLERR, "Pair style mesocnt does not recognize this units style"); + funit = eunit * ang_inv; + + // potential variables + sig = sig_ang * ang; + r = r_ang * ang; + rsq = r * r; + d = 2.0 * r; + d_ang = 2.0 * r_ang; + rc = 3.0 * sig; + cutoff = rc + d; + cutoffsq = cutoff * cutoff; + cutoff_ang = cutoff * ang_inv; + cutoffsq_ang = cutoff_ang * cutoff_ang; + comega = 0.275 * (1.0 - 1.0 / (1.0 + 0.59 * r_ang)); + ctheta = 0.35 + 0.0226 * (r_ang - 6.785); + + // compute spline coefficients + spline_coeff(uinf_data, uinf_coeff, delh_uinf, uinf_points); + spline_coeff(gamma_data, gamma_coeff, delh_gamma, gamma_points); + spline_coeff(phi_data, phi_coeff, delh_phi, delpsi_phi, phi_points); + spline_coeff(usemi_data, usemi_coeff, delh_usemi, delxi_usemi, usemi_points); + + memory->destroy(uinf_data); + memory->destroy(gamma_data); + memory->destroy(phi_data); + memory->destroy(usemi_data); + + int ntypes = atom->ntypes; + for (int i = 1; i <= ntypes; i++) + for (int j = i; j <= ntypes; j++) setflag[i][j] = 1; +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairMesoCNT::init_style() +{ + if (atom->tag_enable == 0) error->all(FLERR, "Pair style mesocnt requires atom IDs"); + if (force->newton_pair == 0) error->all(FLERR, "Pair style mesocnt requires newton pair on"); + + // need a full neighbor list + + neighbor->add_request(this, NeighConst::REQ_FULL); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairMesoCNT::init_one(int /* i */, int /* j */) +{ + return cutoff; +} + +/* ---------------------------------------------------------------------- + calculate energy and forces due to mesoscopic LJ potential +------------------------------------------------------------------------- */ + +void PairMesoCNT::mesolj() { int i, j, k, i1, i2, j1, j2; int endflag, endindex; @@ -119,7 +272,6 @@ void PairMesoCNT::compute(int eflag, int vflag) double temp[3][3]; evdwl = 0.0; - ev_init(eflag, vflag); double **x = atom->x; double **f = atom->f; @@ -284,9 +436,9 @@ void PairMesoCNT::compute(int eflag, int vflag) if (endflag) { sub3(qe, p, delqe); - cross3(delqe, m, t3); - scale3(fend, t3); - add3(t3, torque, torque); + cross3(delqe, m, t3); + scale3(fend, t3); + add3(t3, torque, torque); } cross3(torque, m, ftorque); @@ -390,149 +542,6 @@ void PairMesoCNT::compute(int eflag, int vflag) } } } - - if (vflag_fdotr) virial_fdotr_compute(); -} - -/* ---------------------------------------------------------------------- */ - -void PairMesoCNT::allocate() -{ - allocated = 1; - int ntypes = atom->ntypes; - int np1 = ntypes + 1; - int init_size = 1; - - memory->create(cutsq, np1, np1, "pair:cutsq"); - memory->create(setflag, np1, np1, "pair:setflag"); - for (int i = 1; i <= ntypes; i++) - for (int j = i; j <= ntypes; j++) setflag[i][j] = 0; - - memory->create(end_types, nend_types, "pair:end_types"); - - memory->create(uinf_coeff, uinf_points, 4, "pair:uinf_coeff"); - memory->create(gamma_coeff, gamma_points, 4, "pair:gamma_coeff"); - memory->create(phi_coeff, phi_points, phi_points, 4, 4, "pair:phi_coeff"); - memory->create(usemi_coeff, usemi_points, usemi_points, 4, 4, "pair:usemi_coeff"); - - memory->create(numchainlist, init_size, "pair:numchainlist"); - memory->create(nchainlist, init_size, init_size, "pair:nchainlist"); - memory->create(endlist, init_size, init_size, "pair:endlist"); - memory->create(chainlist, init_size, init_size, init_size, "pair:chainlist"); - - memory->create(w, init_size, "pair:w"); - memory->create(wnode, init_size, "pair:wnode"); - memory->create(dq_w, init_size, 3, "pair:dq_w"); - memory->create(q1_dq_w, init_size, 3, 3, "pair:q1_dq_w"); - memory->create(q2_dq_w, init_size, 3, 3, "pair:q2_dq_w"); - - memory->create(param, 7, "pair:param"); - - memory->create(flocal, 2, 3, "pair:flocal"); - memory->create(fglobal, 4, 3, "pair:fglobal"); - memory->create(basis, 3, 3, "pair:basis"); -} - -/* ---------------------------------------------------------------------- - global settings -------------------------------------------------------------------------- */ - -void PairMesoCNT::settings(int narg, char ** /* arg */) -{ - if (narg != 0) error->all(FLERR, "Illegal pair_style command"); -} - -/* ---------------------------------------------------------------------- - set coeffs for one or more type pairs -------------------------------------------------------------------------- */ - -void PairMesoCNT::coeff(int narg, char **arg) -{ - if (narg < 4) error->all(FLERR, "Incorrect args for pair coefficients"); - read_file(arg[2]); - - nend_types = narg - 3; - - if (!allocated) allocate(); - - // end atom types - for (int i = 3; i < narg; i++) - end_types[i - 3] = utils::inumeric(FLERR, arg[i], false, lmp); - - // units, eV to energy unit conversion - ang = force->angstrom; - ang_inv = 1.0 / ang; - if (strcmp(update->unit_style, "lj") == 0) - error->all(FLERR, "Pair style mesocnt does not support lj units"); - else if (strcmp(update->unit_style, "real") == 0) - eunit = 23.06054966; - else if (strcmp(update->unit_style, "metal") == 0) - eunit = 1.0; - else if (strcmp(update->unit_style, "si") == 0) - eunit = 1.6021765e-19; - else if (strcmp(update->unit_style, "cgs") == 0) - eunit = 1.6021765e-12; - else if (strcmp(update->unit_style, "electron") == 0) - eunit = 3.674932248e-2; - else if (strcmp(update->unit_style, "micro") == 0) - eunit = 1.6021765e-4; - else if (strcmp(update->unit_style, "nano") == 0) - eunit = 1.6021765e2; - else - error->all(FLERR, "Pair style mesocnt does not recognize this units style"); - funit = eunit * ang_inv; - - // potential variables - sig = sig_ang * ang; - r = r_ang * ang; - rsq = r * r; - d = 2.0 * r; - d_ang = 2.0 * r_ang; - rc = 3.0 * sig; - cutoff = rc + d; - cutoffsq = cutoff * cutoff; - cutoff_ang = cutoff * ang_inv; - cutoffsq_ang = cutoff_ang * cutoff_ang; - comega = 0.275 * (1.0 - 1.0 / (1.0 + 0.59 * r_ang)); - ctheta = 0.35 + 0.0226 * (r_ang - 6.785); - - // compute spline coefficients - spline_coeff(uinf_data, uinf_coeff, delh_uinf, uinf_points); - spline_coeff(gamma_data, gamma_coeff, delh_gamma, gamma_points); - spline_coeff(phi_data, phi_coeff, delh_phi, delpsi_phi, phi_points); - spline_coeff(usemi_data, usemi_coeff, delh_usemi, delxi_usemi, usemi_points); - - memory->destroy(uinf_data); - memory->destroy(gamma_data); - memory->destroy(phi_data); - memory->destroy(usemi_data); - - int ntypes = atom->ntypes; - for (int i = 1; i <= ntypes; i++) - for (int j = i; j <= ntypes; j++) setflag[i][j] = 1; -} - -/* ---------------------------------------------------------------------- - init specific to this pair style -------------------------------------------------------------------------- */ - -void PairMesoCNT::init_style() -{ - if (atom->tag_enable == 0) error->all(FLERR, "Pair style mesocnt requires atom IDs"); - if (force->newton_pair == 0) error->all(FLERR, "Pair style mesocnt requires newton pair on"); - - // need a full neighbor list - - neighbor->add_request(this, NeighConst::REQ_FULL); -} - -/* ---------------------------------------------------------------------- - init for one type pair i,j and corresponding j,i -------------------------------------------------------------------------- */ - -double PairMesoCNT::init_one(int /* i */, int /* j */) -{ - return cutoff; } /* ---------------------------------------------------------------------- @@ -792,9 +801,12 @@ void PairMesoCNT::chain_split(int *redlist, int numred, int *nchain, int **chain int cstart = chain[j][0]; int cend = chain[j][nchain[j] - 1]; - if (match_end(atom->type[cstart])) end[j] = 1; - else if (match_end(atom->type[cend])) end[j] = 2; - else end[j] = 0; + if (match_end(atom->type[cstart])) + end[j] = 1; + else if (match_end(atom->type[cend])) + end[j] = 2; + else + end[j] = 0; } } @@ -912,8 +924,8 @@ void PairMesoCNT::read_file(const char *file) MPI_Bcast(uinf_data, uinf_points, MPI_DOUBLE, 0, world); MPI_Bcast(gamma_data, gamma_points, MPI_DOUBLE, 0, world); - MPI_Bcast(&phi_data[0][0], phi_points*phi_points, MPI_DOUBLE, 0, world); - MPI_Bcast(&usemi_data[0][0], usemi_points*usemi_points, MPI_DOUBLE, 0, world); + MPI_Bcast(&phi_data[0][0], phi_points * phi_points, MPI_DOUBLE, 0, world); + MPI_Bcast(&usemi_data[0][0], usemi_points * usemi_points, MPI_DOUBLE, 0, world); } /* ---------------------------------------------------------------------- @@ -1234,7 +1246,8 @@ void PairMesoCNT::spline_coeff(double **data, double ****coeff, double dx, doubl cubic spline evaluation ------------------------------------------------------------------------- */ -inline double PairMesoCNT::spline(double x, double xstart, double dx, double **coeff, int coeff_size) +inline double PairMesoCNT::spline(double x, double xstart, double dx, double **coeff, + int coeff_size) { int i = ceil((x - xstart) / dx); @@ -1262,7 +1275,8 @@ inline double PairMesoCNT::spline(double x, double xstart, double dx, double **c cubic spline derivative ------------------------------------------------------------------------- */ -inline double PairMesoCNT::dspline(double x, double xstart, double dx, double **coeff, int coeff_size) +inline double PairMesoCNT::dspline(double x, double xstart, double dx, double **coeff, + int coeff_size) { int i = ceil((x - xstart) / dx); @@ -1290,8 +1304,8 @@ inline double PairMesoCNT::dspline(double x, double xstart, double dx, double ** bicubic spline evaluation ------------------------------------------------------------------------- */ -inline double PairMesoCNT::spline(double x, double y, double xstart, double ystart, double dx, double dy, - double ****coeff, int coeff_size) +inline double PairMesoCNT::spline(double x, double y, double xstart, double ystart, double dx, + double dy, double ****coeff, int coeff_size) { int i = ceil((x - xstart) / dx); int j = ceil((y - ystart) / dy); @@ -1337,8 +1351,8 @@ inline double PairMesoCNT::spline(double x, double y, double xstart, double ysta bicubic spline partial x derivative ------------------------------------------------------------------------- */ -inline double PairMesoCNT::dxspline(double x, double y, double xstart, double ystart, double dx, double dy, - double ****coeff, int coeff_size) +inline double PairMesoCNT::dxspline(double x, double y, double xstart, double ystart, double dx, + double dy, double ****coeff, int coeff_size) { int i = ceil((x - xstart) / dx); int j = ceil((y - ystart) / dy); @@ -1382,8 +1396,8 @@ inline double PairMesoCNT::dxspline(double x, double y, double xstart, double ys bicubic spline partial y derivative ------------------------------------------------------------------------- */ -inline double PairMesoCNT::dyspline(double x, double y, double xstart, double ystart, double dx, double dy, - double ****coeff, int coeff_size) +inline double PairMesoCNT::dyspline(double x, double y, double xstart, double ystart, double dx, + double dy, double ****coeff, int coeff_size) { int i = ceil((x - xstart) / dx); int j = ceil((y - ystart) / dy); @@ -1522,8 +1536,9 @@ void PairMesoCNT::geometry(const double *r1, const double *r2, const double *p1, weight for substitute CNT chain ------------------------------------------------------------------------- */ -inline void PairMesoCNT::weight(const double *r1, const double *r2, const double *p1, const double *p2, - double &w, double *dr1_w, double *dr2_w, double *dp1_w, double *dp2_w) +inline void PairMesoCNT::weight(const double *r1, const double *r2, const double *p1, + const double *p2, double &w, double *dr1_w, double *dr2_w, + double *dp1_w, double *dp2_w) { double dr, dp, rhoc, rhomin, rho, frac, arg, factor; double r[3], p[3]; diff --git a/src/MESONT/pair_mesocnt.h b/src/MESONT/pair_mesocnt.h index 329150d7b0..70445587e7 100644 --- a/src/MESONT/pair_mesocnt.h +++ b/src/MESONT/pair_mesocnt.h @@ -67,8 +67,7 @@ class PairMesoCNT : public Pair { void sort(int *, int); void read_file(const char *); void read_data(PotentialFileReader &, double *, double &, double &, int); - void read_data(PotentialFileReader &, double **, double &, double &, double &, double &, - int); + void read_data(PotentialFileReader &, double **, double &, double &, double &, double &, int); void spline_coeff(double *, double **, double, int); void spline_coeff(double **, double ****, double, double, int); @@ -79,11 +78,13 @@ class PairMesoCNT : public Pair { void finf(const double *, double &, double **); void fsemi(const double *, double &, double &, double **); + void mesolj(); + // inlined functions for efficiency - - inline void weight(const double *, const double *, const double *, const double *, double &, double *, - double *, double *, double *); - + + inline void weight(const double *, const double *, const double *, const double *, double &, + double *, double *, double *, double *); + inline double spline(double, double, double, double **, int); inline double dspline(double, double, double, double **, int); inline double spline(double, double, double, double, double, double, double ****, int);