fix: try to solve #87 for non-scala components #713
Conversation
|
It ought to be noted that even without considering the force on the fitting group, the force applying on the main group is wrong in the (*atoms)[ia].apply_force(FQ[i] * dq0_2[i]);I think at least we have to rotate the force back to the original frame: ag_rot = atoms->rot.inverse().matrix();
const auto f_ia = FQ[0] * dq0_2[0] + FQ[1] * dq0_2[1] + FQ[2] * dq0_2[2] + FQ[3] * dq0_2[3];
(*atoms)[ia].apply_force(ag_rot * f_ia);@fhh2626 Maybe this is the reason why you found some scripted CVs based on |
|
After more testing I conclude that the approach in this PR yields the correct result. I export the step 0 positions in #!/usr/bin/env python3
import torch
torch.set_printoptions(precision=10)
# function to extract grad
def set_grad(var):
def hook(grad):
var.grad = grad
return hook
main_group = torch.tensor(
[[ 7.00603894, -0.47702751, -0.68594361],
[ 6.1849732 , -0.13302646, 3.0389999 ],
[ 5.2349905 , 3.60400329, 2.87702477],
[ 2.6549758 , 3.21405895, -0.05499306],
[ 1.11002766, 0.02001201, 1.68902446],
[ 0.79599642, 1.8699884 , 5.09001061],
[-1.44499986, 4.62799966, 3.52899942],
[-3.74298744, 1.86798795, 2.08600524],
[-3.74097344, 0.09904095, 5.46300172],
[-4.90700712, 3.24699298, 7.21803419]], requires_grad=True, dtype=torch.float64)
ref_main_group = torch.tensor(
[[6.3061, -1.592, -3.585],
[5.4031, -0.784, 0.045] ,
[3.8061, 2.804, -0.129] ,
[1.2441, 1.799, -2.726] ,
[0.5861, -1.803, -1.516],
[-0.2169, -0.926, 2.133],
[-2.1679, 2.13, 1.258] ,
[-4.4689, 0.434, -1.263],
[-4.8649, -2.25, 1.46] ,
[-5.6269, 0.188, 4.323]], dtype=torch.float64)
fitting_group = torch.tensor(
[[ 8.42800577, 0.0420079 , -0.74403081],
[ 7.00603894, -0.47702751, -0.68594361],
[ 7.03995611, -2.00300543, -1.04193958],
[ 6.232 , -0.252 , 0.617 ],
[ 5.034 , -0.207 , 0.633 ],
[ 6.91001276, -0.26600248, 1.77300131],
[ 6.1849732 , -0.13302646, 3.0389999 ],
[ 7.20098994, -0.31695137, 4.30201369],
[ 5.3 , 1.147 , 3.253 ],
[ 4.096 , 1.082 , 3.6 ],
[ 5.87400947, 2.34899962, 2.93599796],
[ 5.2349905 , 3.60400329, 2.87702477],
[ 6.31703081, 4.65098894, 2.4770467 ],
[ 3.911 , 3.676 , 2.035 ],
[ 2.936 , 4.121 , 2.584 ],
[ 3.81797897, 3.24301039, 0.78301317],
[ 2.6549758 , 3.21405895, -0.05499306],
[ 3.04500477, 2.74098549, -1.55998623],
[ 1.656 , 2.181 , 0.554 ],
[ 0.423 , 2.412 , 0.581 ],
[ 2.00097516, 0.98201007, 1.06700298],
[ 1.11002766, 0.02001201, 1.68902446],
[ 2.01697998, -1.21900134, 2.02197084],
[ 0.49 , 0.565 , 3.013 ],
[-0.693 , 0.401 , 3.207 ],
[ 1.25802926, 1.29900304, 3.82599577],
[ 0.79599642, 1.8699884 , 5.09001061],
[ 2.074998 , 2.57305516, 5.62095664],
[-0.345 , 2.92 , 4.99 ],
[-1.187 , 3.009 , 5.894 ],
[-0.37699068, 3.70699752, 3.88398801],
[-1.44499986, 4.62799966, 3.52899942],
[-0.78001145, 5.37397741, 2.40196061],
[-2.7 , 3.899 , 2.954 ],
[-3.847 , 4.392 , 2.883 ],
[-2.54397388, 2.62898609, 2.51399134],
[-3.74298744, 1.86798795, 2.08600524],
[-3.40099076, 0.78497807, 1.02698123],
[-4.356 , 1.13 , 3.214 ],
[-5.528 , 0.96 , 3.354 ],
[-3.49404539, 0.63298789, 4.18601013],
[-3.74097344, 0.09904095, 5.46300172],
[-2.46600202, -0.44904453, 6.11800127],
[-4.453 , 1.029 , 6.423 ],
[-5.33 , 0.583 , 7.172 ],
[-6.288 , 3.553 , 6.688 ],
[-6.941 , 4.323 , 7.42 ],
[-6.77997934, 3.20500922, 5.49202067],
[-4.07700722, 2.28599541, 6.46500632],
[-4.90700712, 3.24699298, 7.21803419],
[-4.04901756, 4.48202806, 7.26299742]], requires_grad=True, dtype=torch.float64)
ref_fitting_group = torch.tensor(
[[5.69469, -1.80418, -5.22382] ,
[6.37169, -1.08118, -4.10382] ,
[6.61169, 0.365824, -4.58082] ,
[5.62269, -1.09518, -2.77982] ,
[4.39969, -0.986176, -2.72882] ,
[6.36569, -1.21018, -1.66282] ,
[5.83469, -1.37018, -0.318824] ,
[7.00969, -1.45018, 0.678176] ,
[4.84269, -0.296176, 0.126176] ,
[3.77369, -0.596176, 0.648176] ,
[5.15669, 0.991824, -0.113824] ,
[4.27869, 2.10682, 0.190176] ,
[5.03669, 3.42282, -0.0748235] ,
[2.96669, 2.09282, -0.596824] ,
[1.89969, 2.39182, -0.0678235] ,
[3.01469, 1.69682, -1.88382] ,
[1.84869, 1.57082, -2.73482] ,
[2.30269, 1.40782, -4.19882] ,
[0.962686, 0.399824, -2.31982] ,
[-0.255314, 0.522824, -2.22882] ,
[1.57569, -0.758176, -2.00882] ,
[0.900686, -1.92218, -1.46882] ,
[1.91969, -3.07218, -1.33582] ,
[0.223686, -1.65518, -0.123824] ,
[-0.930314, -2.01818, 0.0891765] ,
[0.918686, -0.964176, 0.800176] ,
[0.372686, -0.514176, 2.06618] ,
[1.49069, 0.164824, 2.88418] ,
[-0.814314, 0.435824, 1.92118] ,
[-1.84631, 0.257824, 2.56118] ,
[-0.707314, 1.45482, 1.04518] ,
[-1.78431, 2.37982, 0.745176] ,
[-1.24331, 3.48882, -0.178824] ,
[-3.00931, 1.72082, 0.109176] ,
[-4.14731, 2.02882, 0.455176] ,
[-2.80331, 0.769824, -0.821824] ,
[-3.86731, -0.00917647, -1.42282],
[-3.28731, -0.820176, -2.59882] ,
[-4.57231, -0.946176, -0.440824] ,
[-5.78931, -1.10018, -0.475824] ,
[-3.81431, -1.58718, 0.469176] ,
[-4.35131, -2.44718, 1.50418] ,
[-3.19331, -3.25818, 2.12018] ,
[-5.10031, -1.70118, 2.61018] ,
[-6.15131, -2.13818, 3.07218] ,
[-6.33131, 1.19882, 3.62918] ,
[-7.08531, 1.74482, 4.44018] ,
[-6.45431, 1.44382, 2.30618] ,
[-4.56431, -0.557176, 3.07518] ,
[-5.19131, 0.240824, 4.11318] ,
[-4.10131, 1.05882, 4.83518]], dtype=torch.float64)
# fitting_group = atom_positions
# main_group = atom_positions[[0, 1, 2, 3]]
def get_optimal_rotation(fitting_atoms: torch.tensor, reference_atoms: torch.tensor):
fitting_atoms_centered = fitting_atoms - fitting_atoms.mean(dim=0)
reference_atoms_centered = reference_atoms - reference_atoms.mean(dim=0)
mat_R = torch.linalg.matmul(fitting_atoms_centered.T, reference_atoms_centered)
F00 = mat_R[0][0] + mat_R[1][1] + mat_R[2][2]
F01 = mat_R[1][2] - mat_R[2][1]
F02 = mat_R[2][0] - mat_R[0][2]
F03 = mat_R[0][1] - mat_R[1][0]
F10 = F01
F11 = mat_R[0][0] - mat_R[1][1] - mat_R[2][2]
F12 = mat_R[0][1] + mat_R[1][0]
F13 = mat_R[0][2] + mat_R[2][0]
F20 = F02
F21 = F12
F22 = -mat_R[0][0] + mat_R[1][1] - mat_R[2][2]
F23 = mat_R[1][2] + mat_R[2][1]
F30 = F03
F31 = F13
F32 = F23
F33 = -mat_R[0][0] - mat_R[1][1] + mat_R[2][2]
row0 = torch.stack((F00, F01, F02, F03))
row1 = torch.stack((F10, F11, F12, F13))
row2 = torch.stack((F20, F21, F22, F23))
row3 = torch.stack((F30, F31, F32, F33))
F = torch.stack((row0, row1, row2, row3))
w, v = torch.linalg.eigh(F)
q = v[:, -1]
R00 = q[0] * q[0] + q[1] * q[1] - q[2] * q[2] - q[3] * q[3]
R01 = 2.0 * (q[1] * q[2] - q[0] * q[3])
R02 = 2.0 * (q[1] * q[3] + q[0] * q[2])
R10 = 2.0 * (q[1] * q[2] + q[0] * q[3])
R11 = q[0] * q[0] - q[1] * q[1] + q[2] * q[2] - q[3] * q[3]
R12 = 2.0 * (q[2] * q[3] - q[0] * q[1])
R20 = 2.0 * (q[1] * q[3] - q[0] * q[2])
R21 = 2.0 * (q[2] * q[3] + q[0] * q[1])
R22 = q[0] * q[0] - q[1] * q[1] - q[2] * q[2] + q[3] * q[3]
row0 = torch.stack((R00, R01, R02))
row1 = torch.stack((R10, R11, R12))
row2 = torch.stack((R20, R21, R22))
R = torch.stack((row0, row1, row2))
return w, v, R
def test():
# Using all atoms as the fitting group
w, v, R = get_optimal_rotation(fitting_group, ref_fitting_group)
main_group_centered = main_group - main_group.mean(dim=0)
ref_main_group_centered = ref_main_group - ref_main_group.mean(dim=0)
# Rotate the main group
main_group_r = (R @ main_group_centered.T).T
# main_group_r.register_hook(set_grad(main_group_r))
w2, v2, R2 = get_optimal_rotation(ref_main_group_centered, main_group_r)
q = v2[:, -1]
ref_quat = torch.tensor([1.0, 0.0, 0.0, 0.0], dtype=torch.float64)
q_inner_ref_quat = torch.sum(ref_quat * q)
if (q_inner_ref_quat > 0):
final_q = q
else:
final_q = -1.0 * q
print('Orientation:\n', final_q)
# Force on q
force = torch.tensor([3.40697e-06, -8.24548e-05, -5.71708e-05, 5.95895e-05], dtype=torch.float64)
# Apply force
# torch.sum(force * final_q).backward()
(force[0] * final_q[0]).backward(retain_graph=True)
(force[1] * final_q[1]).backward(retain_graph=True)
(force[2] * final_q[2]).backward(retain_graph=True)
(force[3] * final_q[3]).backward(retain_graph=True)
print('main group force:\n', main_group.grad)
# print('main group (rotated) force:\n', main_group_r.grad)
print('fitting group force:\n', fitting_group.grad)
if __name__ == '__main__':
test()Running the script above I get: Dumping the forces on the main group atoms shows: Fitting group forces: |
|
The latest commit adds the force calculation of the fitting group in It seems I still need to do the same for |
I think I have solved the issue in |
|
Some tests like |
Issue identified by @HanatoK (#713 (comment))
I thought about it for some time, but I think you are right here and opened a PR containing your suggested changes (see #715), so that we can get the tests to run. (Somehow a configuration in your fork seems to prevent it?) Prior to the
That seems at least possible. Besides you guys with Chris, Emad's and Mahmoud's labs have also used the specific combination of orientation in a rotated frame. I'll let the relevant people in each lab know. |
Issue identified by @HanatoK (#713 (comment))
Issue identified by @HanatoK (#713 (comment))
Issue identified by @HanatoK (#713 (comment))
Issue identified by @HanatoK (#713 (comment))
Issue identified by @HanatoK (#713 (comment))
Issue identified by @HanatoK (#713 (comment))
Aha, this takes me back to 2017. Yes, this is probably why Jerome finally implemented independent Euler angles. |
The problem was that the scripted colvars based on orientation never work. They are accurate during equilibration but the forces are wrong when combined with eABF and MtD. |
|
In commit 8bca52f, I have adapted the above python script, which can read the forces on the quaternion, the positions of the main and fitting groups from the .colvars.traj file, calculate the forces on the main and fitting groups by pytorch, and then compare them with the Colvars forces from the debug log output. Running I have no idea about how to implement the debug gradients for forces, and it seems using pytorch is the most straightforward way to verify if the calculation is correct. |
|
The test fails (see https://github.com/HanatoK/colvars/actions/runs/11002205050/job/30548634244#step:11:2287) because the max line length of |
@HanatoK I just pushed an updated version of the "CentOS9-devel" container, which contains spiff from your fork. Please get in touch with @jhenin to merge the changes, and consider adding the spiff repo under the Colvars org? |
The NAMD CI job uses CentOS 7 (https://github.com/HanatoK/colvars/actions/runs/11002205050/job/30801228083#step:1:44). How could I change it to CentOS9-devel? I have made a PR to @jhenin's repo: |
|
The test data in |
We have already had the algorithm for computing the gradients of the fitting group from the gradients of the main group. In cases that the gradients of the main group are not available, we can "intercept" the forces applied to the main group, and then use the same algorithm to convert the forces applied to the main group with the rotation to the forces applied to the fitting group.
…c_fit_forces_impl
This commit introduces the group_force_object that unifies the exertion
of atomic forces in an atom group. If a colvarcomp requires to apply
forces individually on atoms in an atom group, now it can do
```
auto ag_force = atoms->get_group_force_object();
for (size_t ia = 0; ia < atoms->size(); ia++) {
const f_ia = ...;
ag_force.set_atom_force(ia, f_ia);
}
```
If the optimal fitting is enabled for the atom group, group_force_object
will "intercept" the forces applied to the main group, and then calculate
the forces required on the fitting group by calling calc_fit_forces
internally, and then apply the forces on destruction.
This should better match what the existing apply_force does.
This commit dumps the Cartesian coordinates of each frame, which can be read by the compare_with_pytorch.py. compare_with_pytorch.py then calculates the forces on the main and the fitting groups by symbolic gradients, and compares them with the Colvars debug output.
This commit aims to fix the incorrect optimization introduced in commit adebdf6.
For the time being this field is only used for storing the main group forces of non-scalable components with fitting groups.
I agree with you that we should have some approaches to verify the force applied. However, I am not sure how to implement the real fit gradients (not the current propagating "fit forces" way) under the Colvars framework without some fundamental changes.
And also the
Initially I want to do the same thing in this PR, but I wouldn't like to break too many tests, so I chose to only apply forces on the fitting group when |
HanatoK
changed the title
All changes are limited to the copy of the library in `src/external/colvars`. The following is a list of relevant pull requests (bugfixes only) in the library's repository: - 728 Fix undefined behavior when getting the current working directory from std::filesystem Colvars/colvars#728 (@giacomofiorin) - 724 Fix gradients and metric functions of distanceDir Colvars/colvars#724 (@giacomofiorin) - 715 Add missing rotation in orientation component Colvars/colvars#715 (@giacomofiorin) - 713 fix: try to solve #87 for non-scala components Colvars/colvars#713 (@HanatoK) - 706 BUGFIX for Segmentation fault in colvarbias_meta::calc_energy() with useGrids off Colvars/colvars#706 (@alphataubio) - 694 More robust condition to decide when biases run on thread 0 Colvars/colvars#694 (@giacomofiorin) - 675 Fix initialization of histogram output files and move it to the right place Colvars/colvars#675 (@giacomofiorin) Authors: @alphataubio, @giacomofiorin, @HanatoK
Issue identified by @HanatoK (#713 (comment))
- 759 min_image fix; Saves long runs from crashes; Colvars/colvars#759 (@PolyachenkoYA) - 728 Fix undefined behavior when getting the current working directory from std::filesystem Colvars/colvars#728 (@giacomofiorin) - 724 Fix gradients and metric functions of distanceDir Colvars/colvars#724 (@giacomofiorin) - 715 Add missing rotation in orientation component Colvars/colvars#715 (@giacomofiorin) - 713 fix: try to solve lammps#87 for non-scala components Colvars/colvars#713 (@HanatoK) - 706 BUGFIX for Segmentation fault in colvarbias_meta::calc_energy() with useGrids off Colvars/colvars#706 (@alphataubio) - 701 Do not try accessing LAMMPS proxy object before allocating it Colvars/colvars#701 (@giacomofiorin) Authors: @alphataubio, @giacomofiorin, @HanatoK, @PolyachenkoYA
We have already had the algorithm for computing the gradients of the fitting group from the gradients of the main group. In cases that the gradients of the main group are not available, we can "intercept" the forces applied to the atoms of the main group, and then use the same algorithm to convert these forces (with the rotation), to the forces applied to the atoms of the fitting group.
I have tried with a custom test system from NANMA and the applied forces (by setting
COLVARS_DEBUGto true) seems to be correct and consistent with what I get from pytorch autograd, but still need to find a way to verify with the a test example like000_orientation-fitgroup_harmonic-ori-fixed.The python script using pytorch is:
I am not sure if this is the shortest way to fix #87. Anyway, if this approach works there are still a bunch of things to clean up. @giacomofiorin and @jhenin Any ideas?