diff --git a/doc/src/fix_spring_self.rst b/doc/src/fix_spring_self.rst index 4453fd61c5..f8354b17d7 100644 --- a/doc/src/fix_spring_self.rst +++ b/doc/src/fix_spring_self.rst @@ -13,7 +13,7 @@ Syntax * ID, group-ID are documented in :doc:`fix ` command * spring/self = style name of this fix command -* K = spring constant (force/distance units) +* K = spring constant (force/distance units), can be a variable (see below) * dir = xyz, xy, xz, yz, x, y, or z (optional, default: xyz) Examples @@ -22,6 +22,7 @@ Examples .. code-block:: LAMMPS fix tether boundary-atoms spring/self 10.0 + fix var all spring/self v_kvar fix zrest move spring/self 10.0 z Description @@ -42,6 +43,22 @@ directions, but it can be limited to the xy-, xz-, yz-plane and the x-, y-, or z-direction, thus restraining the atoms to a line or a plane, respectively. +The force constant *k* can be specified as an equal-style or atom-style +:doc:`variable `. If the value is a variable, it should be specified +as v_name, where name is the variable name. In this case, the variable +will be evaluated each time step, and its value(s) will be used as +force constant for the spring force. + +Equal-style variables can specify formulas with various mathematical +functions and include :doc:`thermo_style ` command +keywords for the simulation box parameters, time step, and elapsed time. +Thus, it is easy to specify a time-dependent force field. + +Atom-style variables can specify the same formulas as equal-style +variables but can also include per-atom values, such as atom +coordinates. Thus, it is easy to specify a spatially-dependent force +field with optional time-dependence as well. + Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" diff --git a/src/fix_spring_self.cpp b/src/fix_spring_self.cpp index df00a2ba8c..bc59ff9987 100644 --- a/src/fix_spring_self.cpp +++ b/src/fix_spring_self.cpp @@ -21,20 +21,23 @@ #include "atom.h" #include "domain.h" #include "error.h" +#include "input.h" #include "memory.h" #include "respa.h" #include "update.h" +#include "variable.h" #include using namespace LAMMPS_NS; using namespace FixConst; +enum { NONE, CONSTANT, EQUAL, ATOM }; + /* ---------------------------------------------------------------------- */ FixSpringSelf::FixSpringSelf(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), - xoriginal(nullptr) + Fix(lmp, narg, arg), xoriginal(nullptr), kstr(nullptr) { if ((narg < 4) || (narg > 5)) error->all(FLERR,"Illegal fix spring/self command"); @@ -46,8 +49,14 @@ FixSpringSelf::FixSpringSelf(LAMMPS *lmp, int narg, char **arg) : energy_global_flag = 1; respa_level_support = 1; - k = utils::numeric(FLERR,arg[3],false,lmp); - if (k <= 0.0) error->all(FLERR,"Illegal fix spring/self command"); + if (utils::strmatch(arg[3], "^v_")) { + kstr = utils::strdup(arg[3] + 2); + kstyle = NONE; + } else { + k = utils::numeric(FLERR,arg[3],false,lmp); + kstyle = CONSTANT; + if (k <= 0.0) error->all(FLERR,"Illegal force constatnt for fix spring/self command"); + } xflag = yflag = zflag = 1; @@ -123,6 +132,25 @@ int FixSpringSelf::setmask() void FixSpringSelf::init() { + // check variable + + if (kstr) { + kvar = input->variable->find(kstr); + if (kvar < 0) error->all(FLERR, "Variable {} for fix spring/self does not exist", kstr); + if (input->variable->equalstyle(kvar)) + kstyle = EQUAL; + else if (input->variable->atomstyle(kvar)) + kstyle = ATOM; + else + error->all(FLERR, "Variable {} for fix spring/self is invalid style", kstr); + } + + + if ((kstyle == ATOM) && (atom->nmax > maxatom)) { + maxatom = atom->nmax; + memory->destroy(kval); + memory->create(kval, maxatom, "sprint/self:kval"); + } if (utils::strmatch(update->integrate_style,"^respa")) { ilevel_respa = (dynamic_cast(update->integrate))->nlevels-1; if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa); @@ -162,24 +190,55 @@ void FixSpringSelf::post_force(int /*vflag*/) double dx,dy,dz; double unwrap[3]; + // reallocate kval array if necessary + + if ((kstyle == ATOM) && (atom->nmax > maxatom)) { + maxatom = atom->nmax; + memory->destroy(kval); + memory->create(kval, maxatom, "sprint/self:kval"); + } + espring = 0.0; - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - domain->unmap(x[i],image[i],unwrap); - dx = unwrap[0] - xoriginal[i][0]; - dy = unwrap[1] - xoriginal[i][1]; - dz = unwrap[2] - xoriginal[i][2]; - if (!xflag) dx = 0.0; - if (!yflag) dy = 0.0; - if (!zflag) dz = 0.0; - f[i][0] -= k*dx; - f[i][1] -= k*dy; - f[i][2] -= k*dz; - espring += k * (dx*dx + dy*dy + dz*dz); - } + if ((kstyle == CONSTANT) || (kstyle == EQUAL)) { + // update k if equal style variable + if (kstyle == EQUAL) k = input->variable->compute_equal(kvar); + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - xoriginal[i][0]; + dy = unwrap[1] - xoriginal[i][1]; + dz = unwrap[2] - xoriginal[i][2]; + if (!xflag) dx = 0.0; + if (!yflag) dy = 0.0; + if (!zflag) dz = 0.0; + f[i][0] -= k*dx; + f[i][1] -= k*dy; + f[i][2] -= k*dz; + espring += k * (dx*dx + dy*dy + dz*dz); + } + espring *= 0.5; + } else { + // update kval for kstyle == ATOM + input->variable->compute_atom(kvar, igroup, kval, 1, 0); + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - xoriginal[i][0]; + dy = unwrap[1] - xoriginal[i][1]; + dz = unwrap[2] - xoriginal[i][2]; + if (!xflag) dx = 0.0; + if (!yflag) dy = 0.0; + if (!zflag) dz = 0.0; + f[i][0] -= kval[i]*dx; + f[i][1] -= kval[i]*dy; + f[i][2] -= kval[i]*dz; + espring += kval[i] * (dx*dx + dy*dy + dz*dz); + } - espring *= 0.5; + espring *= 0.5; + } } /* ---------------------------------------------------------------------- */ @@ -213,7 +272,8 @@ double FixSpringSelf::compute_scalar() double FixSpringSelf::memory_usage() { - double bytes = (double)atom->nmax*3 * sizeof(double); + double bytes = (double)atom->nmax*4 * sizeof(double); + if (kstyle == ATOM) bytes += (double)atom->nmax * sizeof(double); return bytes; } diff --git a/src/fix_spring_self.h b/src/fix_spring_self.h index f13f2be918..e3220b8cfd 100644 --- a/src/fix_spring_self.h +++ b/src/fix_spring_self.h @@ -50,8 +50,11 @@ class FixSpringSelf : public Fix { protected: double k, espring; double **xoriginal; // original coords of atoms + char *kstr; // name of variable for K + double *kval; // per-atom variable values for K + int kvar, kstyle; int xflag, yflag, zflag; - int ilevel_respa; + int ilevel_respa, maxatom; }; } // namespace LAMMPS_NS diff --git a/unittest/force-styles/tests/fix-timestep-spring_self_atom.yaml b/unittest/force-styles/tests/fix-timestep-spring_self_atom.yaml new file mode 100644 index 0000000000..e8b35e6455 --- /dev/null +++ b/unittest/force-styles/tests/fix-timestep-spring_self_atom.yaml @@ -0,0 +1,77 @@ +--- +lammps_version: 17 Feb 2022 +date_generated: Fri Mar 18 22:18:01 2022 +epsilon: 5e-14 +skip_tests: kokkos_omp +prerequisites: ! | + atom full + fix spring/self +pre_commands: ! | + variable kvar atom 10.0 +post_commands: ! | + fix move all nve + fix test solute spring/self v_kvar xyz +input_file: in.fourmol +natoms: 29 +global_scalar: 0.12623705370750438 +run_pos: ! |2 + 1 -2.7045669792379945e-01 2.4912140072031601e+00 -1.6695897908630153e-01 + 2 3.1003572362014503e-01 2.9612290242130319e+00 -8.5466689099875615e-01 + 3 -7.0398410505917419e-01 1.2305522448803927e+00 -6.2777452858703953e-01 + 4 -1.5818137373512378e+00 1.4837442978868092e+00 -1.2538665277734968e+00 + 5 -9.0719266809604937e-01 9.2652365891304722e-01 3.9954299416556610e-01 + 6 2.4832271436421161e-01 2.8312320769893828e-01 -1.2314303818391423e+00 + 7 3.4143518373063453e-01 -2.2645187306940717e-02 -2.5292229406165907e+00 + 8 1.1743540824610306e+00 -4.8863162985819381e-01 -6.3783828665914544e-01 + 9 1.3800508883119953e+00 -2.5274565574966718e-01 2.8353436538531862e-01 + 10 2.0510726036138696e+00 -1.4604027090169940e+00 -9.8323554549077119e-01 + 11 1.7878052123899795e+00 -1.9921835931655048e+00 -1.8890566136917575e+00 + 12 3.0063000763288268e+00 -4.9013323763786637e-01 -1.6231890083142151e+00 + 13 4.0515352007695906e+00 -8.9202569828585798e-01 -1.6399995011867139e+00 + 14 2.6066952690925067e+00 -4.1788633645045253e-01 -2.6633949696742012e+00 + 15 2.9695337662435719e+00 5.5423141568538492e-01 -1.2342076641871542e+00 + 16 2.6747001492056977e+00 -2.4124097322472577e+00 -2.3429048072365732e-02 + 17 2.2153591049025310e+00 -2.0897997214660862e+00 1.1963106355359285e+00 + 18 2.1369701704115056e+00 3.0158507413628213e+00 -3.5179348337213843e+00 + 19 1.5355837136087338e+00 2.6255292355375399e+00 -4.2353987779878857e+00 + 20 2.7727573005678758e+00 3.6923910449610102e+00 -3.9330842459133470e+00 + 21 4.9040128073205524e+00 -4.0752348172959030e+00 -3.6210314709893159e+00 + 22 4.3582355554440877e+00 -4.2126119427287101e+00 -4.4612844196314150e+00 + 23 5.7439382849307670e+00 -3.5821957939275060e+00 -3.8766361295935892e+00 + 24 2.0689243582422914e+00 3.1513346907271247e+00 3.1550389754829422e+00 + 25 1.3045351331492516e+00 3.2665125705842941e+00 2.5111855257434352e+00 + 26 2.5809237402711318e+00 4.0117602605482858e+00 3.2212060529089945e+00 + 27 -1.9611343130357277e+00 -4.3563411931359841e+00 2.1098293115523528e+00 + 28 -2.7473562684513411e+00 -4.0200819932379330e+00 1.5830052163433952e+00 + 29 -1.3126000191359812e+00 -3.5962518039482934e+00 2.2746342468737817e+00 +run_vel: ! |2 + 1 8.1685220941861477e-03 1.6512578512542727e-02 4.7892799147935001e-03 + 2 5.4427456394786321e-03 5.1693257879352533e-03 -1.4414043022813649e-03 + 3 -8.2272458036248362e-03 -1.2923813188884230e-02 -4.0970749471144546e-03 + 4 -3.7660861920462349e-03 -6.5659911420830365e-03 -1.1120922532834726e-03 + 5 -1.1012635909013241e-02 -9.8847866026321157e-03 -2.8391869073674538e-03 + 6 -3.9665990411620729e-02 4.6803722380071487e-02 3.7135522426802389e-02 + 7 9.1016763589152859e-04 -1.0126055720737583e-02 -5.1556610019025714e-02 + 8 7.9043585267658430e-03 -3.3496064544244345e-03 3.4549326598010660e-02 + 9 1.5620907286754389e-03 3.7378245105921431e-03 1.5036774253075934e-02 + 10 2.9193799040059056e-02 -2.9242248165535563e-02 -1.5014281912567770e-02 + 11 -4.7797459644718264e-03 -3.7436196398511232e-03 -2.3410499103477603e-03 + 12 2.2686069875175316e-03 -3.4732729502899497e-04 -3.0627334265471650e-03 + 13 2.7456854188010020e-03 5.8081889921879817e-03 -7.9308949311655092e-04 + 14 3.5223319737918667e-03 -5.7842699330258648e-03 -3.9396805101296825e-03 + 15 -1.8475459117759364e-03 -5.8469790281561471e-03 6.2849983323582511e-03 + 16 1.8676069228413028e-02 -1.3258381729410438e-02 -4.5625616778429308e-02 + 17 -1.2893668780819389e-02 9.7505325833410258e-03 3.7288200735675299e-02 + 18 -8.0065794869105819e-04 -8.6270473288011819e-04 -1.4483040693746142e-03 + 19 1.2452390836051499e-03 -2.5061097119616180e-03 7.2998631010316650e-03 + 20 3.5930060229538143e-03 3.6938860309035470e-03 3.2322732687958995e-03 + 21 -1.4689220366910704e-03 -2.7352129796142532e-04 7.0581624168175334e-04 + 22 -7.0694199254520765e-03 -4.2577148925037030e-03 2.8079117611209205e-04 + 23 6.0446963117617505e-03 -1.4000131614895336e-03 2.5819754846773601e-03 + 24 3.1926367911308686e-04 -9.9445664741642462e-04 1.4999996978363057e-04 + 25 1.3789754526895179e-04 -4.4335894884219599e-03 -8.1808136698604454e-04 + 26 2.0485904035342870e-03 2.7813358633902757e-03 4.3245727149365584e-03 + 27 4.5604120291626942e-04 -1.0305523027244966e-03 2.1188058375789067e-04 + 28 -6.2544520861839200e-03 1.4127711176141612e-03 -1.8429821884806304e-03 + 29 6.4110631535703737e-04 3.1273432719578029e-03 3.7253671105604122e-03 +... diff --git a/unittest/force-styles/tests/fix-timestep-spring_self_equal.yaml b/unittest/force-styles/tests/fix-timestep-spring_self_equal.yaml new file mode 100644 index 0000000000..382cc2cc73 --- /dev/null +++ b/unittest/force-styles/tests/fix-timestep-spring_self_equal.yaml @@ -0,0 +1,77 @@ +--- +lammps_version: 17 Feb 2022 +date_generated: Fri Mar 18 22:18:01 2022 +epsilon: 5e-14 +skip_tests: kokkos_omp +prerequisites: ! | + atom full + fix spring/self +pre_commands: ! | + variable kvar equal 10.0 +post_commands: ! | + fix move all nve + fix test solute spring/self v_kvar xyz +input_file: in.fourmol +natoms: 29 +global_scalar: 0.12623705370750438 +run_pos: ! |2 + 1 -2.7045669792379945e-01 2.4912140072031601e+00 -1.6695897908630153e-01 + 2 3.1003572362014503e-01 2.9612290242130319e+00 -8.5466689099875615e-01 + 3 -7.0398410505917419e-01 1.2305522448803927e+00 -6.2777452858703953e-01 + 4 -1.5818137373512378e+00 1.4837442978868092e+00 -1.2538665277734968e+00 + 5 -9.0719266809604937e-01 9.2652365891304722e-01 3.9954299416556610e-01 + 6 2.4832271436421161e-01 2.8312320769893828e-01 -1.2314303818391423e+00 + 7 3.4143518373063453e-01 -2.2645187306940717e-02 -2.5292229406165907e+00 + 8 1.1743540824610306e+00 -4.8863162985819381e-01 -6.3783828665914544e-01 + 9 1.3800508883119953e+00 -2.5274565574966718e-01 2.8353436538531862e-01 + 10 2.0510726036138696e+00 -1.4604027090169940e+00 -9.8323554549077119e-01 + 11 1.7878052123899795e+00 -1.9921835931655048e+00 -1.8890566136917575e+00 + 12 3.0063000763288268e+00 -4.9013323763786637e-01 -1.6231890083142151e+00 + 13 4.0515352007695906e+00 -8.9202569828585798e-01 -1.6399995011867139e+00 + 14 2.6066952690925067e+00 -4.1788633645045253e-01 -2.6633949696742012e+00 + 15 2.9695337662435719e+00 5.5423141568538492e-01 -1.2342076641871542e+00 + 16 2.6747001492056977e+00 -2.4124097322472577e+00 -2.3429048072365732e-02 + 17 2.2153591049025310e+00 -2.0897997214660862e+00 1.1963106355359285e+00 + 18 2.1369701704115056e+00 3.0158507413628213e+00 -3.5179348337213843e+00 + 19 1.5355837136087338e+00 2.6255292355375399e+00 -4.2353987779878857e+00 + 20 2.7727573005678758e+00 3.6923910449610102e+00 -3.9330842459133470e+00 + 21 4.9040128073205524e+00 -4.0752348172959030e+00 -3.6210314709893159e+00 + 22 4.3582355554440877e+00 -4.2126119427287101e+00 -4.4612844196314150e+00 + 23 5.7439382849307670e+00 -3.5821957939275060e+00 -3.8766361295935892e+00 + 24 2.0689243582422914e+00 3.1513346907271247e+00 3.1550389754829422e+00 + 25 1.3045351331492516e+00 3.2665125705842941e+00 2.5111855257434352e+00 + 26 2.5809237402711318e+00 4.0117602605482858e+00 3.2212060529089945e+00 + 27 -1.9611343130357277e+00 -4.3563411931359841e+00 2.1098293115523528e+00 + 28 -2.7473562684513411e+00 -4.0200819932379330e+00 1.5830052163433952e+00 + 29 -1.3126000191359812e+00 -3.5962518039482934e+00 2.2746342468737817e+00 +run_vel: ! |2 + 1 8.1685220941861477e-03 1.6512578512542727e-02 4.7892799147935001e-03 + 2 5.4427456394786321e-03 5.1693257879352533e-03 -1.4414043022813649e-03 + 3 -8.2272458036248362e-03 -1.2923813188884230e-02 -4.0970749471144546e-03 + 4 -3.7660861920462349e-03 -6.5659911420830365e-03 -1.1120922532834726e-03 + 5 -1.1012635909013241e-02 -9.8847866026321157e-03 -2.8391869073674538e-03 + 6 -3.9665990411620729e-02 4.6803722380071487e-02 3.7135522426802389e-02 + 7 9.1016763589152859e-04 -1.0126055720737583e-02 -5.1556610019025714e-02 + 8 7.9043585267658430e-03 -3.3496064544244345e-03 3.4549326598010660e-02 + 9 1.5620907286754389e-03 3.7378245105921431e-03 1.5036774253075934e-02 + 10 2.9193799040059056e-02 -2.9242248165535563e-02 -1.5014281912567770e-02 + 11 -4.7797459644718264e-03 -3.7436196398511232e-03 -2.3410499103477603e-03 + 12 2.2686069875175316e-03 -3.4732729502899497e-04 -3.0627334265471650e-03 + 13 2.7456854188010020e-03 5.8081889921879817e-03 -7.9308949311655092e-04 + 14 3.5223319737918667e-03 -5.7842699330258648e-03 -3.9396805101296825e-03 + 15 -1.8475459117759364e-03 -5.8469790281561471e-03 6.2849983323582511e-03 + 16 1.8676069228413028e-02 -1.3258381729410438e-02 -4.5625616778429308e-02 + 17 -1.2893668780819389e-02 9.7505325833410258e-03 3.7288200735675299e-02 + 18 -8.0065794869105819e-04 -8.6270473288011819e-04 -1.4483040693746142e-03 + 19 1.2452390836051499e-03 -2.5061097119616180e-03 7.2998631010316650e-03 + 20 3.5930060229538143e-03 3.6938860309035470e-03 3.2322732687958995e-03 + 21 -1.4689220366910704e-03 -2.7352129796142532e-04 7.0581624168175334e-04 + 22 -7.0694199254520765e-03 -4.2577148925037030e-03 2.8079117611209205e-04 + 23 6.0446963117617505e-03 -1.4000131614895336e-03 2.5819754846773601e-03 + 24 3.1926367911308686e-04 -9.9445664741642462e-04 1.4999996978363057e-04 + 25 1.3789754526895179e-04 -4.4335894884219599e-03 -8.1808136698604454e-04 + 26 2.0485904035342870e-03 2.7813358633902757e-03 4.3245727149365584e-03 + 27 4.5604120291626942e-04 -1.0305523027244966e-03 2.1188058375789067e-04 + 28 -6.2544520861839200e-03 1.4127711176141612e-03 -1.8429821884806304e-03 + 29 6.4110631535703737e-04 3.1273432719578029e-03 3.7253671105604122e-03 +...