leeping / geomeTRIC

Geometry optimization code that includes the TRIC coordinate system
https://geometric.readthedocs.io/
Other
158 stars 68 forks source link

Breaking certain bonds and getting invidual atoms as fragments with 3 degress of freedom #178

Open annulen opened 1 year ago

annulen commented 1 year ago

How can I prevent geomeTRIC from creating specific bonds?

For example, there is a structure in attached file: PMDA_chair_dimer2_in.xyz.txt I'd like to break bonds N(8)—H(19) and N(27)—H(40), so that respective protons become freely floating fragments, each having 3 translational coordinates.

In other words, I'd like to get H(19) and H(40) optimized by their cartesian coordinates, and the rest of both molecules to be handled with regular TRIC coordinates.

This might seem a silly idea for this simple case, but I'm actually struggling with optimization of clusters with larger numbers of molecules of the same kind, e.g. 9 and 12. I believe that separating protons that are involved in H-bond systems would decouple their positions from rotations of "host" molecules thus making the whole coordinate system less coupled and thus improving convergence.

I've tried to achieve this with specially prepared PDB file PMDA_chair_dimer2_in.pdb.txt and used following command

geometric-optimize --ase-class=xtb.ase.calculator.XTB --ase-kwargs='{"method":"GFN2-xTB"}' --engine ase --hessian=file:PMDA_chair_dimer2_631.hess --pdb PMDA_chair_dimer2_in.pdb PMDA_chair_dimer2_in.xyz

However, geomeTRIC seems to ignore bonding information and uses fragments 1-20 and 21-40. Note that ".txt" extensions were added by me specifically to allow GitHub to upload these files. Also, I'm using GFN2 here purely for testing reasons and intend to switch to ORCA (via ASE) for the real work.

In case that it's not currently possible I would be glad to get any pointers how could I implement this feature myself.

annulen commented 1 year ago

I've also tried --frag=yes but it didn't affect anything

leeping commented 1 year ago

The "disconnection" of fragments is done by using a different residue ID for each fragment you don't want bonded. Could you try it?

The CONECT records are not used in the coordinate system that geomeTRIC builds. I understand this can be confusing, and the default behavior could be changed such that if CONECT records are present, geomeTRIC will use them and not attempt to build its own set of bonds. However, I understand there are also situations where CONECT records are used to specify only a portion of bonds in the system (corresponding to a ligand, for example), so I'm not sure if this is a universally useful design. Another option would be to print a warning that CONECT records are not used in internal coordinate construction.