The coordinates should not change on copying the molecule (I just tested
this).
However, I'm guessing that you are comparing two floating point numbers
using "if x==y". If you search around the web, you will see that you should
this doesn't work for a variety of reasons. Instead, you should compare two
non-integers using "if (x-y)<delta" (where delta is 0.0001 or whatever).
An alternative method for your problem would be to search for the isolated
molecule in the complex. See
http://openbabel.org/api/2.3/group__substructure.shtml and
CompileMoleculeQuery. I think (!) these should work from Python okay...
- Noel
On 14 May 2012 04:49, inviso <[email protected]> wrote:
> As a disclaimer, I'm an amazingly poor excuse for a chemist. In this
> case, I'm simply a programming resource in a research project so
> please forgive any amateur mistakes or abuses of terminology. The
> project is being done in Python using pybel.
>
> As part of our process, we will be doing a scf calculation in NWChem
> of both a metal ligand complex and the host portion (the ligand in
> this case). I'm getting a file with the guest metal and another file
> with a series of metal ligand complexes, but I am not provided the
> host portion directly. The goal is to extract the metal portion from
> the complex and then generate NWChem input files.
>
> For example, I might have RhCl3 as the metal and C9H21Cl3N3Rh as the
> complex. I'm reading those two out of .xyz files so I have coordinate
> data as well. I would like to end up with C9H21N3. Assuming that the
> ligand might also contain Cl, it's probably not sufficient to just
> randomly remove a couple Cl and a Rh. My somewhat naive approach was
> to find atoms that shared both atomic number and coordinates. Is that
> an accurate approach?
>
> When doing that, I discovered that the coordinates of the data read in
> from file would change when I made a copy of the molecule. I need to
> hold onto the complex as well as extract/calculate the host so I
> copied the molecule first by doing something like copy_of_complex =
> pybel.Molecule(complex_molecule). I then iterated over the guest
> atoms and removed the matching atoms in the copy of the complex.
> Because the coordinates of the copy did not match what I read in from
> the file, it couldn't find the matching atoms. Done against the
> original metal ligand complex read from file, it found matching atoms.
> In the example code/output below, note that the coordinates of the
> copy are different for the first atom.
>
> Is this occurring because openbabel is storing the coordinates in
> floating point numbers instead of "decimal" style numbers? If I
> remember correctly from looking at the API earlier, they are stored as
> doubles. Any way around this? Or is this unexpected behavior?
>
> Also, is my entire approach flawed? Should I be trying to split the
> complex prior to loading it into openbabel? Should I be reading the
> complex in twice as it seems to be consistent on read? Is there a
> better method available in openbabel to accomplish this?
>
> Finally, one other oddity I've noticed. Even when removing the atoms
> from the original complex, the formula attribute of the pybel molecule
> does not update even tthough the atoms are indeed gone. I'm removing
> them using complex_molecule.OBMol.DeleteAtom(atom_from_complex_to_remove,
> True). In the output below, note that the formula does not change,
> but the number of atoms show that 4 have been removed (as expected).
>
> Thanks in advance for any insights you can provide,
>
> Greg Fortune
>
> --------------
> Code
> --------------
> import pybel
>
> def molecule_diff(minuend, subtrahend):
> for sub_atom in subtrahend.atoms:
> for min_atom in minuend.atoms:
> if(sub_atom.atomicnum == min_atom.atomicnum and
> sub_atom.coords == min_atom.coords):
> t = minuend.OBMol.DeleteAtom(min_atom.OBAtom, True)
>
> def print_atom(atom, indent='\t'):
> print(indent + '%s: %s'% (pybel.ob.etab.GetSymbol(atom.atomicnum),
> str(atom.coords)))
>
> guest = [mol for mol in pybel.readfile('xyz', 'guest.xyz')][0]
> comp = [mol for mol in pybel.readfile('xyz', 'amine.xyz')][0]
> host = pybel.Molecule(comp)
>
> print('guest: ' + guest.formula)
> print_atom(guest.atoms[0], '')
> print('complex: ' + comp.formula)
> print_atom(comp.atoms[0], '')
> print('host: ' + host.formula)
> print_atom(host.atoms[0], '')
>
>
> molecule_diff(comp, guest)
>
> print(comp.formula) #This should give back the formula without RhCl3,
> but it gives the full formula of the original complex
> print(len(comp.atoms)) #The original total is 37 atoms so 33 makes
> sense assuming we removed the 4 atoms from the guest.
>
>
> --------------
> Output
> --------------
> guest: Cl3Rh
> Rh: (0.06029, -2.22815, 0.03786)
> complex: C9H21Cl3N3Rh
> Rh: (0.06029, -2.22815, 0.03786)
> host: C9H21Cl3N3Rh
> Rh: (0.0603, -2.2281, 0.0379)
>
> C9H21Cl3N3Rh
> 33
>
>
> ------------------------------------------------------------------------------
> Live Security Virtual Conference
> Exclusive live event will cover all the ways today's security and
> threat landscape has changed and how IT managers can respond. Discussions
> will include endpoint security, mobile security and the latest in malware
> threats. http://www.accelacomm.com/jaw/sfrnl04242012/114/50122263/
> _______________________________________________
> OpenBabel-discuss mailing list
> [email protected]
> https://lists.sourceforge.net/lists/listinfo/openbabel-discuss
>
------------------------------------------------------------------------------
Live Security Virtual Conference
Exclusive live event will cover all the ways today's security and
threat landscape has changed and how IT managers can respond. Discussions
will include endpoint security, mobile security and the latest in malware
threats. http://www.accelacomm.com/jaw/sfrnl04242012/114/50122263/
_______________________________________________
OpenBabel-discuss mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/openbabel-discuss
