From ac156a043f72b4976b991d21d35b3e7615b631ae Mon Sep 17 00:00:00 2001
From: Gabriel Studer <gabriel.studer@unibas.ch>
Date: Wed, 22 Nov 2023 13:16:12 +0100
Subject: [PATCH] Revert "avoid atan2(0.0, 0.0) in Dihedral::SetAngleICS"
This reverts commit 9153648f1b773a1c1186d0a0ac29e78b659cd1f8.
The intention of this commit was to fix potential issues when operating on a
dihedral connecting 4 atoms where at least 3 atoms are exactly on a line.
Dihedrals are ill defined in this case but "expected behaviour"
(whatever that means) is tested as a unit test.
A fail was observed in two cases. First case that triggered commit
9153648f1b773a1c1186d0a0ac29e78b659cd1f8, second case that showed that the
problem is not solved. Therefore reverting to original state in order to find
different solution.
---
modules/mol/base/src/impl/dihedral.cc | 41 +++++++++++++--------------
1 file changed, 19 insertions(+), 22 deletions(-)
diff --git a/modules/mol/base/src/impl/dihedral.cc b/modules/mol/base/src/impl/dihedral.cc
index 2c8c003c4..ee26230be 100644
--- a/modules/mol/base/src/impl/dihedral.cc
+++ b/modules/mol/base/src/impl/dihedral.cc
@@ -125,29 +125,26 @@ void Dihedral::SetAngleICS(Real angle, bool update_other) {
// not all connectors are living in the same coordinate system. We first
// have to bring all of them into the local coordinate system of a3.
if (a2->GetPrimaryConnector() && a2->GetPrimaryConnector()->GetFirst()==a1) {
- // equivalent to geom::Transpose(conn2->GetLocalRot()) * geom::Vec3(0,0,-1)
- v1=geom::Vec3(0, 0,-1)*conn2->GetLocalRot();
- ConnectorImplP c=GetConnector(a4, a3);
- v2=c->GetDir();
- Real phi1 = 0.0;
- if(std::abs(v1[1]) + std::abs(v1[0]) > Real(1e-6)) phi1=atan2(v1[1], v1[0]);
- geom::Vec3 n2=vec_for_angle(v2, angle+phi1);
- c->SetDir(n2);
- // check if we have to update the other connectors
- const ConnectorImplList& sc=a3->GetSecondaryConnectors();
- if (update_other && sc.size()>1) {
- Real phi2=(angle-phi1);
- if(std::abs(v2[1]) + std::abs(v2[0]) > Real(1e-6)) phi2-=atan2(v2[1], v2[0]);
- for (ConnectorImplList::const_iterator i=sc.begin(); i!=sc.end(); ++i) {
- if (*i==c)
- continue;
- v3=(*i)->GetDir();
- Real phi3=phi2;
- if(std::abs(v3[1]) + std::abs(v3[0]) > Real(1e-6)) phi3+=atan2(v3[1], v3[0]);
- geom::Vec3 n=vec_for_angle(v3, phi3);
- (*i)->SetDir(n);
+ v1=geom::Vec3(0, 0,-1)*conn2->GetLocalRot();
+ v2=geom::Vec3(0, 0, 1);
+ ConnectorImplP c=GetConnector(a4, a3);
+ v3=c->GetDir();
+ Real phi1=atan2(v1[1], v1[0]);
+ geom::Vec3 n3=vec_for_angle(v3, angle+phi1);
+ c->SetDir(n3);
+ // check if we have to update the other connectors
+ const ConnectorImplList& sc=a3->GetSecondaryConnectors();
+ if (update_other && sc.size()>1) {
+ Real phi2=(angle-phi1)-atan2(v3[1], v3[0]);
+ for (ConnectorImplList::const_iterator i=sc.begin(); i!=sc.end(); ++i) {
+ if (*i==c)
+ continue;
+ geom::Vec3 v=(*i)->GetDir();
+ Real phi4=phi2+atan2(v[1], v[0]);
+ geom::Vec3 n=vec_for_angle(v, phi4);
+ (*i)->SetDir(n);
+ }
}
- }
}
}
--
GitLab