GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it
#include"typedef.h"1#include"idem.h"2#include"longtools.h"3#include"getput.h"4#include"voronoi.h"5#include"bravais.h"6#include"datei.h"7#include"matrix.h"8#include"orbit.h"9#include"tools.h"10#include"reduction.h"1112extern int INFO_LEVEL;13/*****************************************************************************14@15@-----------------------------------------------------------------------------16@ FILE: bravais_catalog.c17@-----------------------------------------------------------------------------18@19*****************************************************************************/2021symbol_out *read_symbol_from_string(symb)22char *symb;23{24char string[80],25slash, *str ;26char f[80];27char *fn;28char *dat;29int i, j, k, l, m, n, p, q, x, groesser;30int len;31int no;32int breite;33char merk[15];34char merk1[10];35int index;36char konst[MAXDIM][80];37int konst_dim;38int komp[MAXDIM];39int artgleich[MAXDIM];40int dim = 0, konstit = 0;41int zerleg[MAXDIM][5];42bravais_TYP **grps;43symbol_out *erg;44matrix_TYP *In;4546fn = (char *) malloc(1024 * sizeof(char));4748for (i=0;i<80;i++) string[i] = 0;4950for( i=0; i<MAXDIM; i++)51komp[i] = 0;52for(i=0; i<MAXDIM; i++)53for(j=0; j<5;j++)54zerleg[i][j] = 0;5556/* inserted tilman 05/08/97 */57sprintf(string,"%s",symb);58right_order(string);5960str = string;61len = strlen(str);62while( len > 0 )63{64i = strcspn(str, ";,");65if(i>0)66{67sscanf(str, "%d", &konst_dim);68dim += konst_dim;69if(dim > MAXDIM)70{71printf("dimension of a.d. bravais-group is to big\n");72exit(3);73}74zerleg[konstit][0] = konst_dim;75itoa(konst_dim, konst[konstit]);76k = strcspn(str, ";");77if(k == i)78zerleg[konstit][3] = 1;79j = strcspn(str, "-");80l = strcspn(str, "'");81if(j<i)82{83strcat(konst[konstit], "-");84str = str+j+1;85len = len-j-1;86sscanf(str, "%d", &index);87zerleg[konstit][1] = index;88itoa(index, merk);89strcat(konst[konstit], merk);90memset(merk,'\0',strlen(merk));91if(l<i)92{93strcat(konst[konstit], "'");94zerleg[konstit][2] = 1;95}96i = i-j-1;97}98konstit++;99}100str = str+i+1;101len = len-i-1;102}103104erg = (symbol_out *) malloc(sizeof(symbol_out));105erg->grp = (bravais_TYP *) calloc(1, sizeof(bravais_TYP));106erg->grp->dim = dim;107108/*--------------------------------------------------------------------*\109| swap the atoms into the right order |110\*--------------------------------------------------------------------*/111for(i=0; i<MAXDIM-1; i++)112{113if(zerleg[i][0] != zerleg[i+1][0] || zerleg[i][1] != zerleg[i+1][1] || zerleg[i][2] != zerleg[i+1][2])114zerleg[i][3] = 1;115}116if(zerleg[MAXDIM-1][0] != 0)117zerleg[MAXDIM-1][3] = 1;118i=0;119while(i<MAXDIM)120{121while( i< MAXDIM && zerleg[i][0] == 0)122i++;123k = 1;124while(i<MAXDIM && zerleg[i][3] == 0 && zerleg[i][0] != 0)125{ i++; k++;}126if(i != MAXDIM)127zerleg[i][4] = k;128i++;129}130for(i=0; i<MAXDIM; i++)131{132for(j=i+1; j<MAXDIM; j++)133{134groesser = 1;135if(zerleg[i][4] == 0 && zerleg[j][4] != 0)136groesser = 0;137for(k=0; k<3 && zerleg[j][k] == zerleg[i][k]; k++);138if(k == 0 && zerleg[i][k] < zerleg[j][k])139groesser = 0;140if(k == 1 && zerleg[i][1] > zerleg[j][1])141groesser = 0;142if(k == 2 && zerleg[i][2] < zerleg[j][2])143groesser = 0;144if(k == 3 && zerleg[i][4] < zerleg[j][4])145groesser = 0;146if(groesser == 0)147{148for(k=0; k<5; k++)149{150x = zerleg[i][k]; zerleg[i][k] = zerleg[j][k]; zerleg[j][k] = x;151}152}153}154}155/************************156for(i=0; i<konstit; i++)157{158for(j=0; j<5; j++)159printf("%d ", zerleg[i][j]);160printf("\n");161}162*************************/163164konstit = 0;165for(i=0; i<MAXDIM; i++)166{167if(zerleg[i][4] != 0)168{169itoa(zerleg[i][0], konst[konstit]);170if(zerleg[i][1] != 0)171{172strcat(konst[konstit], "-");173itoa(zerleg[i][1], merk);174strcat(konst[konstit], merk);175memset(merk,'\0',strlen(merk));176if(zerleg[i][2] != 0)177strcat(konst[konstit], "\'");178}179konstit++;180}181}182183/*--------------------------------------------------------------------*\184| read the atoms |185\*--------------------------------------------------------------------*/186grps = (bravais_TYP **) malloc(konstit *sizeof(bravais_TYP *));187dat = ATOMS;188/**************189dat = TOPDIR "/lib/atoms/";190f = (char **) malloc(konstit *sizeof(char *));191***************/192for(i=0; i<konstit; i++)193{194strcpy(f, dat);195itoa(zerleg[i][0], merk);196strcat(f, merk);197if(zerleg[i][1] != 0)198{199strcat(f, "-");200itoa(zerleg[i][1], merk);201strcat(f, merk);202memset(merk,'\0',strlen(merk));203}204if(zerleg[i][2] != 0)205{206strcat(f, "\'");207}208if (INFO_LEVEL & 4) fprintf(stderr,"file-name: %s\n",f);209grps[i] = get_bravais(f);210for(j=0; j<grps[i]->form_no; j++)211Check_mat(grps[i]->form[j]);212}213214/*--------------------------------------------------------------------*\215| calculate the generators of erg->grp |216\*--------------------------------------------------------------------*/217erg->grp->gen_no = 0;218erg->grp->form_no = 0;219erg->grp->zentr_no = 0;220erg->grp->normal_no = 0;221for(i=0; i<konstit; i++)222erg->grp->gen_no += grps[i]->gen_no;223erg->grp->gen = (matrix_TYP **)malloc(erg->grp->gen_no *sizeof(matrix_TYP));224for(i=0; i<erg->grp->gen_no; i++)225erg->grp->gen[i] = init_mat(erg->grp->dim, erg->grp->dim, "");226j=0;227p=0;228for(i=0; i<konstit; i++)229{230for(k=0; k<grps[i]->gen_no; k++)231{232for(l=0; l<zerleg[i][4] * grps[i]->dim; l += grps[i]->dim)233{234for(m=0; m<grps[i]->dim; m++)235for(n=0; n<grps[i]->dim; n++)236erg->grp->gen[p+k]->array.SZ[j+l+m][j+l+n]=grps[i]->gen[k]->array.SZ[m][n];237}238for(m=0; m<j; m++)239erg->grp->gen[p+k]->array.SZ[m][m] = 1;240for(m=( j + zerleg[i][4] * grps[i]->dim); m<erg->grp->dim; m++)241erg->grp->gen[p+k]->array.SZ[m][m] = 1;242}243j = j + zerleg[i][4] * grps[i]->dim;244p = p + grps[i]->gen_no;245}246247/*--------------------------------------------------------------------*\248| calculate the invariant forms of erg->grp |249\*--------------------------------------------------------------------*/250for(i=0; i<konstit; i++)251{252for(j=0; j<grps[i]->form_no; j++)253{254if(grps[i]->form[j]->flags.Symmetric == TRUE)255erg->grp->form_no += (zerleg[i][4] * (zerleg[i][4] +1)/2);256if(grps[i]->form[j]->flags.Symmetric == FALSE)257erg->grp->form_no += (zerleg[i][4] * (zerleg[i][4] -1)/2);258}259}260erg->grp->form = (matrix_TYP **)malloc(erg->grp->form_no *sizeof(matrix_TYP *));261for(i=0; i<erg->grp->form_no; i++)262erg->grp->form[i] = init_mat(erg->grp->dim, erg->grp->dim, "");263j=0;264q=0;265for(i=0; i<konstit; i++)266{267for(k=0; k<grps[i]->form_no; k++)268{269for(l=j; l< j+zerleg[i][4] * grps[i]->dim; l += grps[i]->dim)270{271if(grps[i]->form[k]->flags.Symmetric == TRUE)272{273for(n=0; n<grps[i]->dim; n++)274{275for(p=0; p<grps[i]->dim; p++)276{277erg->grp->form[q]->array.SZ[l+n][l+p] = grps[i]->form[k]->array.SZ[n][p];278}279}280q++;281}282for(m=(l+grps[i]->dim); m<(j+zerleg[i][4] * grps[i]->dim); m +=grps[i]->dim)283{284for(n=0; n<grps[i]->dim; n++)285{286for(p=0; p<grps[i]->dim; p++)287{288erg->grp->form[q]->array.SZ[l+n][m+p] = grps[i]->form[k]->array.SZ[n][p];289erg->grp->form[q]->array.SZ[m+p][l+n] = grps[i]->form[k]->array.SZ[n][p];290}291}292q++;293}294}295}296j = j + zerleg[i][4] * grps[i]->dim;297}298/* assure that the full integral lattice is passed */299long_rein_formspace(erg->grp->form,erg->grp->form_no,1);300301/*--------------------------------------------------------------------*\302| calculate the centralizer of the group |303\*--------------------------------------------------------------------*/304for(i=0; i<konstit; i++)305{306erg->grp->cen_no = erg->grp->cen_no + grps[i]->cen_no;307if(zerleg[i][4] >= 2)308erg->grp->cen_no = erg->grp->cen_no + grps[i]->zentr_no + 1;309if(zerleg[i][4] > 2)310erg->grp->cen_no ++;311}312erg->grp->cen = (matrix_TYP **)malloc(erg->grp->cen_no *sizeof(matrix_TYP));313for(i=0; i<erg->grp->cen_no; i++)314erg->grp->cen[i] = init_mat(erg->grp->dim, erg->grp->dim, "");315j=0;316no = 0;317for(i=0; i<konstit; i++)318{319for(k=0; k<grps[i]->cen_no; k++)320{321for(l=0; l<j; l++)322erg->grp->cen[no]->array.SZ[l][l] = 1;323for(l=j+grps[i]->dim; l<erg->grp->dim; l++)324erg->grp->cen[no]->array.SZ[l][l] = 1;325for(l=0; l<grps[i]->dim; l++)326for(m=0; m<grps[i]->dim; m++)327erg->grp->cen[no]->array.SZ[j+l][j+m] = grps[i]->cen[k]->array.SZ[l][m];328no++;329}330if(zerleg[i][4] >=2)331{332for(k=0; k<grps[i]->zentr_no; k++)333{334for(l=0; l<j; l++)335erg->grp->cen[no]->array.SZ[l][l] = 1;336for(l=j; l<j+grps[i]->dim; l++)337erg->grp->cen[no]->array.SZ[l][l] = -1;338for(l=j+grps[i]->dim; l<erg->grp->dim; l++)339erg->grp->cen[no]->array.SZ[l][l] = 1;340for(l=0; l<grps[i]->dim; l++)341for(m=0; m<grps[i]->dim; m++)342erg->grp->cen[no]->array.SZ[j+l+grps[i]->dim][j+m] = grps[i]->zentr[k]->array.SZ[l][m];343no++;344}345for(l=0; l<j; l++)346erg->grp->cen[no]->array.SZ[l][l] = 1;347for(l=j+(2 * grps[i]->dim); l<erg->grp->dim; l++)348erg->grp->cen[no]->array.SZ[l][l] = 1;349for(l=0; l<grps[i]->dim; l++)350{351erg->grp->cen[no]->array.SZ[j+l+grps[i]->dim][j+l] = 1;352erg->grp->cen[no]->array.SZ[j+l][j+l+grps[i]->dim] = 1;353}354no++;355}356if(zerleg[i][4] > 2)357{358for(l=0; l<j; l++)359erg->grp->cen[no]->array.SZ[l][l] = 1;360for(l=j+(zerleg[i][4] * grps[i]->dim); l<erg->grp->dim; l++)361erg->grp->cen[no]->array.SZ[l][l] = 1;362for(l=0; l< ((zerleg[i][4] -1) * grps[i]->dim); l++)363erg->grp->cen[no]->array.SZ[l+j+grps[i]->dim][l+j] = 1;364for(l=0; l<grps[i]->dim; l++)365erg->grp->cen[no]->array.SZ[j+l][l+j+((zerleg[i][4]-1) * grps[i]->dim)] = 1;366no++;367}368369j = j + zerleg[i][4] * grps[i]->dim;370}371372/*--------------------------------------------------------------------*\373| calculate the normalizer of the group |374\*--------------------------------------------------------------------*/375376for(i=0; i<konstit; i++)377{378artgleich[i] = 0;379erg->grp->normal_no += grps[i]->normal_no;380if(i+1 < konstit)381{382if(zerleg[i][0] == zerleg[i+1][0] && zerleg[i][1] == zerleg[i+1][1] &&383zerleg[i][2] == zerleg[i+1][2] && zerleg[i][3] == zerleg[i+1][3] &&384zerleg[i][4] == zerleg[i+1][4])385{386artgleich[i] = TRUE;387erg->grp->normal_no++;388}389}390}391392if (erg->grp->normal_no > 0)393erg->grp->normal = (matrix_TYP **) malloc(erg->grp->normal_no394* sizeof(matrix_TYP));395else396erg->grp->normal = NULL;397398for(i=0; i<erg->grp->normal_no; i++)399erg->grp->normal[i] = init_mat(erg->grp->dim, erg->grp->dim, "");400j=0;401no = 0;402for(i=0; i<konstit; i++)403{404for(k=0; k<grps[i]->normal_no; k++)405{406for(l=0; l<j;l++)407erg->grp->normal[no]->array.SZ[l][l] = 1;408for(l=(j+zerleg[i][4] * grps[i]->dim); l<erg->grp->dim; l++)409erg->grp->normal[no]->array.SZ[l][l] = 1;410for(l=j; l< (j+zerleg[i][4] * grps[i]->dim); l += grps[i]->dim)411{412for(m=0; m<grps[i]->dim; m++)413for(p=0; p<grps[i]->dim; p++)414erg->grp->normal[no]->array.SZ[m+l][p+l] = grps[i]->normal[k]->array.SZ[m][p];415}416no++;417}418if(artgleich[i] == TRUE)419{420breite = grps[i]->dim * zerleg[i][4];421for(l=0; l<j; l++)422erg->grp->normal[no]->array.SZ[l][l] = 1;423for(l=(j + 2 * breite); l<erg->grp->dim; l++)424erg->grp->normal[no]->array.SZ[l][l] = 1;425for(l=j; l< (j+breite); l++)426{427erg->grp->normal[no]->array.SZ[l][l+breite] = 1;428erg->grp->normal[no]->array.SZ[l+breite][l] = 1;429}430no++;431}432j = j + zerleg[i][4] * grps[i]->dim;433}434435/*--------------------------------------------------------------------*\436| if no additional normalizer or centraliser nessecarry, put identity |437\*--------------------------------------------------------------------*/438if(erg->grp->normal_no == 0 || erg->grp->cen_no == 0)439{440if(erg->grp->normal_no == 0)441{442erg->grp->normal_no = 1;443erg->grp->normal = (matrix_TYP **) malloc(1 *sizeof(matrix_TYP *));444erg->grp->normal[0] = einheitsmatrix(erg->grp->dim);445}446if(erg->grp->cen_no == 0)447{448erg->grp->cen_no = 1;449erg->grp->cen = (matrix_TYP **) malloc(1 *sizeof(matrix_TYP *));450erg->grp->cen[0] = einheitsmatrix(erg->grp->dim);451}452}453454/*--------------------------------------------------------------------*\455| calculate the order of erg->grp |456\*--------------------------------------------------------------------*/457for(j=0; j<100; j++)458erg->grp->divisors[j] = 0;459erg->grp->order = 1;460for(i=0; i<konstit; i++)461{462for(j=0; j<100; j++)463erg->grp->divisors[j] += grps[i]->divisors[j];464erg->grp->order *= grps[i]->order;465}466467/*--------------------------------------------------------------------*\468| find file where erg->grp->zentr are stored |469\*--------------------------------------------------------------------*/470strcpy(fn, TABLEDIM);471/*********************************472strcpy(fn, TOPDIR "/lib/dim");473*********************************/474itoa(erg->grp->dim, merk);475strcat(fn, merk);476strcat(fn, "/");477for(i=0; i<konstit; i++)478{479itoa(zerleg[i][0], merk);480if(zerleg[i][1] != 0)481{482strcat(merk, "-");483itoa(zerleg[i][1], merk1);484strcat(merk, merk1);485}486if(zerleg[i][2] != 0)487strcat(merk, "\'");488for(j=0; j<zerleg[i][4]; j++)489{490strcat(fn, merk);491if(j != (zerleg[i][4] -1))492strcat(fn, ",");493else494{495if(i!= (konstit -1))496strcat(fn, ";");497}498}499}500erg->fn = fn;501502if (INFO_LEVEL & 4) printf("%s\n", erg->fn);503504/* clear the atoms */505for (i=0;i<konstit;i++) free_bravais(grps[i]);506free(grps);507508return(erg);509}510511/****************************************************************************512@513@------------------------------------------------------------------------------514@515@ bravais_TYP *catalog_number(bravais_TYP *G,char *symb,516@ matrix_TYP **TR,int *almost,int *zclass)517@518@ The function searches for a Bravais group Z-equivalent to G in the519@ catalog. It will return this group, a transformation matrix via TR[0],520@ and the coordinates of the group in the catalog via almost[0], zclass[0].521@522@ It will return a transfromation matrix via TR[0], and the position523@ in the catalog via almost[0] and zclass[0].524@525@ bravais_TYP *G : The group in question. Its order must be given.526@ char *symb : The symbol of the group. It can be calculated via symbol(..)527@ matrix_TYP **TR: pointer for the transformation matrix which transforms528@ the given group G to the group returned via konj_bravais,529@ ie. TR[0] * G * TR[0]^-1 = group returned.530@ int *almost : the position of the almost decomposable group in the531@ catalog is returned via this pointer.532@ int *zclass : 2 coordinate of the group in the catalog.533@534@------------------------------------------------------------------------------535@536*****************************************************************************/537bravais_TYP *catalog_number(bravais_TYP *G,char *symb,538matrix_TYP **TR,int *almost,int *zclass)539{540541bravais_TYP *T = NULL,542*Gneu, /* the bravais group G written in a better basis */543*Gtr, /* transposed of Gneu */544*H, /* H is the pointer where a group read from the545catalog is stored */546*Htr; /* transposed of H */547548symbol_out *S;549550matrix_TYP *X, /* strictly speeking not needed, just to551avoid writting TR[0] all the time */552*F, /* a positive definite form for reducing G */553*B, /* a new, good basis for G */554*BI; /* the inverse of B */555556char *file;557558int i,559anz_gperfect=0;560561voronoi_TYP **gp=NULL; /* holds the voronoi data for Gneu. Therefore562we will only calculate it once */563564if (G->dim > MAXDIM){565fprintf(stderr,"This program does only work up to dimension %d\n",566MAXDIM);567exit(3);568}569570/* deal with the v_4 part in a different function */571/* it returns T != NULL iff the group was isomorphic to V_4572or C_2 = <-I_n>. In this case everything has been done */573/* T = catalog_number_v4(G,symb,TR,almost,zclass); */574T = NULL;575576if (T==NULL){577578/* inserted tilman 6/08/97, choose a better basis */579B = init_mat(G->dim,G->dim,"1");580F = rform(G->gen,G->gen_no,B,101);581BI = pair_red(F,B);582free_mat(BI);583free_mat(F);584BI = tr_pose(B);585free_mat(B);586B = BI;587BI = long_mat_inv(B);588Gneu = konj_bravais(G,BI);589590/* choose a good basis for the space of fixed forms of Gneu */591long_rein_formspace(Gneu->form,Gneu->form_no,1);592593/* initialize */594almost[0] = 0;595S = read_symbol_from_string(symb);596get_zentr(S);597Gtr = tr_bravais(Gneu,1,FALSE);598599while (T == NULL && S != NULL){600almost[0]++;601zclass[0] = -1;602while (T == NULL && (zclass[0] < S->grp->zentr_no )){603if (zclass[0] == -1){604H = S->grp;605606/* throw away the normalizer & centralizer, it only hinders607calculation */608for (i=0;i<H->cen_no;i++) free_mat(H->cen[i]);609if (H->cen != NULL && H->cen_no > 0)610free(H->cen);611H->cen_no = 0;612H->cen = NULL;613for (i=0;i<H->normal_no;i++) free_mat(H->normal[i]);614if (H->normal != NULL && H->normal_no > 0)615free(H->normal);616H->normal_no = 0;617H->normal = NULL;618}619else{620H = Z_class(S->grp,S->grp->zentr[zclass[0]]);621}622623if (Gneu->order == 0 ||624(Gneu->order == H->order)){625Htr = tr_bravais(H,1,FALSE);626X = is_z_equivalent_datei(Gneu,Gtr,H,Htr,&gp,&anz_gperfect);627free_bravais(Htr);628}629else{630X = NULL;631}632633if (X == NULL && zclass[0] != -1){634free_bravais(H);635}636else if (X!=NULL) {637T = H;638TR[0] = X;639}640zclass[0]++;641}642643if (zclass[0] > 0 || T == NULL) free_bravais(S->grp);644file = S->fn;645free(S);646647/* get the next almost decomposable group if necessary */648if (T == NULL){649650/* there might be the case that we dealt with all groups in the651catalog, but we didn't find an appropriate one */652if (file == NULL){653fprintf(stderr,"An error occured: This bravais group is not\n");654fprintf(stderr,"in the catalog.\n");655fprintf(stderr,"Please report this to\n");656fprintf(stderr," [email protected]\n");657fprintf(stderr,"immediately, with a copy of your input file and\n");658fprintf(stderr,"this message. (And the bravais group echoed to\n");659fprintf(stderr,"stdout now)\n");660put_bravais(G,NULL,"ERROR IN: catalog_number");661exit(3);662}663S = get_symbol(file);664}665666if (file != NULL) free(file);667}668669/* we started with -1, so we have to correct this "error" */670zclass[0]++;671672/* free the voronoi data for Gneu */673for (i=0;i<anz_gperfect;i++){674clear_voronoi(gp[i]);675free(gp[i]);676}677free(gp);678679/* adjust the transformations matrix */680/* i.e replace TR[0] by TR[0]*BI */681mat_muleq(TR[0],BI);682683free_bravais(Gtr);684free_bravais(Gneu);685free_mat(B);686free_mat(BI);687}688689return T;690691}692693694