pemsley/coot/correlations

PyObject *map_to_model_correlation_stats_per_residue_range_py(
    int imol,                                  # Model molecule number
    const std::string &chain_id,               # Chain identifier
    int imol_map,                              # Map molecule number  
    unsigned int n_residue_per_residue_range,  # Number of residues per range (typically 1)
    short int exclude_NOC_flag                 # Exclude N, O, C atoms (1=yes, 0=no)
)

# Get per-residue correlation stats for chain A
# Use exclude_NOC_flag=1 to get valid side-chain correlations
stats = map_to_model_correlation_stats_per_residue_range_py(
    imol=0,           # Model molecule
    chain_id="A",     # Chain A
    imol_map=1,       # Map molecule
    n_residue_per_residue_range=1,  # Per-residue statistics
    exclude_NOC_flag=1  # MUST be 1 for valid side-chain correlations
)
all_atom = stats[0]      # All-atom correlations
sidechain = stats[1]     # Side-chain correlations
# Find residues with poor all-atom correlation
for res in all_atom:
    resno = res[0][1]
    n_points = res[1][0]
    corr = res[1][1]
    if corr < 0.7:
        print(f"Poor fit: residue {resno}, correlation={corr:.3f}")

PyObject *map_to_model_correlation_py(
    int imol,
    PyObject *residue_specs,        # List of residue specs to evaluate
    PyObject *neighb_residue_specs, # Neighboring residues to exclude from grid
    unsigned short int atom_mask_mode,  # Which atoms to include (see below)
    int imol_map
)

# Evaluate side-chain fit for residues 40-44
residue_specs = [[chain_id, res_no, ins_code] for res_no in range(40, 45)]
correlation = map_to_model_correlation_py(
    imol=1,
    residue_specs=residue_specs,
    neighb_residue_specs=[],  # No neighbors to exclude
    atom_mask_mode=2,  # Side-chain atoms only
    imol_map=2
)
print(f"Side-chain correlation: {correlation}")

PyObject *map_to_model_correlation_stats_py(
    int imol,
    PyObject *residue_specs,
    PyObject *neighb_residue_specs,
    unsigned short int atom_mask_mode,
    int imol_map
)

PyObject *map_to_model_correlation_per_residue_py(
    int imol,
    PyObject *residue_specs,
    unsigned short int atom_mask_mode,
    int imol_map
)

# Get per-residue correlations for a chain
residues = get_residues_in_chain_py(imol=1, chain_id="A")
correlations = map_to_model_correlation_per_residue_py(
    imol=1,
    residue_specs=residues,
    atom_mask_mode=0,  # All atoms
    imol_map=2
)
# Find worst 10 residues
worst_10 = sorted(correlations, key=lambda x: x[1])[:10]
for spec, corr in worst_10:
    print(f"Residue {spec}: correlation = {corr}")

PyObject *all_molecule_ramachandran_score_py(int imol)

rama_data = coot.all_molecule_ramachandran_score_py(0)
per_res = rama_data[5]  # index 5 = per-residue list (confirmed from C++ source)
# Filter for chain A outliers
chain_a_outliers = [(r[1][1], r[1][2], r[2], r[3][1])
                    for r in per_res
                    if isinstance(r, list) and r[1][0] == "A" and r[2] < 0.02]
# Sort worst first
chain_a_outliers.sort(key=lambda x: x[2])
for resno, ins, prob, resname in chain_a_outliers:
    print(f"A/{resno} {resname}  prob={prob:.6f}  *** OUTLIER ***")
# Find single worst outlier across whole molecule
worst = min(per_res, key=lambda x: x[2] if isinstance(x, list) else 999)
chain, resno, ins = worst[1]
print(f"Worst: {chain}/{resno}, prob={worst[2]:.6f}")

PyObject *all_molecule_rotamer_score_py(int imol)

score, n_residues = all_molecule_rotamer_score_py(1)
print(f"Overall rotamer score: {score} for {n_residues} residues")

PyObject *rotamer_graphs_py(int imol)

rotamers = rotamer_graphs_py(1)
# Find worst rotamer (excluding missing atoms)
valid_rotamers = [r for r in rotamers if r[3] > 0]
worst = min(valid_rotamers, key=lambda x: x[3])
chain, resno, score = worst[0], worst[1], worst[3]
print(f"Worst rotamer: {chain} {resno}, score = {score}%")
# Go to worst rotamer
coot.set_go_to_atom_chain_residue_atom_name(chain, resno, 'CA')

void deviant_geometry(int imol)

# 1. Load tutorial data
load_tutorial_model_and_data()
# 2. Run comprehensive validation
imol = 1  # Model molecule
imol_map = 2  # Map molecule
# Ramachandran validation
rama_data = coot.all_molecule_ramachandran_score_py(imol)
per_res = rama_data[5]  # index 5 = per-residue list (confirmed from C++ source)
worst_rama = min(per_res, key=lambda x: x[2] if isinstance(x, list) else 999)
print(f"Worst Ramachandran: {worst_rama[1]}, score={worst_rama[2]:.6f}")
# Rotamer validation  
rotamers = rotamer_graphs_py(imol)
worst_rot = min([r for r in rotamers if r[3] > 0], key=lambda x: x[3])
print(f"Worst rotamer: {worst_rot[0]} {worst_rot[1]}, score={worst_rot[3]}%")
# Density fit validation
fit_stats = map_to_model_correlation_stats_per_residue_range_py(
    imol, "A", imol_map, 1, 0
)
# Geometry validation
deviant_geometry(imol)

float auto_fit_best_rotamer(
    int imol_coords,
    const char *chain_id,
    int resno,
    const char *insertion_code,
    const char *altloc,
    int imol_map,
    int clash_flag,        # 1 to check clashes, 0 to ignore
    float lowest_probability  # Minimum acceptable probability (e.g., 0.01)
)

# Fix worst rotamer
new_score = auto_fit_best_rotamer(
    imol_coords=1,
    chain_id='A',
    resno=91,
    insertion_code='',
    altloc='',
    imol_map=2,
    clash_flag=1,  # Check for clashes
    lowest_probability=0.01
)
print(f"New rotamer score: {new_score}%")

void refine_zone(
    int imol,
    const char *chain_id,
    int resno_start,
    int resno_end,
    const char *altconf
)

# Refine 5 residues around residue 42
refine_zone(1, 'A', 40, 44, '')
accept_regularizement()  # Accept the refinement

void pepflip(
    int imol,
    const char *chain_id,
    int resno,
    const char *ins_code,
    const char *altconf
)

# Flip peptide at residue 41
pepflip(1, 'A', 41, '', '')
refine_zone(1, 'A', 40, 44, '')
accept_regularizement()