I'd be surprised if superimposed atoms are a problem because of the
energy penalties: as the gradients are equal for the two atoms, they
never get split, if the optimisation engine doesn't randomise things a
little.
I'd think the GetHyb() values should probably match the values UFF
uses for coordination, as you propose: 5, 6, 7 for PF5, SF6, IF7
respectively.
I think the "sad teacher" standard you hint at is probably good: deal
with examples high school teachers might want to demonstrate to their
students and be unhappy about if they failed, ignore the rest as too
exotic.
So these are the cases I'm thinking about:
- SO3 should probably be planar, and a special-case if() in
forcefielduff.c: GetCoordination() seems to Just Work.
- PF5, SF6 work with something like the patch above
- S2F10, as built by builder.cpp, post-patch, has the F-S-S-F
dihedral as multiples of 90° rather than 45°, so it lacks the
45°-twist. I'm not sure whether that's considered part of VSEPR, but
I'm tempted to say if you want such positioning, use UFF (Avogadro
always does that. Does anyone use builder.cpp without optimisation?).
- IF7. This relies on code in forcefielduff.cpp that's currently
commented out, and no other force field results in a pentagonal
bipyramid. If I read the comments in forcefielduff.cpp correctly, we
can hope to get this working if builder.cpp chooses the right
directions, and UFF optimisation then fixes bond length, but it might
be too difficult to go from randomly positioned atoms to the right
structure. Experimentation seems to confirm this, though I feel
uncomfortable with the polynomial in cosT I've chosen for the
pentagonal structure (symmetries of 2, 3, and 4 result in pretty
integer coefficients, while cos(2 * M_PI * 1.0/5.0) doesn't)
- XeF4 gets built tetrahedally. UFF seems to have a minimum both for
the (correct, I think) planar structure and the (incorrect, I think)
structure with VX, VY, VZ, and -VX bond directions.
- XeOF4 fails to have the Xe and F all in the same plane both with
and without optimization. Presumably this needs the GetHyb() == 4
case.
- ReH9--
- A hypothetical P2F8 or I2F12 molecule would run into problems: the
UFF_AXIAL_ATOM and UFF_CENTRAL_ATOM tags are per-atom, but in these
examples, both atoms would get both tags, and the code gets confused.
I'd draw the "sad teacher" line between IF7 and XeF4.
So, sorry if this is a bit confusing. I think my current status is:
builder: PF5, SF6 work with patch; IF7 doesn't yet work; S2F10
semi-works (no twist); P2F8 is linear, but no twist; I2F12 doesn't
work.
UFF: PF5, SF6 work; IF7 works with patch; S2F10 works; *P2F8 isn't linear.
I'll try getting IF7 to work in builder, and then try to send you a
clean patch that's applicable in full, if you're okay with that
strategy.
I'm not sure anything like the P2F8/I2F12 case actually exists. I
suppose it's possible to set UFF_AXIAL_ATOM on the bond rather than
the atom - can the same bond be axial wrt its begin atom and
equatorial wrt its end atom?
On Wed, Mar 17, 2010 at 1:27 AM, Geoffrey Hutchison
<[email protected]> wrote:
>> such as "FS(F)(F)(F)(F)F": two of the fluorine atoms are superposed,
> ...
>> (as well as a missed opportunity for a warning in avogadro): the
>
> I think there are a few things going on here. I thought I got all of these
> bugs, but you clearly found a few more. ;-) In Open Babel, the force field
> implementations need a better penalty on superimposed atoms, particularly for
> the angular terms. This is particularly true for UFF, where people are more
> likely to impose strange valences.
>
>> Does it make sense to have new GetHyb() values, and if so, are 4 and 5
>> the right choices?
>
> Yes, it does make sense to have new hybridization values above 3. However,
> I'd use "5" for pentavalent (e.g., PF5), and "6" for octahedral, with "4"
> reserved for square planar (i.e., Pt, Pd, Cu, etc.). It gets harder above
> this, since heavy elements like U can take very high valences. My guess is
> for these cases, it's probably OK for GetHyb() to return 0 as "unknown" and
> position atoms randomly. If someone cares about chiral uranium complexes,
> they probably already have the 3D coordinates for it.
>
> I'll admit that hybridization above 3 has long been on my todo list, so I'm a
> bit embarrassed but also happy that you got it first.
>
>> Are there hypervalent molecules with double bonds
>> which are not tetrahedral in shape (my understanding is H3PO4 is
>> tetrahedral around the P atom)? (Again going on wikipedia, sulfur
>> trioxide might be planar, but the UFF optimisation disagrees. This is
>> potentially a bug in the UFF implementation).
>
> Xenon isn't tetrahedral and XeOF4 is a semi-common VSEPR example. There are
> probably some octahedral examples (e.g. SOF4 perhaps?).
>
> I'll have to check the SO3 bug you filed. To be honest, classic UFF doesn't
> do well with hypervalent molecules, but I knew teachers would want the basic
> VSEPR rules enforced. So the Open Babel implementation has some heuristics to
> override the UFF typing rules. There's also some code in there which attempts
> to do pentavalent species correctly.
>
> If you're willing to help test out some of the common and less-common VSEPR
> examples using Avogadro and OB, I'm definitely interested in cleaning up any
> remaining bugs.
>
> I guess one other issue is that there's still code in
> atom.cpp:OBAtom::GetNewBondVector that should be switched to use the new
> OBBuilder code. This would improve issues with implicit hydrogen atoms and
> hypervalent molecules.
>
> Thanks very much,
> -Geoff
Index: src/builder.cpp
===================================================================
--- src/builder.cpp (revision 3642)
+++ src/builder.cpp (working copy)
@@ -143,11 +143,13 @@
vector3 OBBuilder::GetNewBondVector(OBAtom *atom, double length)
{
- vector3 bond1, bond2, bond3, v1, v2, newbond;
+ vector3 bond1, bond2, bond3, bond4, bond5, v1, v2, newbond;
bond1 = VZero;
bond2 = VZero;
bond3 = VZero;
+ bond4 = VZero;
+ bond5 = VZero;
if (atom == NULL)
return VZero;
@@ -226,222 +228,377 @@
newbond = bond1 - v2 * tan(DEG_TO_RAD*60.0);
else if (atom->GetHyb() == 3)
newbond = bond1 - v2 * tan(DEG_TO_RAD*70.5);
-
+ else if (atom->GetHyb() == 4) {
+ /* the first two atoms are the axial ones; the third, fourth, and fifth atom are equatorial */
+ newbond = bond1;
+ } else if (atom->GetHyb() == 5) {
+ newbond = bond1 - v2 * tan(DEG_TO_RAD*90.0);
+ }
+
newbond = newbond.normalize();
newbond *= length;
newbond += atom->GetVector();
return newbond;
}
-
- //
- // \ \
- // X ---> X--*
- // / /
- //
- if (atom->GetValence() == 2) {
- FOR_NBORS_OF_ATOM (nbr, atom) {
- if (bond1 == VZero)
- bond1 = atom->GetVector() - nbr->GetVector();
- else
- bond2 = atom->GetVector() - nbr->GetVector();
- }
+
+ //
+ // \ \
+ // X ---> X--*
+ // / /
+ //
+ if (atom->GetValence() == 2) {
+ FOR_NBORS_OF_ATOM (nbr, atom) {
+ if (bond1 == VZero)
+ bond1 = atom->GetVector() - nbr->GetVector();
+ else if (bond2 == VZero)
+ bond2 = atom->GetVector() - nbr->GetVector();
+ }
- bond1 = bond1.normalize();
- bond2 = bond2.normalize();
- v1 = bond1 + bond2;
- v1 = v1.normalize();
-
- if (atom->GetHyb() == 2)
- newbond = v1;
- if (atom->GetHyb() == 3) {
- v2 = cross(bond1, bond2);
- //v1 = v1.normalize();
- newbond = v2 + v1 * tan(DEG_TO_RAD*35.25);
+ bond1 = bond1.normalize();
+ bond2 = bond2.normalize();
+ v1 = bond1 + bond2;
+ v1 = v1.normalize();
+
+ if (atom->GetHyb() == 2)
+ newbond = v1;
+ if (atom->GetHyb() == 3) {
+ v2 = cross(bond1, bond2);
+ //v1 = v1.normalize();
+ newbond = v2 + v1 * tan(DEG_TO_RAD*35.25);
+ }
+ if (atom->GetHyb() == 4) {
+ /* add the first equatorial atom, orthogonally to bond1 (and bond2 = -bond1) */
+ /* is atom order correct? I don't think it matters, but I might have to ask a chemist
+ * whether PClF4 would be more likely to have an equatorial or axial Cl-P bond */
+ v1 = cross(bond1, VY);
+ if (v1 == VZero) // This corner-case happened to me, where bond1 was -VY (Noel)
+ v1 = cross(bond1, VX);
+ v1 = v1.normalize();
+ newbond = v1;
+ }
+ if (atom->GetHyb() == 5) {
+ v2 = cross(bond1, bond2);
+ newbond = v2;
+ }
+
+ newbond = newbond.normalize();
+ newbond *= length;
+ newbond += atom->GetVector();
+ return newbond;
}
- newbond = newbond.normalize();
- newbond *= length;
- newbond += atom->GetVector();
- return newbond;
- }
-
- //
- // \ \
- // --X ---> --X--*
- // / /
- //
- if (atom->GetValence() == 3) {
- FOR_NBORS_OF_ATOM (nbr, atom) {
- if (bond1 == VZero)
- bond1 = atom->GetVector() - nbr->GetVector();
- else if (bond2 == VZero)
- bond2 = atom->GetVector() - nbr->GetVector();
- else
- bond3 = atom->GetVector() - nbr->GetVector();
+ /* UFF:
+ * b lg dg o y
+ * b - 45 30 45 30
+ * lg - 45 0 45
+ * dg - 45 30
+ * o - 45
+ * y -
+
+ * 94s:
+ * b lg dg o y
+ * b - 34 34 34 34
+ * lg - 48 21 48
+ * dg - 48 21
+ * o - 48
+ * y -
+
+ //
+ // \ \
+ // --X ---> --X--*
+ // / /
+ //
+ */
+ if (atom->GetValence() == 3) {
+ if (atom->GetHyb() == 3) {
+ FOR_NBORS_OF_ATOM (nbr, atom) {
+ if (bond1 == VZero)
+ bond1 = atom->GetVector() - nbr->GetVector();
+ else if (bond2 == VZero)
+ bond2 = atom->GetVector() - nbr->GetVector();
+ else if (bond3 == VZero)
+ bond3 = atom->GetVector() - nbr->GetVector();
+ }
+
+ bond1 = bond1.normalize();
+ bond2 = bond2.normalize();
+ bond3 = bond3.normalize();
+ newbond = bond1 + bond2 + bond3;
+ newbond = newbond.normalize();
+ newbond *= length;
+ newbond += atom->GetVector();
+ return newbond;
+ }
+
+ if (atom->GetHyb() == 4) {
+ FOR_NBORS_OF_ATOM (nbr, atom) {
+ if (bond1 == VZero)
+ bond1 = atom->GetVector() - nbr->GetVector();
+ else if (bond2 == VZero)
+ bond2 = atom->GetVector() - nbr->GetVector();
+ else if (bond3 == VZero)
+ bond3 = atom->GetVector() - nbr->GetVector();
+ }
+
+ bond1 = bond1.normalize();
+ bond2 = bond2.normalize();
+ bond3 = bond3.normalize();
+
+ v1 = cross(bond1, bond3);
+ v1 = v1.normalize();
+
+ newbond = v1 + tan(DEG_TO_RAD*30.0) * bond3;
+ newbond = newbond.normalize();
+ newbond *= length;
+ newbond += atom->GetVector();
+ return newbond;
+
+ }
+
+ if (atom->GetHyb() == 5) {
+ FOR_NBORS_OF_ATOM (nbr, atom) {
+ if (bond1 == VZero)
+ bond1 = atom->GetVector() - nbr->GetVector();
+ else if (bond2 == VZero)
+ bond2 = atom->GetVector() - nbr->GetVector();
+ else if (bond3 == VZero)
+ bond3 = atom->GetVector() - nbr->GetVector();
+ }
+
+ /* the way things work, newbond is equal to bond1, but will show up at -bond1 next time around */
+ newbond = bond1;
+ newbond = newbond.normalize();
+ newbond *= length;
+ newbond += atom->GetVector();
+ return newbond;
+ }
}
-
- bond1 = bond1.normalize();
- bond2 = bond2.normalize();
- bond3 = bond3.normalize();
- newbond = bond1 + bond2 + bond3;
- newbond = newbond.normalize();
- newbond *= length;
- newbond += atom->GetVector();
- return newbond;
- }
- } else {
- ////////////
- // 2D //
- ////////////
- OBBondIterator i;
+ if (atom->GetValence() == 4) {
+ if (atom->GetHyb() == 5) {
+ FOR_NBORS_OF_ATOM (nbr, atom) {
+ if (bond1 == VZero)
+ bond1 = atom->GetVector() - nbr->GetVector();
+ else if (bond2 == VZero)
+ bond2 = atom->GetVector() - nbr->GetVector();
+ else if (bond3 == VZero)
+ bond3 = atom->GetVector() - nbr->GetVector();
+ else if (bond4 == VZero)
+ bond4 = atom->GetVector() - nbr->GetVector();
+ }
- //
- // a ---> a---*
- //
- if (atom->GetValence() == 0) {
- newbond = atom->GetVector() + VX * length;
- // Check that the vector is still finite before returning
- if (!isfinite(newbond.x()) || !isfinite(newbond.y()))
- newbond.Set(0.0, 0.0, 0.0);
- return newbond;
- }
+ newbond = bond2;
+ newbond = newbond.normalize();
+ newbond *= length;
+ newbond += atom->GetVector();
+ return newbond;
+ }
- // hyb * = 1 //
- // ^^^^^^^^^ //
- // //
- // (a-1)--a ---> (a-1)--a--* angle(a-1, a, *) = 180 //
- // //
- // hyb * = 2 //
- // ^^^^^^^^^ //
- // make sure we place the new atom trans to a-2 (if there is an a-2 atom) //
- // //
- // (a-2) (a-2) //
- // \ \ //
- // (a-1)==a ---> (a-1)==a angle(a-1, a, *) = 120 //
- // \ //
- // * //
- // hyb * = 3 //
- // ^^^^^^^^^ //
- // make sure we place the new atom trans to a-2 (if there is an a-2 atom) //
- // //
- // (a-2) (a-2) //
- // \ \ //
- // (a-1)--a ---> (a-1)--a angle(a-1, a, *) = 109 //
- // \ //
- // * //
- if (atom->GetValence() == 1) {
- OBAtom *nbr = atom->BeginNbrAtom(i);
- if (!nbr)
- return VZero;
- bond1 = atom->GetVector() - nbr->GetVector(); // bond (a-1)--a
+ if (atom->GetHyb() == 4) {
+ FOR_NBORS_OF_ATOM (nbr, atom) {
+ if (bond1 == VZero)
+ bond1 = atom->GetVector() - nbr->GetVector();
+ else if (bond2 == VZero)
+ bond2 = atom->GetVector() - nbr->GetVector();
+ else if (bond3 == VZero)
+ bond3 = atom->GetVector() - nbr->GetVector();
+ else if (bond4 == VZero)
+ bond4 = atom->GetVector() - nbr->GetVector();
+ }
+
+ bond1 = bond1.normalize();
+ bond2 = bond2.normalize();
+ bond3 = bond3.normalize();
+ bond4 = bond4.normalize();
- for (OBAtom *nbr2 = nbr->BeginNbrAtom(i); nbr2; nbr2 = nbr->NextNbrAtom(i))
- if (nbr2 != atom)
- bond2 = nbr->GetVector() - nbr2->GetVector(); // bond (a-2)--(a-1)
+ v1 = cross(bond1, bond3);
+ v1 = v1.normalize();
- int hyb = atom->GetHyb();
- if (hyb == 1)
- newbond = bond1;
- else if (hyb == 2 || hyb == 3) {
- matrix3x3 m;
- m.RotAboutAxisByAngle(VZ, 60.0);
- newbond = m*bond1;
+ newbond = -v1 + tan(DEG_TO_RAD*30.0) * bond3;
+ newbond = newbond.normalize();
+ newbond *= length;
+ newbond += atom->GetVector();
+ return newbond;
+
+ }
+
}
- newbond.normalize();
- newbond *= length;
- newbond += atom->GetVector();
- return newbond;
- } // GetValence() == 1
+
+ if (atom->GetValence() == 5) {
+ if (atom->GetHyb() == 5) {
+ FOR_NBORS_OF_ATOM (nbr, atom) {
+ if (bond1 == VZero)
+ bond1 = atom->GetVector() - nbr->GetVector();
+ else if (bond2 == VZero)
+ bond2 = atom->GetVector() - nbr->GetVector();
+ else if (bond3 == VZero)
+ bond3 = atom->GetVector() - nbr->GetVector();
+ else if (bond4 == VZero)
+ bond4 = atom->GetVector() - nbr->GetVector();
+ else if (bond5 == VZero)
+ bond5 = atom->GetVector() - nbr->GetVector();
+ }
+
+ newbond = bond3;
+ newbond = newbond.normalize();
+ newbond *= length;
+ newbond += atom->GetVector();
+ return newbond;
+ }
+ }
+ } else {
+ ////////////
+ // 2D //
+ ////////////
+ OBBondIterator i;
+
+ //
+ // a ---> a---*
+ //
+ if (atom->GetValence() == 0) {
+ newbond = atom->GetVector() + VX * length;
+ // Check that the vector is still finite before returning
+ if (!isfinite(newbond.x()) || !isfinite(newbond.y()))
+ newbond.Set(0.0, 0.0, 0.0);
+ return newbond;
+ }
+
+ // hyb * = 1 //
+ // ^^^^^^^^^ //
+ // //
+ // (a-1)--a ---> (a-1)--a--* angle(a-1, a, *) = 180 //
+ // //
+ // hyb * = 2 //
+ // ^^^^^^^^^ //
+ // make sure we place the new atom trans to a-2 (if there is an a-2 atom) //
+ // //
+ // (a-2) (a-2) //
+ // \ \ //
+ // (a-1)==a ---> (a-1)==a angle(a-1, a, *) = 120 //
+ // \ //
+ // * //
+ // hyb * = 3 //
+ // ^^^^^^^^^ //
+ // make sure we place the new atom trans to a-2 (if there is an a-2 atom) //
+ // //
+ // (a-2) (a-2) //
+ // \ \ //
+ // (a-1)--a ---> (a-1)--a angle(a-1, a, *) = 109 //
+ // \ //
+ // * //
+ if (atom->GetValence() == 1) {
+ OBAtom *nbr = atom->BeginNbrAtom(i);
+ if (!nbr)
+ return VZero;
+ bond1 = atom->GetVector() - nbr->GetVector(); // bond (a-1)--a
+
+ for (OBAtom *nbr2 = nbr->BeginNbrAtom(i); nbr2; nbr2 = nbr->NextNbrAtom(i))
+ if (nbr2 != atom)
+ bond2 = nbr->GetVector() - nbr2->GetVector(); // bond (a-2)--(a-1)
+
+ int hyb = atom->GetHyb();
+ if (hyb == 1)
+ newbond = bond1;
+ else if (hyb == 2 || hyb == 3 || hyb == 0) {
+ matrix3x3 m;
+ m.RotAboutAxisByAngle(VZ, 60.0);
+ newbond = m*bond1;
+ }
+ newbond.normalize();
+ newbond *= length;
+ newbond += atom->GetVector();
+ return newbond;
+ } // GetValence() == 1
// //
// \ \ //
// X ---> X--* //
// / / //
// //
- if (atom->GetValence() == 2) {
- for (OBAtom *nbr = atom->BeginNbrAtom(i); nbr; nbr = atom->NextNbrAtom(i)) {
- if (bond1 == VZero)
- bond1 = atom->GetVector() - nbr->GetVector();
- else
- bond2 = atom->GetVector() - nbr->GetVector();
+ if (atom->GetValence() == 2) {
+ for (OBAtom *nbr = atom->BeginNbrAtom(i); nbr; nbr = atom->NextNbrAtom(i)) {
+ if (bond1 == VZero)
+ bond1 = atom->GetVector() - nbr->GetVector();
+ else
+ bond2 = atom->GetVector() - nbr->GetVector();
+ }
+ bond1.normalize();
+ bond2.normalize();
+ newbond = bond1 + bond2;
+ newbond.normalize();
+ newbond *= length;
+ newbond += atom->GetVector();
+ return newbond;
}
- bond1.normalize();
- bond2.normalize();
- newbond = bond1 + bond2;
- newbond.normalize();
- newbond *= length;
- newbond += atom->GetVector();
- return newbond;
- }
- // //
- // \ \ //
- // --X ---> --X--* //
- // / / //
- // //
- if (atom->GetValence() == 3) {
- OBStereoFacade stereoFacade((OBMol*)atom->GetParent());
- if (stereoFacade.HasTetrahedralStereo(atom->GetId())) {
- OBBond *hash = 0;
- OBBond *wedge = 0;
- vector<OBBond*> plane;
- for (OBAtom *nbr = atom->BeginNbrAtom(i); nbr; nbr = atom->NextNbrAtom(i)) {
- OBBond *bond = atom->GetBond(nbr);
-
- if (bond->IsWedge()) {
- if (atom == bond->GetBeginAtom())
- wedge = bond;
- else
- hash = bond;
- } else
- if (bond->IsHash()) {
- if (atom == bond->GetBeginAtom())
- hash = bond;
- else
- wedge = bond;
- } else
- plane.push_back(bond);
- }
-
- if (wedge && !plane.empty()) {
- bond2 = atom->GetVector() - wedge->GetNbrAtom(atom)->GetVector();
- bond3 = atom->GetVector() - plane[0]->GetNbrAtom(atom)->GetVector();
- } else if (hash && !plane.empty()) {
- bond2 = atom->GetVector() - hash->GetNbrAtom(atom)->GetVector();
- bond3 = atom->GetVector() - plane[0]->GetNbrAtom(atom)->GetVector();
- } else if (plane.size() >= 2) {
- bond2 = atom->GetVector() - plane[0]->GetNbrAtom(atom)->GetVector();
- bond3 = atom->GetVector() - plane[1]->GetNbrAtom(atom)->GetVector();
- } else if (hash && wedge) {
- bond2 = atom->GetVector() - wedge->GetNbrAtom(atom)->GetVector();
- bond3 = atom->GetVector() - hash->GetNbrAtom(atom)->GetVector();
- }
- } else {
- for (OBAtom *nbr = atom->BeginNbrAtom(i); nbr; nbr = atom->NextNbrAtom(i)) {
- if (bond1 == VZero)
- bond1 = atom->GetVector() - nbr->GetVector();
- else if (bond2 == VZero)
- bond2 = atom->GetVector() - nbr->GetVector();
- else
- bond3 = atom->GetVector() - nbr->GetVector();
- }
+ // //
+ // \ \ //
+ // --X ---> --X--* //
+ // / / //
+ // //
+ if (atom->GetValence() == 3) {
+ OBStereoFacade stereoFacade((OBMol*)atom->GetParent());
+ if (stereoFacade.HasTetrahedralStereo(atom->GetId())) {
+ OBBond *hash = 0;
+ OBBond *wedge = 0;
+ vector<OBBond*> plane;
+ for (OBAtom *nbr = atom->BeginNbrAtom(i); nbr; nbr = atom->NextNbrAtom(i)) {
+ OBBond *bond = atom->GetBond(nbr);
+
+ if (bond->IsWedge()) {
+ if (atom == bond->GetBeginAtom())
+ wedge = bond;
+ else
+ hash = bond;
+ } else
+ if (bond->IsHash()) {
+ if (atom == bond->GetBeginAtom())
+ hash = bond;
+ else
+ wedge = bond;
+ } else
+ plane.push_back(bond);
+ }
+
+ if (wedge && !plane.empty()) {
+ bond2 = atom->GetVector() - wedge->GetNbrAtom(atom)->GetVector();
+ bond3 = atom->GetVector() - plane[0]->GetNbrAtom(atom)->GetVector();
+ } else if (hash && !plane.empty()) {
+ bond2 = atom->GetVector() - hash->GetNbrAtom(atom)->GetVector();
+ bond3 = atom->GetVector() - plane[0]->GetNbrAtom(atom)->GetVector();
+ } else if (plane.size() >= 2) {
+ bond2 = atom->GetVector() - plane[0]->GetNbrAtom(atom)->GetVector();
+ bond3 = atom->GetVector() - plane[1]->GetNbrAtom(atom)->GetVector();
+ } else if (hash && wedge) {
+ bond2 = atom->GetVector() - wedge->GetNbrAtom(atom)->GetVector();
+ bond3 = atom->GetVector() - hash->GetNbrAtom(atom)->GetVector();
+ }
+ } else {
+ for (OBAtom *nbr = atom->BeginNbrAtom(i); nbr; nbr = atom->NextNbrAtom(i)) {
+ if (bond1 == VZero)
+ bond1 = atom->GetVector() - nbr->GetVector();
+ else if (bond2 == VZero)
+ bond2 = atom->GetVector() - nbr->GetVector();
+ else
+ bond3 = atom->GetVector() - nbr->GetVector();
+ }
+ }
+
+ bond2.normalize();
+ bond3.normalize();
+ newbond = -(bond2 + bond3);
+ newbond.normalize();
+ newbond *= length;
+ newbond += atom->GetVector();
+ return newbond;
}
+
+ }
- bond2.normalize();
- bond3.normalize();
- newbond = -(bond2 + bond3);
- newbond.normalize();
- newbond *= length;
- newbond += atom->GetVector();
- return newbond;
- }
-
+ return VZero; //previously undefined
}
- return VZero; //previously undefined
-}
// The OBMol mol contains both the molecule to which we want to connect the
// fragment and the fragment itself. The fragment containing b will be
@@ -722,7 +879,6 @@
// }
}
-
// First we find the most complex fragments in our molecule. Once we have a match,
// vfrag is set for all the atoms in the fragment. A second match (smaller, more
// simple part of the 1st match) is ignored.
@@ -823,7 +979,7 @@
// set coordinates for atoms
OBAtom *atom = workMol.GetAtom(*k);
atom->SetVector(i->second[counter]);
- }
+ }
}
// add the bonds for the fragment
@@ -892,7 +1048,7 @@
vdone.SetBitOn(a->GetIdx());
- // place the atom the atom
+ // place the atom
OBAtom *atom = workMol.GetAtom(a->GetIdx());
atom->SetVector(molvec);
Index: src/forcefields/forcefielduff.cpp
===================================================================
--- src/forcefields/forcefielduff.cpp (revision 3642)
+++ src/forcefields/forcefielduff.cpp (working copy)
@@ -144,6 +144,14 @@
// Rappe form: (1 - cos 4*theta) -- minima at 0, 360 (bad...)
energy = ka * (1.0 + cosT)*cosT*cosT;
break;
+ case 7: // IF7.
+ /* theta = 1/5 * 2 pi. cosT = .30901699
+ * theta = 2/5 * 2 pi. cosT = -.80901699
+ * theta = 3/5 * 2 pi. cosT = -.80901699
+ * theta = 4/5 * 2 pi. cosT = .30901699
+ */
+ energy = ka * c1 * (cosT - .30901699) * (cosT - .30906199) * (cosT + .80901699) * (cosT + .8091699);
+ break;
default: // general (sp3) coordination
energy = ka*(c0 + c1*cosT + c2*(2.0*cosT*cosT - 1.0)); // use cos 2t = (2cos^2 - 1)
}
@@ -162,6 +170,12 @@
case 6: // octahedral
dE = -ka * cosT * (2.0 + 3.0 * cosT) * sinT;
break;
+ case 7: // pentagonal bipyramidal
+ dE =
+ c1 * -ka * (2 * sinT * (cosT - .30906199) * (cosT + .80901699) * (cosT + .8091699) +
+ 2 * sinT * (cosT - .30901699) * (cosT - .30906199) * (cosT + .8091699));
+ //dE = -ka * c1 * sin(5*theta) * 5;
+ break;
default: // general (sp3) coordination
dE = -ka * (c1*sinT + 2.0 * c2*sin(2.0 * theta));
}
@@ -693,15 +707,25 @@
// PR#2971473, thanks to Philipp Rumpf <[email protected]>
coordination = 2; // i.e., sp2
}
+ /* planar coordination of hexavalent molecules.*/
+ if (lonePairs == 0 && b->GetValence() == 3 && b->BOSum() == 6) {
+ coordination = 2;
+ }
+ if (lonePairs == 0 && b->GetValence() == 7) {
+ coordination = 7;
+ }
+ // Check to see if coordination is really correct
+ // if not (e.g., 5- or 7- or 8-coord...)
+ // then create approximate angle bending terms
} else {
coordination = ipar; // coordination of central atom
// Check to see if coordination is really correct
// if not (e.g., 5- or 7- or 8-coord...)
// then create approximate angle bending terms
- if (b->GetValence() > 4) {
- coordination = b->GetValence();
- }
}
+ if (b->GetValence() > 4) {
+ coordination = b->GetValence();
+ }
return coordination;
}
@@ -777,6 +801,46 @@
secondLargestNbr->SetData(label);
} // end work for 5-coordinate angles
+ if (GetCoordination(&*atom, parameterB->_ipar[0]) == 7) {
+ // First, find the two largest neighbors
+ OBAtom *largestNbr, *current, *secondLargestNbr = 0;
+ double largestRadius;
+ OBBondIterator i;
+ largestNbr = atom->BeginNbrAtom(i);
+ // work out the radius
+ parameterA = GetParameterUFF(largestNbr->GetType(), _ffparams);
+ largestRadius = parameterA->_dpar[0];
+
+ for (current = atom->NextNbrAtom(i); current; current = atom->NextNbrAtom(i)) {
+ parameterA = GetParameterUFF(current->GetType(), _ffparams);
+ if (parameterA->_dpar[0] > largestRadius) {
+ // New largest neighbor
+ secondLargestNbr = largestNbr;
+ largestRadius = parameterA->_dpar[0];
+ largestNbr = current;
+ }
+ if (secondLargestNbr == 0) {
+ // save this atom
+ secondLargestNbr = current;
+ }
+ }
+
+ // OK, now we tag the central atom
+ OBPairData *label = new OBPairData;
+ label->SetAttribute("UFF_CENTRAL_ATOM");
+ label->SetValue("True"); // doesn't really matter
+ atom->SetData(label);
+ // And tag the axial substituents
+ label = new OBPairData;
+ label->SetAttribute("UFF_AXIAL_ATOM");
+ label->SetValue("True");
+ largestNbr->SetData(label);
+ label = new OBPairData;
+ label->SetAttribute("UFF_AXIAL_ATOM");
+ label->SetValue("True");
+ secondLargestNbr->SetData(label);
+
+ }
} // end loop through atoms
//
@@ -898,10 +962,11 @@
}
//double currentTheta;
- if (coordination >= 7) {
+ if (coordination > 7) {
// large coordination sphere (e.g., [ReH9]-2 or [Ce(NO3)6]-2)
// just resort to using VDW 1-3 interactions to push atoms into place
// there's not much else we can do without real parameters
+ // wait, what? We can handle IF7 fine.
if (SetupVDWCalculation(a, c, vdwcalc)) {
_vdwcalculations.push_back(vdwcalc);
}
@@ -910,8 +975,64 @@
// The downside is that we can't easily handle lone pairs.
continue;
+ } else if (coordination == 7) { // pentagonal bipyramidal
+ // This section is commented out.
+ // If you want to try to tackle 7-coordinate species, give this a try
+ // and change the first conditional above to >7, not >=7
+ // This doesn't work so well because it's hard to classify between
+ // axial-equatorial (90 degrees) and proximal equatorial (~72 degrees).
+ // You'll probably have to do something like the pre-evaluation of 5-coordinate atoms at the top of this method (this is what I did)
+ //
+ // } else if (coordination == 7) { // pentagonal bipyramidal
+ double currentTheta;
+ currentTheta = a->GetAngle(&*b, &*c);
+
+ anglecalc.c0 = 1.0;
+ if (0) {
+ if (currentTheta >= 155.0) { // axial ligands = linear
+ anglecalc.coord = 1; // like sp
+ anglecalc.theta0 = 180.0;
+ anglecalc.c1 = 1.0;
+ } else if (currentTheta < 155.0 && currentTheta >= 110.0) { // distal equatorial
+ anglecalc.coord = 7; // like sp3
+ anglecalc.theta0 = 144.0;
+ anglecalc.c1 = 1.0;
+ } else if (currentTheta < 110.0 && currentTheta >= 85.0) { // axial-equatorial
+ anglecalc.coord = 4; // like sq. planar or octahedral
+ anglecalc.theta0 = 90.0;
+ anglecalc.c1 = 1.0;
+ } else if (currentTheta < 85.0) { // proximal equatorial
+ anglecalc.coord = 7; // general case (i.e., like sp3)
+ anglecalc.theta0 = 72.0;
+ anglecalc.c1 = 1.0;
+ }
+ anglecalc.c2 = 0.0;
+
+ // // Also add a VDW 1-3 interaction to distort slightly
+ // if (SetupVDWCalculation(a, c, vdwcalc)) {
+ // _vdwcalculations.push_back(vdwcalc);
+ // }
+ } else {
+ if (b->HasData("UFF_CENTRAL_ATOM")
+ && a->HasData("UFF_AXIAL_ATOM")
+ && c->HasData("UFF_AXIAL_ATOM")) { // axial ligands = linear
+ anglecalc.coord = 1; // like sp
+ anglecalc.theta0 = 180.0;
+ anglecalc.c1 = 1.0;
+ } else if ( (a->HasData("UFF_AXIAL_ATOM") && !c->HasData("UFF_AXIAL_ATOM"))
+ || (c->HasData("UFF_AXIAL_ATOM") && !a->HasData("UFF_AXIAL_ATOM")) ) { // axial-equatorial ligands
+ anglecalc.coord = 4; // like sq. planar or octahedral
+ anglecalc.theta0 = 90.0;
+ anglecalc.c1 = 1.0;
+ } else { // equatorial - equatorial
+ anglecalc.coord = 7; // unlike anything else, as theta0 is ignored.
+ anglecalc.theta0 = (currentTheta > 108.0 ? 144.0 : 72.0);
+ anglecalc.c1 = 1.0;
+ }
+ anglecalc.c2 = 0.0;
+ }
+
} else if (coordination == 5) { // trigonal bipyramidal
-
anglecalc.c0 = 1.0;
// We've already done some of our work above -- look for axial markings
if (b->HasData("UFF_CENTRAL_ATOM")
@@ -932,43 +1053,6 @@
}
anglecalc.c2 = 0.0;
}
- /*
- This section is commented out.
- If you want to try to tackle 7-coordinate species, give this a try
- and change the first conditional above to >7, not >=7
- This doesn't work so well because it's hard to classify between
- axial-equatorial (90 degrees) and proximal equatorial (~72 degrees).
- You'll probably have to do something like the pre-evaluation of 5-coordinate atoms at the top of this method
- Honestly, it's just easier to use VDW pushing to get the job done.
-
- } else if (coordination == 7) { // pentagonal bipyramidal
- currentTheta = a->GetAngle(&*b, &*c);
-
- anglecalc.c0 = 1.0;
- if (currentTheta >= 155.0) { // axial ligands = linear
- anglecalc.coord = 1; // like sp
- anglecalc.theta0 = 180.0;
- anglecalc.c1 = 1.0;
- } else if (currentTheta < 155.0 && currentTheta >= 110.0) { // distal equatorial
- anglecalc.coord = 3; // like sp3
- anglecalc.theta0 = 144.0;
- anglecalc.c1 = -1.0;
- } else if (currentTheta < 110.0 && currentTheta >= 85.0) { // axial-equatorial
- anglecalc.coord = 4; // like sq. planar or octahedral
- anglecalc.theta0 = 90.0;
- anglecalc.c1 = 1.0;
- } else if (currentTheta < 85.0) { // proximal equatorial
- anglecalc.coord = 3; // general case (i.e., like sp3)
- anglecalc.theta0 = 72.0;
- anglecalc.c1 = -1.0;
- }
- anglecalc.c2 = 0.0;
-
- // Also add a VDW 1-3 interaction to distort slightly
- if (SetupVDWCalculation(a, c, vdwcalc)) {
- _vdwcalculations.push_back(vdwcalc);
- }
- */
else { // normal coordination: sp, sp2, sp3, square planar, octahedral
anglecalc.coord = coordination;
anglecalc.theta0 = parameterB->_dpar[1];
@@ -982,7 +1066,6 @@
anglecalc.cosT0 = cos(anglecalc.theta0 * DEG_TO_RAD);
anglecalc.zi = parameterA->_dpar[5];
anglecalc.zk = parameterC->_dpar[5];
-
// Precompute the force constant
bondPtr = _mol.GetBond(a,b);
bondorder = bondPtr->GetBondOrder();
Index: data/atomtyp.txt
===================================================================
--- data/atomtyp.txt (revision 3642)
+++ data/atomtyp.txt (working copy)
@@ -1,7 +1,8 @@
##############################################################################
# #
-# Open Babel file: atomtyp.txt # # #
+# Open Babel file: atomtyp.txt #
# #
+# #
# Copyright (c) 1998-2001 by OpenEye Scientific Software, Inc. #
# Some portions Copyright (c) 2001-2008 Geoffrey R. Hutchison #
# Part of the Open Babel package, under the GNU General Public License (GPL)#
@@ -20,7 +21,9 @@
# #
##############################################################################
-INTHYB [D4] 3 #any 4 valent atom
+INTHYB [D4] 3 #any 4-valent atom
+INTHYB [D5] 4 #any 5-valent atom
+INTHYB [D6] 5 #any 6-valent atom
INTHYB [C] 3 #sp3 carbon
INTHYB [c,$(C=*)] 2 #sp2 carbon
INTHYB [$([#6]([#8D1])[#8D1])] 2 #sp2 carbon
@@ -38,8 +41,11 @@
INTHYB [P] 3 #sp3 phosphorus
INTHYB [#15;$([PD1]=*)] 2 #sp2 phosphorus
+INTHYB [PD5] 4 #sp3d phosphorus, as in PF5
+INTHYB [Pv5] 4 #sp3d phosphorus, as in H3PO4
INTHYB [S] 3 #sp3 sulfur
INTHYB [#16;s,$([SD1]=*)] 2 #sp2 sulfur
+INTHYB [SD6] 5 #sp3d2 sulfur, as in SF6
INTHYB [B] 2 #sp2 boron
INTHYB [BD4] 3 #sp3 boron
@@ -220,6 +226,10 @@
######################## Add Extra Definitions Here ##########################
+#INTHYB [U] 3
+#INTHYB [W] 3
+INTHYB [Mo] 5 # for development
+INTHYB [Cr] 4 # for development
############################# End Extra Definitions ##########################
Index: data/UFF.prm
===================================================================
--- data/UFF.prm (revision 3642)
+++ data/UFF.prm (working copy)
@@ -144,7 +144,7 @@
atom [#102] No6+3
atom [#103] Lw6+3
-# Atom r1 theta0 x1 D1 zeta Z1 Vi Uj Xi Hard Radius
+# Atom r1 theta0 x1 D1 zeta Z1 Vi Uj Xi Hard Radius
param H_ 0.354 180 2.886 0.044 12 0.712 0 0 4.528 6.9452 0.371
param H_b 0.46 83.5 2.886 0.044 12 0.712 0 0 4.528 6.9452 0.371
param He4+4 0.849 90 2.362 0.056 15.24 0.098 0 0 9.66 14.92 1.3
------------------------------------------------------------------------------
Download Intel® Parallel Studio Eval
Try the new software tools for yourself. Speed compiling, find bugs
proactively, and fine-tune applications for parallel performance.
See why Intel Parallel Studio got high marks during beta.
http://p.sf.net/sfu/intel-sw-dev
_______________________________________________
OpenBabel-Devel mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/openbabel-devel
