ok, so using old posts to this discussion group I think I figured this out. Sorry to make you waste your time.

My problem was that I was building a neutral tree, but then optimizing the likelihood function, so it changed the branch lengths (duh.) So if I fix the branch lengths, the unoptimized likelihood is what I want.

I then build a second tree, constraining the branch lengths by a variable R. And get the likelihood of this constrained model, and output the relative constraint.

Here's what I have, can you confirm that this is doing what I think it's doing?

global TvTs:= 0.205; /* the global transition/transversion ratio */

HKYM =  {{*,TvTs*t,t,TvTs*t}
    {TvTs*t,*,TvTs*t,t}
    {t,TvTs*t,*,TvTs*t}
    {TvTs*t,t,TvTs*t,*}};


eqf = {{0.32}{0.18}{0.18}{0.32}};


Model  HKY = (HKYM,eqf,1);

ACCEPT_ROOTED_TREES = 1;

t    = 1; 

/* it is important to set the branch length parameter to 1, so that the rate matrix is evaluated properly for the function,
otherwise the scaling factor will be of by a factor of 1/t */

scalingB = computeScalingFactorB (HKYM, eqf);

Tree L_1 = (((Scer:0.23,Spar:0.15):0.121206,Smik:0.32):0.12,Skud:0.329292,Sbay:0.44);
UseModel(HKY);
Tree L_2 = (((Scer:0.23,Spar:0.15):0.121206,Smik:0.32):0.12,Skud:0.329292,Sbay:0.44);


fprintf (stdout, "\nBranch scaling factor computed to be ", scalingB,"\n\n");

/* get all branch names and branch lengths*/

branchNames  = BranchName  (L_1,-1);
branchLengths = BranchLength (L_1,-1);
global R = 1;
for (k = 0; k < Columns (branchNames)-1; k=k+1)
{
ExecuteCommands ("L_1."+branchNames[k]+".t=L_1."+branchNames[k]+".t/scalingB;");
}
branchLengths2 = BranchLength (L_1,-1);
for(k=0; k < Columns(branchNames)-1; k=k+1)
{
ExecuteCommands ("L_2."+branchNames[k]+".t := R * "+branchLengths2[k]+";");

}

characters = {{"A","C","G","T"}
    {"1","","",""}};

DataSet B = ReadDataFile("temp.fa");
DataSetFilter partition1 = CreateFilter(B,1);

LikelihoodFunction neutral_LF = (partition1,L_1);
LikelihoodFunction constrained_LF = (partition1,L_2);
     
fprintf(stdout, "\nunoptimized\n", neutral_LF,"\n");
fprintf(stdout, "\nunoptimized constrained\n", constrained_LF,"\n");

Optimize(res,neutral_LF);
Optimize(res2,constrained_LF);

fprintf(stdout, "\nneutral optimized\n",neutral_LF);
fprintf(stdout, "\nconstrained opt.\n",constrained_LF);
rescaled_R = R/(L_2.Scer.t/BranchLength(L_2,0));
fprintf(stdout, "\nrescaled R = ", rescaled_R);