[lammps-users] questions about the source code fix_spring.cpp

Dear LAMMPS developers and users,

I have several questions about the source code fix_spring.cpp.

(1) fix_spring.cpp can output the constrained spring force acting on a group through ftotal_all[0,1,2], which is summed from vector ftotal[0,1,2] , which in turn is got from the following codes (for fix spring tether):

ftotal[0] = ftotal[1] = ftotal[2] = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
massfrac = mass[type[i]]/masstotal;
f[i][0] -= fxmassfrac;
f[i][1] -= fy
f[i][2] -= fzmassfrac;
ftotal[0] -= fx
ftotal[1] -= fymassfrac;
ftotal[2] -= fz

I think that sum of massfrac over all atoms in a group should be equal to the masstotal of that group, so ftotal[0] should be equal to -fx, ftotal[1] should be equal to -fy, and ftotal[2] should be equal to -fz. So why not output -fx, -fy, -fz directly (ftotal[0] = -fx, etc.), but bother to do the loop?

I modified fix_spring.cpp to output (-fx,-fy,-fz) directly through ftotal_all[3,4,5] and compaired with ftoal_all[0,1,2]. Interestingly I found that ftotal_all[3,4,5] are not equal to ftotal_all[0,1,2], respectively, but exactly 4 times bigger. what’s going on here?

Please see the attached files for what I did.

(2) According to what I understand, LAMMPS doesn’t add the spring energy from “fix spring” to the total energy. If I’m right, why not??

Just before I’m going to send this email, another user Sebastian Fritsch sent an email to this list to ask a similar question about “fix spring/self”.

(3) If we compare the procedures FixSpring::setup, FixSpring::post_force_respa with corresponding procedures in fix_lineforce.cpp (FixLineForce::setup, FixLineForce::post_force_respa), we would find that they are different. For example,

void FixSpring::post_force_respa(int vflag, int ilevel, int iloop)
if (ilevel == nlevels_respa-1) post_force(vflag);

void FixLineForce::post_force_respa(int vflag, int ilevel, int iloop)

Can anyone explain a little bit to me why they need to do differently?

Thank you in advance.

Best regards,


fix_spring_test.cpp (7.57 KB)

fix_spring_test.h (1.26 KB)

re (1) - you are right that the total force could be accumulated outside
the loop, just as -fx,-fy,-fz. I'll change it since it's simpler -
your factor of
4 error must be b/c you did something else wrong.

re (2) - only some fixes have their energy added to the total E (if
the user requests it via fix_modify) - fix spring/self does b/c it's
useful for minimization - I don't think people often use fix spring
in that was - more for PMF calculations - so no one has added this,
but it could be done

re (3) - lineforce is a constraint which needs to applied at all levels
of rRESPA if the atoms are to stay on the line, the spring force is
an added force that should just be added once per timestep, so
its done only at one level of rRESPA


Hi, Steve,

Thanks for your reply. I saw the patch you released. But I have several more questions/suggestions:

(1) I think there is an error in your new released fix_spring.cpp, from line 224 to line 226:
ftotal[0] = -fx;
ftotal[1] = -fy;
ftotal[2] = -fz;
To output the spring force for “spring couple”, ftotal[0] (ftotal[1], ftotal[2]) should be equal to fx (fy,fz), NOT -fx (-fy, -fz). Because (fx,fy,fz) is the spring force acting on group1, (-fx,-fy,-fz) is the spring force acting on group2. The documents of “fix spring” says that “This fix computes a 3-vector of forces, …, In the case of the couple style, it is the force on the fix group (group-ID) or the negative of the force on the 2nd group (group-ID2).” So we should output (fx,fy,fz). Maybe you copied these three lines from “spring tether” but forgot to change the signs.

(2) Since ftotal[0] (ftotal[1], ftotal[2]) are directly from spring force (fx,fy,fz) which are not summed from local atoms, do we still need to do the “MPI_Allreduce” in the procedure FixSpring::compute_vector ?

