From 03ea88c21674daf5a79a2b32019d01e420fa1148 Mon Sep 17 00:00:00 2001 From: monagocs Date: Sat, 1 Aug 2026 14:26:22 +0200 Subject: [PATCH] Fix harmonic force for fix spring/chunk --- src/fix_spring_chunk.cpp | 16 +-- unittest/commands/CMakeLists.txt | 4 + unittest/commands/test_fix_spring_chunk.cpp | 136 ++++++++++++++++++ .../tests/fix-timestep-spring_chunk.yaml | 132 ++++++++--------- 4 files changed, 212 insertions(+), 76 deletions(-) create mode 100644 unittest/commands/test_fix_spring_chunk.cpp diff --git a/src/fix_spring_chunk.cpp b/src/fix_spring_chunk.cpp index 144e8afc693..a6cb6db4f8a 100644 --- a/src/fix_spring_chunk.cpp +++ b/src/fix_spring_chunk.cpp @@ -24,14 +24,11 @@ #include "respa.h" #include "update.h" -#include #include using namespace LAMMPS_NS; using namespace FixConst; -static constexpr double SMALL = 1.0e-10; - /* ---------------------------------------------------------------------- */ FixSpringChunk::FixSpringChunk(LAMMPS *lmp, int narg, char **arg) : @@ -138,7 +135,7 @@ void FixSpringChunk::min_setup(int vflag) void FixSpringChunk::post_force(int /*vflag*/) { int i,m; - double dx,dy,dz,r; + double dx,dy,dz,rsq; // check if first time cchunk will be queried via ccom // if so, lock idchunk for as long as this fix is in place @@ -178,14 +175,13 @@ void FixSpringChunk::post_force(int /*vflag*/) dx = com[m][0] - com0[m][0]; dy = com[m][1] - com0[m][1]; dz = com[m][2] - com0[m][2]; - r = sqrt(dx*dx + dy*dy + dz*dz); - r = MAX(r,SMALL); + rsq = dx*dx + dy*dy + dz*dz; if (masstotal[m] != 0.0) { - fcom[m][0] = k_spring*dx/r / masstotal[m]; - fcom[m][1] = k_spring*dy/r / masstotal[m]; - fcom[m][2] = k_spring*dz/r / masstotal[m]; - esprings += 0.5*k_spring*r*r; + fcom[m][0] = k_spring*dx / masstotal[m]; + fcom[m][1] = k_spring*dy / masstotal[m]; + fcom[m][2] = k_spring*dz / masstotal[m]; + esprings += 0.5*k_spring*rsq; } } diff --git a/unittest/commands/CMakeLists.txt b/unittest/commands/CMakeLists.txt index 473d85f706a..ba91262bdda 100644 --- a/unittest/commands/CMakeLists.txt +++ b/unittest/commands/CMakeLists.txt @@ -44,6 +44,10 @@ add_executable(test_lattice_region test_lattice_region.cpp) target_link_libraries(test_lattice_region PRIVATE lammps GTest::GMock) add_test(NAME LatticeRegion COMMAND test_lattice_region) +add_executable(test_fix_spring_chunk test_fix_spring_chunk.cpp) +target_link_libraries(test_fix_spring_chunk PRIVATE lammps GTest::GMock) +add_test(NAME FixSpringChunk COMMAND test_fix_spring_chunk) + add_executable(test_kspace_auto_slab test_kspace_auto_slab.cpp) target_compile_definitions(test_kspace_auto_slab PRIVATE -DTEST_INPUT_FOLDER=${CMAKE_CURRENT_SOURCE_DIR}) target_link_libraries(test_kspace_auto_slab PRIVATE lammps GTest::GMock) diff --git a/unittest/commands/test_fix_spring_chunk.cpp b/unittest/commands/test_fix_spring_chunk.cpp new file mode 100644 index 00000000000..d9b613b735a --- /dev/null +++ b/unittest/commands/test_fix_spring_chunk.cpp @@ -0,0 +1,136 @@ +/* ---------------------------------------------------------------------- + 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. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "../testing/core.h" + +#include "atom.h" +#include "fix.h" +#include "modify.h" +#include "utils.h" + +#include "gmock/gmock.h" +#include "gtest/gtest.h" + +#include +#include +#include +#include + +// whether to print verbose output (i.e. not capturing LAMMPS screen output). +bool verbose = false; + +namespace LAMMPS_NS { + +class FixSpringChunkTest : public LAMMPSTest { +protected: + void SetUp() override + { + testbinary = "FixSpringChunkTest"; + LAMMPSTest::SetUp(); + + BEGIN_HIDE_OUTPUT(); + command("units lj"); + command("atom_style atomic"); + command("boundary f f f"); + command("region box block -2 2 -2 2 -2 2"); + command("create_box 1 box"); + command("create_atoms 1 single 0.0 0.0 0.0 units box"); + command("mass 1 1.0"); + command("pair_style zero 1.0"); + command("pair_coeff * *"); + command("compute chunks all chunk/atom type"); + command("compute centers all com/chunk chunks"); + command("fix integrate all nve"); + command("fix spring all spring/chunk 2.0 chunks centers"); + command("run 0 post no"); + END_HIDE_OUTPUT(); + } +}; + +TEST_F(FixSpringChunkTest, HarmonicForceAndEnergy) +{ + ASSERT_EQ(lmp->atom->nlocal, 1); + + auto *spring = lmp->modify->get_fix_by_id("spring"); + ASSERT_NE(spring, nullptr); + + auto **force = lmp->atom->f; + EXPECT_DOUBLE_EQ(spring->compute_scalar(), 0.0); + EXPECT_DOUBLE_EQ(force[0][0], 0.0); + EXPECT_DOUBLE_EQ(force[0][1], 0.0); + EXPECT_DOUBLE_EQ(force[0][2], 0.0); + + BEGIN_HIDE_OUTPUT(); + command("displace_atoms all move 0.25 0.0 0.0 units box"); + command("run 0 post no"); + END_HIDE_OUTPUT(); + + EXPECT_DOUBLE_EQ(spring->compute_scalar(), 0.0625); + EXPECT_DOUBLE_EQ(force[0][0], -0.5); + EXPECT_DOUBLE_EQ(force[0][1], 0.0); + EXPECT_DOUBLE_EQ(force[0][2], 0.0); + + // Independently verify that the applied force is the negative numerical + // derivative of the energy returned by the fix. + + constexpr double delta = 1.0e-4; + + BEGIN_HIDE_OUTPUT(); + command("displace_atoms all move 0.0001 0.0 0.0 units box"); + command("run 0 post no"); + END_HIDE_OUTPUT(); + const double eplus = spring->compute_scalar(); + + BEGIN_HIDE_OUTPUT(); + command("displace_atoms all move -0.0002 0.0 0.0 units box"); + command("run 0 post no"); + END_HIDE_OUTPUT(); + const double eminus = spring->compute_scalar(); + + BEGIN_HIDE_OUTPUT(); + command("displace_atoms all move 0.0001 0.0 0.0 units box"); + command("run 0 post no"); + END_HIDE_OUTPUT(); + + const double numerical_force = -(eplus-eminus) / (2.0*delta); + EXPECT_NEAR(force[0][0], numerical_force, 1.0e-12); + + BEGIN_HIDE_OUTPUT(); + command("displace_atoms all move 0.25 0.0 0.0 units box"); + command("run 0 post no"); + END_HIDE_OUTPUT(); + + EXPECT_DOUBLE_EQ(spring->compute_scalar(), 0.25); + EXPECT_DOUBLE_EQ(force[0][0], -1.0); +} + +} // namespace LAMMPS_NS + +int main(int argc, char **argv) +{ + MPI_Init(&argc, &argv); + ::testing::InitGoogleMock(&argc, argv); + + // handle arguments passed via environment variable + if (const char *var = getenv("TEST_ARGS")) { + std::vector env = LAMMPS_NS::utils::split_words(var); + for (auto arg : env) { + if (arg == "-v") verbose = true; + } + } + + if ((argc > 1) && (strcmp(argv[1], "-v") == 0)) verbose = true; + + int rv = RUN_ALL_TESTS(); + MPI_Finalize(); + return rv; +} diff --git a/unittest/force-styles/tests/fix-timestep-spring_chunk.yaml b/unittest/force-styles/tests/fix-timestep-spring_chunk.yaml index 1ad2f5693e4..bf9fa4d2695 100644 --- a/unittest/force-styles/tests/fix-timestep-spring_chunk.yaml +++ b/unittest/force-styles/tests/fix-timestep-spring_chunk.yaml @@ -1,13 +1,13 @@ --- -lammps_version: 17 Feb 2022 -date_generated: Fri Mar 18 22:18:01 2022 +lammps_version: 4 Jul 2026 +date_generated: Sat Aug 1 13:58:57 2026 epsilon: 1e-11 skip_tests: -prerequisites: ! | +prerequisites: | atom full fix spring/chunk -pre_commands: ! "" -post_commands: ! | +pre_commands: "" +post_commands: | compute mol all chunk/atom molecule compute com all com/chunk mol fix move all nve @@ -15,65 +15,65 @@ post_commands: ! | input_file: in.fourmol input_coeffs: coeffs.fourmol natoms: 29 -global_scalar: 0.0015064306233739035 -run_pos: ! |2 - 1 -2.7039466507239462e-01 2.4911394617880092e+00 -1.6690737231056296e-01 - 2 3.1010122768328879e-01 2.9611589338277575e+00 -8.5461248366340781e-01 - 3 -7.0392460911452404e-01 1.2304744899143512e+00 -6.2772412657984222e-01 - 4 -1.5817550043018951e+00 1.4836642574190744e+00 -1.2538199315725116e+00 - 5 -9.0713670028527926e-01 9.2644450405390477e-01 3.9959325076039032e-01 - 6 2.4836774058316266e-01 2.8306202887482079e-01 -1.2313790211702345e+00 - 7 3.4149618856823633e-01 -2.2723103032076725e-02 -2.5291782002728289e+00 - 8 1.1743069518257749e+00 -4.8858533080502903e-01 -6.3786027586713800e-01 - 9 1.3799953447164823e+00 -2.5269321030545128e-01 2.8350792470878511e-01 - 10 2.0510195983770423e+00 -1.4603526091581769e+00 -9.8326956153692791e-01 - 11 1.7877461018297200e+00 -1.9921323697343518e+00 -1.8890923030745981e+00 - 12 3.0062436175009810e+00 -4.9007955039546353e-01 -1.6232218687538209e+00 - 13 4.0514832085934431e+00 -8.9196616161398612e-01 -1.6400326104660659e+00 - 14 2.6066392471730087e+00 -4.1783858526565049e-01 -2.6634324182974152e+00 - 15 2.9694716311437053e+00 5.5428008601350487e-01 -1.2342342595622477e+00 - 16 2.6746458832799949e+00 -2.4123579520309462e+00 -2.3467804719930124e-02 - 17 2.2153006912112221e+00 -2.0897445643946391e+00 1.1962830221078402e+00 - 18 2.1368790264089581e+00 3.0159754143740933e+00 -3.5181977155490975e+00 - 19 1.5354925695842756e+00 2.6256539085382049e+00 -4.2356616598027150e+00 - 20 2.7726661565415887e+00 3.6925157179587416e+00 -3.9333471277222354e+00 - 21 4.9042004233221554e+00 -4.0750508982058031e+00 -3.6211861688914588e+00 - 22 4.3584231714404913e+00 -4.2124280236247813e+00 -4.4614391175232928e+00 - 23 5.7441259009277346e+00 -3.5820118748237770e+00 -3.8767908274855909e+00 - 24 2.0687885826955674e+00 3.1515550738671596e+00 3.1548778739329562e+00 - 25 1.3043993576001649e+00 3.2667329537232512e+00 2.5110244241852970e+00 - 26 2.5807879647318708e+00 4.0119806436999808e+00 3.2210449513704411e+00 - 27 -1.9608878136387959e+00 -4.3563688177752127e+00 2.1096520197822728e+00 - 28 -2.7471097690679747e+00 -4.0201096178908076e+00 1.5828279245795154e+00 - 29 -1.3123535197516709e+00 -3.5962794285999458e+00 2.2744569551087350e+00 -run_vel: ! |2 - 1 8.2391792369945933e-03 1.6438210431911154e-02 4.8466396546764535e-03 - 2 5.5187515357448692e-03 5.1009727439991102e-03 -1.3808762794824425e-03 - 3 -8.1613173851465357e-03 -1.3004672634959346e-02 -4.0419997879360741e-03 - 4 -3.7013110268827817e-03 -6.6504700541051183e-03 -1.0620290427688151e-03 - 5 -1.0953343579757559e-02 -9.9688939579476148e-03 -2.7846596077833229e-03 - 6 -3.9627794721568313e-02 4.6754811634285347e-02 3.7192208638621563e-02 - 7 9.7890622933070946e-04 -1.0206782168930003e-02 -5.1512311389882055e-02 - 8 7.8618361417075872e-03 -3.3078882607816048e-03 3.4531446917267020e-02 - 9 1.5032844351316340e-03 3.7927192074966782e-03 1.5010708889964666e-02 - 10 2.9140729691997380e-02 -2.9194095046793215e-02 -1.5055379032590480e-02 - 11 -4.8448308259422456e-03 -3.6920754803050761e-03 -2.3835363308950717e-03 - 12 2.2084299651528635e-03 -2.9168788114909022e-04 -3.1012020517987618e-03 - 13 2.6919570760867940e-03 5.8731594041239775e-03 -8.3179740851101734e-04 - 14 3.4634010415991240e-03 -5.7379466155056524e-03 -3.9849657546930702e-03 - 15 -1.9160117674839707e-03 -5.7994201845185190e-03 6.2567260808431150e-03 - 16 1.8620286902083130e-02 -1.3206416893569875e-02 -4.5675777789058840e-02 - 17 -1.2957486948358956e-02 9.8088191535149437e-03 3.7259412994009797e-02 - 18 -8.9827950617384468e-04 -7.2910542260445986e-04 -1.7299645349958738e-03 - 19 1.1476174779710083e-03 -2.3725104247018665e-03 7.0182026640168271e-03 - 20 3.4953844133924503e-03 3.8274853118860971e-03 2.9506128446975953e-03 - 21 -1.2679016851845878e-03 -7.6471662408769499e-05 5.4006461308053711e-04 - 22 -6.8683995853632369e-03 -4.0606652263220243e-03 1.1503957004042528e-04 - 23 6.2457166531648027e-03 -1.2029634957433778e-03 2.4162238784361132e-03 - 24 1.7379490178209356e-04 -7.5832030593805120e-04 -2.2597269496325511e-05 - 25 -7.5712376373799337e-06 -4.1974531497480138e-03 -9.9067862512451340e-04 - 26 1.9031216426784775e-03 3.0174722303355750e-03 4.1519755012089715e-03 - 27 7.2014620941965822e-04 -1.0601592030480235e-03 2.1924430684776849e-05 - 28 -5.9903471097165883e-03 1.3831641870693575e-03 -2.0329383278098810e-03 - 29 9.0521129402145969e-04 3.0977363443755040e-03 3.5354109686445011e-03 +global_scalar: 0.001556817966196982 +run_pos: |2 + 1 -2.7045540846604299e-01 2.4912157862348101e+00 -1.6695836395837596e-01 + 2 3.1004048502624637e-01 2.9612352587757633e+00 -8.5466347641229912e-01 + 3 -7.0398532474247566e-01 1.2305507912694424e+00 -6.2777511548740617e-01 + 4 -1.5818157443647520e+00 1.4837405775632138e+00 -1.2538709296439474e+00 + 5 -9.0719744743695763e-01 9.2652083451857281e-01 3.9954225884099936e-01 + 6 2.4831738142821377e-01 2.8313002125795311e-01 -1.2314231878364386e+00 + 7 3.4143546566302785e-01 -2.2646755380746321e-02 -2.5292289877380845e+00 + 8 1.1743551020649519e+00 -4.8863217581806268e-01 -6.3783440464568997e-01 + 9 1.3800522909313353e+00 -2.5274709144373564e-01 2.8353977579385270e-01 + 10 2.0510763901548641e+00 -1.4604062553547221e+00 -9.8323753403920122e-01 + 11 1.7878030624098862e+00 -1.9921862084629869e+00 -1.8890603279409439e+00 + 12 3.0063005719051885e+00 -4.9013338614027419e-01 -1.6231898939166425e+00 + 13 4.0515401638892241e+00 -8.9201999723724390e-01 -1.6400006361702331e+00 + 14 2.6066962025242657e+00 -4.1789242082588229e-01 -2.6634004440571353e+00 + 15 2.9695285865411347e+00 5.5422625048157759e-01 -1.2342022853566925e+00 + 16 2.6747028374936130e+00 -2.4124117866278949e+00 -2.3435829329288980e-02 + 17 2.2153576464983131e+00 -2.0897983998615786e+00 1.1963149962703019e+00 + 18 2.1369701076348355e+00 3.0158508272588955e+00 -3.5179350148222159e+00 + 19 1.5355836508319900e+00 2.6255293214334094e+00 -4.2353989590886165e+00 + 20 2.7727572377911294e+00 3.6923911308568576e+00 -3.9330844270140575e+00 + 21 4.9040130701421862e+00 -4.0752345596616282e+00 -3.6210316876987383e+00 + 22 4.3582358182658352e+00 -4.2126116850945250e+00 -4.4612846363409622e+00 + 23 5.7439385477525091e+00 -3.5821955362933253e+00 -3.8766363463031381e+00 + 24 2.0689242261060365e+00 3.1513349052148687e+00 3.1550388187028293e+00 + 25 1.3045350010129848e+00 3.2665127850720510e+00 2.5111853689632930e+00 + 26 2.5809236081349072e+00 4.0117604750360618e+00 3.2212058961289500e+00 + 27 -1.9611341682504686e+00 -4.3563412093646310e+00 2.1098292074167424e+00 + 28 -2.7473561236660933e+00 -4.0200820094665959e+00 1.5830051122077704e+00 + 29 -1.3125998743507372e+00 -3.5962518201769562e+00 2.2746341427381584e+00 +run_vel: |2 + 1 8.1709479039938568e-03 1.6516026247317694e-02 4.7905268784887945e-03 + 2 5.4505228277743666e-03 5.1787900100627874e-03 -1.4369927027955295e-03 + 3 -8.2294558804281541e-03 -1.2926931464787019e-02 -4.0981176568022949e-03 + 4 -3.7695307846928750e-03 -6.5726691241188981e-03 -1.1181635673188559e-03 + 5 -1.1021587505504849e-02 -9.8910580437927033e-03 -2.8407733353776484e-03 + 6 -3.9676324744187856e-02 4.6816710120698067e-02 3.7148769724226455e-02 + 7 9.1071290577234422e-04 -1.0128904362930896e-02 -5.1567951965237516e-02 + 8 7.9062497393305089e-03 -3.3505281013558553e-03 3.4556953174143826e-02 + 9 1.5641664460075350e-03 3.7367761440903341e-03 1.5047243734689530e-02 + 10 2.9201196054412935e-02 -2.9249357880886114e-02 -1.5018243129380812e-02 + 11 -4.7838474367245037e-03 -3.7479170590153797e-03 -2.3465759352144103e-03 + 12 2.2693939204034663e-03 -3.4752009968797563e-04 -3.0642425526585717e-03 + 13 2.7529227778053038e-03 5.8173276051690630e-03 -7.9484005847234079e-04 + 14 3.5243669696240418e-03 -5.7937781148821223e-03 -3.9480086353338977e-03 + 15 -1.8550456317789994e-03 -5.8552515507128275e-03 6.2936829961238536e-03 + 16 1.8681248709748240e-02 -1.3262244764225677e-02 -4.5638816979057750e-02 + 17 -1.2896521248452229e-02 9.7529879699703899e-03 3.7296369843040272e-02 + 18 -8.0075358612373112e-04 -8.6257381156824311e-04 -1.4485800590071114e-03 + 19 1.2451434459399820e-03 -2.5059787912943996e-03 7.2995871117151534e-03 + 20 3.5929103852771417e-03 3.6940169515049443e-03 3.2319972795414936e-03 + 21 -1.4685215472155209e-03 -2.7312873240596604e-04 7.0548601589149706e-04 + 22 -7.0690194356153148e-03 -4.2573223272425577e-03 2.8046094991073119e-04 + 23 6.0450968015870162e-03 -1.3996205962346756e-03 2.5816452584692855e-03 + 24 3.1906234420715307e-04 -9.9412980711776266e-04 1.4976109294202351e-04 + 25 1.3769621032698454e-04 -4.4332626480807777e-03 -8.1832024392057058e-04 + 26 2.0483890687257840e-03 2.7816627037877371e-03 4.3243338383028062e-03 + 27 4.5626181486638998e-04 -1.0305770378506492e-03 2.1172190907291279e-04 + 28 -6.2542314742682062e-03 1.4127463824377769e-03 -1.8431408632110830e-03 + 29 6.4132692726218230e-04 3.1273185367841167e-03 3.7252084358326302e-03 ...