QChASM / AaronTools.py

Python tools for automating routine tasks encountered when running quantum chemistry computations.
https://aarontools.readthedocs.io/en/latest/
GNU General Public License v3.0
43 stars 8 forks source link

I failed map ligands in a purely organic system. #19

Closed Senppoa closed 9 months ago

Senppoa commented 12 months ago

Hello, dear developers,

I am trying to use AaronTools.py to map ligands in a purely organic system (The TS template system I have tried is "Vinyl_boronic_acid_addition" reaction R/ts1.xyz in the code folder of AARON) it cannot give me a correct result. Results only the ligand which I want to add in the TS structure. So I tried to used perl version AaronTools, below is the perl program reporting error.

微信截图_20231010195010

So how could I specify the ligand atoms in the .xyz file? Is there any other way to map ligands by AaronTools python script in this purely organic reaction?

ajs99778 commented 12 months ago

It looks like the Python AaronTools is failing to figure out where the ligand stops - it thinks the entire structure is the ligand. For organometallic cases, it would just stop when it finds a metal. I'll need to take a closer look at this to see where exactly it's going wrong.

I did find a bug in the mapLigand.py script that caused it to basically not use the -c/--center flag, and I've fixed that. The mapLigand.py script that is now on GitHub should work for organic systems as long as you use the -c or --center flag.

That being said, AaronTools should be able to figure out the ligand in this case without being told what the center is, and I will be looking into a fix for that.

You could also manually specify which atom is the center and which atoms are the ligand atoms in the comment line of the file. I used this:

L:7-11,17-22,28-29,37-42,48-53,55-56,61-66 C:67

to specify the ligand atoms (L:) and the boron center (C:). It's a little bit of a mess because the ligand atoms are not grouped together in the file. For completeness, the command I used with this comment was mapLigand.py ts1.xyz -l 28,29=AcAc -o new.xyz

I think the Perl version is a bit less sophisticated when it comes to figuring out the ligand. If I recall correctly, all atoms that come after the center are assumed to be ligand atoms, so you would need to reorder the atoms and specify the center with C: in the comment.

Senppoa commented 11 months ago

Thanks a lot! I have successfully mapped ligand in this purely organic system by command line operation you provided of AaronTools.

However, this approach does not seem to be able to achieve both ligand replacement and substrate modification at the same time due to changes in atomic number, I do this by python API of AaronTools before. So, if possible, I would suggest adding parameters to specify the central atom and ligands of the system for the map_ligand function in Geometry object.

ajs99778 commented 11 months ago

For the time being, you can run detect_components prior to map_ligand

example:

from AaronTools.geometry import Geometry

geom = Geometry("ts1.xyz")

# set the boron atom as the center
# alternatively, use "67"
geom.detect_components("B") 
geom.map_ligand("AcAc", "28,29")
geom.substitute("CF3", "68")

geom.write(outfile="ts1-AcAc-CF3.xyz")
Senppoa commented 11 months ago

Thanks Dr. Schaefer! It works well.