pyqcprot N < 2

[adw62]

Found a strange error for a specific system: example.zip

I'm prepping the systems with the API like so:

#prep
config = Config()
config.ligand_net_charge = 1
L1 = Ligand('ligand.mol2', config=config)
#L1.antechamber_prepare_mol2()

#Use TIES20 to build all hybrid moleules
for flu_analogue in new_mols:
    if flu_analogue != 'l_H14':
        continue

    config.workdir = flu_analogue
    config.md_engine = 'openmm'

    L2 = Ligand(flu_analogue+'.mol2', config=config)
    #L2.antechamber_prepare_mol2()

    pair = Pair(L1, L2, ligand_net_charge=1, config=config)
    pair.make_atom_names_unique()

    hybrid = pair.superimpose()

which results in:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-7-44456e3772f6> in <module>
     15 
     16 
---> 17     hybrid = pair.superimpose()
     18 
     19     # save meta data

~/miniconda3/envs/TIES_loc/lib/python3.7/site-packages/ties/pair.py in superimpose(self, **kwargs)
    174                                          starting_node_pairs=starting_node_pairs,
    175                                          parmed_ligA=parmed_ligA, parmed_ligZ=parmed_ligZ,
--> 176                                          starting_pair_seed=self.config.superimposition_starting_pair)
    177 
    178         assert len(suptops) == 1

~/miniconda3/envs/TIES_loc/lib/python3.7/site-packages/ties/topology_superimposer.py in superimpose_topologies(top1_nodes, top2_nodes, pair_charge_atol, use_charges, use_coords, starting_node_pairs, force_mismatch, disjoint_components, net_charge_filter, net_charge_threshold, redistribute_charges_over_unmatched, parmed_ligA, parmed_ligZ, align_molecules, partial_rings_allowed, ignore_charges_completely, ignore_bond_types, ignore_coords, use_general_type, use_only_element, check_atom_names_unique, starting_pairs_heuristics, starting_pair_seed)
   3084                                       use_general_type=use_general_type,
   3085                                       starting_pairs_heuristics=starting_pairs_heuristics,
-> 3086                                       starting_pair_seed=starting_pair_seed)
   3087     if not suptops:
   3088         raise Exception('Error: Did not find a single superimposition state.'

~/miniconda3/envs/TIES_loc/lib/python3.7/site-packages/ties/topology_superimposer.py in _superimpose_topologies(top1_nodes, top2_nodes, mda1_nodes, mda2_nodes, starting_node_pairs, ignore_coords, use_general_type, starting_pairs_heuristics, starting_pair_seed)
   3623                                     suptop=suptop,
   3624                                     ignore_coords=ignore_coords,
-> 3625                                     use_element_type=use_general_type)
   3626         if candidate_suptop is None or len(candidate_suptop) == 0:
   3627             # there is no overlap, ignore this case

~/miniconda3/envs/TIES_loc/lib/python3.7/site-packages/ties/topology_superimposer.py in _overlay(n1, n2, parent_n1, parent_n2, bond_types, suptop, ignore_coords, use_element_type)
   2962     # and each which have its own solution, which in turn have to be merged
   2963     for atom_type in combinations:
-> 2964         all_solutions.append(solve_one_combination(atom_type, ignore_coords))
   2965 
   2966     assert len(all_solutions) > 1

~/miniconda3/envs/TIES_loc/lib/python3.7/site-packages/ties/topology_superimposer.py in solve_one_combination(one_atom_species, ignore_coords)
   2715             candidates = [list(v.values())[0] for v in atoms.values()]
   2716             largest_candidates = get_largest(candidates)
-> 2717             best = extract_best_suptop(largest_candidates, ignore_coords)
   2718             return best
   2719 

~/miniconda3/envs/TIES_loc/lib/python3.7/site-packages/ties/topology_superimposer.py in extract_best_suptop(suptops, ignore_coords)
   3311     rmsds = []
   3312     for suptop in candidates:
-> 3313         rmsd = suptop.align_ligands_using_mcs()
   3314         rmsds.append(rmsd)
   3315 

~/miniconda3/envs/TIES_loc/lib/python3.7/site-packages/ties/topology_superimposer.py in align_ligands_using_mcs(self, overwrite_original)
    779         #    return self.rmsd()
    780 
--> 781         rmsd, (rotation_matrix, ligA_mcs_com, ligZ_mcs_com) = superimpose_coordinates(mcs_ligA_coords, mcs_ligZ_coords)
    782 
    783         if not overwrite_original:

~/miniconda3/envs/TIES_loc/lib/python3.7/site-packages/ties/transformation.py in superimpose_coordinates(ref, mob)
    152 
    153     # Calculate rmsd and rotation matrix
--> 154     rmsd = pyqcprot.CalcRMSDRotationalMatrix(ref_origin, mob_origin, N, rotation, None)
    155 
    156     # reshape so that it can be used directly on all coordinates

pyqcprot.pyx in pyqcprot.CalcRMSDRotationalMatrix()

TypeError: must be real number, not NoneType

I have tracked this down to this call of align_ligands_using_mcs():

https://github.com/UCL-CCS/TIES20/blob/fbc77f197f6fd38338730cbc64fa2e066bbfa22e/ties/topology_superimposer.py#L3306-L3310

which then fails here:

https://github.com/UCL-CCS/TIES20/blob/master/ties/topology_superimposer.py#L776-L777

It seems to fail for the case when N <= 2 and by adding this code before the superimpose_coordinates call the problem goes away:

if len(mcs_ligA_coords) <= 2:
         # pyqcprot does not like two or less atoms?
         return self.rmsd()

Why does pyqcprot not like N=<2? What is causing N to be =<2?

[bieniekmateusz]

I gyess you need at least 3 points in order to obtain a meaningful 3D superimposition. 2 points could rotate about the other axis. I'll add an extra check underneath as well. One way our of this is to calculate rmsd() as is, without the superimposition.

[adw62]

That sounds like it makes sense. 'One way our of this is to calculate rmsd() as is, without the superimposition.' Is that the same as just returning self.rmsd() or should something else be run?

[adw62]

This has been addresses here #301