<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def _real(coords1, charges1, coords2, charges2, rcut, alpha, box):
 """Calculate ewald real part. Box has to be a cuboidal box you should transform any other box shape to a cuboidal box before using this. """

    n = coords1.shape[0]
    m = coords2.shape[0]
    # Unit vectors
    a = box[0]
    b = box[1]
    c = box[2]
    
    
    # This is helpful to add the correct number of boxes
    l_max = int(np.ceil(2.0 * rcut / np.min(np.trace(box))))
    result = np.zeros(n)
    
    
    for i in range(n):
        q_i = charges1[i]
        r_i = coords1[i]

        for j in range(m):
            q_j = charges2[j]
            r_j = coords2[j]

            for l_i in range(-l_max, l_max + 1):
                for l_j in range(-l_max, l_max + 1):
                    for l_k in range(-l_max, l_max + 1):
                        nv = l_i * a + l_j * b + l_k * c
                        r_j_n = r_j + nv
                        
                        r_ij = _dist(r_i, r_j_n)
                        if r_ij < 1e-10 or r_ij > rcut:
                            continue
                        value = q_i * q_j * math.erfc(alpha * r_ij) / r_ij
                        result[i] += value

    return result