def METHOD_NAME(mol, mints, verbose=1):
    """Computes array of the EFP contribution to the potential felt by
    QM atoms due to permanent EFP moments. Used for SCF procedure.
    Parameters
    ----------
    mol : :py:class:`psi4.core.Molecule`
        Source of quantum mechanical atoms. As its `EFP` member data, contains
        a :py:class:`pylibefp.core.efp` object that is the source and computer
        of EFP fragments.
    mints : `psi4.core.MintsHelper`
        Integral computer.
    verbose : int, optional
        Whether to print out multipole coordinates and values. 0: no printing.
        1: print charges and dipoles. 2: additionally print quadrupoles and octupoles.
    Returns
    -------
    ndarray
        (nbf, nbf) EFP charge through octupole contribution to the potential
    """
    # get composition counts from libefp
    efpobj = mol.EFP
    nfr = efpobj.get_frag_count()
    natoms = efpobj.get_frag_atom_count()
    # get multipoles count, pos'n, values from libefp
    #   charge + dipoles + quadrupoles + octupoles = 20
    nmp = efpobj.get_multipole_count()
    xyz_mp = np.asarray(efpobj.get_multipole_coordinates(verbose=verbose)).reshape(nmp, 3)
    val_mp = np.asarray(efpobj.get_multipole_values(verbose=verbose)).reshape(nmp, 20)
    #                    0  X  Y  Z  XX   YY   ZZ   XY   XZ   YZ
    prefacs = np.array([ 1, 1, 1, 1, 1/3, 1/3, 1/3, 2/3, 2/3, 2/3,
        1/15, 1/15, 1/15, 3/15, 3/15, 3/15, 3/15, 3/15, 3/15, 6/15])
    #   XXX   YYY   ZZZ   XXY   XXZ   XYY   YYZ   XZZ   YZZ   XYZ
    # EFP permanent moment contribution to the Fock Matrix
    nbf = mints.basisset().nbf()
    V2 = np.zeros((nbf, nbf))
    # Cartesian basis one-electron EFP perturbation
    efp_ints = np.zeros((20, nbf, nbf))
    for imp in range(nmp):
        origin = xyz_mp[imp]
        # get EFP multipole integrals from Psi4
        p4_efp_ints = mints.ao_efp_multipole_potential(origin=origin)
        for pole in range(20):
            efp_ints[pole] = np.asarray(p4_efp_ints[pole])
        # add frag atom Z into multipole charge (when pos'n of atom matches mp)
        for ifr in range(nfr):
            atoms = efpobj.get_frag_atoms(ifr)
            for iat in range(natoms[ifr]):
                xyz_atom = [atoms[iat]['x'], atoms[iat]['y'], atoms[iat]['z']]
                if np.allclose(xyz_atom, origin, atol=1e-10):
                    val_mp[imp, 0] += atoms[iat]['Z']
        # scale multipole integrals by multipole magnitudes. result goes into V
        for pole in range(20):
            efp_ints[pole] *= -prefacs[pole] * val_mp[imp, pole]
            V2 += efp_ints[pole]
    return V2