What is the best way to make a residue-specific SASA constraint?

I know you said remodel/relax application and not script, but the SASA selection is available in scripting.

https://www.rosettacommons.org/docs/latest/scripting_documentation/RosettaScripts/TaskOperations/taskoperations_pages/LayerDesignOperation

The bad news is that there are few cascading effects to deal with:

The <Layer> elements go within <RESIDUE_SELECTORS> element. In the <MoveMap> of <FastRelax> (or similar mover)  there is written that you can add <ResidueSelector>. This in my hands is not true. No idea why the documentation says you can, but you get a "Error: Element 'ResidueSelector': This element is not expected. Expected is one of ( Jump, Chain, Span )." if you do. So you have to use a movemap factory, which can.
So you have to add:

<MOVE_MAP_FACTORIES>
        <MoveMapFactory name="moove" bb="0" chi="0">
            <Backbone residue_selector="layerz" />
            <Chi residue_selector="layerz" />
        </MoveMapFactory>
</MOVE_MAP_FACTORIES>

Which is called by <FastRelax name="chillax" repeats="1" movemap_factory="moove" cartesian="1" />.

Also, cartesian="1" is needed or the constraint gets ignored. But to get cartesian mode to work you have to add -beta_nov16_cart to your run. 

$ROSETTAPATH/source/bin/rosetta_scripts.static.macosclangrelease -s something.pdb -parser:protocol yourscript.xml -beta_nov16_cart

Sorry to double post but I remembered that I had a snippets of code that could something similar to what you asked* and then I re-read the question and realised you know a priori what the phenylalanine is....

Rosetta Relax has an option `-in:file:movemap`, where you can specify explicitly which residues you want to move both backbone (bb) or sidechain (chi). A movemap is also an element in Rosetta scripts and in Pyrosetta.  The movemap file will have to be made external and not in Rosetta.
Note that the values are in pose numbering and not PDB and using the 23A kind of notation does not work IIRC.
Pyrosetta has a really nice function called pdb2pose that does these conversion for you. Here is a chunk of code as an example (same thing as using the application).

import pyrosetta

pyrosetta.init()

pose = pyrosetta.pose_from_pdb("my_protein.pdb")
pdb2pose = pose.pdb_info().pdb2pose

scorefxn = pyrosetta.get_fa_scorefxn()
# or whatever.
# scorefxn = pyrosetta.rosetta.core.scoring.ScoreFunctionFactory.create_score_function("ref2015_cart.wts")

relax = pyrosetta.rosetta.protocols.relax.FastRelax(scorefxn, 3)
movemap = pyrosetta.MoveMap()
for chain, resi in targets:
    posed = pdb2pose(chain=chain, res=resi)
        print(chain, res, '-->', posed)
    movemap.set_bb(posed, True)
    movemap.set_chi(posed, True)
relax.set_movemap(movemap)
relax.apply(pose)
print('Score', scorefxn(pose))
pose.dump_pdb('output.pdb')

*) Not that it matters but the code I was going to post was a python script using pymol that I quickly adapted:

import pymol2

#Wikipedia Tien et al. 2013 (emp.)
maxASA = {'A': 121.0, 'R': 265.0, 'N': 187.0, 'D': 187.0, 'C': 148.0, 'E': 214.0, 'Q': 214.0, 'G': 97.0, 'H': 216.0, 'I': 195.0, 'L': 191.0, 'K': 230.0, 'M': 203.0, 'F': 228.0, 'P': 154.0, 'S': 143.0, 'T': 163.0, 'W': 264.0, 'Y': 255.0, 'V': 165.0}
    

with pymol2.PyMOL() as pymol:
    pymol.cmd.set('dot_solvent', 1)
    pymol.cmd.set('dot_density', 3)

    def is_surface(atom) -> bool:
        sasa = pymol.cmd.get_area(f'resi {atom.resi} and chain {atom.chain} and name CA')
        rsa = sasa/maxASA[self.mutation.from_residue]
        return True if rsa < 0.2 else False #classic
         
    pymol.cmd.load('my_PDB')
    atoms = pymol.cmd.get_model(f"resn PHE")
    for atom in atoms.atom:
         if is_surface(atom):
              resi = int(re.match('\d+', resi).group())
              targets.append((resi, atom.chain))

print(targets)

N