if (force_flag == 0) {
force_flag = 1;
return ftotal_all[n];

Please pardon me that I don’t know what the procedure MPI_Allreduce exactly do. I cannot find the source code defining it. I guess it sums some quantities (in an array) related to the local atoms over different processors. But here (fx,fy,fz) are calculated from the center(s) of mass of group(s), which have been got from group->xcm(igroup,masstotal,xcm). If we don’t call MPI_Allreduce in FixSpring::compute_vector, instead we just assign ftotal[n] directly to ftotal_all[n] ( or just return ftotal[n] directly), will the result be the same or different???

(3) I think we should have an option to add the energy of the spring, as it is part of the “total” potential energy of the “total” system. I read source codes of several other MD packages. All of them add the spring energy to the total energy. Fortunately LAMMPS is flexible. We can add some energy contributions if we want, or don’t add if we don’t want, just like “fix spring/self”. I think at least we should also give such an option to “fix spring”.

I tried to add this option by myself. Please see the attached files fix_spring_suggest.cpp & fix_spring_suggest.h. I added

e_spring = 0.5 * k_spring * dr * dr;

in procedures FixSpring::spring_tether & FixSpring::spring_couple, and added the following procedure:

double FixSpring::compute_scalar()
double all;
return all;

Did I do it correctly? Is it enough? Again, just like question (2), do we need to call MPI_Allreduce? Can we just return e_spring directly?

(4) If you compare my fix_spring_suggest.cpp with the original fix_spring.cpp, you can find that I added two lines in procedures FixSpring::spring_tether & FixSpring::spring_couple:

if (dr >= 0.0) fsign = 1.0;
else fsign = -1.0;

And add the following line to FixSpring::compute_vector:

ftotal_all[3] = fsign*sqrt(ftotal_all[0]*ftotal_all[0]+ftotal_all[1]*ftotal_all[1]+ftotal_all[2]*ftotal_all[2]);

In other words, the original fix_spring.cpp just output the spring force as a vector (fx,fy,fz), but I think it’s better to output also f = fsign * sqrt(fx^2+fy^2+fz^2), where fsign equals to 1.0 or -1.0 depending on if dr >=0.0 or dr < 0.0.

Why do I suggest this modification? I think if we do constrained MD to calculate the pmf by using “fix spring”, we need to average this f = fsign * sqrt(fx^2+fy^2+fz^2) to get the correct constrained force ( and correct pmf). Without it, we can only average fx, fy, fz first (to get ave_fx, ave_fy, ave_fz), and then get sqrt(ave_fx^2+ave_fy^2+ave_fz^2). But this one would give the wrong constrained force, and lead to wrong pmf.

Let’s consider the following case:

Suppose there is a weird strong potential in a system, which repels molecule A away from the origin point (0,0,0). In another sentence, no matter where the position (x,y,z) of the molecule A is, it always feel an strong force pushing it away from (0,0,0). Now let’s use “fix spring” to fix it at a specific distance (e.g., 5.0 angstroms) from (0,0,0):

fix pull ligand spring tether 50.0 0.0 0.0 0.0 5.0

Since we suppose the strong potential likes to push molecule A away from (0,0,0), most of the times the distance between the real position of molecule A and the origin point should be greater than 5.0 angstrom.

Let’s suppose that one time the c.o.m. of the chosen molecule A is located at (6.0, 0.0, 0.0), so dx = 6.0 - 0.0 = 6.0, dy=dz=0.0 - 0.0 = 0.0, r = sqrt(6.0^2+0.0^2+0.0^2) = 6.0, dr = r - r0 = 6.0 - 5.0 = 1.0, and fx = k_smddxdr/r = 50.0, fy=fz=0.0

Another time the chose molecule A moves to (-6.0, 0.0, 0.0). so dx = -6.0, dy=dy=0.0, r = 6.0, dr = 1.0, HOWEVER, fx = k_smddxdr/r = -50.0.

If we output the spring force for these two moments, the original fix_spring.cpp would give:

#fx fy fz
50.0 0.0 0.0
-50.0 0.0 0.0

If we analyze this output file, we would get the following:

ave_fx = (50.0 + (-50.0)) / 2.0 = 0.0; ave_fy = 0.0; ave_fz = 0.0;
f_constrain = sqrt(ave_fx^2+ave_fy^2+ave_fz^2) = 0.0

Since f_constrain equals to zero, we may conclude that the molecule A is in equlibrium when it’s 5.0 angstroms away from the origin point (0.0, 0.0, 0.0). But we know it’s obviously wrong, because we suppose a strong potential always push molecule A away from origin point.

With the modification I did, the fix_spring_suggest.cpp would output the following results for these two specific moments:

fx fy fz f

50.0 0.0 0.0 50.0
-50.0 0.0 0.0 50.0

f_constrain = ave_f = (50.0 + 50.0) / 2.0 = 50.0

Yes, the f_constrain is not equal to zero, the postions which are 5.0 angstrom away from the origin point are not equalibrium positions for molecule A. This result agrees with our assumption.

Although the above example seems abnormal, I think it does show something.

Please let me know if I have made stupid mistakes or asked stupid questions in above discussion.


Xibing He

I’m sorry I forgot to attach the two file fix_spring_suggest.cpp & fix_spring_suggest.h in my last email.

Xibing He

fix_spring_suggest.cpp (7.22 KB)

fix_spring_suggest.h (1.24 KB)

Finally got around to this one. Here are my responses:

(1) yes, this was a sign error
(2) yes, this should not be summed in parallel for the "tether" option
(3) agreed, that we should add the spring energy to the system
     energy and enable it to be part of energy minimization

All of this is in a patch I just posted (18Jan09)

(4) this sounds good, but I want to run it by Paul, who
     understands PMF better than I

Paul - can you take a look at (4) below and see if you agree
we should add a 4th output field that has the sign value
included in the total force?


Hi, Steve,

Thanks for your responses. I read the patch you released. Here are my further opinions:

(1) Like the “tether” option, the “couple” option also should NOT sum the spring force in parallel. We should totally remove the MPI_Allreduce and force_flag in FixSpring::compute_vector(), and ouput ftotal[n] directly. Please check the attached fix_spring_suggest.cpp & fix_spring_suggest.h for what I think it should be done.

(2) Previously I had a weird problem with “fix spring”. One of my job ran well with 7 processors or less, but crashed with 8 processors or more. Dr. Axel Kohlmeyer helped me figured out the reason. In spring_tether() & spring_couple(), the variable r is the denominator when calculate fx,fy,fz. Occasionaly r is zero or underoverflow. To prevent the problem, we can add the following lines:

#define SMALL 1.e-20

// avoid division by zero or under/overflows
if (r < SMALL) r=SMALL;

Please see the attached fix_spring_suggest.cpp & fix_spring_suggest.h for details.

(3) As to the my suggestion of adding a 4th output field that has a sign value included in the total force, I think it’s useful for some cases, although no always useful. But anyway, it doesn’t hurt if we add it, because you can ignore it if you don’t need it.



fix_spring_suggest.cpp (7.44 KB)

fix_spring_suggest.h (1.27 KB)

ok - I patched this one more time - hopefully I got
it all right this time! See the 21Jan09 patch.

I agree with your (1), and added your (2) as a safety
check on the case when r = 0. Also added a 4th
output field (signed force) as you discussed in an
earlier email.

Let me know if all works as expected.


But how to calculate the PMF with this command?

Any suggestion?