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) {
MPI_Allreduce(ftotal,ftotal_all,3,MPI_DOUBLE,MPI_SUM,world);
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;
MPI_Allreduce(&e_spring,&all,1,MPI_DOUBLE,MPI_SUM,world);
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.
Regards,
Xibing He
UPenn