import math


def ligand_rmsd(pose1, pose2):
    """Calculate RMSD of two ligands

    Args:
        pose1 (pose): Pose 1
        pose2 (pose): Pose 2

    Returns:
        float: RMSD score
    """
    # define ligand residue id of pose 1
    ref_res_num = 0
    for j in range(1, pose1.size() + 1):
        if not pose1.residue(j).is_protein():
            ref_res_num = j
            break

    pose1_rsd = pose1.conformation().residue(ref_res_num)

    inp_res_num = 0

    for j in range(1, pose2.size() + 1):
        if not pose2.residue(j).is_protein():
            inp_res_num = j
            break

    pose2_rsd = pose2.conformation().residue(inp_res_num)

    j = 1
    dist_sum = 0
    for i in range(1, pose1_rsd.nheavyatoms() + 1):
        if j <= pose2_rsd.nheavyatoms():
            x_dist = (pose1_rsd.atom(i).xyz()[
                      0] - pose2_rsd.atom(i).xyz()[0]) ** 2
            y_dist = (pose1_rsd.atom(i).xyz()[
                      1] - pose2_rsd.atom(i).xyz()[1]) ** 2
            z_dist = (pose1_rsd.atom(i).xyz()[
                      2] - pose2_rsd.atom(i).xyz()[2]) ** 2
            dist_sum += x_dist + y_dist + z_dist
        j += 1

    rmsd = math.sqrt(dist_sum / pose1_rsd.nheavyatoms())

    return rmsd