source: branches/nameserver/GDE/PHYML20130708/phyml/src/init.c

Last change on this file was 10307, checked in by aboeckma, 12 years ago

added most recent version of phyml
File size: 158.0 KB
Line 
1/*
2
3PhyML:  a program that  computes maximum likelihood phylogenies from
4DNA or AA homologous sequences.
5
6Copyright (C) Stephane Guindon. Oct 2003 onward.
7
8All parts of the source except where indicated are distributed under
9the GNU public licence. See http://www.opensource.org for details.
10
11*/
12
13#include "init.h"
14
15void Init_Eigen_Struct(eigen *this)
16{
17  this->next = NULL;
18  this->prev = NULL;
19}
20
21//////////////////////////////////////////////////////////////
22//////////////////////////////////////////////////////////////
23
24void Init_Scalar_Dbl(scalar_dbl *p)
25{
26  p->v     = -1.;
27  p->onoff = ON;
28  p->next  = NULL;
29  p->prev  = NULL;
30}
31
32//////////////////////////////////////////////////////////////
33//////////////////////////////////////////////////////////////
34
35void Init_Scalar_Int(scalar_int *p)
36{
37  p->v    = -1.;
38  p->next = NULL;
39  p->prev = NULL;
40}
41
42//////////////////////////////////////////////////////////////
43//////////////////////////////////////////////////////////////
44
45void Init_Vect_Dbl(int len, vect_dbl *p)
46{
47  p->len  = len;
48  p->next = NULL;
49  p->prev = NULL;
50  p->v    = NULL;
51}
52
53//////////////////////////////////////////////////////////////
54//////////////////////////////////////////////////////////////
55
56void Init_Vect_Int(int len, vect_int *p)
57{
58  p->len  = len;
59  p->next = NULL;
60  p->prev = NULL;
61  p->v    = NULL;
62}
63
64//////////////////////////////////////////////////////////////
65//////////////////////////////////////////////////////////////
66
67void Init_String(t_string *ts)
68{
69  ts->len  = -1.;
70  ts->next = NULL;
71  ts->prev = NULL;
72}
73
74//////////////////////////////////////////////////////////////
75//////////////////////////////////////////////////////////////
76
77void Init_Triplet_Struct(triplet *t)
78{
79  t->next = NULL;
80  t->prev = NULL;
81}
82
83//////////////////////////////////////////////////////////////
84//////////////////////////////////////////////////////////////
85
86void Init_Efrq(t_efrq *f)
87{
88  f->next = NULL;
89  f->prev = NULL;
90}
91
92//////////////////////////////////////////////////////////////
93//////////////////////////////////////////////////////////////
94
95void Init_Tree(t_tree *tree, int n_otu)
96{
97  tree->n_otu                     = n_otu;
98  tree->mat                       = NULL;
99  tree->n_root                    = NULL;
100  tree->e_root                    = NULL;
101  tree->ps_tree                   = NULL;
102  tree->short_l                   = NULL;
103  tree->mutmap                    = NULL;
104  tree->next                      = NULL;
105  tree->prev                      = NULL;
106  tree->next                      = NULL;
107  tree->prev                      = NULL;
108  tree->mixt_tree                 = NULL;
109  tree->geo                       = NULL;
110
111  tree->is_mixt_tree              = NO; 
112  tree->tree_num                  = 0;
113  tree->depth_curr_path           = 0;
114  tree->has_bip                   = NO;
115  tree->n_moves                   = 0;
116  tree->n_improvements            = 0;
117  tree->bl_from_node_stamps       = 0;
118  tree->lock_topo                 = NO;
119  tree->ps_page_number            = 0;
120  tree->init_lnL                  = UNLIKELY;
121  tree->best_lnL                  = UNLIKELY;
122  tree->old_lnL                   = UNLIKELY;
123  tree->c_lnL                     = UNLIKELY;
124  tree->sum_min_sum_scale         = .0;
125  tree->n_swap                    = 0;
126  tree->best_pars                 = 1E+5;
127  tree->n_pattern                 = -1;
128  tree->n_root_pos                = -1.;
129  tree->write_labels              = YES;
130  tree->write_br_lens             = YES;
131  tree->print_boot_val            = 0;
132  tree->print_alrt_val            = 0;
133  tree->num_curr_branch_available = 0;
134  tree->tip_order_score           = .0;
135  tree->write_tax_names           = YES;
136  tree->update_alias_subpatt      = NO;
137  tree->bl_ndigits                = 8;
138  tree->n_short_l                 = 100;
139  tree->norm_scale                = 0.0;
140  tree->br_len_recorded           = NO;
141  tree->max_spr_depth             = 0;
142  tree->apply_lk_scaling          = YES;
143  tree->dp                        = 0;
144}
145
146//////////////////////////////////////////////////////////////
147//////////////////////////////////////////////////////////////
148
149
150void Init_Edge_Light(t_edge *b, int num)
151{
152  b->num                  = num;
153  b->bip_score            = 0;
154  b->dist_btw_edges       = .0;
155  b->topo_dist_btw_edges  = 0;
156  b->has_zero_br_len      = NO;
157  b->n_jumps              = 0;
158  b->l_var->v             = -1.;
159  b->does_exist           = YES;
160  b->l->v                 = -1.;
161  b->bin_cod_num          = -1.;
162  b->l->onoff             = ON;
163
164  b->next                 = NULL;
165  b->prev                 = NULL;
166  b->next_mixt            = NULL;
167  b->prev_mixt            = NULL;
168
169  b->p_lk_left            = NULL;
170  b->p_lk_rght            = NULL;
171  b->p_lk_loc_left        = NULL;
172  b->p_lk_loc_rght        = NULL;
173  b->Pij_rr               = NULL;
174  b->labels               = NULL;
175}
176
177//////////////////////////////////////////////////////////////
178//////////////////////////////////////////////////////////////
179
180
181void Init_Node_Light(t_node *n, int num)
182{
183  n->num                    = num;
184  n->tax                    = -1;
185  n->dist_to_root           = .0;
186  n->common                 = 1;
187  n->ext_node               = NULL;
188  n->name                   = NULL;
189  n->ori_name               = NULL;
190  n->c_seq                  = NULL;
191  n->c_seq_anc              = NULL;
192  n->y_rank                 = 0.;
193  n->y_rank_ori             = 0.;
194  n->y_rank_max             = 0.;
195  n->y_rank_min             = 0.;
196  n->anc                    = NULL;
197  n->rank                   = 0;
198  n->match_node             = NULL;
199  n->id_rank                = 0;
200  n->next                   = NULL;
201  n->prev                   = NULL;
202  /* n->next                  = NULL; */
203  /* n->prev                 = NULL; */
204}
205
206//////////////////////////////////////////////////////////////
207//////////////////////////////////////////////////////////////
208
209
210void Init_NNI(nni *a_nni)
211{
212  a_nni->left         = NULL;
213  a_nni->rght         = NULL;
214  a_nni->b            = NULL;
215  a_nni->init_l       = -1.;
216  a_nni->init_lk      = .0;
217  a_nni->score        = +1.0;
218  a_nni->best_l       = -1.;
219  a_nni->swap_node_v1 = NULL;
220  a_nni->swap_node_v2 = NULL;
221  a_nni->swap_node_v3 = NULL;
222  a_nni->swap_node_v4 = NULL;
223  a_nni->lk0          = UNLIKELY;
224  a_nni->lk1          = UNLIKELY;
225  a_nni->lk2          = UNLIKELY;
226  a_nni->l0           = -1.0;
227  a_nni->l1           = -1.0;
228  a_nni->l2           = -1.0;
229}
230
231//////////////////////////////////////////////////////////////
232//////////////////////////////////////////////////////////////
233
234void Init_Nexus_Format(nexcom **com)
235{
236
237  /*****************************/
238
239  strcpy(com[0]->name,"dimensions");
240  com[0]->nparm = 2;
241  com[0]->nxt_token_t = NEXUS_PARM;
242  com[0]->cur_token_t = NEXUS_COM;
243
244  com[0]->parm[0] = Make_Nexus_Parm();
245  strcpy(com[0]->parm[0]->name,"ntax");
246  com[0]->parm[0]->fp = Read_Nexus_Dimensions;
247  com[0]->parm[0]->com = com[0];
248  com[0]->parm[0]->nxt_token_t = NEXUS_EQUAL;
249  com[0]->parm[0]->cur_token_t = NEXUS_PARM;
250
251  com[0]->parm[1] = Make_Nexus_Parm();
252  strcpy(com[0]->parm[1]->name,"nchar");
253  com[0]->parm[1]->fp = Read_Nexus_Dimensions;
254  com[0]->parm[1]->com = com[0];
255  com[0]->parm[1]->nxt_token_t = NEXUS_EQUAL;
256  com[0]->parm[1]->cur_token_t = NEXUS_PARM;
257
258  /*****************************/
259
260  strcpy(com[1]->name,"format");
261  com[1]->nparm = 11;
262  com[1]->nxt_token_t = NEXUS_PARM;
263  com[1]->cur_token_t = NEXUS_COM;
264
265  com[1]->parm[0] = Make_Nexus_Parm();
266  strcpy(com[1]->parm[0]->name,"datatype");
267  com[1]->parm[0]->fp = Read_Nexus_Format;
268  com[1]->parm[0]->com = com[1];
269  com[1]->parm[0]->nxt_token_t = NEXUS_EQUAL;
270  com[1]->parm[0]->cur_token_t = NEXUS_PARM;
271
272  com[1]->parm[1] = Make_Nexus_Parm();
273  strcpy(com[1]->parm[1]->name,"respectcase");
274  com[1]->parm[1]->fp = Read_Nexus_Format;
275  com[1]->parm[1]->com = com[1];
276  com[1]->parm[1]->nxt_token_t = NEXUS_PARM;
277  com[1]->parm[1]->cur_token_t = NEXUS_VALUE;
278
279  com[1]->parm[2] = Make_Nexus_Parm();
280  strcpy(com[1]->parm[2]->name,"missing");
281  com[1]->parm[2]->fp = Read_Nexus_Format;
282  com[1]->parm[2]->com = com[1];
283  com[1]->parm[2]->nxt_token_t = NEXUS_EQUAL;
284  com[1]->parm[2]->cur_token_t = NEXUS_PARM;
285
286  com[1]->parm[3] = Make_Nexus_Parm();
287  strcpy(com[1]->parm[3]->name,"gap");
288  com[1]->parm[3]->fp = Read_Nexus_Format;
289  com[1]->parm[3]->com = com[1];
290  com[1]->parm[3]->nxt_token_t = NEXUS_EQUAL;
291  com[1]->parm[3]->cur_token_t = NEXUS_PARM;
292
293  com[1]->parm[4] = Make_Nexus_Parm();
294  strcpy(com[1]->parm[4]->name,"symbols");
295  com[1]->parm[4]->fp = Read_Nexus_Format;
296  com[1]->parm[4]->com = com[1];
297  com[1]->parm[4]->nxt_token_t = NEXUS_EQUAL;
298  com[1]->parm[4]->cur_token_t = NEXUS_PARM;
299
300  com[1]->parm[5] = Make_Nexus_Parm();
301  strcpy(com[1]->parm[5]->name,"equate");
302  com[1]->parm[5]->fp = Read_Nexus_Format;
303  com[1]->parm[5]->com = com[1];
304  com[1]->parm[5]->nxt_token_t = NEXUS_EQUAL;
305  com[1]->parm[5]->cur_token_t = NEXUS_PARM;
306
307  com[1]->parm[6] = Make_Nexus_Parm();
308  strcpy(com[1]->parm[6]->name,"matchchar");
309  com[1]->parm[6]->fp = Read_Nexus_Format;
310  com[1]->parm[6]->com = com[1];
311  com[1]->parm[6]->nxt_token_t = NEXUS_EQUAL;
312  com[1]->parm[6]->cur_token_t = NEXUS_PARM;
313
314  com[1]->parm[7] = Make_Nexus_Parm();
315  strcpy(com[1]->parm[7]->name,"transpose");
316  com[1]->parm[7]->fp = Read_Nexus_Format;
317  com[1]->parm[7]->com = com[1];
318  com[1]->parm[7]->nxt_token_t = NEXUS_PARM;
319  com[1]->parm[7]->cur_token_t = NEXUS_VALUE;
320
321  com[1]->parm[8] = Make_Nexus_Parm();
322  strcpy(com[1]->parm[8]->name,"interleave");
323  com[1]->parm[8]->fp = Read_Nexus_Format;
324  com[1]->parm[8]->com = com[1];
325  com[1]->parm[8]->nxt_token_t = NEXUS_PARM;
326  com[1]->parm[8]->cur_token_t = NEXUS_VALUE;
327
328  com[1]->parm[9] = Make_Nexus_Parm();
329  strcpy(com[1]->parm[9]->name,"items");
330  com[1]->parm[9]->fp = Read_Nexus_Format;
331  com[1]->parm[9]->com = com[1];
332  com[1]->parm[9]->nxt_token_t = NEXUS_EQUAL;
333  com[1]->parm[9]->cur_token_t = NEXUS_PARM;
334
335  com[1]->parm[10] = Make_Nexus_Parm();
336  strcpy(com[1]->parm[10]->name,"statesformat");
337  com[1]->parm[10]->fp = Read_Nexus_Format;
338  com[1]->parm[10]->com = com[1];
339  com[1]->parm[10]->nxt_token_t = NEXUS_EQUAL;
340  com[1]->parm[10]->cur_token_t = NEXUS_PARM;
341
342  /*****************************/
343
344  strcpy(com[2]->name,"eliminate");
345  com[2]->nparm = 0;
346  com[2]->nxt_token_t = NEXUS_VALUE;
347  com[2]->cur_token_t = NEXUS_COM;
348
349  /*****************************/
350
351  strcpy(com[3]->name,"taxlabels");
352  com[3]->nparm = 0;
353  com[3]->nxt_token_t = -1;
354  com[3]->cur_token_t = -1;
355 
356 /*****************************/
357
358  strcpy(com[4]->name,"charstatelabels");
359  com[4]->nparm = 0;
360  com[4]->nxt_token_t = -1;
361  com[4]->cur_token_t = -1;
362
363  /*****************************/
364
365  strcpy(com[5]->name,"charlabels");
366  com[5]->nparm = 0;
367  com[5]->nxt_token_t = -1;
368  com[5]->cur_token_t = -1;
369
370  /*****************************/
371
372  strcpy(com[6]->name,"statelabels");
373  com[6]->nparm = 0;
374  com[6]->nxt_token_t = -1;
375  com[6]->cur_token_t = -1;
376
377  /*****************************/
378
379  strcpy(com[7]->name,"matrix");
380  com[7]->nparm = 1;
381  com[7]->nxt_token_t = NEXUS_COM;
382  com[7]->cur_token_t = NEXUS_VALUE; /* This will allow us to skip directly
383                                        to the matrix reading function */
384
385  com[7]->parm[0] = Make_Nexus_Parm();
386  strcpy(com[7]->parm[0]->name,"matrix");
387  com[7]->parm[0]->fp = Read_Nexus_Matrix;
388  com[7]->parm[0]->com = com[7];
389  com[7]->parm[0]->nxt_token_t = NEXUS_COM;
390  com[7]->parm[0]->cur_token_t = -1; 
391
392  /*****************************/
393
394  strcpy(com[8]->name,"begin");
395  com[8]->nparm = 3;
396
397  com[8]->nxt_token_t = NEXUS_PARM;
398  com[8]->cur_token_t = NEXUS_COM;
399
400  com[8]->parm[0] = Make_Nexus_Parm();
401  strcpy(com[8]->parm[0]->name,"data");
402  com[8]->parm[0]->fp = Read_Nexus_Begin;
403  com[8]->parm[0]->com = com[8];
404  com[8]->parm[0]->nxt_token_t = NEXUS_COM;
405  com[8]->parm[0]->cur_token_t = NEXUS_PARM;
406
407
408  com[8]->parm[1] = Make_Nexus_Parm();
409  strcpy(com[8]->parm[1]->name,"trees");
410  com[8]->parm[1]->fp = Read_Nexus_Begin;
411  com[8]->parm[1]->com = com[8];
412  com[8]->parm[1]->nxt_token_t = NEXUS_COM;
413  com[8]->parm[1]->cur_token_t = NEXUS_PARM;
414
415
416  com[8]->parm[2] = Make_Nexus_Parm();
417  strcpy(com[8]->parm[2]->name,"taxa");
418  com[8]->parm[2]->fp = Read_Nexus_Taxa;
419  com[8]->parm[2]->com = com[8];
420  com[8]->parm[2]->nxt_token_t = NEXUS_COM;
421  com[8]->parm[2]->cur_token_t = NEXUS_VALUE; 
422
423
424  /*****************************/
425
426  strcpy(com[9]->name,"end");
427  com[9]->nparm = 0;
428  com[9]->nxt_token_t = -1;
429  com[9]->cur_token_t = -1;
430 
431  /*****************************/
432
433  strcpy(com[10]->name,"translate");
434  com[10]->nparm = 1;
435  com[10]->nxt_token_t = NEXUS_COM;
436  com[10]->cur_token_t = NEXUS_VALUE;
437
438  com[10]->parm[0] = Make_Nexus_Parm();
439  strcpy(com[10]->parm[0]->name,"translate");
440  com[10]->parm[0]->fp = Read_Nexus_Translate;
441  com[10]->parm[0]->com = com[10];
442  com[10]->parm[0]->nxt_token_t = NEXUS_COM;
443  com[10]->parm[0]->cur_token_t = -1; 
444
445  /*****************************/
446
447  strcpy(com[11]->name,"tree");
448  com[11]->nparm = 1;
449  com[11]->nxt_token_t = NEXUS_COM;
450  com[11]->cur_token_t = NEXUS_VALUE;
451
452  com[11]->parm[0] = Make_Nexus_Parm();
453  strcpy(com[11]->parm[0]->name,"tree");
454  com[11]->parm[0]->fp = Read_Nexus_Tree;
455  com[11]->parm[0]->com = com[11];
456  com[11]->parm[0]->nxt_token_t = -1;
457  com[11]->parm[0]->cur_token_t = -1; 
458
459
460  /*****************************/
461
462}
463
464//////////////////////////////////////////////////////////////
465//////////////////////////////////////////////////////////////
466
467void Init_Mat(matrix *mat, calign *data)
468{
469  int i;
470
471  mat->n_otu = data->n_otu;
472  mat->r = mat->n_otu;
473  mat->curr_int = mat->n_otu;
474  mat->method = 1;
475
476  For(i,data->n_otu)
477    {
478      strcpy(mat->name[i],data->c_seq[i]->name);
479      mat->on_off[i] = 1;
480    }
481}
482
483//////////////////////////////////////////////////////////////
484//////////////////////////////////////////////////////////////
485
486void Set_Defaults_Input(option* io)
487{
488  io->fp_in_align                = NULL;
489  io->fp_in_tree                 = NULL;
490  io->fp_in_constraint_tree      = NULL;
491  io->fp_out_tree                = NULL;
492  io->fp_out_trees               = NULL;
493  io->fp_out_boot_tree           = NULL;
494  io->fp_out_boot_stats          = NULL;
495  io->fp_out_stats               = NULL;
496  io->long_tax_names             = NULL;
497  io->short_tax_names            = NULL;
498  io->lon                        = NULL;
499  io->lat                        = NULL;
500  io->z_scores                   = NULL;
501  io->cstr_tree                  = NULL;
502  io->next                       = NULL;
503  io->prev                       = NULL;
504  io->tree                       = NULL;
505  io->mod                        = NULL;
506  strcpy(io->nt_or_cd,"nucleotides");
507  io->n_data_sets                = 1;
508  io->interleaved                = 1;
509  io->in_tree                    = 0;
510  io->out_tree_file_open_mode    = 1;
511  io->out_stats_file_open_mode   = 1;
512  io->init_len                   = -1;
513  io->n_otu                      = -1;
514  io->n_data_set_asked           = -1;
515  io->print_boot_trees           = 1;
516  io->n_part                     = 1;
517  io->ratio_test                 = ABAYES;
518  io->multigene                  = 0;
519  io->config_multigene           = 0;
520  io->curr_interface             = 0;
521  io->r_seed                     = -1;
522  io->collapse_boot              = 0;
523  io->random_boot_seq_order      = 1;
524  io->print_trace                = 0;
525  io->print_site_lnl             = 0;
526  io->m4_model                   = NO;
527  io->rm_ambigu                  = 0;
528  io->append_run_ID              = 0;
529  io->quiet                      = 0;
530  io->datatype                   = NT;
531  io->colalias                   = YES;
532  io->data_file_format           = PHYLIP;
533  io->tree_file_format           = PHYLIP;
534  io->boot_prog_every            = 20;
535  io->mem_question               = YES;
536  io->do_alias_subpatt           = NO;
537  io->lk_approx                  = EXACT;
538  io->codpos                     = -1;
539  io->mutmap                     = NO;
540  io->state_len                  = 1;
541
542  MCMC_Init_MCMC_Struct(NULL,io,io->mcmc);
543  RATES_Init_Rate_Struct(io->rates,NULL,-1);
544  io->rates->model               = GUINDON;
545}
546
547//////////////////////////////////////////////////////////////
548//////////////////////////////////////////////////////////////
549
550void Init_Rmat(t_rmat *rmat)
551{
552  rmat->n_diff_rr = 1;
553}
554
555//////////////////////////////////////////////////////////////
556//////////////////////////////////////////////////////////////
557
558void Set_Defaults_Model(t_mod *mod)
559{
560  strcpy(mod->modelname->s,"HKY85");
561  strcpy(mod->custom_mod_string->s,"000000");
562  mod->next                    = NULL;
563  mod->prev                    = NULL;
564  mod->next_mixt               = NULL;
565  mod->prev_mixt               = NULL;
566  mod->r_mat                   = NULL;
567  mod->e_frq                   = NULL;
568  mod->whichmodel              = HKY85;
569  mod->ras->n_catg             = 4;
570  mod->n_mixt_classes          = 0;
571  mod->mod_num                 = 0;
572  mod->update_eigen            = NO;
573  mod->is_mixt_mod             = NO;
574  mod->ras->normalise_rr       = YES;
575
576  mod->kappa->v                = 4.0;
577  mod->ras->alpha->v           = 1.0;
578  mod->lambda->v               = 1.0;
579  mod->ras->pinvar->v          = 0.0;
580  mod->l_var_sigma             = 1.E-2;
581  mod->e_frq_weight->v         = 1.0;
582  mod->r_mat_weight->v         = 1.0;
583
584  mod->bootstrap               = 0;
585  mod->ras->invar              = NO;
586  mod->ns                      = 4;
587  mod->use_m4mod               = NO;
588  mod->ras->gamma_median       = NO;
589  mod->m4mod                   = NULL;
590 
591  /* mod->r_mat->n_diff_rr        = 0; */
592  /* mod->r_mat->rr               = NULL; */
593  /* mod->r_mat->rr_val           = NULL; */
594  /* mod->r_mat->n_rr_per_cat     = NULL; */
595
596  mod->io                       = NULL;
597  mod->log_l                    = NO;
598  mod->ras->free_mixt_rates     = NO;
599  mod->gamma_mgf_bl             = NO;
600  mod->br_len_multiplier->v     = 1.0;
601
602 
603  mod->ras->parent_class_number = 0;
604  mod->ras->init_r_proba        = YES;
605  mod->ras->init_rr             = YES;
606
607
608#if !(defined PHYTIME || defined SERGEII)
609  mod->l_min = 1.E-8;
610  mod->l_max = 100.0;
611#else
612  mod->l_min = 1.E-8;
613  mod->l_max = 2.0;
614#endif
615
616  mod->l_var_min  = mod->l_min;
617  mod->l_var_max  = mod->l_max;
618
619
620}
621
622//////////////////////////////////////////////////////////////
623//////////////////////////////////////////////////////////////
624
625void Set_Defaults_Optimiz(t_opt *s_opt)
626{
627  s_opt->print                   = YES;
628  s_opt->last_opt                = YES;
629  s_opt->opt_subst_param         = YES;
630  s_opt->opt_alpha               = YES;
631  s_opt->opt_kappa               = YES;
632  s_opt->opt_bl                  = YES;
633  s_opt->opt_lambda              = NO;
634  s_opt->opt_pinvar              = NO;
635  s_opt->opt_cov_delta           = NO;
636  s_opt->opt_cov_alpha           = NO;
637  s_opt->opt_cov_free_rates      = NO;
638  s_opt->opt_rr                  = NO;
639  s_opt->init_lk                 = UNLIKELY;
640  s_opt->n_it_max                = 1000;
641  s_opt->opt_topo                = YES;
642  s_opt->topo_search             = NNI_MOVE;
643  s_opt->random_input_tree       = 0;
644  s_opt->n_rand_starts           = 5;
645  s_opt->brent_it_max            = BRENT_IT_MAX;
646  s_opt->steph_spr               = YES;
647  s_opt->user_state_freq         = NO;
648  s_opt->min_diff_lk_local       = 1.E-04;
649  s_opt->min_diff_lk_global      = 1.E-03;
650  s_opt->min_diff_lk_move        = 1.E-02;
651  s_opt->p_moves_to_examine      = 0.15;
652  s_opt->fast_nni                = NO;
653  s_opt->greedy                  = NO;
654  s_opt->general_pars            = NO;
655  s_opt->tree_size_mult          = 1;
656  s_opt->opt_five_branch         = YES;
657
658  s_opt->pars_thresh             = 5;
659
660  s_opt->hybrid_thresh             = NO;
661  s_opt->quickdirty                = NO;
662  s_opt->spr_pars                  = YES;
663  s_opt->spr_lnL                   = NO;
664  s_opt->min_depth_path            = 0;
665  s_opt->max_depth_path            = 20;
666  s_opt->deepest_path              = 20;
667  s_opt->max_delta_lnL_spr         = 50.;
668  s_opt->br_len_in_spr             = 10;
669  s_opt->opt_free_mixt_rates       = YES;
670  s_opt->constrained_br_len        = NO;
671  s_opt->opt_gamma_br_len          = NO;
672  s_opt->first_opt_free_mixt_rates = YES;
673
674  s_opt->wim_n_rgrft             = -1;
675  s_opt->wim_n_globl             = -1;
676  s_opt->wim_max_dist            = -1;
677  s_opt->wim_n_optim             = -1;
678  s_opt->wim_n_best              = -1;
679  s_opt->wim_inside_opt          =  0;
680
681  s_opt->opt_rmat_weight         = NO;
682  s_opt->opt_efrq_weight         = NO;
683
684  s_opt->skip_tree_traversal     = NO;
685  s_opt->serial_free_rates       = YES;
686
687  s_opt->curr_opt_free_rates     = NO;
688}
689
690//////////////////////////////////////////////////////////////
691//////////////////////////////////////////////////////////////
692
693void XML_Init_Attribute(xml_attr *attr)
694{
695}
696
697//////////////////////////////////////////////////////////////
698//////////////////////////////////////////////////////////////
699
700void XML_Init_Node(xml_node *parent, xml_node *new_node, char *name)
701{
702  if(name) strcpy(new_node->name,name);
703
704  new_node->parent   = parent ? parent : NULL;
705  new_node->next     = NULL;
706  new_node->prev     = NULL;
707  new_node->child    = NULL;
708  new_node->ds->obj  = NULL;
709  new_node->ds->next = NULL;
710 
711  if(parent)
712    { 
713      if(!parent->child)
714        {
715          parent->child = new_node;
716        }
717      else
718        {
719          xml_node *node = parent->child;
720          while(node->next) node = node->next;
721          node->next = new_node;
722          new_node->prev = node;
723        }
724    }
725 
726  new_node->attr = NULL;
727}
728
729//////////////////////////////////////////////////////////////
730//////////////////////////////////////////////////////////////
731
732void RATES_Init_Rate_Struct(t_rate *rates, t_rate *existing_rates, int n_otu)
733{
734  int i;
735
736  if(existing_rates && (existing_rates->model != -1))
737    {
738      rates->model = existing_rates->model;
739    }
740  else
741    {
742      rates->model = NONE;
743    }
744
745
746  if(rates->model == NONE)
747    rates->model_log_rates = NO;
748  else if(rates->model == THORNE)
749    rates->model_log_rates = YES;
750  else if(rates->model == GUINDON)
751    rates->model_log_rates = YES;
752  else if(rates->model == GAMMA)
753    rates->model_log_rates = NO;
754  else if(rates->model == STRICTCLOCK)
755    rates->model_log_rates = NO;
756  else
757    {
758      PhyML_Printf("\n== Please initialize model properly.");
759      PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
760      Warn_And_Exit("");
761    }
762
763
764  rates->met_within_gibbs = NO;
765  rates->c_lnL_rates      = UNLIKELY;
766  rates->c_lnL_jps        = UNLIKELY;
767  rates->adjust_rates     = 0;
768  rates->use_rates        = 1;
769  rates->lexp             = 1.E-3;
770  rates->norm_fact        = 1.0;
771  rates->inflate_var      = 1.0;
772  rates->nd_t_recorded    = NO;
773  rates->br_r_recorded    = NO;
774
775  rates->birth_rate       = 1.E-2;
776  rates->birth_rate_min   = 1.E-6;
777  rates->birth_rate_max   = 1.E+1;
778
779  if(rates->model_log_rates == YES)
780    {
781      rates->max_rate  =  LOG(10.);
782      rates->min_rate  = -LOG(10.);
783      /* rates->max_rate  =  MDBL_MAX; */
784      /* rates->min_rate  = -MDBL_MAX; */
785    }
786  else
787    {
788      rates->max_rate  = 10.0;
789      rates->min_rate  = 0.0;
790    }
791  /* rates->max_rate         = 6.0; */
792  /* rates->min_rate         = 0.0; */
793
794
795  rates->clock_r       = 1.E-4;
796  rates->min_clock     = 1.E-8;
797  rates->max_clock     = 1.E-3;
798
799  /* rates->clock_r       = 3.E-4; */
800  /* rates->max_clock     = 1.E-3; */
801  /* rates->min_clock     = 1.E-5; */
802
803  if(rates->model != GAMMA)
804    {
805      rates->nu            = 1.E-3;
806      rates->min_nu        = 0.0;
807      rates->max_nu        = 1.0;
808    }
809  else
810    {
811      rates->nu            = 1.E-1;
812      rates->min_nu        = 1.E-2;
813      rates->max_nu        = 100.0;
814    }
815
816  /* rates->max_nu        = 2.0; */
817
818  /* rates->nu            = 1.E-4; */
819  /* rates->max_nu        = 1.E-1; */
820  /* rates->min_nu        = 1.E-5; */
821
822  rates->min_dt        = 0.0;
823
824  rates->step_rate     = 1.E-4;
825  rates->approx        = 1;
826  rates->bl_from_rt    = 0;
827
828  rates->update_mean_l = NO;
829  rates->update_cov_l  = NO;
830
831  rates->p_max         = 0.01;
832
833  rates->true_tree_size = 0.0;
834
835  if(n_otu > 0)
836    {
837      For(i,(2*n_otu-2)*(2*n_otu-2)) rates->cov_l[i] = 0.0;
838     
839      For(i,2*n_otu-2) 
840        {
841          rates->n_jps[i]  =  -1;
842          rates->t_jps[i]  =  -1;
843          rates->mean_r[i] = 1.0;     
844          rates->mean_l[i] = 0.0;     
845        }
846     
847      For(i,2*n_otu-1) 
848        {
849          if(rates->model_log_rates == YES)
850            {
851              rates->nd_r[i]   = 0.0;
852              rates->br_r[i]   = 0.0;
853            }
854          else
855            {
856              rates->nd_r[i]   = 1.0;
857              rates->br_r[i]   = 1.0;
858            }
859
860          rates->mean_t[i] = 0.0;
861          rates->nd_t[i]   = 0.0;
862          rates->true_t[i] = 0.0;
863          if(i < n_otu)
864            {
865              rates->t_has_prior[i] = YES;
866              rates->t_prior_max[i] = 0.0;
867              rates->t_prior_min[i] = 0.0;
868            }
869          else
870            {
871              rates->t_has_prior[i] = NO;
872              rates->t_prior_max[i] =  BIG;
873              rates->t_prior_min[i] = -BIG;
874            }
875
876          rates->br_do_updt[i] = YES;
877          rates->has_survived[i] = NO;
878         
879          rates->t_ranked[i] = i;
880        }
881    }
882
883  rates->calib = NULL;
884 
885}
886
887//////////////////////////////////////////////////////////////
888//////////////////////////////////////////////////////////////
889
890void Init_One_Spr(t_spr *a_spr)
891{
892  a_spr->lnL             = UNLIKELY;
893  a_spr->pars            = 1E+5;
894  a_spr->depth_path      = 0;
895  a_spr->dist            = 0;
896  a_spr->init_target_l   = -1.;
897  a_spr->init_target_v   = -1.;
898  a_spr->l0              = -1.;
899  a_spr->l1              = -1.;
900  a_spr->l2              = -1.;
901  a_spr->v0              = -1.;
902  a_spr->v1              = -1.;
903  a_spr->v2              = -1.;
904  a_spr->n_link          = NULL;
905  a_spr->n_opp_to_link   = NULL;
906  a_spr->b_opp_to_link   = NULL;
907  a_spr->b_target        = NULL;
908  a_spr->b_init_target   = NULL;
909  a_spr->next            = NULL;
910  a_spr->prev            = NULL;
911  a_spr->next           = NULL;
912  a_spr->prev          = NULL;
913}
914
915//////////////////////////////////////////////////////////////
916//////////////////////////////////////////////////////////////
917
918void Init_Model(calign *data, t_mod *mod, option *io)
919{
920  int i,j;
921  phydbl sum,aux;
922  int result;
923  phydbl *dr, *di, *space;
924
925  mod->ns = io->mod->ns;
926
927  if(io->datatype == GENERIC) mod->whichmodel = JC69;
928
929  /* if(!mod->ras->invar) For(i,data->crunch_len) data->invar[i] = 0; */
930
931  dr      = (phydbl *)mCalloc(  mod->ns,sizeof(phydbl));
932  di      = (phydbl *)mCalloc(  mod->ns,sizeof(phydbl));
933  space   = (phydbl *)mCalloc(2*mod->ns,sizeof(phydbl));
934 
935  if(mod->log_l == YES)
936    {
937      mod->l_min = LOG(mod->l_min);
938      mod->l_max = LOG(mod->l_max);
939    }
940
941  if(mod->ras->init_r_proba == YES)
942    {
943      For(i,mod->ras->n_catg) mod->ras->gamma_r_proba->v[i]          = (phydbl)1./mod->ras->n_catg;
944      For(i,mod->ras->n_catg) mod->ras->gamma_r_proba_unscaled->v[i] = (phydbl)(i+1);
945    }
946  else
947    {
948      mod->ras->gamma_r_proba_unscaled->v[0] = mod->ras->gamma_r_proba->v[0];
949      for(i=1;i<mod->ras->n_catg;i++)
950        mod->ras->gamma_r_proba_unscaled->v[i] = 
951          mod->ras->gamma_r_proba_unscaled->v[i-1] +
952          mod->ras->gamma_r_proba->v[i];
953    }
954
955  if(mod->ras->init_rr == YES)
956    {
957      if(mod->ras->n_catg > 1)
958        {
959          For(i,mod->ras->n_catg) mod->ras->gamma_rr->v[i]               = (phydbl)i;
960          For(i,mod->ras->n_catg) mod->ras->gamma_rr_unscaled->v[i]      = (phydbl)i;
961        }
962      else
963        {
964          mod->ras->gamma_rr->v[0]          = 1.0;
965          mod->ras->gamma_rr_unscaled->v[0] = 1.0;
966        }
967    }
968
969
970
971
972  mod->br_len_multiplier->v = 1.0;
973
974  For(i,mod->ns) 
975    {
976      mod->e_frq->pi->v[i] = data->b_frq[i];     
977      mod->e_frq->pi_unscaled->v[i] = mod->e_frq->pi->v[i] * 100.;
978    }
979 
980  if(io->datatype == NT)
981    {
982      /* Set the substitution parameters to their default values
983         if they are not fixed by the user */
984      if(mod->s_opt->opt_kappa == YES) 
985        {
986          mod->kappa->v  = 4.0;
987          mod->lambda->v = 1.0;
988        }
989
990      /* if(mod->s_opt->opt_rr) */
991      if(mod->whichmodel == GTR || mod->whichmodel == CUSTOM)
992        {
993          int i;
994
995          For(i,6) 
996            {
997              mod->r_mat->rr->v[i]     = 1.0;
998              mod->r_mat->rr_val->v[i] = 1.0;
999            }
1000        }
1001    }
1002
1003  if(mod->s_opt->opt_alpha)  mod->ras->alpha->v  = 1.0;
1004  if(mod->s_opt->opt_pinvar) mod->ras->pinvar->v = 0.2;
1005 
1006
1007  if(io->datatype == NT) /* Nucleotides */
1008    { 
1009      /* init for nucleotides */
1010      mod->lambda->v    = 1.;
1011      /* mod->update_eigen = YES; */
1012     
1013      if(mod->whichmodel == JC69)
1014        {
1015          mod->e_frq->pi->v[0] = mod->e_frq->pi->v[1] = mod->e_frq->pi->v[2] = mod->e_frq->pi->v[3] = .25;
1016          mod->kappa->v = 1.;
1017          mod->s_opt->opt_state_freq = NO;
1018          mod->s_opt->opt_kappa      = NO;
1019          mod->s_opt->opt_lambda     = NO;
1020          mod->update_eigen          = NO;
1021        }
1022     
1023      if(mod->whichmodel == K80)
1024        {
1025          mod->e_frq->pi->v[0] = mod->e_frq->pi->v[1] = mod->e_frq->pi->v[2] = mod->e_frq->pi->v[3] = .25;
1026          mod->s_opt->opt_state_freq = NO;
1027          mod->s_opt->opt_lambda     = NO;
1028          mod->update_eigen          = NO;
1029        }
1030           
1031      if(mod->whichmodel == F81)
1032        {
1033          mod->kappa->v = 1.;
1034          mod->update_eigen          = NO;
1035        }
1036
1037      if(mod->whichmodel == F84)
1038        {
1039          aux = ((mod->e_frq->pi->v[0]+mod->e_frq->pi->v[2])-(mod->e_frq->pi->v[1]+mod->e_frq->pi->v[3]))/(2.*mod->kappa->v);
1040          mod->lambda->v = ((mod->e_frq->pi->v[1]+mod->e_frq->pi->v[3]) + aux)/((mod->e_frq->pi->v[0]+mod->e_frq->pi->v[2]) - aux); 
1041          mod->update_eigen          = NO;
1042        }
1043
1044      if(mod->whichmodel == TN93)
1045        {
1046          mod->update_eigen          = NO;
1047          if(io->mod->s_opt->opt_kappa) io->mod->s_opt->opt_lambda = 1;
1048        }
1049
1050      if(mod->whichmodel == GTR)
1051        {
1052          mod->kappa->v = 1.;
1053          mod->update_eigen      = YES;
1054          io->mod->s_opt->opt_rr = YES;
1055        }
1056
1057      if(mod->whichmodel == CUSTOM)
1058        {
1059          mod->kappa->v = 1.;
1060          mod->update_eigen = YES;
1061          /*      io->mod->s_opt->opt_rr     = YES; */ /* What if the user decided not to optimise the rates? */
1062        }
1063             
1064      if(mod->whichmodel == GTR)
1065        {                 
1066          mod->custom_mod_string->s[0] = '0';
1067          mod->custom_mod_string->s[1] = '1';
1068          mod->custom_mod_string->s[2] = '2';
1069          mod->custom_mod_string->s[3] = '3';
1070          mod->custom_mod_string->s[4] = '4';
1071          mod->custom_mod_string->s[5] = '5';
1072          Translate_Custom_Mod_String(mod);
1073        }
1074     
1075      if(mod->s_opt->user_state_freq == YES) 
1076        {
1077          For(i,4) 
1078            {
1079              mod->e_frq->pi->v[i] = mod->user_b_freq->v[i];
1080            }         
1081        }
1082
1083      if(((mod->whichmodel == GTR)    || 
1084          (mod->whichmodel == CUSTOM) || 
1085          (mod->whichmodel == HKY85)) &&
1086         (mod->use_m4mod == NO))
1087        {
1088          mod->update_eigen = YES;
1089          Update_Eigen(mod);
1090          mod->update_eigen = NO;
1091        }
1092      /* if((mod->whichmodel != GTR)    &&  */
1093      /*         (mod->whichmodel != CUSTOM) &&  */
1094      /*         (mod->whichmodel != HKY85)) mod->update_eigen = NO; */
1095    }
1096  else if(mod->io->datatype == AA)
1097    { 
1098
1099      /* init for amino-acids */
1100      /* see comments of PMat_Empirical for details */
1101      /* read pi and Q from file */
1102           
1103      /* These initialisations are needed when analysing multiple
1104       * data sets
1105       */
1106      For(i,mod->ns*mod->ns) mod->r_mat->qmat->v[i] = .0;
1107      For(i,mod->ns        ) mod->e_frq->pi->v[i]   = .0;
1108
1109      switch(mod->whichmodel)
1110        {
1111        case DAYHOFF : 
1112          {
1113            Init_Qmat_Dayhoff(mod->r_mat->qmat->v,mod->e_frq->pi->v);
1114            if(mod->s_opt->opt_state_freq)
1115              For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i];
1116            break;
1117          }
1118        case JTT : 
1119          {
1120            Init_Qmat_JTT(mod->r_mat->qmat->v,mod->e_frq->pi->v);
1121            if(mod->s_opt->opt_state_freq)
1122              For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i];
1123            break;
1124          }
1125        case MTREV : 
1126          {
1127            Init_Qmat_MtREV(mod->r_mat->qmat->v,mod->e_frq->pi->v);
1128            if(mod->s_opt->opt_state_freq)
1129              For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i];
1130            break;
1131          }
1132        case LG : 
1133          {
1134            Init_Qmat_LG(mod->r_mat->qmat->v,mod->e_frq->pi->v);
1135            if(mod->s_opt->opt_state_freq)
1136              For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i];         
1137            break;
1138          }
1139        case WAG : 
1140          {
1141            Init_Qmat_WAG(mod->r_mat->qmat->v,mod->e_frq->pi->v);
1142            if(mod->s_opt->opt_state_freq)
1143              For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i];
1144            break;
1145          }
1146        case DCMUT : 
1147          {
1148            Init_Qmat_DCMut(mod->r_mat->qmat->v,mod->e_frq->pi->v);
1149            if(mod->s_opt->opt_state_freq)
1150              For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i];
1151            break;
1152          }
1153        case RTREV : 
1154          {
1155            Init_Qmat_RtREV(mod->r_mat->qmat->v,mod->e_frq->pi->v);
1156            if(mod->s_opt->opt_state_freq)
1157              For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i];
1158            break;
1159          }
1160        case CPREV : 
1161          {
1162            Init_Qmat_CpREV(mod->r_mat->qmat->v,mod->e_frq->pi->v);
1163            if(mod->s_opt->opt_state_freq)
1164              For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i];
1165            break;
1166          }
1167        case VT : 
1168          {
1169            Init_Qmat_VT(mod->r_mat->qmat->v,mod->e_frq->pi->v);
1170            if(mod->s_opt->opt_state_freq)
1171              For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i];
1172            break;
1173          }
1174        case BLOSUM62 : 
1175          {
1176            Init_Qmat_Blosum62(mod->r_mat->qmat->v,mod->e_frq->pi->v);
1177            if(mod->s_opt->opt_state_freq)
1178              For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i];
1179            break;
1180          }
1181        case MTMAM : 
1182          {
1183            Init_Qmat_MtMam(mod->r_mat->qmat->v,mod->e_frq->pi->v);
1184            if(mod->s_opt->opt_state_freq)
1185              For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i];
1186            break;
1187          }
1188        case MTART : 
1189          {
1190            Init_Qmat_MtArt(mod->r_mat->qmat->v,mod->e_frq->pi->v);
1191            if(mod->s_opt->opt_state_freq)
1192              For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i];
1193            break;
1194          }
1195        case HIVW : 
1196          {
1197            Init_Qmat_HIVw(mod->r_mat->qmat->v,mod->e_frq->pi->v);
1198            if(mod->s_opt->opt_state_freq)
1199              For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i];
1200            break;
1201          }
1202        case HIVB : 
1203          {
1204            Init_Qmat_HIVb(mod->r_mat->qmat->v,mod->e_frq->pi->v);
1205            if(mod->s_opt->opt_state_freq)
1206              For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i];
1207            break;
1208          }
1209        case CUSTOMAA :
1210          {
1211            Read_Qmat(mod->r_mat->qmat->v,mod->e_frq->pi->v,mod->fp_aa_rate_mat);
1212/*          Print_Qmat_AA(mod->r_mat->qmat->v,mod->e_frq->pi->v); */
1213            if(mod->s_opt->opt_state_freq) For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i];           
1214            break;
1215          }
1216        default : 
1217          {
1218            Init_Qmat_LG(mod->r_mat->qmat->v,mod->e_frq->pi->v);
1219            if(mod->s_opt->opt_state_freq)
1220              For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i];         
1221            break;
1222          }
1223        }
1224 
1225/*       /\* multiply the nth col of Q by the nth term of pi/100 just as in PAML *\/ */
1226      For(i,mod->ns) For(j,mod->ns) mod->r_mat->qmat->v[i*mod->ns+j] *= mod->e_frq->pi->v[j] / 100.0;
1227     
1228      /* compute diagonal terms of Q and mean rate mr = l/t */
1229      mod->mr->v= .0;
1230      For (i,mod->ns)
1231        {
1232          sum=.0;
1233          For(j, mod->ns) sum += mod->r_mat->qmat->v[i*mod->ns+j];
1234          mod->r_mat->qmat->v[i*mod->ns+i] = -sum;
1235          mod->mr->v += mod->e_frq->pi->v[i] * sum;
1236        }
1237
1238      /* scale imod->nstantaneous rate matrix so that mu=1 */
1239      For (i,mod->ns*mod->ns) mod->r_mat->qmat->v[i] /= mod->mr->v;
1240
1241      /* compute eigenvectors/values */
1242      result = 0;
1243
1244      For(i,mod->ns*mod->ns) mod->r_mat->qmat_buff->v[i] = mod->r_mat->qmat->v[i];
1245
1246      if(!Eigen(1,mod->r_mat->qmat_buff->v,mod->eigen->size,mod->eigen->e_val,
1247                mod->eigen->e_val_im,mod->eigen->r_e_vect,
1248                mod->eigen->r_e_vect_im,mod->eigen->space))
1249        {
1250          /* compute inverse(Vr) into Vi */
1251          For (i,mod->ns*mod->ns) mod->eigen->l_e_vect[i] = mod->eigen->r_e_vect[i];
1252          if(!Matinv(mod->eigen->l_e_vect,mod->eigen->size,mod->eigen->size,YES))
1253            {
1254              PhyML_Printf("\n== Err in file %s at line %d.",__FILE__,__LINE__);
1255              Exit("\n");     
1256            }
1257         
1258          /* compute the diagonal terms of EXP(D) */
1259          For(i,mod->ns) mod->eigen->e_val[i] = (phydbl)EXP(mod->eigen->e_val[i]);
1260        }
1261      else
1262        {
1263          if (result==-1) PhyML_Printf("\n== Eigenvalues/vectors computation does not converge : computation cancelled");
1264          else if (result==1) PhyML_Printf("\n== Complex eigenvalues/vectors : computation cancelled");
1265        }
1266    }
1267  else if(mod->io->datatype == GENERIC)
1268    {
1269      /* Uniform state frequencies */
1270      For(i,mod->ns)  mod->e_frq->pi->v[i] = 1./(phydbl)mod->ns;
1271      mod->kappa->v = 1;
1272      mod->s_opt->opt_state_freq = NO;
1273      mod->s_opt->opt_kappa      = NO;
1274      mod->s_opt->opt_lambda     = NO;
1275      mod->update_eigen          = NO;
1276    }
1277  else
1278    {
1279      PhyML_Printf("\n== Datatype not recognized.\n");
1280      PhyML_Printf("\n== Err. in file %s at line %d\n",__FILE__,__LINE__);
1281      Warn_And_Exit("");
1282    }
1283
1284  if(!mod->use_m4mod) Set_Model_Parameters(mod);     
1285                 
1286  Init_Eigen_Struct(mod->eigen);
1287 
1288  free(dr);free(di);free(space);
1289}
1290
1291//////////////////////////////////////////////////////////////
1292//////////////////////////////////////////////////////////////
1293
1294int Init_Qmat_Dayhoff(phydbl *daa, phydbl *pi)
1295{
1296  /* Dayhoff's model data
1297   * Dayhoff, M.O., Schwartz, R.M., Orcutt, B.C. (1978)
1298   * "A model of evolutionary change in proteins."
1299   * Dayhoff, M.O.(ed.) Atlas of Protein Sequence Structur., Vol5, Suppl3.
1300   * National Biomedical Research Foundation, Washington DC, pp.345-352.
1301   */
1302  int i,j,naa;
1303
1304  naa = 20;
1305
1306/*   PhyML_Printf("\n\n. REMINDER : THIS IS NOT DAYHOFF !!!\n\n"); */
1307 
1308/*   daa[1*20 + 0] = 0.538903;  */
1309/*   daa[2*20 + 0] = 0.412504;  */
1310/*   daa[2*20 + 1] = 0.736081;  */
1311/*   daa[3*20 + 0] = 0.586915;  */
1312/*   daa[3*20 + 1] = 0.108051;  */
1313/*   daa[3*20 + 2] = 5.446642;  */
1314/*   daa[4*20 + 0] = 2.189718;  */
1315/*   daa[4*20 + 1] = 0.830604;  */
1316/*   daa[4*20 + 2] = 0.573426;  */
1317/*   daa[4*20 + 3] = 0.077565;  */
1318/*   daa[5*20 + 0] = 1.08213;  */
1319/*   daa[5*20 + 1] = 2.950693;  */
1320/*   daa[5*20 + 2] = 1.739514;  */
1321/*   daa[5*20 + 3] = 0.559035;  */
1322/*   daa[5*20 + 4] = 0.111314;  */
1323/*   daa[6*20 + 0] = 1.386865;  */
1324/*   daa[6*20 + 1] = 0.358434;  */
1325/*   daa[6*20 + 2] = 0.541447;  */
1326/*   daa[6*20 + 3] = 5.406871;  */
1327/*   daa[6*20 + 4] = 0.003738;  */
1328/*   daa[6*20 + 5] = 4.124043;  */
1329/*   daa[7*20 + 0] = 2.085747;  */
1330/*   daa[7*20 + 1] = 0.453132;  */
1331/*   daa[7*20 + 2] = 1.815844;  */
1332/*   daa[7*20 + 3] = 1.08647;  */
1333/*   daa[7*20 + 4] = 0.526418;  */
1334/*   daa[7*20 + 5] = 0.347693;  */
1335/*   daa[7*20 + 6] = 0.438476;  */
1336/*   daa[8*20 + 0] = 0.407898;  */
1337/*   daa[8*20 + 1] = 2.689322;  */
1338/*   daa[8*20 + 2] = 5.386808;  */
1339/*   daa[8*20 + 3] = 0.884563;  */
1340/*   daa[8*20 + 4] = 0.658583;  */
1341/*   daa[8*20 + 5] = 4.588358;  */
1342/*   daa[8*20 + 6] = 0.386218;  */
1343/*   daa[8*20 + 7] = 0.368668;  */
1344/*   daa[9*20 + 0] = 0.10177;  */
1345/*   daa[9*20 + 1] = 0.104875;  */
1346/*   daa[9*20 + 2] = 0.16239;  */
1347/*   daa[9*20 + 3] = 0.011698;  */
1348/*   daa[9*20 + 4] = 0.253282;  */
1349/*   daa[9*20 + 5] = 0.083872;  */
1350/*   daa[9*20 + 6] = 0.041767;  */
1351/*   daa[9*20 + 7] = 0.009778;  */
1352/*   daa[9*20 + 8] = 0.1042;  */
1353/*   daa[10*20 + 0] = 0.248202;  */
1354/*   daa[10*20 + 1] = 0.375742;  */
1355/*   daa[10*20 + 2] = 0.093863;  */
1356/*   daa[10*20 + 3] = 0.019241;  */
1357/*   daa[10*20 + 4] = 0.572688;  */
1358/*   daa[10*20 + 5] = 0.703538;  */
1359/*   daa[10*20 + 6] = 0.071961;  */
1360/*   daa[10*20 + 7] = 0.03006;  */
1361/*   daa[10*20 + 8] = 0.418222;  */
1362/*   daa[10*20 + 9] = 3.702051;  */
1363/*   daa[11*20 + 0] = 0.652019;  */
1364/*   daa[11*20 + 1] = 5.940478;  */
1365/*   daa[11*20 + 2] = 2.352253;  */
1366/*   daa[11*20 + 3] = 0.231001;  */
1367/*   daa[11*20 + 4] = 0.027995;  */
1368/*   daa[11*20 + 5] = 3.646743;  */
1369/*   daa[11*20 + 6] = 1.507981;  */
1370/*   daa[11*20 + 7] = 0.331175;  */
1371/*   daa[11*20 + 8] = 0.698362;  */
1372/*   daa[11*20 + 9] = 0.140326;  */
1373/*   daa[11*20 + 10] = 0.171396;  */
1374/*   daa[12*20 + 0] = 0.694226;  */
1375/*   daa[12*20 + 1] = 0.419899;  */
1376/*   daa[12*20 + 2] = 0.326927;  */
1377/*   daa[12*20 + 3] = 0.039488;  */
1378/*   daa[12*20 + 4] = 0.844827;  */
1379/*   daa[12*20 + 5] = 1.394214;  */
1380/*   daa[12*20 + 6] = 0.133235;  */
1381/*   daa[12*20 + 7] = 0.085075;  */
1382/*   daa[12*20 + 8] = 0.347092;  */
1383/*   daa[12*20 + 9] = 4.051255;  */
1384/*   daa[12*20 + 10] = 6.650794;  */
1385/*   daa[12*20 + 11] = 0.617549;  */
1386/*   daa[13*20 + 0] = 0.155206;  */
1387/*   daa[13*20 + 1] = 0.057971;  */
1388/*   daa[13*20 + 2] = 0.09816;  */
1389/*   daa[13*20 + 3] = 0.020441;  */
1390/*   daa[13*20 + 4] = 0.904305;  */
1391/*   daa[13*20 + 5] = 0.052719;  */
1392/*   daa[13*20 + 6] = 0.0219;  */
1393/*   daa[13*20 + 7] = 0.046668;  */
1394/*   daa[13*20 + 8] = 0.890005;  */
1395/*   daa[13*20 + 9] = 0.844963;  */
1396/*   daa[13*20 + 10] = 2.348881;  */
1397/*   daa[13*20 + 11] = 0.028372;  */
1398/*   daa[13*20 + 12] = 1.671635;  */
1399/*   daa[14*20 + 0] = 1.433475;  */
1400/*   daa[14*20 + 1] = 0.328393;  */
1401/*   daa[14*20 + 2] = 0.173181;  */
1402/*   daa[14*20 + 3] = 0.431874;  */
1403/*   daa[14*20 + 4] = 0.09902;  */
1404/*   daa[14*20 + 5] = 0.592324;  */
1405/*   daa[14*20 + 6] = 0.488352;  */
1406/*   daa[14*20 + 7] = 0.23865;  */
1407/*   daa[14*20 + 8] = 0.462856;  */
1408/*   daa[14*20 + 9] = 0.057048;  */
1409/*   daa[14*20 + 10] = 0.233532;  */
1410/*   daa[14*20 + 11] = 0.387808;  */
1411/*   daa[14*20 + 12] = 0.096377;  */
1412/*   daa[14*20 + 13] = 0.079912;  */
1413/*   daa[15*20 + 0] = 4.887126;  */
1414/*   daa[15*20 + 1] = 0.883923;  */
1415/*   daa[15*20 + 2] = 4.627163;  */
1416/*   daa[15*20 + 3] = 1.122164;  */
1417/*   daa[15*20 + 4] = 3.186667;  */
1418/*   daa[15*20 + 5] = 1.085947;  */
1419/*   daa[15*20 + 6] = 0.569339;  */
1420/*   daa[15*20 + 7] = 1.993432;  */
1421/*   daa[15*20 + 8] = 0.867972;  */
1422/*   daa[15*20 + 9] = 0.070512;  */
1423/*   daa[15*20 + 10] = 0.163009;  */
1424/*   daa[15*20 + 11] = 0.718913;  */
1425/*   daa[15*20 + 12] = 0.301103;  */
1426/*   daa[15*20 + 13] = 0.32579;  */
1427/*   daa[15*20 + 14] = 1.449582;  */
1428/*   daa[16*20 + 0] = 2.030538;  */
1429/*   daa[16*20 + 1] = 0.639463;  */
1430/*   daa[16*20 + 2] = 2.076294;  */
1431/*   daa[16*20 + 3] = 0.377239;  */
1432/*   daa[16*20 + 4] = 1.42848;  */
1433/*   daa[16*20 + 5] = 0.979403;  */
1434/*   daa[16*20 + 6] = 0.647562;  */
1435/*   daa[16*20 + 7] = 0.145556;  */
1436/*   daa[16*20 + 8] = 0.493329;  */
1437/*   daa[16*20 + 9] = 0.973405;  */
1438/*   daa[16*20 + 10] = 0.271824;  */
1439/*   daa[16*20 + 11] = 1.20033;  */
1440/*   daa[16*20 + 12] = 1.659187;  */
1441/*   daa[16*20 + 13] = 0.1217;  */
1442/*   daa[16*20 + 14] = 0.571399;  */
1443/*   daa[16*20 + 15] = 6.641034;  */
1444/*   daa[17*20 + 0] = 0.131405;  */
1445/*   daa[17*20 + 1] = 0.552911;  */
1446/*   daa[17*20 + 2] = 0.079985;  */
1447/*   daa[17*20 + 3] = 0.060514;  */
1448/*   daa[17*20 + 4] = 0.633662;  */
1449/*   daa[17*20 + 5] = 0.21823;  */
1450/*   daa[17*20 + 6] = 0.074988;  */
1451/*   daa[17*20 + 7] = 0.169114;  */
1452/*   daa[17*20 + 8] = 0.847725;  */
1453/*   daa[17*20 + 9] = 0.10627;  */
1454/*   daa[17*20 + 10] = 0.622044;  */
1455/*   daa[17*20 + 11] = 0.060755;  */
1456/*   daa[17*20 + 12] = 0.719575;  */
1457/*   daa[17*20 + 13] = 3.14824;  */
1458/*   daa[17*20 + 14] = 0.077123;  */
1459/*   daa[17*20 + 15] = 0.276716;  */
1460/*   daa[17*20 + 16] = 0.148883;  */
1461/*   daa[18*20 + 0] = 0.165179;  */
1462/*   daa[18*20 + 1] = 0.224883;  */
1463/*   daa[18*20 + 2] = 0.528334;  */
1464/*   daa[18*20 + 3] = 0.121252;  */
1465/*   daa[18*20 + 4] = 1.174118;  */
1466/*   daa[18*20 + 5] = 0.177062;  */
1467/*   daa[18*20 + 6] = 0.074715;  */
1468/*   daa[18*20 + 7] = 0.042356;  */
1469/*   daa[18*20 + 8] = 5.911187;  */
1470/*   daa[18*20 + 9] = 0.192481;  */
1471/*   daa[18*20 + 10] = 0.321454;  */
1472/*   daa[18*20 + 11] = 0.090556;  */
1473/*   daa[18*20 + 12] = 0.406415;  */
1474/*   daa[18*20 + 13] = 10.908861;  */
1475/*   daa[18*20 + 14] = 0.070752;  */
1476/*   daa[18*20 + 15] = 0.328483;  */
1477/*   daa[18*20 + 16] = 0.181539;  */
1478/*   daa[18*20 + 17] = 3.823886;  */
1479/*   daa[19*20 + 0] = 1.78517;  */
1480/*   daa[19*20 + 1] = 0.166975;  */
1481/*   daa[19*20 + 2] = 0.106482;  */
1482/*   daa[19*20 + 3] = 0.041707;  */
1483/*   daa[19*20 + 4] = 1.876812;  */
1484/*   daa[19*20 + 5] = 0.22421;  */
1485/*   daa[19*20 + 6] = 0.247356;  */
1486/*   daa[19*20 + 7] = 0.06688;  */
1487/*   daa[19*20 + 8] = 0.1436;  */
1488/*   daa[19*20 + 9] = 9.60184;  */
1489/*   daa[19*20 + 10] = 1.599119;  */
1490/*   daa[19*20 + 11] = 0.17319;  */
1491/*   daa[19*20 + 12] = 1.645134;  */
1492/*   daa[19*20 + 13] = 0.438571;  */
1493/*   daa[19*20 + 14] = 0.252486;  */
1494/*   daa[19*20 + 15] = 0.105536;  */
1495/*   daa[19*20 + 16] = 1.789097;  */
1496/*   daa[19*20 + 17] = 0.147951;  */
1497/*   daa[19*20 + 18] = 0.200571; */
1498 
1499 
1500   
1501/*   for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j]; */
1502 
1503 
1504/*   pi[0] = 0.093862;  */
1505/*   pi[1] = 0.05757;  */
1506/*   pi[2] = 0.037017;  */
1507/*   pi[3] = 0.060952;  */
1508/*   pi[4] = 0.011004;  */
1509/*   pi[5] = 0.044631;  */
1510/*   pi[6] = 0.083648;  */
1511/*   pi[7] = 0.053004;  */
1512/*   pi[8] = 0.022528;  */
1513/*   pi[9] = 0.060152;  */
1514/*   pi[10] = 0.092314;  */
1515/*   pi[11] = 0.065886;  */
1516/*   pi[12] = 0.021075;  */
1517/*   pi[13] = 0.033928;  */
1518/*   pi[14] = 0.043333;  */
1519/*   pi[15] = 0.052957;  */
1520/*   pi[16] = 0.053465;  */
1521/*   pi[17] = 0.00928;  */
1522/*   pi[18] = 0.030025;  */
1523/*   pi[19] = 0.073369; */
1524 
1525 
1526  daa[ 1*20+ 0] =   27.00; daa[ 2*20+ 0] =   98.00; daa[ 2*20+ 1] =   32.00; daa[ 3*20+ 0] =  120.00;
1527  daa[ 3*20+ 1] =    0.00; daa[ 3*20+ 2] =  905.00; daa[ 4*20+ 0] =   36.00; daa[ 4*20+ 1] =   23.00;
1528  daa[ 4*20+ 2] =    0.00; daa[ 4*20+ 3] =    0.00; daa[ 5*20+ 0] =   89.00; daa[ 5*20+ 1] =  246.00;
1529  daa[ 5*20+ 2] =  103.00; daa[ 5*20+ 3] =  134.00; daa[ 5*20+ 4] =    0.00; daa[ 6*20+ 0] =  198.00;
1530  daa[ 6*20+ 1] =    1.00; daa[ 6*20+ 2] =  148.00; daa[ 6*20+ 3] = 1153.00; daa[ 6*20+ 4] =    0.00;
1531  daa[ 6*20+ 5] =  716.00; daa[ 7*20+ 0] =  240.00; daa[ 7*20+ 1] =    9.00; daa[ 7*20+ 2] =  139.00;
1532  daa[ 7*20+ 3] =  125.00; daa[ 7*20+ 4] =   11.00; daa[ 7*20+ 5] =   28.00; daa[ 7*20+ 6] =   81.00;
1533  daa[ 8*20+ 0] =   23.00; daa[ 8*20+ 1] =  240.00; daa[ 8*20+ 2] =  535.00; daa[ 8*20+ 3] =   86.00;
1534  daa[ 8*20+ 4] =   28.00; daa[ 8*20+ 5] =  606.00; daa[ 8*20+ 6] =   43.00; daa[ 8*20+ 7] =   10.00;
1535  daa[ 9*20+ 0] =   65.00; daa[ 9*20+ 1] =   64.00; daa[ 9*20+ 2] =   77.00; daa[ 9*20+ 3] =   24.00;
1536  daa[ 9*20+ 4] =   44.00; daa[ 9*20+ 5] =   18.00; daa[ 9*20+ 6] =   61.00; daa[ 9*20+ 7] =    0.00;
1537  daa[ 9*20+ 8] =    7.00; daa[10*20+ 0] =   41.00; daa[10*20+ 1] =   15.00; daa[10*20+ 2] =   34.00;
1538  daa[10*20+ 3] =    0.00; daa[10*20+ 4] =    0.00; daa[10*20+ 5] =   73.00; daa[10*20+ 6] =   11.00;
1539  daa[10*20+ 7] =    7.00; daa[10*20+ 8] =   44.00; daa[10*20+ 9] =  257.00; daa[11*20+ 0] =   26.00;
1540  daa[11*20+ 1] =  464.00; daa[11*20+ 2] =  318.00; daa[11*20+ 3] =   71.00; daa[11*20+ 4] =    0.00;
1541  daa[11*20+ 5] =  153.00; daa[11*20+ 6] =   83.00; daa[11*20+ 7] =   27.00; daa[11*20+ 8] =   26.00;
1542  daa[11*20+ 9] =   46.00; daa[11*20+10] =   18.00; daa[12*20+ 0] =   72.00; daa[12*20+ 1] =   90.00;
1543  daa[12*20+ 2] =    1.00; daa[12*20+ 3] =    0.00; daa[12*20+ 4] =    0.00; daa[12*20+ 5] =  114.00;
1544  daa[12*20+ 6] =   30.00; daa[12*20+ 7] =   17.00; daa[12*20+ 8] =    0.00; daa[12*20+ 9] =  336.00;
1545  daa[12*20+10] =  527.00; daa[12*20+11] =  243.00; daa[13*20+ 0] =   18.00; daa[13*20+ 1] =   14.00;
1546  daa[13*20+ 2] =   14.00; daa[13*20+ 3] =    0.00; daa[13*20+ 4] =    0.00; daa[13*20+ 5] =    0.00;
1547  daa[13*20+ 6] =    0.00; daa[13*20+ 7] =   15.00; daa[13*20+ 8] =   48.00; daa[13*20+ 9] =  196.00;
1548  daa[13*20+10] =  157.00; daa[13*20+11] =    0.00; daa[13*20+12] =   92.00; daa[14*20+ 0] =  250.00;
1549  daa[14*20+ 1] =  103.00; daa[14*20+ 2] =   42.00; daa[14*20+ 3] =   13.00; daa[14*20+ 4] =   19.00;
1550  daa[14*20+ 5] =  153.00; daa[14*20+ 6] =   51.00; daa[14*20+ 7] =   34.00; daa[14*20+ 8] =   94.00;
1551  daa[14*20+ 9] =   12.00; daa[14*20+10] =   32.00; daa[14*20+11] =   33.00; daa[14*20+12] =   17.00;
1552  daa[14*20+13] =   11.00; daa[15*20+ 0] =  409.00; daa[15*20+ 1] =  154.00; daa[15*20+ 2] =  495.00;
1553  daa[15*20+ 3] =   95.00; daa[15*20+ 4] =  161.00; daa[15*20+ 5] =   56.00; daa[15*20+ 6] =   79.00;
1554  daa[15*20+ 7] =  234.00; daa[15*20+ 8] =   35.00; daa[15*20+ 9] =   24.00; daa[15*20+10] =   17.00;
1555  daa[15*20+11] =   96.00; daa[15*20+12] =   62.00; daa[15*20+13] =   46.00; daa[15*20+14] =  245.00;
1556  daa[16*20+ 0] =  371.00; daa[16*20+ 1] =   26.00; daa[16*20+ 2] =  229.00; daa[16*20+ 3] =   66.00;
1557  daa[16*20+ 4] =   16.00; daa[16*20+ 5] =   53.00; daa[16*20+ 6] =   34.00; daa[16*20+ 7] =   30.00;
1558  daa[16*20+ 8] =   22.00; daa[16*20+ 9] =  192.00; daa[16*20+10] =   33.00; daa[16*20+11] =  136.00;
1559  daa[16*20+12] =  104.00; daa[16*20+13] =   13.00; daa[16*20+14] =   78.00; daa[16*20+15] =  550.00;
1560  daa[17*20+ 0] =    0.00; daa[17*20+ 1] =  201.00; daa[17*20+ 2] =   23.00; daa[17*20+ 3] =    0.00;
1561  daa[17*20+ 4] =    0.00; daa[17*20+ 5] =    0.00; daa[17*20+ 6] =    0.00; daa[17*20+ 7] =    0.00;
1562  daa[17*20+ 8] =   27.00; daa[17*20+ 9] =    0.00; daa[17*20+10] =   46.00; daa[17*20+11] =    0.00;
1563  daa[17*20+12] =    0.00; daa[17*20+13] =   76.00; daa[17*20+14] =    0.00; daa[17*20+15] =   75.00;
1564  daa[17*20+16] =    0.00; daa[18*20+ 0] =   24.00; daa[18*20+ 1] =    8.00; daa[18*20+ 2] =   95.00;
1565  daa[18*20+ 3] =    0.00; daa[18*20+ 4] =   96.00; daa[18*20+ 5] =    0.00; daa[18*20+ 6] =   22.00;
1566  daa[18*20+ 7] =    0.00; daa[18*20+ 8] =  127.00; daa[18*20+ 9] =   37.00; daa[18*20+10] =   28.00;
1567  daa[18*20+11] =   13.00; daa[18*20+12] =    0.00; daa[18*20+13] =  698.00; daa[18*20+14] =    0.00;
1568  daa[18*20+15] =   34.00; daa[18*20+16] =   42.00; daa[18*20+17] =   61.00; daa[19*20+ 0] =  208.00;
1569  daa[19*20+ 1] =   24.00; daa[19*20+ 2] =   15.00; daa[19*20+ 3] =   18.00; daa[19*20+ 4] =   49.00;
1570  daa[19*20+ 5] =   35.00; daa[19*20+ 6] =   37.00; daa[19*20+ 7] =   54.00; daa[19*20+ 8] =   44.00;
1571  daa[19*20+ 9] =  889.00; daa[19*20+10] =  175.00; daa[19*20+11] =   10.00; daa[19*20+12] =  258.00;
1572  daa[19*20+13] =   12.00; daa[19*20+14] =   48.00; daa[19*20+15] =   30.00; daa[19*20+16] =  157.00;
1573  daa[19*20+17] =    0.00; daa[19*20+18] =   28.00;
1574 
1575  for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j];
1576 
1577  pi[ 0] = 0.087127; pi[ 1] = 0.040904; pi[ 2] = 0.040432; pi[ 3] = 0.046872;
1578  pi[ 4] = 0.033474; pi[ 5] = 0.038255; pi[ 6] = 0.049530; pi[ 7] = 0.088612;
1579  pi[ 8] = 0.033618; pi[ 9] = 0.036886; pi[10] = 0.085357; pi[11] = 0.080482;
1580  pi[12] = 0.014753; pi[13] = 0.039772; pi[14] = 0.050680; pi[15] = 0.069577;
1581  pi[16] = 0.058542; pi[17] = 0.010494; pi[18] = 0.029916; pi[19] = 0.064718;
1582 
1583 
1584   
1585  return 1;
1586}
1587
1588//////////////////////////////////////////////////////////////
1589//////////////////////////////////////////////////////////////
1590
1591
1592int Init_Qmat_DCMut(phydbl *daa, phydbl *pi)
1593{
1594  /*
1595     DCMut : new implementation based on Dayhoff et al.'s raw data and amino acid mutabilities
1596     C. Kosiol and N. Goldman. (2005),
1597     ``Different versions of the Dayhoff rate matrix'',
1598     Mol. Biol. Evol., 22. 193-199.
1599  */
1600
1601  int i,j,naa;
1602
1603  naa = 20;
1604
1605  daa[ 1*20+ 0] =   26.78280; daa[ 2*20+ 0] =   98.44740; daa[ 2*20+ 1] =   32.70590; daa[ 3*20+ 0] =  119.98050; 
1606  daa[ 3*20+ 1] =    0.00000; daa[ 3*20+ 2] =  893.15150; daa[ 4*20+ 0] =   36.00160; daa[ 4*20+ 1] =   23.23740; 
1607  daa[ 4*20+ 2] =    0.00000; daa[ 4*20+ 3] =    0.00000; daa[ 5*20+ 0] =   88.77530; daa[ 5*20+ 1] =  243.99390; 
1608  daa[ 5*20+ 2] =  102.85090; daa[ 5*20+ 3] =  134.85510; daa[ 5*20+ 4] =    0.00000; daa[ 6*20+ 0] =  196.11670; 
1609  daa[ 6*20+ 1] =    0.00000; daa[ 6*20+ 2] =  149.34090; daa[ 6*20+ 3] = 1138.86590; daa[ 6*20+ 4] =    0.00000; 
1610  daa[ 6*20+ 5] =  708.60220; daa[ 7*20+ 0] =  238.61110; daa[ 7*20+ 1] =    8.77910; daa[ 7*20+ 2] =  138.53520; 
1611  daa[ 7*20+ 3] =  124.09810; daa[ 7*20+ 4] =   10.72780; daa[ 7*20+ 5] =   28.15810; daa[ 7*20+ 6] =   81.19070; 
1612  daa[ 8*20+ 0] =   22.81160; daa[ 8*20+ 1] =  238.31480; daa[ 8*20+ 2] =  529.00240; daa[ 8*20+ 3] =   86.82410; 
1613  daa[ 8*20+ 4] =   28.27290; daa[ 8*20+ 5] =  601.16130; daa[ 8*20+ 6] =   43.94690; daa[ 8*20+ 7] =   10.68020; 
1614  daa[ 9*20+ 0] =   65.34160; daa[ 9*20+ 1] =   63.26290; daa[ 9*20+ 2] =   76.80240; daa[ 9*20+ 3] =   23.92480; 
1615  daa[ 9*20+ 4] =   43.80740; daa[ 9*20+ 5] =   18.03930; daa[ 9*20+ 6] =   60.95260; daa[ 9*20+ 7] =    0.00000; 
1616  daa[ 9*20+ 8] =    7.69810; daa[10*20+ 0] =   40.64310; daa[10*20+ 1] =   15.49240; daa[10*20+ 2] =   34.11130; 
1617  daa[10*20+ 3] =    0.00000; daa[10*20+ 4] =    0.00000; daa[10*20+ 5] =   73.07720; daa[10*20+ 6] =   11.28800; 
1618  daa[10*20+ 7] =    7.15140; daa[10*20+ 8] =   44.35040; daa[10*20+ 9] =  255.66850; daa[11*20+ 0] =   25.86350; 
1619  daa[11*20+ 1] =  461.01240; daa[11*20+ 2] =  314.83710; daa[11*20+ 3] =   71.69130; daa[11*20+ 4] =    0.00000; 
1620  daa[11*20+ 5] =  151.90780; daa[11*20+ 6] =   83.00780; daa[11*20+ 7] =   26.76830; daa[11*20+ 8] =   27.04750; 
1621  daa[11*20+ 9] =   46.08570; daa[11*20+10] =   18.06290; daa[12*20+ 0] =   71.78400; daa[12*20+ 1] =   89.63210; 
1622  daa[12*20+ 2] =    0.00000; daa[12*20+ 3] =    0.00000; daa[12*20+ 4] =    0.00000; daa[12*20+ 5] =  112.74990; 
1623  daa[12*20+ 6] =   30.48030; daa[12*20+ 7] =   17.03720; daa[12*20+ 8] =    0.00000; daa[12*20+ 9] =  333.27320; 
1624  daa[12*20+10] =  523.01150; daa[12*20+11] =  241.17390; daa[13*20+ 0] =   18.36410; daa[13*20+ 1] =   13.69060; 
1625  daa[13*20+ 2] =   13.85030; daa[13*20+ 3] =    0.00000; daa[13*20+ 4] =    0.00000; daa[13*20+ 5] =    0.00000; 
1626  daa[13*20+ 6] =    0.00000; daa[13*20+ 7] =   15.34780; daa[13*20+ 8] =   47.59270; daa[13*20+ 9] =  195.19510; 
1627  daa[13*20+10] =  156.51600; daa[13*20+11] =    0.00000; daa[13*20+12] =   92.18600; daa[14*20+ 0] =  248.59200; 
1628  daa[14*20+ 1] =  102.83130; daa[14*20+ 2] =   41.92440; daa[14*20+ 3] =   13.39400; daa[14*20+ 4] =   18.75500; 
1629  daa[14*20+ 5] =  152.61880; daa[14*20+ 6] =   50.70030; daa[14*20+ 7] =   34.71530; daa[14*20+ 8] =   93.37090; 
1630  daa[14*20+ 9] =   11.91520; daa[14*20+10] =   31.62580; daa[14*20+11] =   33.54190; daa[14*20+12] =   17.02050; 
1631  daa[14*20+13] =   11.05060; daa[15*20+ 0] =  405.18700; daa[15*20+ 1] =  153.15900; daa[15*20+ 2] =  488.58920; 
1632  daa[15*20+ 3] =   95.60970; daa[15*20+ 4] =  159.83560; daa[15*20+ 5] =   56.18280; daa[15*20+ 6] =   79.39990; 
1633  daa[15*20+ 7] =  232.22430; daa[15*20+ 8] =   35.36430; daa[15*20+ 9] =   24.79550; daa[15*20+10] =   17.14320; 
1634  daa[15*20+11] =   95.45570; daa[15*20+12] =   61.99510; daa[15*20+13] =   45.99010; daa[15*20+14] =  242.72020; 
1635  daa[16*20+ 0] =  368.03650; daa[16*20+ 1] =   26.57450; daa[16*20+ 2] =  227.16970; daa[16*20+ 3] =   66.09300; 
1636  daa[16*20+ 4] =   16.23660; daa[16*20+ 5] =   52.56510; daa[16*20+ 6] =   34.01560; daa[16*20+ 7] =   30.66620; 
1637  daa[16*20+ 8] =   22.63330; daa[16*20+ 9] =  190.07390; daa[16*20+10] =   33.10900; daa[16*20+11] =  135.05990; 
1638  daa[16*20+12] =  103.15340; daa[16*20+13] =   13.66550; daa[16*20+14] =   78.28570; daa[16*20+15] =  543.66740; 
1639  daa[17*20+ 0] =    0.00000; daa[17*20+ 1] =  200.13750; daa[17*20+ 2] =   22.49680; daa[17*20+ 3] =    0.00000; 
1640  daa[17*20+ 4] =    0.00000; daa[17*20+ 5] =    0.00000; daa[17*20+ 6] =    0.00000; daa[17*20+ 7] =    0.00000; 
1641  daa[17*20+ 8] =   27.05640; daa[17*20+ 9] =    0.00000; daa[17*20+10] =   46.17760; daa[17*20+11] =    0.00000; 
1642  daa[17*20+12] =    0.00000; daa[17*20+13] =   76.23540; daa[17*20+14] =    0.00000; daa[17*20+15] =   74.08190; 
1643  daa[17*20+16] =    0.00000; daa[18*20+ 0] =   24.41390; daa[18*20+ 1] =    7.80120; daa[18*20+ 2] =   94.69400; 
1644  daa[18*20+ 3] =    0.00000; daa[18*20+ 4] =   95.31640; daa[18*20+ 5] =    0.00000; daa[18*20+ 6] =   21.47170; 
1645  daa[18*20+ 7] =    0.00000; daa[18*20+ 8] =  126.54000; daa[18*20+ 9] =   37.48340; daa[18*20+10] =   28.65720; 
1646  daa[18*20+11] =   13.21420; daa[18*20+12] =    0.00000; daa[18*20+13] =  695.26290; daa[18*20+14] =    0.00000; 
1647  daa[18*20+15] =   33.62890; daa[18*20+16] =   41.78390; daa[18*20+17] =   60.80700; daa[19*20+ 0] =  205.95640; 
1648  daa[19*20+ 1] =   24.03680; daa[19*20+ 2] =   15.80670; daa[19*20+ 3] =   17.83160; daa[19*20+ 4] =   48.46780; 
1649  daa[19*20+ 5] =   34.69830; daa[19*20+ 6] =   36.72500; daa[19*20+ 7] =   53.81650; daa[19*20+ 8] =   43.87150; 
1650  daa[19*20+ 9] =  881.00380; daa[19*20+10] =  174.51560; daa[19*20+11] =   10.38500; daa[19*20+12] =  256.59550; 
1651  daa[19*20+13] =   12.36060; daa[19*20+14] =   48.50260; daa[19*20+15] =   30.38360; daa[19*20+16] =  156.19970; 
1652  daa[19*20+17] =    0.00000; daa[19*20+18] =   27.93790; 
1653 
1654  for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j];
1655
1656  pi[ 0] = 0.087127; pi[ 1] = 0.040904; pi[ 2] = 0.040432; pi[ 3] = 0.046872; 
1657  pi[ 4] = 0.033474; pi[ 5] = 0.038255; pi[ 6] = 0.049530; pi[ 7] = 0.088612;
1658  pi[ 8] = 0.033619; pi[ 9] = 0.036886; pi[10] = 0.085357; pi[11] = 0.080481;
1659  pi[12] = 0.014753; pi[13] = 0.039772; pi[14] = 0.050680; pi[15] = 0.069577;
1660  pi[16] = 0.058542; pi[17] = 0.010494; pi[18] = 0.029916; pi[19] = 0.064718; 
1661
1662  return 1;
1663}
1664
1665//////////////////////////////////////////////////////////////
1666//////////////////////////////////////////////////////////////
1667
1668
1669int Init_Qmat_MtArt(phydbl *daa, phydbl *pi) // Added by Federico Abascal
1670{
1671    /*
1672       Federico Abascal, April 2005 (c).
1673       
1674       This model has been derived from 36 artropoda mitochondrial genomes.
1675
1676       Each gene of the given species was aligned individually. Then, alignments of the whole set
1677       of 13 genes where concatenated and passed through GBlocks (Castresana, 2000, in JME) with
1678       parameters and output:
1679
1680            Minimum Number Of Sequences For A Conserved Position: 20
1681            Minimum Number Of Sequences For A Flanking Position: 32
1682            Maximum Number Of Contiguous Nonconserved Positions: 8
1683            Minimum Length Of A Block: 10
1684            Allowed Gap Positions: With Half
1685            Use Similarity Matrices: Yes
1686
1687            Flank positions of the 40 selected block(s)
1688            Flanks: [6  22]  [26  44]  [61  70]  [77  143]  [145  185]  [208  236]  [309  640] 
1689            [644  802]  [831  941]  [956  966]  [973  1062]  [1085  1339]  [1343  1702] 
1690            [1754  1831]  [1840  1911]  [1916  1987]  [2011  2038]  [2097  2118]  [2125  2143] 
1691            [2179  2215]  [2243  2268]  [2277  2288]  [2333  2347]  [2476  2518]  [2539  2558] 
1692            [2600  2613]  [2637  2672]  [2738  2759]  [2784  2839]  [2882  2924]  [2948  3097] 
1693            [3113  3123]  [3210  3235]  [3239  3322]  [3348  3392]  [3406  3526]  [3588  3617] 
1694            [3660  3692]  [3803  3830]  [3909  3928] 
1695
1696            New number of positions in MtArt-strict.phy.fasta-gb: <b> 2664 </b> (67% of the original 3933 positions)
1697       
1698       The species included in the analysis were:
1699            Harpiosquilla harpax          [NCBI_TaxID 287944]   
1700            Ixodes uriae                  [NCBI_TaxID 59655]     
1701            Heptathela hangzhouensis      [NCBI_TaxID 216259]   
1702            Triops longicaudatus          [NCBI_TaxID 58777]     
1703            Gryllotalpa orientalis        [NCBI_TaxID 213494]   
1704            lepidopsocid RS-2001          [NCBI_TaxID 159971]   
1705            Locusta migratoria            [NCBI_TaxID 7004]     
1706            Drosophila yakuba             [NCBI_TaxID 7245]     
1707            Ostrinia furnacalis           [NCBI_TaxID 93504]     
1708            Megabalanus volcano           [NCBI_TaxID 266495]   
1709            Periplaneta fuliginosa        [NCBI_TaxID 36977]     
1710            Thermobia domestica           [NCBI_TaxID 89055]     
1711            Aleurochiton aceris           [NCBI_TaxID 266942]   
1712            Schizaphis graminum           [NCBI_TaxID 13262]     
1713            Pteronarcys princeps          [NCBI_TaxID 285953]   
1714            Aleurodicus dugesii           [NCBI_TaxID 30099]     
1715            Pollicipes polymerus          [NCBI_TaxID 36137]     
1716            Gomphiocephalus hodgsoni      [NCBI_TaxID 221270]   
1717            Habronattus oregonensis       [NCBI_TaxID 130930]   
1718            Speleonectes tulumensis       [NCBI_TaxID 84346]     
1719            Hutchinsoniella macracantha   [NCBI_TaxID 84335]     
1720            Haemaphysalis flava           [NCBI_TaxID 181088]   
1721            Scutigera coleoptrata         [NCBI_TaxID 29022]     
1722            Vargula hilgendorfii          [NCBI_TaxID 6674]     
1723            Tricholepidion gertschi       [NCBI_TaxID 89825]     
1724            Varroa destructor             [NCBI_TaxID 109461]   
1725            Bombyx mandarina              [NCBI_TaxID 7092]     
1726            Thyropygus sp.                [NCBI_TaxID 174155]   
1727            Tribolium castaneum           [NCBI_TaxID 7070]     
1728            Pagurus longicarpus           [NCBI_TaxID 111067]   
1729            Limulus polyphemus            [NCBI_TaxID 6850]     
1730            Tetrodontophora bielanensis   [NCBI_TaxID 48717]     
1731            Penaeus monodon               [NCBI_TaxID 6687]     
1732            Daphnia pulex                 [NCBI_TaxID 6669]     
1733            Apis mellifera                [NCBI_TaxID 7469]     
1734            Anopheles gambiae             [NCBI_TaxID 7165]     
1735             
1736       The topology used for inferring the model was:
1737            (((Daph_pulex,Trio_longi),((((((Aleu_aceri,Aleu_duges),Schi_grami),lepi_RS_20),
1738            ((((Ostr_furna,Bomb_manda),(Dros_yakub,Anop_gambi)),Apis_melli),Trib_casta)),
1739            ((Gryl_orien,Locu_migra),(Pter_princ,Peri_fulig))),(Tric_gerts,Ther_domes)),
1740            (Scut_coleo,Thyr_sp),Varg_hilge,Hutc_macra,((((Ixod_uriae,Haem_flava),Varr_destr),
1741            (Habr_orego,Hept_hangz)),Limu_polyp),(Poll_polym,Mega_volca),(Gomp_hodgs,Tetr_biela),
1742            ((Pagu_longi,Pena_monod),Harp_harpa),Spel_tulum));
1743           
1744            Note this is not the ML topoLOGy but the consensus one (based on morphoLOGical data,
1745            phyLOGenetic reconstruction using nuclear genes, etc). Where relationships are
1746            not clear, a polytomy was introduced (it contains quite a lot of polytomies!).
1747       
1748       The model was estimated using (the great and helpful) Ziheng Yang's Paml software package.
1749       A four-categorized gamma distribution was used to account for heterogeneity (alpha
1750       was estimated to be 0.47821). Sites with ambiguity data were taken into account.
1751       
1752       If you would like the data related to this matrix, please contact fabascal@uvigo.es.
1753       Federico Abascal (c)2005.
1754    */
1755
1756    int i,j,naa;
1757    naa = 20;
1758    daa[1*20+ 0] = 0.2;     daa[2*20+ 0] = 0.2;     daa[2*20+ 1] = 0.2;     daa[3*20+ 0] = 0.6;
1759    daa[3*20+ 1] = 4.3;     daa[3*20+ 2] = 500.2;   daa[4*20+ 0] = 253.5;   daa[4*20+ 1] = 35.5;
1760    daa[4*20+ 2] = 98.2;    daa[4*20+ 3] = 10.6;    daa[5*20+ 0] = 0.2;     daa[5*20+ 1] = 154.0;
1761    daa[5*20+ 2] = 261.8;   daa[5*20+ 3] = 0.2;     daa[5*20+ 4] = 0.2;     daa[6*20+ 0] = 0.2;
1762    daa[6*20+ 1] = 0.2;     daa[6*20+ 2] = 183.0;   daa[6*20+ 3] = 861.8;   daa[6*20+ 4] = 0.2;
1763    daa[6*20+ 5] = 261.6;   daa[7*20+ 0] = 199.8;   daa[7*20+ 1] = 0.2;     daa[7*20+ 2] = 120.5;
1764    daa[7*20+ 3] = 12.5;    daa[7*20+ 4] = 80.5;    daa[7*20+ 5] = 2.6;     daa[7*20+ 6] = 43.9;
1765    daa[8*20+ 0] = 0.2;     daa[8*20+ 1] = 41.3;    daa[8*20+ 2] = 179.5;   daa[8*20+ 3] = 0.2;
1766    daa[8*20+ 4] = 12.4;    daa[8*20+ 5] = 313.5;   daa[8*20+ 6] = 15.2;    daa[8*20+ 7] = 0.2;
1767    daa[9*20+ 0] = 25.7;    daa[9*20+ 1] = 1.8;     daa[9*20+ 2] = 21.3;    daa[9*20+ 3] = 6.6;
1768    daa[9*20+ 4] = 63.0;    daa[9*20+ 5] = 10.5;    daa[9*20+ 6] = 6.8;     daa[9*20+ 7] = 2.7;
1769    daa[9*20+ 8] = 0.2;     daa[10*20+ 0] = 3.7;    daa[10*20+ 1] = 1.8;    daa[10*20+ 2] = 12.6;
1770    daa[10*20+ 3] = 1.2;    daa[10*20+ 4] = 78.7;   daa[10*20+ 5] = 16.3;   daa[10*20+ 6] = 1.7;
1771    daa[10*20+ 7] = 1.4;    daa[10*20+ 8] = 5.5;    daa[10*20+ 9] = 514.5;  daa[11*20+ 0] = 0.2;
1772    daa[11*20+ 1] = 208.6;  daa[11*20+ 2] = 467.3;  daa[11*20+ 3] = 1.7;    daa[11*20+ 4] = 0.2;
1773    daa[11*20+ 5] = 349.3;  daa[11*20+ 6] = 106.3;  daa[11*20+ 7] = 0.2;    daa[11*20+ 8] = 0.2;
1774    daa[11*20+ 9] = 3.5;    daa[11*20+ 10] = 3.8;   daa[12*20+ 0] = 120.6;  daa[12*20+ 1] = 5.2;
1775    daa[12*20+ 2] = 78.8;   daa[12*20+ 3] = 0.2;    daa[12*20+ 4] = 312.3;  daa[12*20+ 5] = 67.3;
1776    daa[12*20+ 6] = 0.2;    daa[12*20+ 7] = 55.7;   daa[12*20+ 8] = 0.2;    daa[12*20+ 9] = 514.8;
1777    daa[12*20+ 10] = 885.5; daa[12*20+ 11] = 105.6; daa[13*20+ 0] = 13.1;   daa[13*20+ 1] = 4.7;
1778    daa[13*20+ 2] = 19.7;   daa[13*20+ 3] = 0.2;    daa[13*20+ 4] = 184.1;  daa[13*20+ 5] = 0.2;
1779    daa[13*20+ 6] = 0.2;    daa[13*20+ 7] = 0.8;    daa[13*20+ 8] = 13.8;   daa[13*20+ 9] = 117.9;
1780    daa[13*20+ 10] = 262.6; daa[13*20+ 11] = 10.7;  daa[13*20+ 12] = 321.6; daa[14*20+ 0] = 49.3;
1781    daa[14*20+ 1] = 0.2;    daa[14*20+ 2] = 16.5;   daa[14*20+ 3] = 0.2;    daa[14*20+ 4] = 0.2;
1782    daa[14*20+ 5] = 39.3;   daa[14*20+ 6] = 7.9;    daa[14*20+ 7] = 0.2;    daa[14*20+ 8] = 0.8;
1783    daa[14*20+ 9] = 0.2;    daa[14*20+ 10] = 12.2;  daa[14*20+ 11] = 16.8;  daa[14*20+ 12] = 5.3;
1784    daa[14*20+ 13] = 14.6;  daa[15*20+ 0] = 673.0;  daa[15*20+ 1] = 2.7;    daa[15*20+ 2] = 398.4;
1785    daa[15*20+ 3] = 44.4;   daa[15*20+ 4] = 664.2;  daa[15*20+ 5] = 52.4;   daa[15*20+ 6] = 31.5;
1786    daa[15*20+ 7] = 226.0;  daa[15*20+ 8] = 10.6;   daa[15*20+ 9] = 7.2;    daa[15*20+ 10] = 8.2;
1787    daa[15*20+ 11] = 144.2; daa[15*20+ 12] = 111.7; daa[15*20+ 13] = 36.1;  daa[15*20+ 14] = 86.5;
1788    daa[16*20+ 0] = 243.9;  daa[16*20+ 1] = 0.2;    daa[16*20+ 2] = 165.9;  daa[16*20+ 3] = 0.2;
1789    daa[16*20+ 4] = 182.8;  daa[16*20+ 5] = 43.7;   daa[16*20+ 6] = 43.4;   daa[16*20+ 7] = 0.2;
1790    daa[16*20+ 8] = 18.6;   daa[16*20+ 9] = 203.7;  daa[16*20+ 10] = 47.8;  daa[16*20+ 11] = 69.5;
1791    daa[16*20+ 12] = 288.6; daa[16*20+ 13] = 13.5;  daa[16*20+ 14] = 46.8;  daa[16*20+ 15] = 660.4;
1792    daa[17*20+ 0] = 0.2;    daa[17*20+ 1] = 0.2;    daa[17*20+ 2] = 7.7;    daa[17*20+ 3] = 0.2;
1793    daa[17*20+ 4] = 21.6;   daa[17*20+ 5] = 6.7;    daa[17*20+ 6] = 11.0;   daa[17*20+ 7] = 1.9;
1794    daa[17*20+ 8] = 0.2;    daa[17*20+ 9] = 0.2;    daa[17*20+ 10] = 21.1;  daa[17*20+ 11] = 16.0;
1795    daa[17*20+ 12] = 70.7;  daa[17*20+ 13] = 53.7;  daa[17*20+ 14] = 0.2;   daa[17*20+ 15] = 2.4;
1796    daa[17*20+ 16] = 0.2;   daa[18*20+ 0] = 1.2;    daa[18*20+ 1] = 3.9;    daa[18*20+ 2] = 251.2;
1797    daa[18*20+ 3] = 0.2;    daa[18*20+ 4] = 72.0;   daa[18*20+ 5] = 86.7;   daa[18*20+ 6] = 7.7;
1798    daa[18*20+ 7] = 8.6;    daa[18*20+ 8] = 191.4;  daa[18*20+ 9] = 12.3;   daa[18*20+ 10] = 19.8;
1799    daa[18*20+ 11] = 117.1; daa[18*20+ 12] = 70.9;  daa[18*20+ 13] = 791.6; daa[18*20+ 14] = 18.4;
1800    daa[18*20+ 15] = 30.5;  daa[18*20+ 16] = 46.0;  daa[18*20+ 17] = 37.7;  daa[19*20+ 0] = 339.9;
1801    daa[19*20+ 1] = 0.2;    daa[19*20+ 2] = 22.6;   daa[19*20+ 3] = 0.2;    daa[19*20+ 4] = 350.4;
1802    daa[19*20+ 5] = 0.2;    daa[19*20+ 6] = 13.6;   daa[19*20+ 7] = 2.6;    daa[19*20+ 8] = 0.2;
1803    daa[19*20+ 9] = 1854.5; daa[19*20+ 10] = 84.7;  daa[19*20+ 11] = 26.1;  daa[19*20+ 12] = 281.3;
1804    daa[19*20+ 13] = 51.9;  daa[19*20+ 14] = 31.7;  daa[19*20+ 15] = 60.6;  daa[19*20+ 16] = 544.1;
1805    daa[19*20+ 17] = 0.2;   daa[19*20+ 18] = 1.6;
1806       
1807/*  MtArt.old: esta es la MtArt que hice con 26 secuencias (2 outgroups) con una topoLOGia incorrecta
1808    daa[1*20+ 0] = 0.2;     daa[2*20+ 0] = 0.2;     daa[2*20+ 1] = 8.0;     daa[3*20+ 0] = 0.2;
1809    daa[3*20+ 1] = 0.2;     daa[3*20+ 2] = 441.7;   daa[4*20+ 0] = 287.9;   daa[4*20+ 1] = 48.4;
1810    daa[4*20+ 2] = 82.4;    daa[4*20+ 3] = 0.2;     daa[5*20+ 0] = 0.2;     daa[5*20+ 1] = 149.9;
1811    daa[5*20+ 2] = 278.6;   daa[5*20+ 3] = 0.2;     daa[5*20+ 4] = 21.7;    daa[6*20+ 0] = 6.6;
1812    daa[6*20+ 1] = 0.2;     daa[6*20+ 2] = 213.9;   daa[6*20+ 3] = 760.8;   daa[6*20+ 4] = 0.2;
1813    daa[6*20+ 5] = 292.9;   daa[7*20+ 0] = 228.2;   daa[7*20+ 1] = 0.2;     daa[7*20+ 2] = 97.1;
1814    daa[7*20+ 3] = 10.4;    daa[7*20+ 4] = 98.4;    daa[7*20+ 5] = 4.0;     daa[7*20+ 6] = 48.7;
1815    daa[8*20+ 0] = 0.2;     daa[8*20+ 1] = 56.7;    daa[8*20+ 2] = 156.4;   daa[8*20+ 3] = 24.5;
1816    daa[8*20+ 4] = 15.5;    daa[8*20+ 5] = 328.6;   daa[8*20+ 6] = 7.0;     daa[8*20+ 7] = 8.4;
1817    daa[9*20+ 0] = 26.4;    daa[9*20+ 1] = 1.6;     daa[9*20+ 2] = 40.1;    daa[9*20+ 3] = 0.2;
1818    daa[9*20+ 4] = 22.1;    daa[9*20+ 5] = 13.8;    daa[9*20+ 6] = 0.2;     daa[9*20+ 7] = 3.6;
1819    daa[9*20+ 8] = 0.2;     daa[10*20+ 0] = 3.4;    daa[10*20+ 1] = 0.6;    daa[10*20+ 2] = 13.8;
1820    daa[10*20+ 3] = 0.7;    daa[10*20+ 4] = 76.9;   daa[10*20+ 5] = 12.1;   daa[10*20+ 6] = 5.4;
1821    daa[10*20+ 7] = 2.5;    daa[10*20+ 8] = 2.9;    daa[10*20+ 9] = 542.6;  daa[11*20+ 0] = 0.2;
1822    daa[11*20+ 1] = 240.2;  daa[11*20+ 2] = 602.8;  daa[11*20+ 3] = 35.5;   daa[11*20+ 4] = 0.2;
1823    daa[11*20+ 5] = 357.6;  daa[11*20+ 6] = 62.6;   daa[11*20+ 7] = 0.2;    daa[11*20+ 8] = 3.3;
1824    daa[11*20+ 9] = 0.2;    daa[11*20+ 10] = 17.5;  daa[12*20+ 0] = 119.0;  daa[12*20+ 1] = 0.2;
1825    daa[12*20+ 2] = 91.4;   daa[12*20+ 3] = 6.4;    daa[12*20+ 4] = 332.3;  daa[12*20+ 5] = 65.4;
1826    daa[12*20+ 6] = 0.2;    daa[12*20+ 7] = 60.4;   daa[12*20+ 8] = 2.4;    daa[12*20+ 9] = 492.5;
1827    daa[12*20+ 10] = 815.8; daa[12*20+ 11] = 67.3;  daa[13*20+ 0] = 8.2;    daa[13*20+ 1] = 6.4;
1828    daa[13*20+ 2] = 31.5;   daa[13*20+ 3] = 3.4;    daa[13*20+ 4] = 174.4;  daa[13*20+ 5] = 5.7;
1829    daa[13*20+ 6] = 5.7;    daa[13*20+ 7] = 2.1;    daa[13*20+ 8] = 11.0;   daa[13*20+ 9] = 94.4;
1830    daa[13*20+ 10] = 243.3; daa[13*20+ 11] = 12.3;  daa[13*20+ 12] = 357.8; daa[14*20+ 0] = 62.5;
1831    daa[14*20+ 1] = 0.4;    daa[14*20+ 2] = 17.5;   daa[14*20+ 3] = 0.2;    daa[14*20+ 4] = 0.2;
1832    daa[14*20+ 5] = 48.6;   daa[14*20+ 6] = 17.7;   daa[14*20+ 7] = 2.7;    daa[14*20+ 8] = 0.2;
1833    daa[14*20+ 9] = 0.2;    daa[14*20+ 10] = 11.2;  daa[14*20+ 11] = 21.7;  daa[14*20+ 12] = 5.2;
1834    daa[14*20+ 13] = 12.6;  daa[15*20+ 0] = 659.0;  daa[15*20+ 1] = 5.2;    daa[15*20+ 2] = 469.8;
1835    daa[15*20+ 3] = 52.3;   daa[15*20+ 4] = 570.7;  daa[15*20+ 5] = 47.8;   daa[15*20+ 6] = 37.3;
1836    daa[15*20+ 7] = 227.8;  daa[15*20+ 8] = 12.7;   daa[15*20+ 9] = 12.3;   daa[15*20+ 10] = 7.4;
1837    daa[15*20+ 11] = 189.0; daa[15*20+ 12] = 155.3; daa[15*20+ 13] = 43.8;  daa[15*20+ 14] = 103.4;
1838    daa[16*20+ 0] = 276.4;  daa[16*20+ 1] = 1.6;    daa[16*20+ 2] = 175.6;  daa[16*20+ 3] = 0.2;
1839    daa[16*20+ 4] = 96.2;   daa[16*20+ 5] = 71.4;   daa[16*20+ 6] = 37.4;   daa[16*20+ 7] = 0.2;
1840    daa[16*20+ 8] = 14.2;   daa[16*20+ 9] = 212.5;  daa[16*20+ 10] = 38.5;  daa[16*20+ 11] = 97.4;
1841    daa[16*20+ 12] = 254.7; daa[16*20+ 13] = 2.1;   daa[16*20+ 14] = 41.6;  daa[16*20+ 15] = 670.6;
1842    daa[17*20+ 0] = 6.2;    daa[17*20+ 1] = 0.2;    daa[17*20+ 2] = 0.2;    daa[17*20+ 3] = 5.6;
1843    daa[17*20+ 4] = 0.2;    daa[17*20+ 5] = 0.2;    daa[17*20+ 6] = 3.1;    daa[17*20+ 7] = 0.4;
1844    daa[17*20+ 8] = 0.2;    daa[17*20+ 9] = 15.2;   daa[17*20+ 10] = 11.5;  daa[17*20+ 11] = 32.6;
1845    daa[17*20+ 12] = 82.4;  daa[17*20+ 13] = 81.9;  daa[17*20+ 14] = 0.2;   daa[17*20+ 15] = 9.7;
1846    daa[17*20+ 16] = 0.2;   daa[18*20+ 0] = 1.6;    daa[18*20+ 1] = 7.7;    daa[18*20+ 2] = 242.5;
1847    daa[18*20+ 3] = 0.2;    daa[18*20+ 4] = 88.0;   daa[18*20+ 5] = 93.1;   daa[18*20+ 6] = 0.2;
1848    daa[18*20+ 7] = 6.0;    daa[18*20+ 8] = 113.7;  daa[18*20+ 9] = 22.1;   daa[18*20+ 10] = 17.2;
1849    daa[18*20+ 11] = 138.5; daa[18*20+ 12] = 37.6;  daa[18*20+ 13] = 770.2; daa[18*20+ 14] = 5.3;
1850    daa[18*20+ 15] = 25.0;  daa[18*20+ 16] = 55.5;  daa[18*20+ 17] = 69.3;  daa[19*20+ 0] = 307.8;
1851    daa[19*20+ 1] = 2.2;    daa[19*20+ 2] = 6.9;    daa[19*20+ 3] = 0.2;    daa[19*20+ 4] = 405.7;
1852    daa[19*20+ 5] = 0.8;    daa[19*20+ 6] = 20.2;   daa[19*20+ 7] = 5.7;    daa[19*20+ 8] = 0.2;
1853    daa[19*20+ 9] = 1687.9; daa[19*20+ 10] = 49.4;  daa[19*20+ 11] = 23.4;  daa[19*20+ 12] = 329.9;
1854    daa[19*20+ 13] = 86.3;  daa[19*20+ 14] = 27.3;  daa[19*20+ 15] = 95.0;  daa[19*20+ 16] = 443.0;
1855    daa[19*20+ 17] = 2.4;   daa[19*20+ 18] = 0.2;
18563*/   
1857    for (i=0; i<naa; i++) for (j=0; j<i; j++) daa[j*naa+i] = daa[i*naa+j];
1858
1859    pi[0] = 0.054116;       pi[1] = 0.018227;       pi[2] = 0.039903;       pi[3] = 0.020160;       pi[4] = 0.009709;
1860    pi[5] = 0.018781;       pi[6] = 0.024289;       pi[7] = 0.068183;       pi[8] = 0.024518;       pi[9] = 0.092639;
1861    pi[10] = 0.148658;      pi[11] = 0.021718;      pi[12] = 0.061453;      pi[13] = 0.088668;      pi[14] = 0.041826;
1862    pi[15] = 0.091030;      pi[16] = 0.049194;      pi[17] = 0.029786;      pi[18] = 0.039443;      pi[19] = 0.057701;
1863/*  pi[0] = 0.055505;       pi[1] = 0.018320;       pi[2] = 0.036614;       pi[3] = 0.019517;    pi[4] = 0.009878;
1864    pi[5] = 0.018732;       pi[6] = 0.021591;       pi[7] = 0.068414;       pi[8] = 0.023638;    pi[9] = 0.092558;
1865    pi[10] = 0.154697;      pi[11] = 0.020168;      pi[12] = 0.062033;      pi[13] = 0.089433;   pi[14] = 0.040642;
1866    pi[15] = 0.090298;      pi[16] = 0.050134;      pi[17] = 0.027693;      pi[18] = 0.038448;   pi[19] = 0.061687;
1867*/   
1868    return 1;
1869}
1870
1871//////////////////////////////////////////////////////////////
1872//////////////////////////////////////////////////////////////
1873
1874
1875int Init_Qmat_HIVb(phydbl *daa, phydbl *pi) //added by FEDE
1876{
1877    /*
1878      Nickle DC, Heath L, Jensen MA, Gilbert PB, Mullins JI, Kosakovsky Pond SL.
1879      HIV-Specific Probabilistic Models of Protein Evolution.
1880      PLoS ONE. 2007 Jun 6;2:e503.
1881
1882          [thanks to Sergei L. Kosakovsky]
1883         
1884          Translated from HYPHY to Phyml format by Federico Abascal.
1885    */
1886
1887    int i,j,naa;
1888    naa = 20;
1889
1890    daa[1*20+0]= 0.307507;        daa[2*20+0]= 0.005;           daa[2*20+1]= 0.295543;        daa[3*20+0]= 1.45504;         
1891    daa[3*20+1]= 0.005;           daa[3*20+2]= 17.6612;         daa[4*20+0]= 0.123758;        daa[4*20+1]= 0.351721;       
1892    daa[4*20+2]= 0.0860642;       daa[4*20+3]= 0.005;           daa[5*20+0]= 0.0551128;       daa[5*20+1]= 3.4215;         
1893    daa[5*20+2]= 0.672052;        daa[5*20+3]= 0.005;           daa[5*20+4]= 0.005;           daa[6*20+0]= 1.48135;         
1894    daa[6*20+1]= 0.0749218;       daa[6*20+2]= 0.0792633;       daa[6*20+3]= 10.5872;         daa[6*20+4]= 0.005;           
1895    daa[6*20+5]= 2.5602;          daa[7*20+0]= 2.13536;         daa[7*20+1]= 3.65345;         daa[7*20+2]= 0.323401;       
1896    daa[7*20+3]= 2.83806;         daa[7*20+4]= 0.897871;        daa[7*20+5]= 0.0619137;       daa[7*20+6]= 3.92775;         
1897    daa[8*20+0]= 0.0847613;       daa[8*20+1]= 9.04044;         daa[8*20+2]= 7.64585;         daa[8*20+3]= 1.9169;         
1898    daa[8*20+4]= 0.240073;        daa[8*20+5]= 7.05545;         daa[8*20+6]= 0.11974;         daa[8*20+7]= 0.005;           
1899    daa[9*20+0]= 0.005;           daa[9*20+1]= 0.677289;        daa[9*20+2]= 0.680565;        daa[9*20+3]= 0.0176792;       
1900    daa[9*20+4]= 0.005;           daa[9*20+5]= 0.005;           daa[9*20+6]= 0.00609079;      daa[9*20+7]= 0.005;           
1901    daa[9*20+8]= 0.103111;        daa[10*20+0]= 0.215256;       daa[10*20+1]= 0.701427;       daa[10*20+2]= 0.005;         
1902    daa[10*20+3]= 0.00876048;     daa[10*20+4]= 0.129777;       daa[10*20+5]= 1.49456;        daa[10*20+6]= 0.005;         
1903    daa[10*20+7]= 0.005;          daa[10*20+8]= 1.74171;        daa[10*20+9]= 5.95879;        daa[11*20+0]= 0.005;         
1904    daa[11*20+1]= 20.45;          daa[11*20+2]= 7.90443;        daa[11*20+3]= 0.005;          daa[11*20+4]= 0.005;         
1905    daa[11*20+5]= 6.54737;        daa[11*20+6]= 4.61482;        daa[11*20+7]= 0.521705;       daa[11*20+8]= 0.005;         
1906    daa[11*20+9]= 0.322319;       daa[11*20+10]= 0.0814995;     daa[12*20+0]= 0.0186643;      daa[12*20+1]= 2.51394;       
1907    daa[12*20+2]= 0.005;          daa[12*20+3]= 0.005;          daa[12*20+4]= 0.005;          daa[12*20+5]= 0.303676;       
1908    daa[12*20+6]= 0.175789;       daa[12*20+7]= 0.005;          daa[12*20+8]= 0.005;          daa[12*20+9]= 11.2065;       
1909    daa[12*20+10]= 5.31961;       daa[12*20+11]= 1.28246;       daa[13*20+0]= 0.0141269;      daa[13*20+1]= 0.005;         
1910    daa[13*20+2]= 0.005;          daa[13*20+3]= 0.005;          daa[13*20+4]= 9.29815;        daa[13*20+5]= 0.005;         
1911    daa[13*20+6]= 0.005;          daa[13*20+7]= 0.291561;       daa[13*20+8]= 0.145558;       daa[13*20+9]= 3.39836;       
1912    daa[13*20+10]= 8.52484;       daa[13*20+11]= 0.0342658;     daa[13*20+12]= 0.188025;      daa[14*20+0]= 2.12217;       
1913    daa[14*20+1]= 1.28355;        daa[14*20+2]= 0.00739578;     daa[14*20+3]= 0.0342658;      daa[14*20+4]= 0.005;         
1914    daa[14*20+5]= 4.47211;        daa[14*20+6]= 0.0120226;      daa[14*20+7]= 0.005;          daa[14*20+8]= 2.45318;       
1915    daa[14*20+9]= 0.0410593;      daa[14*20+10]= 2.07757;       daa[14*20+11]= 0.0313862;     daa[14*20+12]= 0.005;         
1916    daa[14*20+13]= 0.005;         daa[15*20+0]= 2.46633;        daa[15*20+1]= 3.4791;         daa[15*20+2]= 13.1447;       
1917    daa[15*20+3]= 0.52823;        daa[15*20+4]= 4.69314;        daa[15*20+5]= 0.116311;       daa[15*20+6]= 0.005;         
1918    daa[15*20+7]= 4.38041;        daa[15*20+8]= 0.382747;       daa[15*20+9]= 1.21803;        daa[15*20+10]= 0.927656;     
1919    daa[15*20+11]= 0.504111;      daa[15*20+12]= 0.005;         daa[15*20+13]= 0.956472;      daa[15*20+14]= 5.37762;       
1920    daa[16*20+0]= 15.9183;        daa[16*20+1]= 2.86868;        daa[16*20+2]= 6.88667;        daa[16*20+3]= 0.274724;       
1921    daa[16*20+4]= 0.739969;       daa[16*20+5]= 0.243589;       daa[16*20+6]= 0.289774;       daa[16*20+7]= 0.369615;       
1922    daa[16*20+8]= 0.711594;       daa[16*20+9]= 8.61217;        daa[16*20+10]= 0.0437673;     daa[16*20+11]= 4.67142;       
1923    daa[16*20+12]= 4.94026;       daa[16*20+13]= 0.0141269;     daa[16*20+14]= 2.01417;       daa[16*20+15]= 8.93107;       
1924    daa[17*20+0]= 0.005;          daa[17*20+1]= 0.991338;       daa[17*20+2]= 0.005;          daa[17*20+3]= 0.005;         
1925    daa[17*20+4]= 2.63277;        daa[17*20+5]= 0.026656;       daa[17*20+6]= 0.005;          daa[17*20+7]= 1.21674;       
1926    daa[17*20+8]= 0.0695179;      daa[17*20+9]= 0.005;          daa[17*20+10]= 0.748843;      daa[17*20+11]= 0.005;         
1927    daa[17*20+12]= 0.089078;      daa[17*20+13]= 0.829343;      daa[17*20+14]= 0.0444506;     daa[17*20+15]= 0.0248728;     
1928    daa[17*20+16]= 0.005;         daa[18*20+0]= 0.005;          daa[18*20+1]= 0.00991826;     daa[18*20+2]= 1.76417;       
1929    daa[18*20+3]= 0.674653;       daa[18*20+4]= 7.57932;        daa[18*20+5]= 0.113033;       daa[18*20+6]= 0.0792633;     
1930    daa[18*20+7]= 0.005;          daa[18*20+8]= 18.6943;        daa[18*20+9]= 0.148168;       daa[18*20+10]= 0.111986;     
1931    daa[18*20+11]= 0.005;         daa[18*20+12]= 0.005;         daa[18*20+13]= 15.34;         daa[18*20+14]= 0.0304381;     
1932    daa[18*20+15]= 0.648024;      daa[18*20+16]= 0.105652;      daa[18*20+17]= 1.28022;       daa[19*20+0]= 7.61428;       
1933    daa[19*20+1]= 0.0812454;      daa[19*20+2]= 0.026656;       daa[19*20+3]= 1.04793;        daa[19*20+4]= 0.420027;       
1934    daa[19*20+5]= 0.0209153;      daa[19*20+6]= 1.02847;        daa[19*20+7]= 0.953155;       daa[19*20+8]= 0.005;         
1935    daa[19*20+9]= 17.7389;        daa[19*20+10]= 1.41036;       daa[19*20+11]= 0.265829;      daa[19*20+12]= 6.8532;       
1936    daa[19*20+13]= 0.723274;      daa[19*20+14]= 0.005;         daa[19*20+15]= 0.0749218;     daa[19*20+16]= 0.709226;     
1937    daa[19*20+17]= 0.005;         daa[19*20+18]= 0.0410593;     
1938    for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j];
1939
1940    pi[0]= 0.060490222;           pi[1]= 0.066039665;           pi[2]= 0.044127815;           pi[3]= 0.042109048;           
1941    pi[4]= 0.020075899;           pi[5]= 0.053606488;           pi[6]= 0.071567447;           pi[7]= 0.072308239;           
1942    pi[8]= 0.022293943;           pi[9]= 0.069730629;           pi[10]= 0.098851122;          pi[11]= 0.056968211;         
1943    pi[12]= 0.019768318;          pi[13]= 0.028809447;          pi[14]= 0.046025282;          pi[15]= 0.05060433;           
1944    pi[16]= 0.053636813;          pi[17]= 0.033011601;          pi[18]= 0.028350243;          pi[19]= 0.061625237;         
1945    return 1;
1946}
1947
1948//////////////////////////////////////////////////////////////
1949//////////////////////////////////////////////////////////////
1950
1951int Init_Qmat_HIVw(phydbl *daa, phydbl *pi)
1952{
1953    /*
1954          Nickle DC, Heath L, Jensen MA, Gilbert PB, Mullins JI, Kosakovsky Pond SL.
1955      HIV-Specific Probabilistic Models of Protein Evolution.
1956      PLoS ONE. 2007 Jun 6;2:e503.
1957
1958          [thanks to Sergei L. Kosakovsky]
1959         
1960          Translated from HYPHY to Phyml format by Federico Abascal.
1961    */
1962
1963    int i,j,naa;
1964    naa = 20;
1965
1966    daa[1*20+0]= 0.0744808;       daa[2*20+0]= 0.617509;        daa[2*20+1]= 0.16024;         daa[3*20+0]= 4.43521;         
1967    daa[3*20+1]= 0.0674539;       daa[3*20+2]= 29.4087;         daa[4*20+0]= 0.167653;        daa[4*20+1]= 2.86364;         
1968    daa[4*20+2]= 0.0604932;       daa[4*20+3]= 0.005;           daa[5*20+0]= 0.005;           daa[5*20+1]= 10.6746;         
1969    daa[5*20+2]= 0.342068;        daa[5*20+3]= 0.005;           daa[5*20+4]= 0.005;           daa[6*20+0]= 5.56325;         
1970    daa[6*20+1]= 0.0251632;       daa[6*20+2]= 0.201526;        daa[6*20+3]= 12.1233;         daa[6*20+4]= 0.005;           
1971    daa[6*20+5]= 3.20656;         daa[7*20+0]= 1.8685;          daa[7*20+1]= 13.4379;         daa[7*20+2]= 0.0604932;       
1972    daa[7*20+3]= 10.3969;         daa[7*20+4]= 0.0489798;       daa[7*20+5]= 0.0604932;       daa[7*20+6]= 14.7801;         
1973    daa[8*20+0]= 0.005;           daa[8*20+1]= 6.84405;         daa[8*20+2]= 8.59876;         daa[8*20+3]= 2.31779;         
1974    daa[8*20+4]= 0.005;           daa[8*20+5]= 18.5465;         daa[8*20+6]= 0.005;           daa[8*20+7]= 0.005;           
1975    daa[9*20+0]= 0.005;           daa[9*20+1]= 1.34069;         daa[9*20+2]= 0.987028;        daa[9*20+3]= 0.145124;       
1976    daa[9*20+4]= 0.005;           daa[9*20+5]= 0.0342252;       daa[9*20+6]= 0.0390512;       daa[9*20+7]= 0.005;           
1977    daa[9*20+8]= 0.005;           daa[10*20+0]= 0.16024;        daa[10*20+1]= 0.586757;       daa[10*20+2]= 0.005;         
1978    daa[10*20+3]= 0.005;          daa[10*20+4]= 0.005;          daa[10*20+5]= 2.89048;        daa[10*20+6]= 0.129839;       
1979    daa[10*20+7]= 0.0489798;      daa[10*20+8]= 1.76382;        daa[10*20+9]= 9.10246;        daa[11*20+0]= 0.592784;       
1980    daa[11*20+1]= 39.8897;        daa[11*20+2]= 10.6655;        daa[11*20+3]= 0.894313;       daa[11*20+4]= 0.005;         
1981    daa[11*20+5]= 13.0705;        daa[11*20+6]= 23.9626;        daa[11*20+7]= 0.279425;       daa[11*20+8]= 0.22406;       
1982    daa[11*20+9]= 0.817481;       daa[11*20+10]= 0.005;         daa[12*20+0]= 0.005;          daa[12*20+1]= 3.28652;       
1983    daa[12*20+2]= 0.201526;       daa[12*20+3]= 0.005;          daa[12*20+4]= 0.005;          daa[12*20+5]= 0.005;         
1984    daa[12*20+6]= 0.005;          daa[12*20+7]= 0.0489798;      daa[12*20+8]= 0.005;          daa[12*20+9]= 17.3064;       
1985    daa[12*20+10]= 11.3839;       daa[12*20+11]= 4.09564;       daa[13*20+0]= 0.597923;       daa[13*20+1]= 0.005;         
1986    daa[13*20+2]= 0.005;          daa[13*20+3]= 0.005;          daa[13*20+4]= 0.362959;       daa[13*20+5]= 0.005;         
1987    daa[13*20+6]= 0.005;          daa[13*20+7]= 0.005;          daa[13*20+8]= 0.005;          daa[13*20+9]= 1.48288;       
1988    daa[13*20+10]= 7.48781;       daa[13*20+11]= 0.005;         daa[13*20+12]= 0.005;         daa[14*20+0]= 1.00981;       
1989    daa[14*20+1]= 0.404723;       daa[14*20+2]= 0.344848;       daa[14*20+3]= 0.005;          daa[14*20+4]= 0.005;         
1990    daa[14*20+5]= 3.04502;        daa[14*20+6]= 0.005;          daa[14*20+7]= 0.005;          daa[14*20+8]= 13.9444;       
1991    daa[14*20+9]= 0.005;          daa[14*20+10]= 9.83095;       daa[14*20+11]= 0.111928;      daa[14*20+12]= 0.005;         
1992    daa[14*20+13]= 0.0342252;     daa[15*20+0]= 8.5942;         daa[15*20+1]= 8.35024;        daa[15*20+2]= 14.5699;       
1993    daa[15*20+3]= 0.427881;       daa[15*20+4]= 1.12195;        daa[15*20+5]= 0.16024;        daa[15*20+6]= 0.005;         
1994    daa[15*20+7]= 6.27966;        daa[15*20+8]= 0.725157;       daa[15*20+9]= 0.740091;       daa[15*20+10]= 6.14396;       
1995    daa[15*20+11]= 0.005;         daa[15*20+12]= 0.392575;      daa[15*20+13]= 4.27939;       daa[15*20+14]= 14.249;       
1996    daa[16*20+0]= 24.1422;        daa[16*20+1]= 0.928203;       daa[16*20+2]= 4.54206;        daa[16*20+3]= 0.630395;       
1997    daa[16*20+4]= 0.005;          daa[16*20+5]= 0.203091;       daa[16*20+6]= 0.458743;       daa[16*20+7]= 0.0489798;     
1998    daa[16*20+8]= 0.95956;        daa[16*20+9]= 9.36345;        daa[16*20+10]= 0.005;         daa[16*20+11]= 4.04802;       
1999    daa[16*20+12]= 7.41313;       daa[16*20+13]= 0.114512;      daa[16*20+14]= 4.33701;       daa[16*20+15]= 6.34079;       
2000    daa[17*20+0]= 0.005;          daa[17*20+1]= 5.96564;        daa[17*20+2]= 0.005;          daa[17*20+3]= 0.005;         
2001    daa[17*20+4]= 5.49894;        daa[17*20+5]= 0.0443298;      daa[17*20+6]= 0.005;          daa[17*20+7]= 2.8258;         
2002    daa[17*20+8]= 0.005;          daa[17*20+9]= 0.005;          daa[17*20+10]= 1.37031;       daa[17*20+11]= 0.005;         
2003    daa[17*20+12]= 0.005;         daa[17*20+13]= 0.005;         daa[17*20+14]= 0.005;         daa[17*20+15]= 1.10156;       
2004    daa[17*20+16]= 0.005;         daa[18*20+0]= 0.005;          daa[18*20+1]= 0.005;          daa[18*20+2]= 5.06475;       
2005    daa[18*20+3]= 2.28154;        daa[18*20+4]= 8.34835;        daa[18*20+5]= 0.005;          daa[18*20+6]= 0.005;         
2006    daa[18*20+7]= 0.005;          daa[18*20+8]= 47.4889;        daa[18*20+9]= 0.114512;       daa[18*20+10]= 0.005;         
2007    daa[18*20+11]= 0.005;         daa[18*20+12]= 0.579198;      daa[18*20+13]= 4.12728;       daa[18*20+14]= 0.005;         
2008    daa[18*20+15]= 0.933142;      daa[18*20+16]= 0.490608;      daa[18*20+17]= 0.005;         daa[19*20+0]= 24.8094;       
2009    daa[19*20+1]= 0.279425;       daa[19*20+2]= 0.0744808;      daa[19*20+3]= 2.91786;        daa[19*20+4]= 0.005;         
2010    daa[19*20+5]= 0.005;          daa[19*20+6]= 2.19952;        daa[19*20+7]= 2.79622;        daa[19*20+8]= 0.827479;       
2011    daa[19*20+9]= 24.8231;        daa[19*20+10]= 2.95344;       daa[19*20+11]= 0.128065;      daa[19*20+12]= 14.7683;       
2012    daa[19*20+13]= 2.28;          daa[19*20+14]= 0.005;         daa[19*20+15]= 0.862637;      daa[19*20+16]= 0.005;         
2013    daa[19*20+17]= 0.005;         daa[19*20+18]= 1.35482;       
2014    for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j];
2015
2016    pi[0]= 0.0377494;             pi[1]= 0.057321;              pi[2]= 0.0891129;             pi[3]= 0.0342034;             
2017    pi[4]= 0.0240105;             pi[5]= 0.0437824;             pi[6]= 0.0618606;             pi[7]= 0.0838496;             
2018    pi[8]= 0.0156076;             pi[9]= 0.0983641;             pi[10]= 0.0577867;            pi[11]= 0.0641682;           
2019    pi[12]= 0.0158419;            pi[13]= 0.0422741;            pi[14]= 0.0458601;            pi[15]= 0.0550846;           
2020    pi[16]= 0.0813774;            pi[17]= 0.019597;             pi[18]= 0.0205847;            pi[19]= 0.0515639;           
2021    return 1;
2022}
2023
2024//////////////////////////////////////////////////////////////
2025//////////////////////////////////////////////////////////////
2026
2027int Init_Qmat_JTT(phydbl *daa, phydbl *pi)
2028{
2029  int i,j,naa;
2030
2031  /* JTT's model data
2032   * D.T.Jones, W.R.Taylor and J.M.Thornton
2033   * "The rapid generation of mutation data matrices from protein sequences"
2034   * CABIOS  vol.8 no.3 1992 pp275-282
2035   */
2036
2037  naa = 20;
2038
2039/*   PhyML_Printf("\n\n. REMINDER : THIS IS NOT JTT !!!\n\n"); */
2040
2041
2042/*   daa[1*20 + 0] = 0.592439084;  */
2043/*   daa[2*20 + 0] = 0.427686326;  */
2044/*   daa[2*20 + 1] = 0.795121783;  */
2045/*   daa[3*20 + 0] = 0.649119246;  */
2046/*   daa[3*20 + 1] = 0.116285502;  */
2047/*   daa[3*20 + 2] = 6.004659556;  */
2048/*   daa[4*20 + 0] = 2.384113318;  */
2049/*   daa[4*20 + 1] = 0.925080814;  */
2050/*   daa[4*20 + 2] = 0.648459384;  */
2051/*   daa[4*20 + 3] = 0.088632079;  */
2052/*   daa[5*20 + 0] = 1.188499461;  */
2053/*   daa[5*20 + 1] = 3.230029905;  */
2054/*   daa[5*20 + 2] = 1.910660418;  */
2055/*   daa[5*20 + 3] = 0.599639374;  */
2056/*   daa[5*20 + 4] = 0.123291416;  */
2057/*   daa[6*20 + 0] = 1.534092487;  */
2058/*   daa[6*20 + 1] = 0.390654981;  */
2059/*   daa[6*20 + 2] = 0.566472743;  */
2060/*   daa[6*20 + 3] = 5.916482256;  */
2061/*   daa[6*20 + 4] = 0.001457992;  */
2062/*   daa[6*20 + 5] = 4.564894849;  */
2063/*   daa[7*20 + 0] = 2.282012114;  */
2064/*   daa[7*20 + 1] = 0.502549742;  */
2065/*   daa[7*20 + 2] = 1.97013036;  */
2066/*   daa[7*20 + 3] = 1.198016556;  */
2067/*   daa[7*20 + 4] = 0.565245125;  */
2068/*   daa[7*20 + 5] = 0.382529852;  */
2069/*   daa[7*20 + 6] = 0.477515495;  */
2070/*   daa[8*20 + 0] = 0.447587408;  */
2071/*   daa[8*20 + 1] = 2.948365572;  */
2072/*   daa[8*20 + 2] = 5.900694893;  */
2073/*   daa[8*20 + 3] = 0.96157692;  */
2074/*   daa[8*20 + 4] = 0.757153271;  */
2075/*   daa[8*20 + 5] = 5.048587596;  */
2076/*   daa[8*20 + 6] = 0.410583482;  */
2077/*   daa[8*20 + 7] = 0.406115605;  */
2078/*   daa[9*20 + 0] = 0.105337934;  */
2079/*   daa[9*20 + 1] = 0.115421183;  */
2080/*   daa[9*20 + 2] = 0.178948502;  */
2081/*   daa[9*20 + 3] = 0.01319564;  */
2082/*   daa[9*20 + 4] = 0.269637325;  */
2083/*   daa[9*20 + 5] = 0.087618072;  */
2084/*   daa[9*20 + 6] = 0.045232757;  */
2085/*   daa[9*20 + 7] = 0.00958642;  */
2086/*   daa[9*20 + 8] = 0.116145048;  */
2087/*   daa[10*20 + 0] = 0.269460438;  */
2088/*   daa[10*20 + 1] = 0.416549984;  */
2089/*   daa[10*20 + 2] = 0.100490196;  */
2090/*   daa[10*20 + 3] = 0.021637571;  */
2091/*   daa[10*20 + 4] = 0.610442865;  */
2092/*   daa[10*20 + 5] = 0.782347175;  */
2093/*   daa[10*20 + 6] = 0.077565847;  */
2094/*   daa[10*20 + 7] = 0.034916545;  */
2095/*   daa[10*20 + 8] = 0.456334781;  */
2096/*   daa[10*20 + 9] = 4.057302736;  */
2097/*   daa[11*20 + 0] = 0.71706921;  */
2098/*   daa[11*20 + 1] = 6.469532749;  */
2099/*   daa[11*20 + 2] = 2.606469511;  */
2100/*   daa[11*20 + 3] = 0.239488567;  */
2101/*   daa[11*20 + 4] = 0.020594132;  */
2102/*   daa[11*20 + 5] = 3.978479772;  */
2103/*   daa[11*20 + 6] = 1.669092348;  */
2104/*   daa[11*20 + 7] = 0.368262509;  */
2105/*   daa[11*20 + 8] = 0.744835535;  */
2106/*   daa[11*20 + 9] = 0.153626461;  */
2107/*   daa[11*20 + 10] = 0.186106976;  */
2108/*   daa[12*20 + 0] = 0.766271778;  */
2109/*   daa[12*20 + 1] = 0.465902609;  */
2110/*   daa[12*20 + 2] = 0.352429433;  */
2111/*   daa[12*20 + 3] = 0.04670109;  */
2112/*   daa[12*20 + 4] = 0.930016402;  */
2113/*   daa[12*20 + 5] = 1.529072358;  */
2114/*   daa[12*20 + 6] = 0.149497081;  */
2115/*   daa[12*20 + 7] = 0.092795029;  */
2116/*   daa[12*20 + 8] = 0.367452788;  */
2117/*   daa[12*20 + 9] = 4.443975824;  */
2118/*   daa[12*20 + 10] = 7.305963709;  */
2119/*   daa[12*20 + 11] = 0.691067831;  */
2120/*   daa[13*20 + 0] = 0.164314785;  */
2121/*   daa[13*20 + 1] = 0.064046051;  */
2122/*   daa[13*20 + 2] = 0.108166343;  */
2123/*   daa[13*20 + 3] = 0.021096631;  */
2124/*   daa[13*20 + 4] = 0.9869692;  */
2125/*   daa[13*20 + 5] = 0.054959009;  */
2126/*   daa[13*20 + 6] = 0.025919221;  */
2127/*   daa[13*20 + 7] = 0.053577857;  */
2128/*   daa[13*20 + 8] = 0.933637734;  */
2129/*   daa[13*20 + 9] = 0.928915383;  */
2130/*   daa[13*20 + 10] = 2.549658015;  */
2131/*   daa[13*20 + 11] = 0.032602723;  */
2132/*   daa[13*20 + 12] = 1.837397986;  */
2133/*   daa[14*20 + 0] = 1.57937466;  */
2134/*   daa[14*20 + 1] = 0.360749538;  */
2135/*   daa[14*20 + 2] = 0.189183366;  */
2136/*   daa[14*20 + 3] = 0.475387263;  */
2137/*   daa[14*20 + 4] = 0.108255878;  */
2138/*   daa[14*20 + 5] = 0.645294956;  */
2139/*   daa[14*20 + 6] = 0.537278362;  */
2140/*   daa[14*20 + 7] = 0.264516985;  */
2141/*   daa[14*20 + 8] = 0.510062949;  */
2142/*   daa[14*20 + 9] = 0.063380048;  */
2143/*   daa[14*20 + 10] = 0.25549658;  */
2144/*   daa[14*20 + 11] = 0.42259906;  */
2145/*   daa[14*20 + 12] = 0.105330979;  */
2146/*   daa[14*20 + 13] = 0.092650543;  */
2147/*   daa[15*20 + 0] = 5.358849732;  */
2148/*   daa[15*20 + 1] = 0.980362668;  */
2149/*   daa[15*20 + 2] = 5.114953674;  */
2150/*   daa[15*20 + 3] = 1.242898123;  */
2151/*   daa[15*20 + 4] = 3.417258976;  */
2152/*   daa[15*20 + 5] = 1.199117804;  */
2153/*   daa[15*20 + 6] = 0.61024061;  */
2154/*   daa[15*20 + 7] = 2.200746452;  */
2155/*   daa[15*20 + 8] = 0.928717085;  */
2156/*   daa[15*20 + 9] = 0.07737779;  */
2157/*   daa[15*20 + 10] = 0.17657022;  */
2158/*   daa[15*20 + 11] = 0.768621405;  */
2159/*   daa[15*20 + 12] = 0.330215244;  */
2160/*   daa[15*20 + 13] = 0.359985962;  */
2161/*   daa[15*20 + 14] = 1.559718375;  */
2162/*   daa[16*20 + 0] = 2.21115657;  */
2163/*   daa[16*20 + 1] = 0.711040813;  */
2164/*   daa[16*20 + 2] = 2.273163111;  */
2165/*   daa[16*20 + 3] = 0.408851733;  */
2166/*   daa[16*20 + 4] = 1.57572444;  */
2167/*   daa[16*20 + 5] = 1.074674746;  */
2168/*   daa[16*20 + 6] = 0.707096016;  */
2169/*   daa[16*20 + 7] = 0.157418108;  */
2170/*   daa[16*20 + 8] = 0.543088926;  */
2171/*   daa[16*20 + 9] = 1.056024176;  */
2172/*   daa[16*20 + 10] = 0.294118282;  */
2173/*   daa[16*20 + 11] = 1.325676171;  */
2174/*   daa[16*20 + 12] = 1.834248079;  */
2175/*   daa[16*20 + 13] = 0.136987103;  */
2176/*   daa[16*20 + 14] = 0.63778855;  */
2177/*   daa[16*20 + 15] = 7.294506488;  */
2178/*   daa[17*20 + 0] = 0.142106331;  */
2179/*   daa[17*20 + 1] = 0.600995696;  */
2180/*   daa[17*20 + 2] = 0.079481793;  */
2181/*   daa[17*20 + 3] = 0.069125482;  */
2182/*   daa[17*20 + 4] = 0.701931839;  */
2183/*   daa[17*20 + 5] = 0.23881661;  */
2184/*   daa[17*20 + 6] = 0.082469159;  */
2185/*   daa[17*20 + 7] = 0.185722676;  */
2186/*   daa[17*20 + 8] = 0.911694299;  */
2187/*   daa[17*20 + 9] = 0.11679729;  */
2188/*   daa[17*20 + 10] = 0.688555894;  */
2189/*   daa[17*20 + 11] = 0.063328042;  */
2190/*   daa[17*20 + 12] = 0.780079225;  */
2191/*   daa[17*20 + 13] = 3.432193724;  */
2192/*   daa[17*20 + 14] = 0.088827331;  */
2193/*   daa[17*20 + 15] = 0.293602516;  */
2194/*   daa[17*20 + 16] = 0.158083832;  */
2195/*   daa[18*20 + 0] = 0.178618354;  */
2196/*   daa[18*20 + 1] = 0.243634289;  */
2197/*   daa[18*20 + 2] = 0.567657832;  */
2198/*   daa[18*20 + 3] = 0.134890583;  */
2199/*   daa[18*20 + 4] = 1.28093676;  */
2200/*   daa[18*20 + 5] = 0.180560506;  */
2201/*   daa[18*20 + 6] = 0.081426157;  */
2202/*   daa[18*20 + 7] = 0.044067219;  */
2203/*   daa[18*20 + 8] = 6.464181222;  */
2204/*   daa[18*20 + 9] = 0.206728215;  */
2205/*   daa[18*20 + 10] = 0.334379787;  */
2206/*   daa[18*20 + 11] = 0.102051407;  */
2207/*   daa[18*20 + 12] = 0.431775879;  */
2208/*   daa[18*20 + 13] = 11.87781125;  */
2209/*   daa[18*20 + 14] = 0.078070175;  */
2210/*   daa[18*20 + 15] = 0.359969114;  */
2211/*   daa[18*20 + 16] = 0.199064371;  */
2212/*   daa[18*20 + 17] = 4.185797039;  */
2213/*   daa[19*20 + 0] = 1.964727121;  */
2214/*   daa[19*20 + 1] = 0.18579405;  */
2215/*   daa[19*20 + 2] = 0.113068304;  */
2216/*   daa[19*20 + 3] = 0.046258503;  */
2217/*   daa[19*20 + 4] = 2.069983598;  */
2218/*   daa[19*20 + 5] = 0.247259847;  */
2219/*   daa[19*20 + 6] = 0.272991021;  */
2220/*   daa[19*20 + 7] = 0.073432734;  */
2221/*   daa[19*20 + 8] = 0.158879333;  */
2222/*   daa[19*20 + 9] = 10.52590329;  */
2223/*   daa[19*20 + 10] = 1.74812825;  */
2224/*   daa[19*20 + 11] = 0.188152512;  */
2225/*   daa[19*20 + 12] = 1.800124087;  */
2226/*   daa[19*20 + 13] = 0.478606732;  */
2227/*   daa[19*20 + 14] = 0.277216066;  */
2228/*   daa[19*20 + 15] = 0.104728903;  */
2229/*   daa[19*20 + 16] = 1.963379491;  */
2230/*   daa[19*20 + 17] = 0.16552242;  */
2231/*   daa[19*20 + 18] = 0.222040517; */
2232   
2233/*   for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j]; */
2234
2235/*   pi[0] = 0.093862;  */
2236/*   pi[1] = 0.05757;  */
2237/*   pi[2] = 0.037017;  */
2238/*   pi[3] = 0.060952;  */
2239/*   pi[4] = 0.011004;  */
2240/*   pi[5] = 0.044631;  */
2241/*   pi[6] = 0.083648;  */
2242/*   pi[7] = 0.053004;  */
2243/*   pi[8] = 0.022528;  */
2244/*   pi[9] = 0.060152;  */
2245/*   pi[10] = 0.092314;  */
2246/*   pi[11] = 0.065886;  */
2247/*   pi[12] = 0.021075;  */
2248/*   pi[13] = 0.033928;  */
2249/*   pi[14] = 0.043333;  */
2250/*   pi[15] = 0.052957;  */
2251/*   pi[16] = 0.053465;  */
2252/*   pi[17] = 0.00928;  */
2253/*   pi[18] = 0.030025;  */
2254/*   pi[19] = 0.073369; */
2255 
2256  daa[ 1*20+ 0] =   58.00; daa[ 2*20+ 0] =   54.00; daa[ 2*20+ 1] =   45.00; daa[ 3*20+ 0] =   81.00;
2257  daa[ 3*20+ 1] =   16.00; daa[ 3*20+ 2] =  528.00; daa[ 4*20+ 0] =   56.00; daa[ 4*20+ 1] =  113.00;
2258  daa[ 4*20+ 2] =   34.00; daa[ 4*20+ 3] =   10.00; daa[ 5*20+ 0] =   57.00; daa[ 5*20+ 1] =  310.00;
2259  daa[ 5*20+ 2] =   86.00; daa[ 5*20+ 3] =   49.00; daa[ 5*20+ 4] =    9.00; daa[ 6*20+ 0] =  105.00;
2260  daa[ 6*20+ 1] =   29.00; daa[ 6*20+ 2] =   58.00; daa[ 6*20+ 3] =  767.00; daa[ 6*20+ 4] =    5.00;
2261  daa[ 6*20+ 5] =  323.00; daa[ 7*20+ 0] =  179.00; daa[ 7*20+ 1] =  137.00; daa[ 7*20+ 2] =   81.00;
2262  daa[ 7*20+ 3] =  130.00; daa[ 7*20+ 4] =   59.00; daa[ 7*20+ 5] =   26.00; daa[ 7*20+ 6] =  119.00;
2263  daa[ 8*20+ 0] =   27.00; daa[ 8*20+ 1] =  328.00; daa[ 8*20+ 2] =  391.00; daa[ 8*20+ 3] =  112.00;
2264  daa[ 8*20+ 4] =   69.00; daa[ 8*20+ 5] =  597.00; daa[ 8*20+ 6] =   26.00; daa[ 8*20+ 7] =   23.00;
2265  daa[ 9*20+ 0] =   36.00; daa[ 9*20+ 1] =   22.00; daa[ 9*20+ 2] =   47.00; daa[ 9*20+ 3] =   11.00;
2266  daa[ 9*20+ 4] =   17.00; daa[ 9*20+ 5] =    9.00; daa[ 9*20+ 6] =   12.00; daa[ 9*20+ 7] =    6.00;
2267  daa[ 9*20+ 8] =   16.00; daa[10*20+ 0] =   30.00; daa[10*20+ 1] =   38.00; daa[10*20+ 2] =   12.00;
2268  daa[10*20+ 3] =    7.00; daa[10*20+ 4] =   23.00; daa[10*20+ 5] =   72.00; daa[10*20+ 6] =    9.00;
2269  daa[10*20+ 7] =    6.00; daa[10*20+ 8] =   56.00; daa[10*20+ 9] =  229.00; daa[11*20+ 0] =   35.00;
2270  daa[11*20+ 1] =  646.00; daa[11*20+ 2] =  263.00; daa[11*20+ 3] =   26.00; daa[11*20+ 4] =    7.00;
2271  daa[11*20+ 5] =  292.00; daa[11*20+ 6] =  181.00; daa[11*20+ 7] =   27.00; daa[11*20+ 8] =   45.00;
2272  daa[11*20+ 9] =   21.00; daa[11*20+10] =   14.00; daa[12*20+ 0] =   54.00; daa[12*20+ 1] =   44.00;
2273  daa[12*20+ 2] =   30.00; daa[12*20+ 3] =   15.00; daa[12*20+ 4] =   31.00; daa[12*20+ 5] =   43.00;
2274  daa[12*20+ 6] =   18.00; daa[12*20+ 7] =   14.00; daa[12*20+ 8] =   33.00; daa[12*20+ 9] =  479.00;
2275  daa[12*20+10] =  388.00; daa[12*20+11] =   65.00; daa[13*20+ 0] =   15.00; daa[13*20+ 1] =    5.00;
2276  daa[13*20+ 2] =   10.00; daa[13*20+ 3] =    4.00; daa[13*20+ 4] =   78.00; daa[13*20+ 5] =    4.00;
2277  daa[13*20+ 6] =    5.00; daa[13*20+ 7] =    5.00; daa[13*20+ 8] =   40.00; daa[13*20+ 9] =   89.00;
2278  daa[13*20+10] =  248.00; daa[13*20+11] =    4.00; daa[13*20+12] =   43.00; daa[14*20+ 0] =  194.00;
2279  daa[14*20+ 1] =   74.00; daa[14*20+ 2] =   15.00; daa[14*20+ 3] =   15.00; daa[14*20+ 4] =   14.00;
2280  daa[14*20+ 5] =  164.00; daa[14*20+ 6] =   18.00; daa[14*20+ 7] =   24.00; daa[14*20+ 8] =  115.00;
2281  daa[14*20+ 9] =   10.00; daa[14*20+10] =  102.00; daa[14*20+11] =   21.00; daa[14*20+12] =   16.00;
2282  daa[14*20+13] =   17.00; daa[15*20+ 0] =  378.00; daa[15*20+ 1] =  101.00; daa[15*20+ 2] =  503.00;
2283  daa[15*20+ 3] =   59.00; daa[15*20+ 4] =  223.00; daa[15*20+ 5] =   53.00; daa[15*20+ 6] =   30.00;
2284  daa[15*20+ 7] =  201.00; daa[15*20+ 8] =   73.00; daa[15*20+ 9] =   40.00; daa[15*20+10] =   59.00;
2285  daa[15*20+11] =   47.00; daa[15*20+12] =   29.00; daa[15*20+13] =   92.00; daa[15*20+14] =  285.00;
2286  daa[16*20+ 0] =  475.00; daa[16*20+ 1] =   64.00; daa[16*20+ 2] =  232.00; daa[16*20+ 3] =   38.00;
2287  daa[16*20+ 4] =   42.00; daa[16*20+ 5] =   51.00; daa[16*20+ 6] =   32.00; daa[16*20+ 7] =   33.00;
2288  daa[16*20+ 8] =   46.00; daa[16*20+ 9] =  245.00; daa[16*20+10] =   25.00; daa[16*20+11] =  103.00;
2289  daa[16*20+12] =  226.00; daa[16*20+13] =   12.00; daa[16*20+14] =  118.00; daa[16*20+15] =  477.00;
2290  daa[17*20+ 0] =    9.00; daa[17*20+ 1] =  126.00; daa[17*20+ 2] =    8.00; daa[17*20+ 3] =    4.00;
2291  daa[17*20+ 4] =  115.00; daa[17*20+ 5] =   18.00; daa[17*20+ 6] =   10.00; daa[17*20+ 7] =   55.00;
2292  daa[17*20+ 8] =    8.00; daa[17*20+ 9] =    9.00; daa[17*20+10] =   52.00; daa[17*20+11] =   10.00;
2293  daa[17*20+12] =   24.00; daa[17*20+13] =   53.00; daa[17*20+14] =    6.00; daa[17*20+15] =   35.00;
2294  daa[17*20+16] =   12.00; daa[18*20+ 0] =   11.00; daa[18*20+ 1] =   20.00; daa[18*20+ 2] =   70.00;
2295  daa[18*20+ 3] =   46.00; daa[18*20+ 4] =  209.00; daa[18*20+ 5] =   24.00; daa[18*20+ 6] =    7.00;
2296  daa[18*20+ 7] =    8.00; daa[18*20+ 8] =  573.00; daa[18*20+ 9] =   32.00; daa[18*20+10] =   24.00;
2297  daa[18*20+11] =    8.00; daa[18*20+12] =   18.00; daa[18*20+13] =  536.00; daa[18*20+14] =   10.00;
2298  daa[18*20+15] =   63.00; daa[18*20+16] =   21.00; daa[18*20+17] =   71.00; daa[19*20+ 0] =  298.00;
2299  daa[19*20+ 1] =   17.00; daa[19*20+ 2] =   16.00; daa[19*20+ 3] =   31.00; daa[19*20+ 4] =   62.00;
2300  daa[19*20+ 5] =   20.00; daa[19*20+ 6] =   45.00; daa[19*20+ 7] =   47.00; daa[19*20+ 8] =   11.00;
2301  daa[19*20+ 9] =  961.00; daa[19*20+10] =  180.00; daa[19*20+11] =   14.00; daa[19*20+12] =  323.00;
2302  daa[19*20+13] =   62.00; daa[19*20+14] =   23.00; daa[19*20+15] =   38.00; daa[19*20+16] =  112.00;
2303  daa[19*20+17] =   25.00; daa[19*20+18] =   16.00;
2304 
2305  for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j];
2306
2307  pi[ 0] = 0.076748; pi[ 1] = 0.051691; pi[ 2] = 0.042645; pi[ 3] = 0.051544;
2308  pi[ 4] = 0.019803; pi[ 5] = 0.040752; pi[ 6] = 0.061830; pi[ 7] = 0.073152;
2309  pi[ 8] = 0.022944; pi[ 9] = 0.053761; pi[10] = 0.091904; pi[11] = 0.058676;
2310  pi[12] = 0.023826; pi[13] = 0.040126; pi[14] = 0.050901; pi[15] = 0.068765;
2311  pi[16] = 0.058565; pi[17] = 0.014261; pi[18] = 0.032102; pi[19] = 0.066005;
2312 
2313 return 1;
2314}
2315
2316//////////////////////////////////////////////////////////////
2317//////////////////////////////////////////////////////////////
2318
2319
2320int Init_Qmat_MtREV(phydbl *daa, phydbl *pi)
2321{
2322  /* J. Adachi and M. Hasegawa, ``Model of amino acid substitution in proteins
2323     encoded by mitochondrial DNA'' J. Mol. Evol. 42, 459 (1996) */
2324
2325  int i,j,naa;
2326
2327  naa = 20;
2328
2329
2330  daa[ 1*20+ 0] =   23.18; daa[ 2*20+ 0] =   26.95; daa[ 2*20+ 1] =   13.24; daa[ 3*20+ 0] =   17.67;
2331  daa[ 3*20+ 1] =    1.90; daa[ 3*20+ 2] =  794.38; daa[ 4*20+ 0] =   59.93; daa[ 4*20+ 1] =  103.33;
2332  daa[ 4*20+ 2] =   58.94; daa[ 4*20+ 3] =    1.90; daa[ 5*20+ 0] =    1.90; daa[ 5*20+ 1] =  220.99;
2333  daa[ 5*20+ 2] =  173.56; daa[ 5*20+ 3] =   55.28; daa[ 5*20+ 4] =   75.24; daa[ 6*20+ 0] =    9.77;
2334  daa[ 6*20+ 1] =    1.90; daa[ 6*20+ 2] =   63.05; daa[ 6*20+ 3] =  583.55; daa[ 6*20+ 4] =    1.90;
2335  daa[ 6*20+ 5] =  313.56; daa[ 7*20+ 0] =  120.71; daa[ 7*20+ 1] =   23.03; daa[ 7*20+ 2] =   53.30;
2336  daa[ 7*20+ 3] =   56.77; daa[ 7*20+ 4] =   30.71; daa[ 7*20+ 5] =    6.75; daa[ 7*20+ 6] =   28.28;
2337  daa[ 8*20+ 0] =   13.90; daa[ 8*20+ 1] =  165.23; daa[ 8*20+ 2] =  496.13; daa[ 8*20+ 3] =  113.99;
2338  daa[ 8*20+ 4] =  141.49; daa[ 8*20+ 5] =  582.40; daa[ 8*20+ 6] =   49.12; daa[ 8*20+ 7] =    1.90;
2339  daa[ 9*20+ 0] =   96.49; daa[ 9*20+ 1] =    1.90; daa[ 9*20+ 2] =   27.10; daa[ 9*20+ 3] =    4.34;
2340  daa[ 9*20+ 4] =   62.73; daa[ 9*20+ 5] =    8.34; daa[ 9*20+ 6] =    3.31; daa[ 9*20+ 7] =    5.98;
2341  daa[ 9*20+ 8] =   12.26; daa[10*20+ 0] =   25.46; daa[10*20+ 1] =   15.58; daa[10*20+ 2] =   15.16;
2342  daa[10*20+ 3] =    1.90; daa[10*20+ 4] =   25.65; daa[10*20+ 5] =   39.70; daa[10*20+ 6] =    1.90;
2343  daa[10*20+ 7] =    2.41; daa[10*20+ 8] =   11.49; daa[10*20+ 9] =  329.09; daa[11*20+ 0] =    8.36;
2344  daa[11*20+ 1] =  141.40; daa[11*20+ 2] =  608.70; daa[11*20+ 3] =    2.31; daa[11*20+ 4] =    1.90;
2345  daa[11*20+ 5] =  465.58; daa[11*20+ 6] =  313.86; daa[11*20+ 7] =   22.73; daa[11*20+ 8] =  127.67;
2346  daa[11*20+ 9] =   19.57; daa[11*20+10] =   14.88; daa[12*20+ 0] =  141.88; daa[12*20+ 1] =    1.90;
2347  daa[12*20+ 2] =   65.41; daa[12*20+ 3] =    1.90; daa[12*20+ 4] =    6.18; daa[12*20+ 5] =   47.37;
2348  daa[12*20+ 6] =    1.90; daa[12*20+ 7] =    1.90; daa[12*20+ 8] =   11.97; daa[12*20+ 9] =  517.98;
2349  daa[12*20+10] =  537.53; daa[12*20+11] =   91.37; daa[13*20+ 0] =    6.37; daa[13*20+ 1] =    4.69;
2350  daa[13*20+ 2] =   15.20; daa[13*20+ 3] =    4.98; daa[13*20+ 4] =   70.80; daa[13*20+ 5] =   19.11;
2351  daa[13*20+ 6] =    2.67; daa[13*20+ 7] =    1.90; daa[13*20+ 8] =   48.16; daa[13*20+ 9] =   84.67;
2352  daa[13*20+10] =  216.06; daa[13*20+11] =    6.44; daa[13*20+12] =   90.82; daa[14*20+ 0] =   54.31;
2353  daa[14*20+ 1] =   23.64; daa[14*20+ 2] =   73.31; daa[14*20+ 3] =   13.43; daa[14*20+ 4] =   31.26;
2354  daa[14*20+ 5] =  137.29; daa[14*20+ 6] =   12.83; daa[14*20+ 7] =    1.90; daa[14*20+ 8] =   60.97;
2355  daa[14*20+ 9] =   20.63; daa[14*20+10] =   40.10; daa[14*20+11] =   50.10; daa[14*20+12] =   18.84;
2356  daa[14*20+13] =   17.31; daa[15*20+ 0] =  387.86; daa[15*20+ 1] =    6.04; daa[15*20+ 2] =  494.39;
2357  daa[15*20+ 3] =   69.02; daa[15*20+ 4] =  277.05; daa[15*20+ 5] =   54.11; daa[15*20+ 6] =   54.71;
2358  daa[15*20+ 7] =  125.93; daa[15*20+ 8] =   77.46; daa[15*20+ 9] =   47.70; daa[15*20+10] =   73.61;
2359  daa[15*20+11] =  105.79; daa[15*20+12] =  111.16; daa[15*20+13] =   64.29; daa[15*20+14] =  169.90;
2360  daa[16*20+ 0] =  480.72; daa[16*20+ 1] =    2.08; daa[16*20+ 2] =  238.46; daa[16*20+ 3] =   28.01;
2361  daa[16*20+ 4] =  179.97; daa[16*20+ 5] =   94.93; daa[16*20+ 6] =   14.82; daa[16*20+ 7] =   11.17;
2362  daa[16*20+ 8] =   44.78; daa[16*20+ 9] =  368.43; daa[16*20+10] =  126.40; daa[16*20+11] =  136.33;
2363  daa[16*20+12] =  528.17; daa[16*20+13] =   33.85; daa[16*20+14] =  128.22; daa[16*20+15] =  597.21;
2364  daa[17*20+ 0] =    1.90; daa[17*20+ 1] =   21.95; daa[17*20+ 2] =   10.68; daa[17*20+ 3] =   19.86;
2365  daa[17*20+ 4] =   33.60; daa[17*20+ 5] =    1.90; daa[17*20+ 6] =    1.90; daa[17*20+ 7] =   10.92;
2366  daa[17*20+ 8] =    7.08; daa[17*20+ 9] =    1.90; daa[17*20+10] =   32.44; daa[17*20+11] =   24.00;
2367  daa[17*20+12] =   21.71; daa[17*20+13] =    7.84; daa[17*20+14] =    4.21; daa[17*20+15] =   38.58;
2368  daa[17*20+16] =    9.99; daa[18*20+ 0] =    6.48; daa[18*20+ 1] =    1.90; daa[18*20+ 2] =  191.36;
2369  daa[18*20+ 3] =   21.21; daa[18*20+ 4] =  254.77; daa[18*20+ 5] =   38.82; daa[18*20+ 6] =   13.12;
2370  daa[18*20+ 7] =    3.21; daa[18*20+ 8] =  670.14; daa[18*20+ 9] =   25.01; daa[18*20+10] =   44.15;
2371  daa[18*20+11] =   51.17; daa[18*20+12] =   39.96; daa[18*20+13] =  465.58; daa[18*20+14] =   16.21;
2372  daa[18*20+15] =   64.92; daa[18*20+16] =   38.73; daa[18*20+17] =   26.25; daa[19*20+ 0] =  195.06;
2373  daa[19*20+ 1] =    7.64; daa[19*20+ 2] =    1.90; daa[19*20+ 3] =    1.90; daa[19*20+ 4] =    1.90;
2374  daa[19*20+ 5] =   19.00; daa[19*20+ 6] =   21.14; daa[19*20+ 7] =    2.53; daa[19*20+ 8] =    1.90;
2375  daa[19*20+ 9] = 1222.94; daa[19*20+10] =   91.67; daa[19*20+11] =    1.90; daa[19*20+12] =  387.54;
2376  daa[19*20+13] =    6.35; daa[19*20+14] =    8.23; daa[19*20+15] =    1.90; daa[19*20+16] =  204.54;
2377  daa[19*20+17] =    5.37; daa[19*20+18] =    1.90;
2378 
2379  for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j];
2380
2381  pi[ 0] = 0.072000; pi[ 1] = 0.019000; pi[ 2] = 0.039000; pi[ 3] = 0.019000;
2382  pi[ 4] = 0.006000; pi[ 5] = 0.025000; pi[ 6] = 0.024000; pi[ 7] = 0.056000;
2383  pi[ 8] = 0.028000; pi[ 9] = 0.088000; pi[10] = 0.169000; pi[11] = 0.023000;
2384  pi[12] = 0.054000; pi[13] = 0.061000; pi[14] = 0.054000; pi[15] = 0.072000;
2385  pi[16] = 0.086000; pi[17] = 0.029000; pi[18] = 0.033000; pi[19] = 0.043000;
2386
2387  return 1;
2388}
2389
2390
2391//////////////////////////////////////////////////////////////
2392//////////////////////////////////////////////////////////////
2393
2394
2395int Init_Qmat_LG(phydbl *daa, phydbl *pi)
2396{
2397  int i,j,naa;
2398
2399  /* LG model
2400   * Si Quang Le & Olivier Gascuel
2401   * 'An improved general amino-acid replacement matrix'
2402   */
2403
2404  naa = 20;
2405
2406  daa[ 1*20+ 0] =   0.425093;  daa[ 2*20+ 0] =   0.276818;  daa[ 2*20+ 1] =   0.751878;  daa[ 3*20+ 0] =   0.395144; 
2407  daa[ 3*20+ 1] =   0.123954;  daa[ 3*20+ 2] =   5.076149;  daa[ 4*20+ 0] =   2.489084;  daa[ 4*20+ 1] =   0.534551; 
2408  daa[ 4*20+ 2] =   0.528768;  daa[ 4*20+ 3] =   0.062556;  daa[ 5*20+ 0] =   0.969894;  daa[ 5*20+ 1] =   2.807908; 
2409  daa[ 5*20+ 2] =   1.695752;  daa[ 5*20+ 3] =   0.523386;  daa[ 5*20+ 4] =   0.084808;  daa[ 6*20+ 0] =   1.038545; 
2410  daa[ 6*20+ 1] =   0.363970;  daa[ 6*20+ 2] =   0.541712;  daa[ 6*20+ 3] =   5.243870;  daa[ 6*20+ 4] =   0.003499; 
2411  daa[ 6*20+ 5] =   4.128591;  daa[ 7*20+ 0] =   2.066040;  daa[ 7*20+ 1] =   0.390192;  daa[ 7*20+ 2] =   1.437645; 
2412  daa[ 7*20+ 3] =   0.844926;  daa[ 7*20+ 4] =   0.569265;  daa[ 7*20+ 5] =   0.267959;  daa[ 7*20+ 6] =   0.348847; 
2413  daa[ 8*20+ 0] =   0.358858;  daa[ 8*20+ 1] =   2.426601;  daa[ 8*20+ 2] =   4.509238;  daa[ 8*20+ 3] =   0.927114; 
2414  daa[ 8*20+ 4] =   0.640543;  daa[ 8*20+ 5] =   4.813505;  daa[ 8*20+ 6] =   0.423881;  daa[ 8*20+ 7] =   0.311484; 
2415  daa[ 9*20+ 0] =   0.149830;  daa[ 9*20+ 1] =   0.126991;  daa[ 9*20+ 2] =   0.191503;  daa[ 9*20+ 3] =   0.010690; 
2416  daa[ 9*20+ 4] =   0.320627;  daa[ 9*20+ 5] =   0.072854;  daa[ 9*20+ 6] =   0.044265;  daa[ 9*20+ 7] =   0.008705; 
2417  daa[ 9*20+ 8] =   0.108882;  daa[10*20+ 0] =   0.395337;  daa[10*20+ 1] =   0.301848;  daa[10*20+ 2] =   0.068427; 
2418  daa[10*20+ 3] =   0.015076;  daa[10*20+ 4] =   0.594007;  daa[10*20+ 5] =   0.582457;  daa[10*20+ 6] =   0.069673; 
2419  daa[10*20+ 7] =   0.044261;  daa[10*20+ 8] =   0.366317;  daa[10*20+ 9] =   4.145067;  daa[11*20+ 0] =   0.536518; 
2420  daa[11*20+ 1] =   6.326067;  daa[11*20+ 2] =   2.145078;  daa[11*20+ 3] =   0.282959;  daa[11*20+ 4] =   0.013266; 
2421  daa[11*20+ 5] =   3.234294;  daa[11*20+ 6] =   1.807177;  daa[11*20+ 7] =   0.296636;  daa[11*20+ 8] =   0.697264; 
2422  daa[11*20+ 9] =   0.159069;  daa[11*20+10] =   0.137500;  daa[12*20+ 0] =   1.124035;  daa[12*20+ 1] =   0.484133; 
2423  daa[12*20+ 2] =   0.371004;  daa[12*20+ 3] =   0.025548;  daa[12*20+ 4] =   0.893680;  daa[12*20+ 5] =   1.672569; 
2424  daa[12*20+ 6] =   0.173735;  daa[12*20+ 7] =   0.139538;  daa[12*20+ 8] =   0.442472;  daa[12*20+ 9] =   4.273607; 
2425  daa[12*20+10] =   6.312358;  daa[12*20+11] =   0.656604;  daa[13*20+ 0] =   0.253701;  daa[13*20+ 1] =   0.052722; 
2426  daa[13*20+ 2] =   0.089525;  daa[13*20+ 3] =   0.017416;  daa[13*20+ 4] =   1.105251;  daa[13*20+ 5] =   0.035855; 
2427  daa[13*20+ 6] =   0.018811;  daa[13*20+ 7] =   0.089586;  daa[13*20+ 8] =   0.682139;  daa[13*20+ 9] =   1.112727; 
2428  daa[13*20+10] =   2.592692;  daa[13*20+11] =   0.023918;  daa[13*20+12] =   1.798853;  daa[14*20+ 0] =   1.177651; 
2429  daa[14*20+ 1] =   0.332533;  daa[14*20+ 2] =   0.161787;  daa[14*20+ 3] =   0.394456;  daa[14*20+ 4] =   0.075382; 
2430  daa[14*20+ 5] =   0.624294;  daa[14*20+ 6] =   0.419409;  daa[14*20+ 7] =   0.196961;  daa[14*20+ 8] =   0.508851; 
2431  daa[14*20+ 9] =   0.078281;  daa[14*20+10] =   0.249060;  daa[14*20+11] =   0.390322;  daa[14*20+12] =   0.099849; 
2432  daa[14*20+13] =   0.094464;  daa[15*20+ 0] =   4.727182;  daa[15*20+ 1] =   0.858151;  daa[15*20+ 2] =   4.008358; 
2433  daa[15*20+ 3] =   1.240275;  daa[15*20+ 4] =   2.784478;  daa[15*20+ 5] =   1.223828;  daa[15*20+ 6] =   0.611973; 
2434  daa[15*20+ 7] =   1.739990;  daa[15*20+ 8] =   0.990012;  daa[15*20+ 9] =   0.064105;  daa[15*20+10] =   0.182287; 
2435  daa[15*20+11] =   0.748683;  daa[15*20+12] =   0.346960;  daa[15*20+13] =   0.361819;  daa[15*20+14] =   1.338132; 
2436  daa[16*20+ 0] =   2.139501;  daa[16*20+ 1] =   0.578987;  daa[16*20+ 2] =   2.000679;  daa[16*20+ 3] =   0.425860; 
2437  daa[16*20+ 4] =   1.143480;  daa[16*20+ 5] =   1.080136;  daa[16*20+ 6] =   0.604545;  daa[16*20+ 7] =   0.129836; 
2438  daa[16*20+ 8] =   0.584262;  daa[16*20+ 9] =   1.033739;  daa[16*20+10] =   0.302936;  daa[16*20+11] =   1.136863; 
2439  daa[16*20+12] =   2.020366;  daa[16*20+13] =   0.165001;  daa[16*20+14] =   0.571468;  daa[16*20+15] =   6.472279; 
2440  daa[17*20+ 0] =   0.180717;  daa[17*20+ 1] =   0.593607;  daa[17*20+ 2] =   0.045376;  daa[17*20+ 3] =   0.029890; 
2441  daa[17*20+ 4] =   0.670128;  daa[17*20+ 5] =   0.236199;  daa[17*20+ 6] =   0.077852;  daa[17*20+ 7] =   0.268491; 
2442  daa[17*20+ 8] =   0.597054;  daa[17*20+ 9] =   0.111660;  daa[17*20+10] =   0.619632;  daa[17*20+11] =   0.049906; 
2443  daa[17*20+12] =   0.696175;  daa[17*20+13] =   2.457121;  daa[17*20+14] =   0.095131;  daa[17*20+15] =   0.248862; 
2444  daa[17*20+16] =   0.140825;  daa[18*20+ 0] =   0.218959;  daa[18*20+ 1] =   0.314440;  daa[18*20+ 2] =   0.612025; 
2445  daa[18*20+ 3] =   0.135107;  daa[18*20+ 4] =   1.165532;  daa[18*20+ 5] =   0.257336;  daa[18*20+ 6] =   0.120037; 
2446  daa[18*20+ 7] =   0.054679;  daa[18*20+ 8] =   5.306834;  daa[18*20+ 9] =   0.232523;  daa[18*20+10] =   0.299648; 
2447  daa[18*20+11] =   0.131932;  daa[18*20+12] =   0.481306;  daa[18*20+13] =   7.803902;  daa[18*20+14] =   0.089613; 
2448  daa[18*20+15] =   0.400547;  daa[18*20+16] =   0.245841;  daa[18*20+17] =   3.151815;  daa[19*20+ 0] =   2.547870; 
2449  daa[19*20+ 1] =   0.170887;  daa[19*20+ 2] =   0.083688;  daa[19*20+ 3] =   0.037967;  daa[19*20+ 4] =   1.959291; 
2450  daa[19*20+ 5] =   0.210332;  daa[19*20+ 6] =   0.245034;  daa[19*20+ 7] =   0.076701;  daa[19*20+ 8] =   0.119013; 
2451  daa[19*20+ 9] =  10.649107;  daa[19*20+10] =   1.702745;  daa[19*20+11] =   0.185202;  daa[19*20+12] =   1.898718; 
2452  daa[19*20+13] =   0.654683;  daa[19*20+14] =   0.296501;  daa[19*20+15] =   0.098369;  daa[19*20+16] =   2.188158; 
2453  daa[19*20+17] =   0.189510;  daa[19*20+18] =   0.249313; 
2454 
2455 
2456  for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j];
2457 
2458  pi[0] = 0.079066; pi[1] = 0.055941; pi[2] = 0.041977; pi[3] = 0.053052; 
2459  pi[4] = 0.012937; pi[5] = 0.040767; pi[6] = 0.071586; pi[7] = 0.057337; 
2460  pi[8] = 0.022355; pi[9] = 0.062157; pi[10] = 0.099081; pi[11] = 0.064600; 
2461  pi[12] = 0.022951; pi[13] = 0.042302; pi[14] = 0.044040; pi[15] = 0.061197; 
2462  pi[16] = 0.053287; pi[17] = 0.012066; pi[18] = 0.034155; pi[19] = 0.069147; 
2463 
2464  return 1;
2465}
2466
2467//////////////////////////////////////////////////////////////
2468//////////////////////////////////////////////////////////////
2469
2470
2471int Init_Qmat_WAG(phydbl *daa, phydbl *pi)
2472{
2473  int i,j,naa;
2474
2475  /* WAG's model data
2476   * Simon Whelan and Nick Goldman
2477   * 'A general empirical model of protein evolution derived from multiple
2478   *  protein families using a maximum-likelihood approach'
2479   * MBE (2001) 18:691-699
2480   */
2481
2482
2483  naa = 20;
2484
2485  daa[ 1*20+ 0] =  55.15710; daa[ 2*20+ 0] =  50.98480; daa[ 2*20+ 1] =  63.53460; 
2486  daa[ 3*20+ 0] =  73.89980; daa[ 3*20+ 1] =  14.73040; daa[ 3*20+ 2] = 542.94200; 
2487  daa[ 4*20+ 0] = 102.70400; daa[ 4*20+ 1] =  52.81910; daa[ 4*20+ 2] =  26.52560; 
2488  daa[ 4*20+ 3] =   3.02949; daa[ 5*20+ 0] =  90.85980; daa[ 5*20+ 1] = 303.55000; 
2489  daa[ 5*20+ 2] = 154.36400; daa[ 5*20+ 3] =  61.67830; daa[ 5*20+ 4] =   9.88179; 
2490  daa[ 6*20+ 0] = 158.28500; daa[ 6*20+ 1] =  43.91570; daa[ 6*20+ 2] =  94.71980; 
2491  daa[ 6*20+ 3] = 617.41600; daa[ 6*20+ 4] =   2.13520; daa[ 6*20+ 5] = 546.94700; 
2492  daa[ 7*20+ 0] = 141.67200; daa[ 7*20+ 1] =  58.46650; daa[ 7*20+ 2] = 112.55600; 
2493  daa[ 7*20+ 3] =  86.55840; daa[ 7*20+ 4] =  30.66740; daa[ 7*20+ 5] =  33.00520; 
2494  daa[ 7*20+ 6] =  56.77170; daa[ 8*20+ 0] =  31.69540; daa[ 8*20+ 1] = 213.71500; 
2495  daa[ 8*20+ 2] = 395.62900; daa[ 8*20+ 3] =  93.06760; daa[ 8*20+ 4] =  24.89720; 
2496  daa[ 8*20+ 5] = 429.41100; daa[ 8*20+ 6] =  57.00250; daa[ 8*20+ 7] =  24.94100; 
2497  daa[ 9*20+ 0] =  19.33350; daa[ 9*20+ 1] =  18.69790; daa[ 9*20+ 2] =  55.42360; 
2498  daa[ 9*20+ 3] =   3.94370; daa[ 9*20+ 4] =  17.01350; daa[ 9*20+ 5] =  11.39170; 
2499  daa[ 9*20+ 6] =  12.73950; daa[ 9*20+ 7] =   3.04501; daa[ 9*20+ 8] =  13.81900; 
2500  daa[10*20+ 0] =  39.79150; daa[10*20+ 1] =  49.76710; daa[10*20+ 2] =  13.15280; 
2501  daa[10*20+ 3] =   8.48047; daa[10*20+ 4] =  38.42870; daa[10*20+ 5] =  86.94890; 
2502  daa[10*20+ 6] =  15.42630; daa[10*20+ 7] =   6.13037; daa[10*20+ 8] =  49.94620; 
2503  daa[10*20+ 9] = 317.09700; daa[11*20+ 0] =  90.62650; daa[11*20+ 1] = 535.14200; 
2504  daa[11*20+ 2] = 301.20100; daa[11*20+ 3] =  47.98550; daa[11*20+ 4] =   7.40339; 
2505  daa[11*20+ 5] = 389.49000; daa[11*20+ 6] = 258.44300; daa[11*20+ 7] =  37.35580; 
2506  daa[11*20+ 8] =  89.04320; daa[11*20+ 9] =  32.38320; daa[11*20+10] =  25.75550; 
2507  daa[12*20+ 0] =  89.34960; daa[12*20+ 1] =  68.31620; daa[12*20+ 2] =  19.82210; 
2508  daa[12*20+ 3] =  10.37540; daa[12*20+ 4] =  39.04820; daa[12*20+ 5] = 154.52600; 
2509  daa[12*20+ 6] =  31.51240; daa[12*20+ 7] =  17.41000; daa[12*20+ 8] =  40.41410; 
2510  daa[12*20+ 9] = 425.74600; daa[12*20+10] = 485.40200; daa[12*20+11] =  93.42760; 
2511  daa[13*20+ 0] =  21.04940; daa[13*20+ 1] =  10.27110; daa[13*20+ 2] =   9.61621; 
2512  daa[13*20+ 3] =   4.67304; daa[13*20+ 4] =  39.80200; daa[13*20+ 5] =   9.99208; 
2513  daa[13*20+ 6] =   8.11339; daa[13*20+ 7] =   4.99310; daa[13*20+ 8] =  67.93710; 
2514  daa[13*20+ 9] = 105.94700; daa[13*20+10] = 211.51700; daa[13*20+11] =   8.88360; 
2515  daa[13*20+12] = 119.06300; daa[14*20+ 0] = 143.85500; daa[14*20+ 1] =  67.94890; 
2516  daa[14*20+ 2] =  19.50810; daa[14*20+ 3] =  42.39840; daa[14*20+ 4] =  10.94040; 
2517  daa[14*20+ 5] =  93.33720; daa[14*20+ 6] =  68.23550; daa[14*20+ 7] =  24.35700; 
2518  daa[14*20+ 8] =  69.61980; daa[14*20+ 9] =   9.99288; daa[14*20+10] =  41.58440; 
2519  daa[14*20+11] =  55.68960; daa[14*20+12] =  17.13290; daa[14*20+13] =  16.14440; 
2520  daa[15*20+ 0] = 337.07900; daa[15*20+ 1] = 122.41900; daa[15*20+ 2] = 397.42300; 
2521  daa[15*20+ 3] = 107.17600; daa[15*20+ 4] = 140.76600; daa[15*20+ 5] = 102.88700; 
2522  daa[15*20+ 6] =  70.49390; daa[15*20+ 7] = 134.18200; daa[15*20+ 8] =  74.01690; 
2523  daa[15*20+ 9] =  31.94400; daa[15*20+10] =  34.47390; daa[15*20+11] =  96.71300; 
2524  daa[15*20+12] =  49.39050; daa[15*20+13] =  54.59310; daa[15*20+14] = 161.32800; 
2525  daa[16*20+ 0] = 212.11100; daa[16*20+ 1] =  55.44130; daa[16*20+ 2] = 203.00600; 
2526  daa[16*20+ 3] =  37.48660; daa[16*20+ 4] =  51.29840; daa[16*20+ 5] =  85.79280; 
2527  daa[16*20+ 6] =  82.27650; daa[16*20+ 7] =  22.58330; daa[16*20+ 8] =  47.33070; 
2528  daa[16*20+ 9] = 145.81600; daa[16*20+10] =  32.66220; daa[16*20+11] = 138.69800; 
2529  daa[16*20+12] = 151.61200; daa[16*20+13] =  17.19030; daa[16*20+14] =  79.53840; 
2530  daa[16*20+15] = 437.80200; daa[17*20+ 0] =  11.31330; daa[17*20+ 1] = 116.39200; 
2531  daa[17*20+ 2] =   7.19167; daa[17*20+ 3] =  12.97670; daa[17*20+ 4] =  71.70700; 
2532  daa[17*20+ 5] =  21.57370; daa[17*20+ 6] =  15.65570; daa[17*20+ 7] =  33.69830; 
2533  daa[17*20+ 8] =  26.25690; daa[17*20+ 9] =  21.24830; daa[17*20+10] =  66.53090; 
2534  daa[17*20+11] =  13.75050; daa[17*20+12] =  51.57060; daa[17*20+13] = 152.96400; 
2535  daa[17*20+14] =  13.94050; daa[17*20+15] =  52.37420; daa[17*20+16] =  11.08640; 
2536  daa[18*20+ 0] =  24.07350; daa[18*20+ 1] =  38.15330; daa[18*20+ 2] = 108.60000; 
2537  daa[18*20+ 3] =  32.57110; daa[18*20+ 4] =  54.38330; daa[18*20+ 5] =  22.77100; 
2538  daa[18*20+ 6] =  19.63030; daa[18*20+ 7] =  10.36040; daa[18*20+ 8] = 387.34400; 
2539  daa[18*20+ 9] =  42.01700; daa[18*20+10] =  39.86180; daa[18*20+11] =  13.32640; 
2540  daa[18*20+12] =  42.84370; daa[18*20+13] = 645.42800; daa[18*20+14] =  21.60460; 
2541  daa[18*20+15] =  78.69930; daa[18*20+16] =  29.11480; daa[18*20+17] = 248.53900; 
2542  daa[19*20+ 0] = 200.60100; daa[19*20+ 1] =  25.18490; daa[19*20+ 2] =  19.62460; 
2543  daa[19*20+ 3] =  15.23350; daa[19*20+ 4] = 100.21400; daa[19*20+ 5] =  30.12810; 
2544  daa[19*20+ 6] =  58.87310; daa[19*20+ 7] =  18.72470; daa[19*20+ 8] =  11.83580; 
2545  daa[19*20+ 9] = 782.13000; daa[19*20+10] = 180.03400; daa[19*20+11] =  30.54340; 
2546  daa[19*20+12] = 205.84500; daa[19*20+13] =  64.98920; daa[19*20+14] =  31.48870; 
2547  daa[19*20+15] =  23.27390; daa[19*20+16] = 138.82300; daa[19*20+17] =  36.53690; 
2548  daa[19*20+18] =  31.47300; 
2549
2550  for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j];
2551 
2552  pi[0] = 0.0866279; pi[1] =  0.043972; pi[2] =  0.0390894; pi[3] =  0.0570451;
2553  pi[4] =  0.0193078; pi[5] =  0.0367281; pi[6] =  0.0580589; pi[7] =  0.0832518;
2554  pi[8] =  0.0244313; pi[9] =  0.048466; pi[10] =  0.086209; pi[11] = 0.0620286;
2555  pi[12] = 0.0195027; pi[13] =  0.0384319; pi[14] =  0.0457631; pi[15] = 0.0695179;
2556  pi[16] =  0.0610127; pi[17] =  0.0143859; pi[18] =  0.0352742; pi[19] =  0.0708956;
2557
2558
2559  return 1;
2560}
2561
2562//////////////////////////////////////////////////////////////
2563//////////////////////////////////////////////////////////////
2564
2565
2566int Init_Qmat_RtREV(phydbl *daa, phydbl *pi)
2567{
2568    /*
2569       This model has been 'translated' from John Huelsenbeck and Fredrik Ronquist
2570       MrBayes program into PHYML format by Federico Abascal. Many thanks to them.
2571    */
2572   
2573    /*
2574    Dimmic M.W., J.S. Rest, D.P. Mindell, and D. Goldstein. 2002. RArtREV:
2575    An amino acid substitution matrix for inference of retrovirus and
2576    reverse transcriptase phyLOGeny. Journal of Molecular Evolution
2577    55: 65-73.
2578    */
2579
2580    int i,j,naa;
2581    naa = 20;
2582
2583    daa[1*20+0]= 34;         daa[2*20+0]= 51;         daa[2*20+1]= 35;         daa[3*20+0]= 10;         
2584    daa[3*20+1]= 30;         daa[3*20+2]= 384;        daa[4*20+0]= 439;        daa[4*20+1]= 92;         
2585    daa[4*20+2]= 128;        daa[4*20+3]= 1;          daa[5*20+0]= 32;         daa[5*20+1]= 221;       
2586    daa[5*20+2]= 236;        daa[5*20+3]= 78;         daa[5*20+4]= 70;         daa[6*20+0]= 81;         
2587    daa[6*20+1]= 10;         daa[6*20+2]= 79;         daa[6*20+3]= 542;        daa[6*20+4]= 1;         
2588    daa[6*20+5]= 372;        daa[7*20+0]= 135;        daa[7*20+1]= 41;         daa[7*20+2]= 94;         
2589    daa[7*20+3]= 61;         daa[7*20+4]= 48;         daa[7*20+5]= 18;         daa[7*20+6]= 70;         
2590    daa[8*20+0]= 30;         daa[8*20+1]= 90;         daa[8*20+2]= 320;        daa[8*20+3]= 91;         
2591    daa[8*20+4]= 124;        daa[8*20+5]= 387;        daa[8*20+6]= 34;         daa[8*20+7]= 68;         
2592    daa[9*20+0]= 1;          daa[9*20+1]= 24;         daa[9*20+2]= 35;         daa[9*20+3]= 1;         
2593    daa[9*20+4]= 104;        daa[9*20+5]= 33;         daa[9*20+6]= 1;          daa[9*20+7]= 1;         
2594    daa[9*20+8]= 34;         daa[10*20+0]= 45;        daa[10*20+1]= 18;        daa[10*20+2]= 15;       
2595    daa[10*20+3]= 5;         daa[10*20+4]= 110;       daa[10*20+5]= 54;        daa[10*20+6]= 21;       
2596    daa[10*20+7]= 3;         daa[10*20+8]= 51;        daa[10*20+9]= 385;       daa[11*20+0]= 38;       
2597    daa[11*20+1]= 593;       daa[11*20+2]= 123;       daa[11*20+3]= 20;        daa[11*20+4]= 16;       
2598    daa[11*20+5]= 309;       daa[11*20+6]= 141;       daa[11*20+7]= 30;        daa[11*20+8]= 76;       
2599    daa[11*20+9]= 34;        daa[11*20+10]= 23;       daa[12*20+0]= 235;       daa[12*20+1]= 57;       
2600    daa[12*20+2]= 1;         daa[12*20+3]= 1;         daa[12*20+4]= 156;       daa[12*20+5]= 158;       
2601    daa[12*20+6]= 1;         daa[12*20+7]= 37;        daa[12*20+8]= 116;       daa[12*20+9]= 375;       
2602    daa[12*20+10]= 581;      daa[12*20+11]= 134;      daa[13*20+0]= 1;         daa[13*20+1]= 7;         
2603    daa[13*20+2]= 49;        daa[13*20+3]= 1;         daa[13*20+4]= 70;        daa[13*20+5]= 1;         
2604    daa[13*20+6]= 1;         daa[13*20+7]= 7;         daa[13*20+8]= 141;       daa[13*20+9]= 64;       
2605    daa[13*20+10]= 179;      daa[13*20+11]= 14;       daa[13*20+12]= 247;      daa[14*20+0]= 97;       
2606    daa[14*20+1]= 24;        daa[14*20+2]= 33;        daa[14*20+3]= 55;        daa[14*20+4]= 1;         
2607    daa[14*20+5]= 68;        daa[14*20+6]= 52;        daa[14*20+7]= 17;        daa[14*20+8]= 44;       
2608    daa[14*20+9]= 10;        daa[14*20+10]= 22;       daa[14*20+11]= 43;       daa[14*20+12]= 1;       
2609    daa[14*20+13]= 11;       daa[15*20+0]= 460;       daa[15*20+1]= 102;       daa[15*20+2]= 294;       
2610    daa[15*20+3]= 136;       daa[15*20+4]= 75;        daa[15*20+5]= 225;       daa[15*20+6]= 95;       
2611    daa[15*20+7]= 152;       daa[15*20+8]= 183;       daa[15*20+9]= 4;         daa[15*20+10]= 24;       
2612    daa[15*20+11]= 77;       daa[15*20+12]= 1;        daa[15*20+13]= 20;       daa[15*20+14]= 134;     
2613    daa[16*20+0]= 258;       daa[16*20+1]= 64;        daa[16*20+2]= 148;       daa[16*20+3]= 55;       
2614    daa[16*20+4]= 117;       daa[16*20+5]= 146;       daa[16*20+6]= 82;        daa[16*20+7]= 7;         
2615    daa[16*20+8]= 49;        daa[16*20+9]= 72;        daa[16*20+10]= 25;       daa[16*20+11]= 110;     
2616    daa[16*20+12]= 131;      daa[16*20+13]= 69;       daa[16*20+14]= 62;       daa[16*20+15]= 671;     
2617    daa[17*20+0]= 5;         daa[17*20+1]= 13;        daa[17*20+2]= 16;        daa[17*20+3]= 1;         
2618    daa[17*20+4]= 55;        daa[17*20+5]= 10;        daa[17*20+6]= 17;        daa[17*20+7]= 23;       
2619    daa[17*20+8]= 48;        daa[17*20+9]= 39;        daa[17*20+10]= 47;       daa[17*20+11]= 6;       
2620    daa[17*20+12]= 111;      daa[17*20+13]= 182;      daa[17*20+14]= 9;        daa[17*20+15]= 14;       
2621    daa[17*20+16]= 1;        daa[18*20+0]= 55;        daa[18*20+1]= 47;        daa[18*20+2]= 28;       
2622    daa[18*20+3]= 1;         daa[18*20+4]= 131;       daa[18*20+5]= 45;        daa[18*20+6]= 1;         
2623    daa[18*20+7]= 21;        daa[18*20+8]= 307;       daa[18*20+9]= 26;        daa[18*20+10]= 64;       
2624    daa[18*20+11]= 1;        daa[18*20+12]= 74;       daa[18*20+13]= 1017;     daa[18*20+14]= 14;       
2625    daa[18*20+15]= 31;       daa[18*20+16]= 34;       daa[18*20+17]= 176;      daa[19*20+0]= 197;       
2626    daa[19*20+1]= 29;        daa[19*20+2]= 21;        daa[19*20+3]= 6;         daa[19*20+4]= 295;       
2627    daa[19*20+5]= 36;        daa[19*20+6]= 35;        daa[19*20+7]= 3;         daa[19*20+8]= 1;         
2628    daa[19*20+9]= 1048;      daa[19*20+10]= 112;      daa[19*20+11]= 19;       daa[19*20+12]= 236;     
2629    daa[19*20+13]= 92;       daa[19*20+14]= 25;       daa[19*20+15]= 39;       daa[19*20+16]= 196;     
2630    daa[19*20+17]= 26;       daa[19*20+18]= 59;       
2631
2632    for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j];
2633
2634    pi[0]= 0.0646;           pi[1]= 0.0453;           pi[2]= 0.0376;           pi[3]= 0.0422;           
2635    pi[4]= 0.0114;           pi[5]= 0.0606;           pi[6]= 0.0607;           pi[7]= 0.0639;           
2636    pi[8]= 0.0273;           pi[9]= 0.0679;           pi[10]= 0.1018;          pi[11]= 0.0751;         
2637    pi[12]= 0.015;           pi[13]= 0.0287;          pi[14]= 0.0681;          pi[15]= 0.0488;         
2638    pi[16]= 0.0622;          pi[17]= 0.0251;          pi[18]= 0.0318;          pi[19]= 0.0619;         
2639    return 1;
2640}
2641
2642//////////////////////////////////////////////////////////////
2643//////////////////////////////////////////////////////////////
2644
2645
2646int Init_Qmat_CpREV(phydbl *daa, phydbl *pi)
2647{
2648    /*
2649       This model has been 'translated' from John Huelsenbeck and Fredrik Ronquist
2650       MrBayes program into PHYML format by Federico Abascal. Many thanks to them.
2651    */
2652
2653    /*
2654      Adachi, J., P. Waddell, W. Martin, and M. Hasegawa. 2000. Plastid         
2655      genome phyLOGeny and a model of amino acid substitution for proteins   
2656      encoded by chloroplast DNA. Journal of Molecular Evolution             
2657      50:348-358.
2658    */
2659
2660    int i,j,naa;
2661    naa = 20;
2662
2663    daa[1*20+0]= 105;        daa[2*20+0]= 227;        daa[2*20+1]= 357;        daa[3*20+0]= 175;       
2664    daa[3*20+1]= 43;         daa[3*20+2]= 4435;       daa[4*20+0]= 669;        daa[4*20+1]= 823;       
2665    daa[4*20+2]= 538;        daa[4*20+3]= 10;         daa[5*20+0]= 157;        daa[5*20+1]= 1745;       
2666    daa[5*20+2]= 768;        daa[5*20+3]= 400;        daa[5*20+4]= 10;         daa[6*20+0]= 499;       
2667    daa[6*20+1]= 152;        daa[6*20+2]= 1055;       daa[6*20+3]= 3691;       daa[6*20+4]= 10;         
2668    daa[6*20+5]= 3122;       daa[7*20+0]= 665;        daa[7*20+1]= 243;        daa[7*20+2]= 653;       
2669    daa[7*20+3]= 431;        daa[7*20+4]= 303;        daa[7*20+5]= 133;        daa[7*20+6]= 379;       
2670    daa[8*20+0]= 66;         daa[8*20+1]= 715;        daa[8*20+2]= 1405;       daa[8*20+3]= 331;       
2671    daa[8*20+4]= 441;        daa[8*20+5]= 1269;       daa[8*20+6]= 162;        daa[8*20+7]= 19;         
2672    daa[9*20+0]= 145;        daa[9*20+1]= 136;        daa[9*20+2]= 168;        daa[9*20+3]= 10;         
2673    daa[9*20+4]= 280;        daa[9*20+5]= 92;         daa[9*20+6]= 148;        daa[9*20+7]= 40;         
2674    daa[9*20+8]= 29;         daa[10*20+0]= 197;       daa[10*20+1]= 203;       daa[10*20+2]= 113;       
2675    daa[10*20+3]= 10;        daa[10*20+4]= 396;       daa[10*20+5]= 286;       daa[10*20+6]= 82;       
2676    daa[10*20+7]= 20;        daa[10*20+8]= 66;        daa[10*20+9]= 1745;      daa[11*20+0]= 236;       
2677    daa[11*20+1]= 4482;      daa[11*20+2]= 2430;      daa[11*20+3]= 412;       daa[11*20+4]= 48;       
2678    daa[11*20+5]= 3313;      daa[11*20+6]= 2629;      daa[11*20+7]= 263;       daa[11*20+8]= 305;       
2679    daa[11*20+9]= 345;       daa[11*20+10]= 218;      daa[12*20+0]= 185;       daa[12*20+1]= 125;       
2680    daa[12*20+2]= 61;        daa[12*20+3]= 47;        daa[12*20+4]= 159;       daa[12*20+5]= 202;       
2681    daa[12*20+6]= 113;       daa[12*20+7]= 21;        daa[12*20+8]= 10;        daa[12*20+9]= 1772;     
2682    daa[12*20+10]= 1351;     daa[12*20+11]= 193;      daa[13*20+0]= 68;        daa[13*20+1]= 53;       
2683    daa[13*20+2]= 97;        daa[13*20+3]= 22;        daa[13*20+4]= 726;       daa[13*20+5]= 10;       
2684    daa[13*20+6]= 145;       daa[13*20+7]= 25;        daa[13*20+8]= 127;       daa[13*20+9]= 454;       
2685    daa[13*20+10]= 1268;     daa[13*20+11]= 72;       daa[13*20+12]= 327;      daa[14*20+0]= 490;       
2686    daa[14*20+1]= 87;        daa[14*20+2]= 173;       daa[14*20+3]= 170;       daa[14*20+4]= 285;       
2687    daa[14*20+5]= 323;       daa[14*20+6]= 185;       daa[14*20+7]= 28;        daa[14*20+8]= 152;       
2688    daa[14*20+9]= 117;       daa[14*20+10]= 219;      daa[14*20+11]= 302;      daa[14*20+12]= 100;     
2689    daa[14*20+13]= 43;       daa[15*20+0]= 2440;      daa[15*20+1]= 385;       daa[15*20+2]= 2085;     
2690    daa[15*20+3]= 590;       daa[15*20+4]= 2331;      daa[15*20+5]= 396;       daa[15*20+6]= 568;       
2691    daa[15*20+7]= 691;       daa[15*20+8]= 303;       daa[15*20+9]= 216;       daa[15*20+10]= 516;     
2692    daa[15*20+11]= 868;      daa[15*20+12]= 93;       daa[15*20+13]= 487;      daa[15*20+14]= 1202;     
2693    daa[16*20+0]= 1340;      daa[16*20+1]= 314;       daa[16*20+2]= 1393;      daa[16*20+3]= 266;       
2694    daa[16*20+4]= 576;       daa[16*20+5]= 241;       daa[16*20+6]= 369;       daa[16*20+7]= 92;       
2695    daa[16*20+8]= 32;        daa[16*20+9]= 1040;      daa[16*20+10]= 156;      daa[16*20+11]= 918;     
2696    daa[16*20+12]= 645;      daa[16*20+13]= 148;      daa[16*20+14]= 260;      daa[16*20+15]= 2151;     
2697    daa[17*20+0]= 14;        daa[17*20+1]= 230;       daa[17*20+2]= 40;        daa[17*20+3]= 18;       
2698    daa[17*20+4]= 435;       daa[17*20+5]= 53;        daa[17*20+6]= 63;        daa[17*20+7]= 82;       
2699    daa[17*20+8]= 69;        daa[17*20+9]= 42;        daa[17*20+10]= 159;      daa[17*20+11]= 10;       
2700    daa[17*20+12]= 86;       daa[17*20+13]= 468;      daa[17*20+14]= 49;       daa[17*20+15]= 73;       
2701    daa[17*20+16]= 29;       daa[18*20+0]= 56;        daa[18*20+1]= 323;       daa[18*20+2]= 754;       
2702    daa[18*20+3]= 281;       daa[18*20+4]= 1466;      daa[18*20+5]= 391;       daa[18*20+6]= 142;       
2703    daa[18*20+7]= 10;        daa[18*20+8]= 1971;      daa[18*20+9]= 89;        daa[18*20+10]= 189;     
2704    daa[18*20+11]= 247;      daa[18*20+12]= 215;      daa[18*20+13]= 2370;     daa[18*20+14]= 97;       
2705    daa[18*20+15]= 522;      daa[18*20+16]= 71;       daa[18*20+17]= 346;      daa[19*20+0]= 968;       
2706    daa[19*20+1]= 92;        daa[19*20+2]= 83;        daa[19*20+3]= 75;        daa[19*20+4]= 592;       
2707    daa[19*20+5]= 54;        daa[19*20+6]= 200;       daa[19*20+7]= 91;        daa[19*20+8]= 25;       
2708    daa[19*20+9]= 4797;      daa[19*20+10]= 865;      daa[19*20+11]= 249;      daa[19*20+12]= 475;     
2709    daa[19*20+13]= 317;      daa[19*20+14]= 122;      daa[19*20+15]= 167;      daa[19*20+16]= 760;     
2710    daa[19*20+17]= 10;       daa[19*20+18]= 119;     
2711
2712    for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j];
2713
2714    pi[0]= 0.076;            pi[1]= 0.062;            pi[2]= 0.041;            pi[3]= 0.037;           
2715    pi[4]= 0.009;            pi[5]= 0.038;            pi[6]= 0.049;            pi[7]= 0.084;           
2716    pi[8]= 0.025;            pi[9]= 0.081;            pi[10]= 0.101;           pi[11]= 0.05;           
2717    pi[12]= 0.022;           pi[13]= 0.051;           pi[14]= 0.043;           pi[15]= 0.062;           
2718    pi[16]= 0.054;           pi[17]= 0.018;           pi[18]= 0.031;           pi[19]= 0.066;           
2719    return 1;
2720}
2721
2722//////////////////////////////////////////////////////////////
2723//////////////////////////////////////////////////////////////
2724
2725
2726int Init_Qmat_VT(phydbl *daa, phydbl *pi)
2727{
2728    /*
2729       This model has been 'translated' from John Huelsenbeck and Fredrik Ronquist
2730       MrBayes program into PHYML format by Federico Abascal. Many thanks to them.
2731    */
2732
2733    /*
2734      Muller, T., and M. Vingron. 2000. Modeling amino acid replacement.         
2735      Journal of Computational Biology 7:761-776.                             
2736    */
2737
2738    int i,j,naa;
2739    naa = 20;
2740
2741    /* daa[1*20+0]= 0.233108;   daa[2*20+0]= 0.199097;   daa[2*20+1]= 0.210797;   daa[3*20+0]= 0.265145;    */
2742    /* daa[3*20+1]= 0.105191;   daa[3*20+2]= 0.883422;   daa[4*20+0]= 0.227333;   daa[4*20+1]= 0.031726;    */
2743    /* daa[4*20+2]= 0.027495;   daa[4*20+3]= 0.010313;   daa[5*20+0]= 0.310084;   daa[5*20+1]= 0.493763;    */
2744    /* daa[5*20+2]= 0.2757;     daa[5*20+3]= 0.205842;   daa[5*20+4]= 0.004315;   daa[6*20+0]= 0.567957;    */
2745    /* daa[6*20+1]= 0.25524;    daa[6*20+2]= 0.270417;   daa[6*20+3]= 1.599461;   daa[6*20+4]= 0.005321;    */
2746    /* daa[6*20+5]= 0.960976;   daa[7*20+0]= 0.876213;   daa[7*20+1]= 0.156945;   daa[7*20+2]= 0.362028;    */
2747    /* daa[7*20+3]= 0.311718;   daa[7*20+4]= 0.050876;   daa[7*20+5]= 0.12866;    daa[7*20+6]= 0.250447;    */
2748    /* daa[8*20+0]= 0.078692;   daa[8*20+1]= 0.213164;   daa[8*20+2]= 0.290006;   daa[8*20+3]= 0.134252;    */
2749    /* daa[8*20+4]= 0.016695;   daa[8*20+5]= 0.315521;   daa[8*20+6]= 0.104458;   daa[8*20+7]= 0.058131;    */
2750    /* daa[9*20+0]= 0.222972;   daa[9*20+1]= 0.08151;    daa[9*20+2]= 0.087225;   daa[9*20+3]= 0.01172;     */
2751    /* daa[9*20+4]= 0.046398;   daa[9*20+5]= 0.054602;   daa[9*20+6]= 0.046589;   daa[9*20+7]= 0.051089;    */
2752    /* daa[9*20+8]= 0.020039;   daa[10*20+0]= 0.42463;   daa[10*20+1]= 0.192364;  daa[10*20+2]= 0.069245;   */
2753    /* daa[10*20+3]= 0.060863;  daa[10*20+4]= 0.091709;  daa[10*20+5]= 0.24353;   daa[10*20+6]= 0.151924;   */
2754    /* daa[10*20+7]= 0.087056;  daa[10*20+8]= 0.103552;  daa[10*20+9]= 2.08989;   daa[11*20+0]= 0.393245;   */
2755    /* daa[11*20+1]= 1.755838;  daa[11*20+2]= 0.50306;   daa[11*20+3]= 0.261101;  daa[11*20+4]= 0.004067;   */
2756    /* daa[11*20+5]= 0.738208;  daa[11*20+6]= 0.88863;   daa[11*20+7]= 0.193243;  daa[11*20+8]= 0.153323;   */
2757    /* daa[11*20+9]= 0.093181;  daa[11*20+10]= 0.201204; daa[12*20+0]= 0.21155;   daa[12*20+1]= 0.08793;    */
2758    /* daa[12*20+2]= 0.05742;   daa[12*20+3]= 0.012182;  daa[12*20+4]= 0.02369;   daa[12*20+5]= 0.120801;   */
2759    /* daa[12*20+6]= 0.058643;  daa[12*20+7]= 0.04656;   daa[12*20+8]= 0.021157;  daa[12*20+9]= 0.493845;   */
2760    /* daa[12*20+10]= 1.105667; daa[12*20+11]= 0.096474; daa[13*20+0]= 0.116646;  daa[13*20+1]= 0.042569;   */
2761    /* daa[13*20+2]= 0.039769;  daa[13*20+3]= 0.016577;  daa[13*20+4]= 0.051127;  daa[13*20+5]= 0.026235;   */
2762    /* daa[13*20+6]= 0.028168;  daa[13*20+7]= 0.050143;  daa[13*20+8]= 0.079807;  daa[13*20+9]= 0.32102;    */
2763    /* daa[13*20+10]= 0.946499; daa[13*20+11]= 0.038261; daa[13*20+12]= 0.173052; daa[14*20+0]= 0.399143;   */
2764    /* daa[14*20+1]= 0.12848;   daa[14*20+2]= 0.083956;  daa[14*20+3]= 0.160063;  daa[14*20+4]= 0.011137;   */
2765    /* daa[14*20+5]= 0.15657;   daa[14*20+6]= 0.205134;  daa[14*20+7]= 0.124492;  daa[14*20+8]= 0.078892;   */
2766    /* daa[14*20+9]= 0.054797;  daa[14*20+10]= 0.169784; daa[14*20+11]= 0.212302; daa[14*20+12]= 0.010363;  */
2767    /* daa[14*20+13]= 0.042564; daa[15*20+0]= 1.817198;  daa[15*20+1]= 0.292327;  daa[15*20+2]= 0.847049;   */
2768    /* daa[15*20+3]= 0.461519;  daa[15*20+4]= 0.17527;   daa[15*20+5]= 0.358017;  daa[15*20+6]= 0.406035;   */
2769    /* daa[15*20+7]= 0.612843;  daa[15*20+8]= 0.167406;  daa[15*20+9]= 0.081567;  daa[15*20+10]= 0.214977;  */
2770    /* daa[15*20+11]= 0.400072; daa[15*20+12]= 0.090515; daa[15*20+13]= 0.138119; daa[15*20+14]= 0.430431;  */
2771    /* daa[16*20+0]= 0.877877;  daa[16*20+1]= 0.204109;  daa[16*20+2]= 0.471268;  daa[16*20+3]= 0.178197;   */
2772    /* daa[16*20+4]= 0.079511;  daa[16*20+5]= 0.248992;  daa[16*20+6]= 0.321028;  daa[16*20+7]= 0.136266;   */
2773    /* daa[16*20+8]= 0.101117;  daa[16*20+9]= 0.376588;  daa[16*20+10]= 0.243227; daa[16*20+11]= 0.446646;  */
2774    /* daa[16*20+12]= 0.184609; daa[16*20+13]= 0.08587;  daa[16*20+14]= 0.207143; daa[16*20+15]= 1.767766;  */
2775    /* daa[17*20+0]= 0.030309;  daa[17*20+1]= 0.046417;  daa[17*20+2]= 0.010459;  daa[17*20+3]= 0.011393;   */
2776    /* daa[17*20+4]= 0.007732;  daa[17*20+5]= 0.021248;  daa[17*20+6]= 0.018844;  daa[17*20+7]= 0.02399;    */
2777    /* daa[17*20+8]= 0.020009;  daa[17*20+9]= 0.034954;  daa[17*20+10]= 0.083439; daa[17*20+11]= 0.023321;  */
2778    /* daa[17*20+12]= 0.022019; daa[17*20+13]= 0.12805;  daa[17*20+14]= 0.014584; daa[17*20+15]= 0.035933;  */
2779    /* daa[17*20+16]= 0.020437; daa[18*20+0]= 0.087061;  daa[18*20+1]= 0.09701;   daa[18*20+2]= 0.093268;   */
2780    /* daa[18*20+3]= 0.051664;  daa[18*20+4]= 0.042823;  daa[18*20+5]= 0.062544;  daa[18*20+6]= 0.0552;     */
2781    /* daa[18*20+7]= 0.037568;  daa[18*20+8]= 0.286027;  daa[18*20+9]= 0.086237;  daa[18*20+10]= 0.189842;  */
2782    /* daa[18*20+11]= 0.068689; daa[18*20+12]= 0.073223; daa[18*20+13]= 0.898663; daa[18*20+14]= 0.032043;  */
2783    /* daa[18*20+15]= 0.121979; daa[18*20+16]= 0.094617; daa[18*20+17]= 0.124746; daa[19*20+0]= 1.230985;   */
2784    /* daa[19*20+1]= 0.113146;  daa[19*20+2]= 0.049824;  daa[19*20+3]= 0.048769;  daa[19*20+4]= 0.163831;   */
2785    /* daa[19*20+5]= 0.112027;  daa[19*20+6]= 0.205868;  daa[19*20+7]= 0.082579;  daa[19*20+8]= 0.068575;   */
2786    /* daa[19*20+9]= 3.65443;   daa[19*20+10]= 1.337571; daa[19*20+11]= 0.144587; daa[19*20+12]= 0.307309;  */
2787    /* daa[19*20+13]= 0.247329; daa[19*20+14]= 0.129315; daa[19*20+15]= 0.1277;   daa[19*20+16]= 0.740372;  */
2788    /* daa[19*20+17]= 0.022134; daa[19*20+18]= 0.125733;  */
2789
2790    /* for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j]; */
2791
2792    /* pi[0]= 0.078837;         pi[1]= 0.051238;         pi[2]= 0.042313;         pi[3]= 0.053066;          */
2793    /* pi[4]= 0.015175;         pi[5]= 0.036713;         pi[6]= 0.061924;         pi[7]= 0.070852;          */
2794    /* pi[8]= 0.023082;         pi[9]= 0.062056;         pi[10]= 0.096371;        pi[11]= 0.057324;         */
2795    /* pi[12]= 0.023771;        pi[13]= 0.043296;        pi[14]= 0.043911;        pi[15]= 0.063403;         */
2796    /* pi[16]= 0.055897;        pi[17]= 0.013272;        pi[18]= 0.034399;        pi[19]= 0.073101;         */
2797    /* return 1; */
2798
2799    /*
2800       This model has been 'translated' from TREE-PUZZLE. Many thanks to  Heiko A. Schmidt,
2801       Korbinian Strimmer, and Arndt von Haeseler.
2802    */
2803
2804    /*
2805      Muller, T., and M. Vingron. 2000. Modeling amino acid replacement.         
2806      Journal of Computational Biology 7:761-776.                             
2807    */
2808
2809    daa[ 0*20+ 1] = 1.2412691067876198;  daa[ 0*20+ 2] = 1.2184237953498958;
2810    daa[ 0*20+ 3] = 1.3759368509441177;  daa[ 0*20+ 4] = 2.4731223087544874;
2811    daa[ 0*20+ 5] = 2.2155167805137470;  daa[ 0*20+ 6] = 2.3379911207495061;
2812    daa[ 0*20+ 7] = 3.3386555146457697;  daa[ 0*20+ 8] = 0.9615841926910841;
2813    daa[ 0*20+ 9] = 0.8908203061925510;  daa[ 0*20+10] = 1.0778497408764076;
2814    daa[ 0*20+11] = 1.4932055816372476;  daa[ 0*20+12] = 1.9006455961717605;
2815    daa[ 0*20+13] = 0.6883439026872615;  daa[ 0*20+14] = 2.7355620089953550;
2816    daa[ 0*20+15] = 6.4208961859142883;  daa[ 0*20+16] = 5.2892514169776437;
2817    daa[ 0*20+17] = 0.5488578478106930;  daa[ 0*20+18] = 0.5411769916657778;
2818    daa[ 0*20+19] = 4.6501894691803214;
2819   
2820    daa[ 1*20+ 2] = 1.5720770753326880;  daa[ 1*20+ 3] = 0.7550654439001206;
2821    daa[ 1*20+ 4] = 1.4414262567428417;  daa[ 1*20+ 5] = 5.5120819705248678;
2822    daa[ 1*20+ 6] = 1.3542404860613146;  daa[ 1*20+ 7] = 1.3121700301622004;
2823    daa[ 1*20+ 8] = 4.9238668283945266;  daa[ 1*20+ 9] = 0.4323005487925516;
2824    daa[ 1*20+10] = 0.8386701149158265;  daa[ 1*20+11] = 10.0173308173660018;
2825    daa[ 1*20+12] = 1.2488638689609959;  daa[ 1*20+13] = 0.4224945197276290;
2826    daa[ 1*20+14] = 1.3091837782420783;  daa[ 1*20+15] = 1.9202994262316166;
2827    daa[ 1*20+16] = 1.3363401740560601;  daa[ 1*20+17] = 1.5170142153962840;
2828    daa[ 1*20+18] = 0.8912614404565405;  daa[ 1*20+19] = 0.7807017855806767;
2829   
2830    daa[ 2*20+ 3] = 7.8584219153689405;  daa[ 2*20+ 4] = 0.9784679122774127;
2831    daa[ 2*20+ 5] = 3.0143201670924822;  daa[ 2*20+ 6] = 2.0093434778398112;
2832    daa[ 2*20+ 7] = 2.4117632898861809;  daa[ 2*20+ 8] = 6.1974384977884114;
2833    daa[ 2*20+ 9] = 0.9179291175331520;  daa[ 2*20+10] = 0.4098311270816011;
2834    daa[ 2*20+11] = 4.4034547578962568;  daa[ 2*20+12] = 0.9378803706165143;
2835    daa[ 2*20+13] = 0.5044944273324311;  daa[ 2*20+14] = 0.7103720531974738;
2836    daa[ 2*20+15] = 6.1234512396801764;  daa[ 2*20+16] = 3.8852506105922231;
2837    daa[ 2*20+17] = 0.1808525752605976;  daa[ 2*20+18] = 1.0894926581511342;
2838    daa[ 2*20+19] = 0.4586061981719967;
2839   
2840    daa[ 3*20+ 4] = 0.2272488448121475;  daa[ 3*20+ 5] = 1.6562495638176040;
2841    daa[ 3*20+ 6] = 9.6883451875685065;  daa[ 3*20+ 7] = 1.9142079025990228;
2842    daa[ 3*20+ 8] = 2.1459640610133781;  daa[ 3*20+ 9] = 0.2161660372725585;
2843    daa[ 3*20+10] = 0.3574207468998517;  daa[ 3*20+11] = 1.4521790561663968;
2844    daa[ 3*20+12] = 0.4075239926000898;  daa[ 3*20+13] = 0.1675129724559251;
2845    daa[ 3*20+14] = 1.0714605979577547;  daa[ 3*20+15] = 2.2161944596741829;
2846    daa[ 3*20+16] = 1.5066839872944762;  daa[ 3*20+17] = 0.2496584188151770;
2847    daa[ 3*20+18] = 0.7447620891784513;  daa[ 3*20+19] = 0.4594535241660911;
2848   
2849    daa[ 4*20+ 5] = 0.4587469126746136;  daa[ 4*20+ 6] = 0.4519167943192672;
2850    daa[ 4*20+ 7] = 1.1034605684472507;  daa[ 4*20+ 8] = 1.5196756759380692;
2851    daa[ 4*20+ 9] = 0.9126668032539315;  daa[ 4*20+10] = 1.4081315998413697;
2852    daa[ 4*20+11] = 0.3371091785647479;  daa[ 4*20+12] = 1.2213054800811556;
2853    daa[ 4*20+13] = 1.6953951980808002;  daa[ 4*20+14] = 0.4326227078645523;
2854    daa[ 4*20+15] = 3.6366815408744255;  daa[ 4*20+16] = 1.7557065205837685;
2855    daa[ 4*20+17] = 1.6275179891253113;  daa[ 4*20+18] = 2.1579775140421025;
2856    daa[ 4*20+19] = 2.2627456996290891;
2857   
2858    daa[ 5*20+ 6] = 6.8124601839937675;  daa[ 5*20+ 7] = 0.8776110594765502;
2859    daa[ 5*20+ 8] = 7.9943228564946525;  daa[ 5*20+ 9] = 0.4882733432879921;
2860    daa[ 5*20+10] = 1.3318097154194044;  daa[ 5*20+11] = 6.0519085243118811;
2861    daa[ 5*20+12] = 1.9106190827629084;  daa[ 5*20+13] = 0.3573432522499545;
2862    daa[ 5*20+14] = 2.3019177728300728;  daa[ 5*20+15] = 2.3193703643237220;
2863    daa[ 5*20+16] = 2.1576510103471440;  daa[ 5*20+17] = 0.8959082681546182;
2864    daa[ 5*20+18] = 0.9183596801412757;  daa[ 5*20+19] = 0.6366932501396869;
2865   
2866    daa[ 6*20+ 7] = 1.3860121390169038;  daa[ 6*20+ 8] = 1.6360079688522375;
2867    daa[ 6*20+ 9] = 0.4035497929633328;  daa[ 6*20+10] = 0.5610717242294755;
2868    daa[ 6*20+11] = 4.3290086529582830;  daa[ 6*20+12] = 0.7471936218068498;
2869    daa[ 6*20+13] = 0.2317194387691585;  daa[ 6*20+14] = 1.5132807416252063;
2870    daa[ 6*20+15] = 1.8273535587773553;  daa[ 6*20+16] = 1.5839981708584689;
2871    daa[ 6*20+17] = 0.4198391148111098;  daa[ 6*20+18] = 0.5818111331782764;
2872    daa[ 6*20+19] = 0.8940572875547330;
2873   
2874    daa[ 7*20+ 8] = 0.8561248973045037;  daa[ 7*20+ 9] = 0.2888075033037488;
2875    daa[ 7*20+10] = 0.3578662395745526;  daa[ 7*20+11] = 0.8945563662345198;
2876    daa[ 7*20+12] = 0.5954812791740037;  daa[ 7*20+13] = 0.3693722640980460;
2877    daa[ 7*20+14] = 0.7744933618134962;  daa[ 7*20+15] = 3.0637776193717610;
2878    daa[ 7*20+16] = 0.7147489676267383;  daa[ 7*20+17] = 0.9349753595598769;
2879    daa[ 7*20+18] = 0.3374467649724478;  daa[ 7*20+19] = 0.6193321034173915;
2880   
2881    daa[ 8*20+ 9] = 0.5787937115407940;  daa[ 8*20+10] = 1.0765007949562073;
2882    daa[ 8*20+11] = 1.8085136096039203;  daa[ 8*20+12] = 1.3808291710019667;
2883    daa[ 8*20+13] = 1.3629765501081097;  daa[ 8*20+14] = 1.8370555852070649;
2884    daa[ 8*20+15] = 1.9699895187387506;  daa[ 8*20+16] = 1.6136654573285647;
2885    daa[ 8*20+17] = 0.6301954684360302;  daa[ 8*20+18] = 7.7587442309146040;
2886    daa[ 8*20+19] = 0.5333220944030346;
2887   
2888    daa[ 9*20+10] = 6.0019110258426362;  daa[ 9*20+11] = 0.6244297525127139;
2889    daa[ 9*20+12] = 6.7597899772045418;  daa[ 9*20+13] = 2.2864286949316077;
2890    daa[ 9*20+14] = 0.4811402387911145;  daa[ 9*20+15] = 0.6047491507504744;
2891    daa[ 9*20+16] = 2.6344778384442731;  daa[ 9*20+17] = 0.5604648274060783;
2892    daa[ 9*20+18] = 0.8626796044156272;  daa[ 9*20+19] = 14.8729334615190609;
2893   
2894    daa[10*20+11] = 0.5642322882556321;  daa[10*20+12] = 8.0327792947421148;
2895    daa[10*20+13] = 4.3611548063555778;  daa[10*20+14] = 1.0084320519837335;
2896    daa[10*20+15] = 0.8953754669269811;  daa[10*20+16] = 1.0192004372506540;
2897    daa[10*20+17] = 1.5183114434679339;  daa[10*20+18] = 1.2452243224541324;
2898    daa[10*20+19] = 3.5458093276667237;
2899   
2900    daa[11*20+12] = 1.7129670976916258;  daa[11*20+13] = 0.3910559903834828;
2901    daa[11*20+14] = 1.3918935593582853;  daa[11*20+15] = 1.9776630140912268;
2902    daa[11*20+16] = 2.5513781312660280;  daa[11*20+17] = 0.5851920879490173;
2903    daa[11*20+18] = 0.7835447533710449;  daa[11*20+19] = 0.7801080335991272;
2904   
2905    daa[12*20+13] = 2.3201373546296349;  daa[12*20+14] = 0.4953193808676289;
2906    daa[12*20+15] = 1.0657482318076852;  daa[12*20+16] = 3.3628488360462363;
2907    daa[12*20+17] = 1.4680478689711018;  daa[12*20+18] = 1.0899165770956820;
2908    daa[12*20+19] = 4.0584577156753401;
2909   
2910    daa[13*20+14] = 0.3746821107962129;  daa[13*20+15] = 1.1079144700606407;
2911    daa[13*20+16] = 0.6882725908872254;  daa[13*20+17] = 3.3448437239772266;
2912    daa[13*20+18] = 10.3848523331334590;  daa[13*20+19] = 1.7039730522675411;
2913   
2914    daa[14*20+15] = 3.5465914843628927;  daa[14*20+16] = 1.9485376673137556;
2915    daa[14*20+17] = 0.4326058001438786;  daa[14*20+18] = 0.4819109019647465;
2916    daa[14*20+19] = 0.5985498912985666;
2917   
2918    daa[15*20+16] = 8.8479984061248178;  daa[15*20+17] = 0.6791126595939816;
2919    daa[15*20+18] = 0.9547229305958682;  daa[15*20+19] = 0.9305232113028208;
2920   
2921    daa[16*20+17] = 0.4514203099376473;  daa[16*20+18] = 0.8564314184691215;
2922    daa[16*20+19] = 3.4242218450865543;
2923   
2924    daa[17*20+18] = 4.5377235790405388;  daa[17*20+19] = 0.5658969249032649;
2925
2926    daa[18*20+19] = 1.0000000000000000;
2927
2928    for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[i*naa+j] = daa[j*naa+i];
2929
2930    pi[ 0]=0.0770764620135024 ; pi[ 1]=0.0500819370772208 ;
2931    pi[ 2]=0.0462377395993731 ; pi[ 3]=0.0537929860758246 ;
2932    pi[ 4]=0.0144533387583345 ; pi[ 5]=0.0408923608974345 ;
2933    pi[ 6]=0.0633579339160905 ; pi[ 7]=0.0655672355884439 ;
2934    pi[ 8]=0.0218802687005936 ; pi[ 9]=0.0591969699027449 ;
2935    pi[10]=0.0976461276528445 ; pi[11]=0.0592079410822730 ;
2936    pi[12]=0.0220695876653368 ; pi[13]=0.0413508521834260 ;
2937    pi[14]=0.0476871596856874 ; pi[15]=0.0707295165111524 ;
2938    pi[16]=0.0567759161524817 ; pi[17]=0.0127019797647213 ;
2939    pi[18]=0.0323746050281867 ; pi[19]=0.0669190817443274 ;
2940
2941    return 1;
2942}
2943
2944//////////////////////////////////////////////////////////////
2945//////////////////////////////////////////////////////////////
2946
2947
2948int Init_Qmat_Blosum62(phydbl *daa, phydbl *pi)
2949{
2950
2951    /*
2952       This model has been 'translated' from John Huelsenbeck and Fredrik Ronquist
2953       MrBayes program into PHYML format by Federico Abascal. Many thanks to them.
2954    */
2955
2956    /*
2957      Henikoff, S., and J. G. Henikoff. 1992. Amino acid substitution           
2958      matrices from protein blocks. Proc. Natl. Acad. Sci., U.S.A.           
2959      89:10915-10919.                                                         
2960    */
2961
2962    int i,j,naa;
2963    naa = 20;
2964
2965    daa[1*20+0]= 0.735790389698;  daa[2*20+0]= 0.485391055466;  daa[2*20+1]= 1.297446705134;  daa[3*20+0]= 0.543161820899; 
2966    daa[3*20+1]= 0.500964408555;  daa[3*20+2]= 3.180100048216;  daa[4*20+0]= 1.45999531047;   daa[4*20+1]= 0.227826574209; 
2967    daa[4*20+2]= 0.397358949897;  daa[4*20+3]= 0.240836614802;  daa[5*20+0]= 1.199705704602;  daa[5*20+1]= 3.020833610064; 
2968    daa[5*20+2]= 1.839216146992;  daa[5*20+3]= 1.190945703396;  daa[5*20+4]= 0.32980150463;   daa[6*20+0]= 1.1709490428;   
2969    daa[6*20+1]= 1.36057419042;   daa[6*20+2]= 1.24048850864;   daa[6*20+3]= 3.761625208368;  daa[6*20+4]= 0.140748891814; 
2970    daa[6*20+5]= 5.528919177928;  daa[7*20+0]= 1.95588357496;   daa[7*20+1]= 0.418763308518;  daa[7*20+2]= 1.355872344485; 
2971    daa[7*20+3]= 0.798473248968;  daa[7*20+4]= 0.418203192284;  daa[7*20+5]= 0.609846305383;  daa[7*20+6]= 0.423579992176; 
2972    daa[8*20+0]= 0.716241444998;  daa[8*20+1]= 1.456141166336;  daa[8*20+2]= 2.414501434208;  daa[8*20+3]= 0.778142664022; 
2973    daa[8*20+4]= 0.354058109831;  daa[8*20+5]= 2.43534113114;   daa[8*20+6]= 1.626891056982;  daa[8*20+7]= 0.539859124954; 
2974    daa[9*20+0]= 0.605899003687;  daa[9*20+1]= 0.232036445142;  daa[9*20+2]= 0.283017326278;  daa[9*20+3]= 0.418555732462; 
2975    daa[9*20+4]= 0.774894022794;  daa[9*20+5]= 0.236202451204;  daa[9*20+6]= 0.186848046932;  daa[9*20+7]= 0.189296292376; 
2976    daa[9*20+8]= 0.252718447885;  daa[10*20+0]= 0.800016530518; daa[10*20+1]= 0.622711669692; daa[10*20+2]= 0.211888159615; 
2977    daa[10*20+3]= 0.218131577594; daa[10*20+4]= 0.831842640142; daa[10*20+5]= 0.580737093181; daa[10*20+6]= 0.372625175087; 
2978    daa[10*20+7]= 0.217721159236; daa[10*20+8]= 0.348072209797; daa[10*20+9]= 3.890963773304; daa[11*20+0]= 1.295201266783; 
2979    daa[11*20+1]= 5.411115141489; daa[11*20+2]= 1.593137043457; daa[11*20+3]= 1.032447924952; daa[11*20+4]= 0.285078800906; 
2980    daa[11*20+5]= 3.945277674515; daa[11*20+6]= 2.802427151679; daa[11*20+7]= 0.752042440303; daa[11*20+8]= 1.022507035889; 
2981    daa[11*20+9]= 0.406193586642; daa[11*20+10]= 0.445570274261;daa[12*20+0]= 1.253758266664; daa[12*20+1]= 0.983692987457; 
2982    daa[12*20+2]= 0.648441278787; daa[12*20+3]= 0.222621897958; daa[12*20+4]= 0.76768882348;  daa[12*20+5]= 2.494896077113; 
2983    daa[12*20+6]= 0.55541539747;  daa[12*20+7]= 0.459436173579; daa[12*20+8]= 0.984311525359; daa[12*20+9]= 3.364797763104; 
2984    daa[12*20+10]= 6.030559379572;daa[12*20+11]= 1.073061184332;daa[13*20+0]= 0.492964679748; daa[13*20+1]= 0.371644693209; 
2985    daa[13*20+2]= 0.354861249223; daa[13*20+3]= 0.281730694207; daa[13*20+4]= 0.441337471187; daa[13*20+5]= 0.14435695975; 
2986    daa[13*20+6]= 0.291409084165; daa[13*20+7]= 0.368166464453; daa[13*20+8]= 0.714533703928; daa[13*20+9]= 1.517359325954; 
2987    daa[13*20+10]= 2.064839703237;daa[13*20+11]= 0.266924750511;daa[13*20+12]= 1.77385516883; daa[14*20+0]= 1.173275900924; 
2988    daa[14*20+1]= 0.448133661718; daa[14*20+2]= 0.494887043702; daa[14*20+3]= 0.730628272998; daa[14*20+4]= 0.356008498769; 
2989    daa[14*20+5]= 0.858570575674; daa[14*20+6]= 0.926563934846; daa[14*20+7]= 0.504086599527; daa[14*20+8]= 0.527007339151; 
2990    daa[14*20+9]= 0.388355409206; daa[14*20+10]= 0.374555687471;daa[14*20+11]= 1.047383450722;daa[14*20+12]= 0.454123625103;
2991    daa[14*20+13]= 0.233597909629;daa[15*20+0]= 4.325092687057; daa[15*20+1]= 1.12278310421;  daa[15*20+2]= 2.904101656456; 
2992    daa[15*20+3]= 1.582754142065; daa[15*20+4]= 1.197188415094; daa[15*20+5]= 1.934870924596; daa[15*20+6]= 1.769893238937; 
2993    daa[15*20+7]= 1.509326253224; daa[15*20+8]= 1.11702976291;  daa[15*20+9]= 0.35754441246;  daa[15*20+10]= 0.352969184527;
2994    daa[15*20+11]= 1.752165917819;daa[15*20+12]= 0.918723415746;daa[15*20+13]= 0.540027644824;daa[15*20+14]= 1.169129577716;
2995    daa[16*20+0]= 1.729178019485; daa[16*20+1]= 0.914665954563; daa[16*20+2]= 1.898173634533; daa[16*20+3]= 0.934187509431; 
2996    daa[16*20+4]= 1.119831358516; daa[16*20+5]= 1.277480294596; daa[16*20+6]= 1.071097236007; daa[16*20+7]= 0.641436011405; 
2997    daa[16*20+8]= 0.585407090225; daa[16*20+9]= 1.17909119726;  daa[16*20+10]= 0.915259857694;daa[16*20+11]= 1.303875200799;
2998    daa[16*20+12]= 1.488548053722;daa[16*20+13]= 0.488206118793;daa[16*20+14]= 1.005451683149;daa[16*20+15]= 5.15155629227; 
2999    daa[17*20+0]= 0.465839367725; daa[17*20+1]= 0.426382310122; daa[17*20+2]= 0.191482046247; daa[17*20+3]= 0.145345046279; 
3000    daa[17*20+4]= 0.527664418872; daa[17*20+5]= 0.758653808642; daa[17*20+6]= 0.407635648938; daa[17*20+7]= 0.508358924638; 
3001    daa[17*20+8]= 0.30124860078;  daa[17*20+9]= 0.34198578754;  daa[17*20+10]= 0.6914746346;  daa[17*20+11]= 0.332243040634;
3002    daa[17*20+12]= 0.888101098152;daa[17*20+13]= 2.074324893497;daa[17*20+14]= 0.252214830027;daa[17*20+15]= 0.387925622098;
3003    daa[17*20+16]= 0.513128126891;daa[18*20+0]= 0.718206697586; daa[18*20+1]= 0.720517441216; daa[18*20+2]= 0.538222519037; 
3004    daa[18*20+3]= 0.261422208965; daa[18*20+4]= 0.470237733696; daa[18*20+5]= 0.95898974285;  daa[18*20+6]= 0.596719300346; 
3005    daa[18*20+7]= 0.308055737035; daa[18*20+8]= 4.218953969389; daa[18*20+9]= 0.674617093228; daa[18*20+10]= 0.811245856323;
3006    daa[18*20+11]= 0.7179934869;  daa[18*20+12]= 0.951682162246;daa[18*20+13]= 6.747260430801;daa[18*20+14]= 0.369405319355;
3007    daa[18*20+15]= 0.796751520761;daa[18*20+16]= 0.801010243199;daa[18*20+17]= 4.054419006558;daa[19*20+0]= 2.187774522005; 
3008    daa[19*20+1]= 0.438388343772; daa[19*20+2]= 0.312858797993; daa[19*20+3]= 0.258129289418; daa[19*20+4]= 1.116352478606; 
3009    daa[19*20+5]= 0.530785790125; daa[19*20+6]= 0.524253846338; daa[19*20+7]= 0.25334079019;  daa[19*20+8]= 0.20155597175; 
3010    daa[19*20+9]= 8.311839405458; daa[19*20+10]= 2.231405688913;daa[19*20+11]= 0.498138475304;daa[19*20+12]= 2.575850755315;
3011    daa[19*20+13]= 0.838119610178;daa[19*20+14]= 0.496908410676;daa[19*20+15]= 0.561925457442;daa[19*20+16]= 2.253074051176;
3012    daa[19*20+17]= 0.266508731426;daa[19*20+18]= 1;             
3013
3014    for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j];
3015
3016    pi[0]= 0.074;                 pi[1]= 0.052;                 pi[2]= 0.045;                 pi[3]= 0.054;                 
3017    pi[4]= 0.025;                 pi[5]= 0.034;                 pi[6]= 0.054;                 pi[7]= 0.074;                 
3018    pi[8]= 0.026;                 pi[9]= 0.068;                 pi[10]= 0.099;                pi[11]= 0.058;               
3019    pi[12]= 0.025;                pi[13]= 0.047;                pi[14]= 0.039;                pi[15]= 0.057;               
3020    pi[16]= 0.051;                pi[17]= 0.013;                pi[18]= 0.032;                pi[19]= 0.073;               
3021
3022    return 1;
3023  }
3024
3025//////////////////////////////////////////////////////////////
3026//////////////////////////////////////////////////////////////
3027
3028
3029int Init_Qmat_MtMam(phydbl *daa, phydbl *pi)
3030{
3031    /*
3032       This model has been 'translated' from Ziheng Yang's PAML program
3033       into PHYML format by Federico Abascal. Many thanks to them.
3034    */
3035
3036    /*
3037      Cao, Y. et al. 1998 Conflict amongst individual mitochondrial
3038      proteins in resolving the phyLOGeny of eutherian orders. Journal
3039      of Molecular Evolution 15:1600-1611.
3040    */
3041
3042    int i,j,naa;
3043    naa = 20;
3044   
3045    daa[1*20+0]= 32;              daa[2*20+0]= 2;    daa[2*20+1]= 4;               daa[3*20+0]= 11;
3046    daa[3*20+1]= 0;               daa[3*20+2]= 864;  daa[4*20+0]= 0;               daa[4*20+1]= 186;
3047    daa[4*20+2]= 0;               daa[4*20+3]= 0;    daa[5*20+0]= 0;               daa[5*20+1]= 246;
3048    daa[5*20+2]= 8;               daa[5*20+3]= 49;   daa[5*20+4]= 0;               daa[6*20+0]= 0;
3049    daa[6*20+1]= 0;               daa[6*20+2]= 0;    daa[6*20+3]= 569;             daa[6*20+4]= 0;
3050    daa[6*20+5]= 274;             daa[7*20+0]= 78;   daa[7*20+1]= 18;              daa[7*20+2]= 47;
3051    daa[7*20+3]= 79;              daa[7*20+4]= 0;    daa[7*20+5]= 0;               daa[7*20+6]= 22;
3052    daa[8*20+0]= 8;               daa[8*20+1]= 232;  daa[8*20+2]= 458;             daa[8*20+3]= 11;
3053    daa[8*20+4]= 305;             daa[8*20+5]= 550;  daa[8*20+6]= 22;              daa[8*20+7]= 0;
3054    daa[9*20+0]= 75;              daa[9*20+1]= 0;    daa[9*20+2]= 19;              daa[9*20+3]= 0;
3055    daa[9*20+4]= 41;              daa[9*20+5]= 0;    daa[9*20+6]= 0;               daa[9*20+7]= 0;
3056    daa[9*20+8]= 0;               daa[10*20+0]= 21;  daa[10*20+1]= 6;              daa[10*20+2]= 0;
3057    daa[10*20+3]= 0;              daa[10*20+4]= 27;  daa[10*20+5]= 20;             daa[10*20+6]= 0;
3058    daa[10*20+7]= 0;              daa[10*20+8]= 26;  daa[10*20+9]= 232;            daa[11*20+0]= 0;
3059    daa[11*20+1]= 50;             daa[11*20+2]= 408; daa[11*20+3]= 0;              daa[11*20+4]= 0;
3060    daa[11*20+5]= 242;            daa[11*20+6]= 215; daa[11*20+7]= 0;              daa[11*20+8]= 0;
3061    daa[11*20+9]= 6;              daa[11*20+10]= 4;  daa[12*20+0]= 76;             daa[12*20+1]= 0;
3062    daa[12*20+2]= 21;             daa[12*20+3]= 0;   daa[12*20+4]= 0;              daa[12*20+5]= 22;
3063    daa[12*20+6]= 0;              daa[12*20+7]= 0;   daa[12*20+8]= 0;              daa[12*20+9]= 378;
3064    daa[12*20+10]= 609;           daa[12*20+11]= 59; daa[13*20+0]= 0;              daa[13*20+1]= 0;
3065    daa[13*20+2]= 6;              daa[13*20+3]= 5;   daa[13*20+4]= 7;              daa[13*20+5]= 0;
3066    daa[13*20+6]= 0;              daa[13*20+7]= 0;   daa[13*20+8]= 0;              daa[13*20+9]= 57;
3067    daa[13*20+10]= 246;           daa[13*20+11]= 0;  daa[13*20+12]= 11;            daa[14*20+0]= 53;
3068    daa[14*20+1]= 9;              daa[14*20+2]= 33;  daa[14*20+3]= 2;              daa[14*20+4]= 0;
3069    daa[14*20+5]= 51;             daa[14*20+6]= 0;   daa[14*20+7]= 0;              daa[14*20+8]= 53;
3070    daa[14*20+9]= 5;              daa[14*20+10]= 43; daa[14*20+11]= 18;            daa[14*20+12]= 0;
3071    daa[14*20+13]= 17;            daa[15*20+0]= 342; daa[15*20+1]= 3;              daa[15*20+2]= 446;
3072    daa[15*20+3]= 16;             daa[15*20+4]= 347; daa[15*20+5]= 30;             daa[15*20+6]= 21;
3073    daa[15*20+7]= 112;            daa[15*20+8]= 20;  daa[15*20+9]= 0;              daa[15*20+10]= 74;
3074    daa[15*20+11]= 65;            daa[15*20+12]= 47; daa[15*20+13]= 90;            daa[15*20+14]= 202;
3075    daa[16*20+0]= 681;            daa[16*20+1]= 0;   daa[16*20+2]= 110;            daa[16*20+3]= 0;
3076    daa[16*20+4]= 114;            daa[16*20+5]= 0;   daa[16*20+6]= 4;              daa[16*20+7]= 0;
3077    daa[16*20+8]= 1;              daa[16*20+9]= 360; daa[16*20+10]= 34;            daa[16*20+11]= 50;
3078    daa[16*20+12]= 691;           daa[16*20+13]= 8;  daa[16*20+14]= 78;            daa[16*20+15]= 614;
3079    daa[17*20+0]= 5;              daa[17*20+1]= 16;  daa[17*20+2]= 6;              daa[17*20+3]= 0;
3080    daa[17*20+4]= 65;             daa[17*20+5]= 0;   daa[17*20+6]= 0;              daa[17*20+7]= 0;
3081    daa[17*20+8]= 0;              daa[17*20+9]= 0;   daa[17*20+10]= 12;            daa[17*20+11]= 0;
3082    daa[17*20+12]= 13;            daa[17*20+13]= 0;  daa[17*20+14]= 7;             daa[17*20+15]= 17;
3083    daa[17*20+16]= 0;             daa[18*20+0]= 0;   daa[18*20+1]= 0;              daa[18*20+2]= 156;
3084    daa[18*20+3]= 0;              daa[18*20+4]= 530; daa[18*20+5]= 54;             daa[18*20+6]= 0;
3085    daa[18*20+7]= 1;              daa[18*20+8]= 1525;daa[18*20+9]= 16;             daa[18*20+10]= 25;
3086    daa[18*20+11]= 67;            daa[18*20+12]= 0;  daa[18*20+13]= 682;           daa[18*20+14]= 8;
3087    daa[18*20+15]= 107;           daa[18*20+16]= 0;  daa[18*20+17]= 14;            daa[19*20+0]= 398;
3088    daa[19*20+1]= 0;              daa[19*20+2]= 0;   daa[19*20+3]= 10;             daa[19*20+4]= 0;
3089    daa[19*20+5]= 33;             daa[19*20+6]= 20;  daa[19*20+7]= 5;              daa[19*20+8]= 0;
3090    daa[19*20+9]= 2220;           daa[19*20+10]= 100;daa[19*20+11]= 0;             daa[19*20+12]= 832;
3091    daa[19*20+13]= 6;             daa[19*20+14]= 0;  daa[19*20+15]= 0;             daa[19*20+16]= 237;
3092    daa[19*20+17]= 0;             daa[19*20+18]= 0;
3093   
3094    for (i=0; i<naa; i++) for (j=0; j<i; j++) daa[j*naa+i] = daa[i*naa+j];
3095   
3096    pi[0]= 0.0692;  pi[1]=  0.0184;  pi[2]= 0.04;    pi[3]= 0.0186;
3097    pi[4]= 0.0065;  pi[5]=  0.0238;  pi[6]= 0.0236;  pi[7]= 0.0557;
3098    pi[8]= 0.0277;  pi[9]=  0.0905;  pi[10]=0.1675;  pi[11]= 0.0221;
3099    pi[12]=0.0561;  pi[13]= 0.0611;  pi[14]=0.0536;  pi[15]= 0.0725;
3100    pi[16]=0.087;   pi[17]= 0.0293;  pi[18]=0.034;   pi[19]= 0.0428;
3101   
3102    return 1;
3103}
3104                         
3105//////////////////////////////////////////////////////////////
3106//////////////////////////////////////////////////////////////
3107
3108
3109void M4_Init_Model(m4 *m4mod, calign *data, t_mod *mod)
3110{
3111  int i,j,ct;
3112  phydbl fq;
3113
3114  if(mod->io->datatype == NT)      m4mod->n_o = 4;
3115  else if(mod->io->datatype == AA) m4mod->n_o = 20;
3116  else 
3117    {
3118      PhyML_Printf("\n== Not implemented yet.");
3119      PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
3120      Warn_And_Exit("");
3121    }
3122
3123  mod->ns = m4mod->n_o * m4mod->n_h;
3124
3125  For(i,m4mod->n_o) m4mod->o_fq[i] = mod->e_frq->pi->v[i]; /*! At that stage, the mod->pi vector as been initialized
3126                                                             under a standard non covarion type of model. Use these
3127                                                             frequencies as they have been set according to the
3128                                                             nucleotide substitution model chosen (e.g., 1/4 for JC69). !*/
3129  For(i,(int)(m4mod->n_h)) m4mod->multipl[i] = 1.;
3130
3131  ct = 0;
3132  For(i,m4mod->n_o-1)
3133    {
3134      for(j=i+1;j<m4mod->n_o;j++)
3135        {
3136          m4mod->o_rr[ct] = MAX(mod->r_mat->qmat->v[i*m4mod->n_o+j],1.E-5);       
3137          ct++;
3138        }
3139    }
3140 
3141  For(i,(int)(m4mod->n_h*(m4mod->n_h-1)/2)) m4mod->h_rr[i] = 1.;
3142  fq = (phydbl)(1./m4mod->n_h);
3143 
3144
3145  if(mod->s_opt->opt_cov_delta) m4mod->delta = 1.0;
3146  if(mod->s_opt->opt_cov_alpha) m4mod->alpha = 1.0;
3147  For(i,m4mod->n_h) m4mod->h_fq[i] = fq;
3148  For(i,m4mod->n_h) m4mod->h_fq_unscaled[i] = 1.0;
3149  For(i,m4mod->n_h) m4mod->multipl[i] = (phydbl)i;
3150  For(i,m4mod->n_h) m4mod->multipl_unscaled[i] = (phydbl)i;
3151
3152  Switch_Eigen(YES,mod);
3153  M4_Update_Qmat(m4mod,mod);
3154}
3155
3156//////////////////////////////////////////////////////////////
3157//////////////////////////////////////////////////////////////
3158
3159
3160//////////////////////////////////////////////////////////////
3161//////////////////////////////////////////////////////////////
3162
3163//////////////////////////////////////////////////////////////
3164//////////////////////////////////////////////////////////////
3165
3166//////////////////////////////////////////////////////////////
3167//////////////////////////////////////////////////////////////
3168
3169//////////////////////////////////////////////////////////////
3170//////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.