the mutation matrix i used was:

Mutationmatrix

0.5447  0.0250  0.0311  0.0298  0.0807  0.0426  0.0497  0.0517  0.0203  0.0211  0.0251  0.0356  0.0328  0.0138  0.0557  0.1110  0.0554  0.0073  0.0148  0.0580  
0.0151  0.5685  0.0275  0.0129  0.0168  0.0571  0.0261  0.0118  0.0351  0.0074  0.0101  0.1103  0.0173  0.0072  0.0120  0.0240  0.0203  0.0121  0.0114  0.0080  
0.0141  0.0207  0.4755  0.0592  0.0135  0.0280  0.0255  0.0254  0.0452  0.0057  0.0050  0.0339  0.0077  0.0050  0.0099  0.0385  0.0295  0.0048  0.0152  0.0061  
0.0202  0.0145  0.0887  0.5687  0.0067  0.0386  0.1123  0.0206  0.0346  0.0028  0.0043  0.0286  0.0050  0.0039  0.0175  0.0380  0.0221  0.0048  0.0091  0.0038  
0.0120  0.0041  0.0044  0.0015  0.5719  0.0018  0.0009  0.0027  0.0031  0.0047  0.0036  0.0019  0.0046  0.0061  0.0016  0.0114  0.0100  0.0085  0.0061  0.0112  
0.0174  0.0389  0.0253  0.0234  0.0050  0.4155  0.0531  0.0077  0.0356  0.0054  0.0119  0.0503  0.0237  0.0050  0.0126  0.0238  0.0190  0.0073  0.0114  0.0089  
0.0404  0.0352  0.0457  0.1346  0.0050  0.1052  0.5182  0.0184  0.0320  0.0062  0.0081  0.0668  0.0132  0.0036  0.0243  0.0395  0.0354  0.0073  0.0106  0.0141  
0.0482  0.0182  0.0523  0.0283  0.0168  0.0176  0.0212  0.7729  0.0153  0.0029  0.0043  0.0221  0.0064  0.0047  0.0120  0.0428  0.0144  0.0097  0.0045  0.0054  
0.0050  0.0143  0.0245  0.0125  0.0050  0.0213  0.0097  0.0040  0.5384  0.0022  0.0046  0.0149  0.0073  0.0149  0.0089  0.0089  0.0074  0.0157  0.0364  0.0037  
0.0169  0.0099  0.0102  0.0033  0.0252  0.0106  0.0061  0.0025  0.0071  0.4988  0.1017  0.0117  0.0901  0.0398  0.0063  0.0102  0.0288  0.0145  0.0140  0.1760  
0.0264  0.0176  0.0116  0.0066  0.0252  0.0304  0.0104  0.0048  0.0198  0.1330  0.5999  0.0201  0.1866  0.0914  0.0162  0.0194  0.0245  0.0242  0.0349  0.0803  
0.0273  0.1405  0.0575  0.0324  0.0101  0.0942  0.0630  0.0182  0.0468  0.0112  0.0147  0.4960  0.0218  0.0064  0.0220  0.0402  0.0317  0.0048  0.0121  0.0114  
0.0090  0.0079  0.0047  0.0020  0.0084  0.0158  0.0044  0.0019  0.0081  0.0307  0.0486  0.0078  0.4201  0.0226  0.0031  0.0071  0.0118  0.0060  0.0091  0.0181  
0.0062  0.0054  0.0050  0.0026  0.0185  0.0055  0.0020  0.0023  0.0275  0.0223  0.0392  0.0037  0.0373  0.6051  0.0037  0.0071  0.0072  0.0713  0.1391  0.0137  
0.0265  0.0095  0.0105  0.0123  0.0050  0.0146  0.0143  0.0061  0.0173  0.0037  0.0073  0.0136  0.0055  0.0039  0.7212  0.0274  0.0122  0.0024  0.0038  0.0073  
0.0743  0.0267  0.0570  0.0375  0.0513  0.0389  0.0325  0.0307  0.0244  0.0085  0.0123  0.0351  0.0173  0.0105  0.0385  0.4331  0.0904  0.0157  0.0106  0.0135  
0.0374  0.0228  0.0441  0.0221  0.0454  0.0313  0.0295  0.0104  0.0203  0.0242  0.0158  0.0279  0.0291  0.0108  0.0173  0.0912  0.5182  0.0048  0.0118  0.0405  
0.0007  0.0021  0.0011  0.0007  0.0059  0.0018  0.0009  0.0011  0.0066  0.0019  0.0024  0.0006  0.0023  0.0163  0.0005  0.0024  0.0007  0.7037  0.0182  0.0019  
0.0049  0.0062  0.0110  0.0044  0.0135  0.0091  0.0043  0.0016  0.0488  0.0057  0.0109  0.0052  0.0109  0.1014  0.0026  0.0052  0.0057  0.0580  0.6111  0.0057  
0.0533  0.0122  0.0124  0.0051  0.0698  0.0201  0.0160  0.0053  0.0137  0.2016  0.0703  0.0136  0.0610  0.0279  0.0141  0.0186  0.0552  0.0169  0.0159  0.5124  

my @PAM1= @mutMatrix;
my $e;   # e is equal to the DBL_EPSILON constant
$e=1;$e/=2 while $e/2+1>1;
my $k = 1;
my $mutRate = sum map{ $freq[ $_ ] * ( 1 - $mutMatrix[ $_ ][ $_ ] ) } 0..19;
my $iter = 0;

my @lnPAM=(); # in case the position in the matrix is 0 
LINE: for (my $i=0; $i<20; $i++){
    for (my $j=0; $j<20; $j++){
        if ($mutMatrix[$i][$j] > 0){
            $lnPAM[$i][$j] = log $mutMatrix[$i][$j];
        }else{
            next LINE;
        }
    }
}


# The following loop goes until the difference reach the computer's
# precision.
while( abs( $mutRate - 0.01) > $e ) {
    ++$iter;
    $k *= $mutRate * 100;
    map $_ /= ( $mutRate * 100), @$_ for @lnPAM;
    @PAM1 = map[ map exp( $_ ), @$_ ], @lnPAM;
    $mutRate = sum map{ $freq[ $_ ] * ( 1 - $PAM1[ $_ ][ $_ ] ) } 0..19;
    printf "After %d iteration(s), k=%f, RateMutation=%f\n", $iter, $k, $mutRate;
}

print "\n";

print "PAM1 Matrix:\n\n";
for (my $i=0; $i<20; $i++){ #prints Pam1 matrix                                                                 print PAM1 matrix
    for (my $j=0; $j<20; $j++){
        printf("%.4f  ", $PAM1[$i][$j]);
    }
    print "\n";
}