###################################################################### ## ISOMERS Save this file as ISOMERS to use it, # # stay in the # ## same directory, get into Maple (by typing: maple ) # ## and then type: read ISOMERS # ## Then follow the instructions given there # ## # ## Written by Doron Zeilberger, Rutgers University , # ## zeilberg at math dot rutgers dot edu. # ###################################################################### print(`First Written: May, 2010: tested for Maple 13 `): print(` Version of July 20, 2014; Thanks to John Ogilvie`): print(): print(`This is ISOMERS, A Maple package to count isomers`): print(`Warning: It is VERY slow!`): print(): print(`The most current version is available on WWW at:`): print(` http://www.math.rutgers.edu/~zeilberg/tokhniot/ISOMERS .`): print(`Please report all bugs to: zeilberg at math dot rutgers dot edu .`): print(): print(`For general help, and a list of the MAIN functions,`): print(` type "ezra();". For specific help type "ezra(procedure_name);" `): #print(`For a list of the supporting functions type: ezra1();`): print(): ezra1:=proc() if args=NULL then print(`The supporting procedures are CH, CNOH, CNOHtrees `): print(` DegSeqsCH, Image1, Images, LaCG, MatToG, NIVecs, NIVecs1, SM, `): print(` UCG, Vecs `): else ezra(args): fi: end: ezraCNO:=proc() if args=NULL then print(`The procedures for the generalized case of`): print(`Carbon, Nitrogen and Oxygen are `): print(`CNOH, CNOHtrees, DegSeqsCNOH, ImagesCNO, UCGcno `): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(` ISOMERS: A Maple package for `): print(`The MAIN procedures are`): print(` CH, CHNO, CHtrees `): elif nargs=1 and args[1]=CH then print(`CH(m,n): all the isomers of a molecule with m Carbons`): print(`and n Hydrogens. `): print(``): print(`The output is the set of isomers. Since the carbon atoms are indistinguishale, each such isomer is an ULABELLED graph`): print(`with m Carbons. However each isomer is represented, for convenience sake, as a Labelled graph`): print(`where the chosen labeling is random, and at the end of the day one has to remove the labels`): print(`Each such isomer graph, let's call it G, is represended as a list of length m where G[i] is the set of`): print(`neighboring vertices(carbons), followed by 1 for a single bond, 2 for a double bond, extra`): print(`For example: A:=CH(3,2); `): print(`outputs`): print(`{[{[2, 2], [3, 2]}, {[1, 2], [3, 1]}, {[1, 2], [2, 1]}], [{[2, 3], [3, 1]}, {[1, 3], [3, 1]}, {[1, 1], [2, 1]}]}`): print(`a set of two elements, meaning that C3H2 has two isomers`): print(`The first one is`): print(`[{[2, 2], [3, 2]}, {[1, 2], [3, 1]}, {[1, 2], [2, 1]}]`): print(`which means that (lableling the three carbons arbitrarily C1,C2,C3`): print(`C1 is connected to C2 by a double bond, and to C3 by a double bond`): print(`C2 is connected to C1 by a double bond (of course we already know that), and C2 is connected to C3 by a single bond`): print(`and hence C2 is also connected to a Hydrogen atom`): print(`Finally, C3 is connected to C1 by a double bond (we already know it, by symmetry), and C3 is connected to C2 by a single bond (again we already know it)`): print(`and hence C3 is connected to a hydrogen atom`): print(`The second isomer in A is`): print(` [{[2, 3], [3, 1]}, {[1, 3], [3, 1]}, {[1, 1], [2, 1]}]}`): print(`which means that (lableling the three carbons arbitrarily C1,C2,C3`): print(`C1 is connected to C2 by a triple bond, and to C3 by a single bond`): print(`C2 is connected to C1 by a triple bond (of course we already know that), and C2 is connected to C3 by a single bond`): print(`Finally, C3 is connected to C1 by a single bond (we already know it, by symmetry), and C3 is connected to C2 by a single bond`): print(`and hence C3 is connected to the two hydrogen atoms`): print(`To see the set of 217 isomers of Benzene (including itself), try:`) : print(`CH(6,6);`): elif nargs=1 and args[1]=CHtrees then print(`CHtrees(m,n): all the isomers of a molecule with m Carbons`): print(`and n Hydrogens that are trees. For example, for the set if isomers`): print(`of Benzene that are trees, type:`): print(`CHtrees(6,6);`): elif nargs=1 and args[1]=CHNO then print(`CHNO(m1,n,m2,m3): all the isomers of a molecule with m1 Carbons`): print(` n Hydrogens `): print(`m2 Nitrogens and m3 Oxygens,`): print(`The convention is the same as for CH (q.v). and the labellings are: first the m1 Carbons, then the m2 Nitrogens, then the `): print(`m3 Oxygens. `): print(`For example, for the set of Nicotine (Warning: It will never finish!`): print(` you would have typed:`): print(`CHNO(10,14,2,0);`): print(`BUT don't!, since it would take many years`): print(`But to see the 1005 isomers of`): print(`C4H4N2, type`): print(`CHNO(4,4,2,0);`): elif nargs=1 and args[1]=CHNOtrees then print(`CHNOtrees(m1,m2,m3,n): `): print(`all the isomers of a molecule with m1 Carbons`): print(`m2 Nitrogens and m3 oxygens, `): print(`and n Hydrogens that are trees. For example, for the set if isomers`): print(`of Benzene that are trees, type:`): print(`CHNOtrees(6,0,0,6);`): elif nargs=1 and args[1]=CNOH then print(`CNOH(m1,m2,m3,n): all the isomers of a molecule with m1 Carbons`): print(`m2 Nitrogens and m3 Oxygens,`): print(`and n Hydrogens. For example, for the set of isomers`): print(`of Benzene, type:`): print(`CNOH(6,0,0,6);`): elif nargs=1 and args[1]=DegSeqsCH then print(`DegSeqsCH(m,n): all the possible degree-sequences of `): print(`a compound with m Carbon and n Hydroden atoms. `): print(`For example, for Benzene, do:`): print(`DegSeqsCH(6,6);`): elif nargs=1 and args[1]=DegSeqsCNOH then print(`DegSeqsCNOH(m1,m2,m3,n): all the possible degree-sequences `): print(`for triples of degree-sequences`): print(`a compound with m1 Carbon, m2 Nitrogen, and m3 oxygen`): print(`and n Hydroden atoms.`): print(`For example, for C4NOH6 do:`); print(`DegSeqsCNOH(4,1,1,6);`): elif nargs=1 and args[1]=Image1 then print(`Image1(G,pi): the images of the graph G`): print(`under the permutation pi`): print(`try Image1([{[2,2]},{[1,2]}],[2,1]);`): elif nargs=1 and args[1]=Images then print(`Images(G): all the images of the graph G, try:`): print(`Images([{[2,2]},{[1,2]}]);`): elif nargs=1 and args[1]=ImagesCNO then print(`ImagesCNO(G,m1,m2,m3): all the images of the graph G, `): print(`with m1 Carbon, m2 Nitrogen, and m3 Oxygen atoms. Try:`): print(`ImagesCNO([{[2,2]},{[1,2]}],1,1,0);`): elif nargs=1 and args[1]=LaCG then print(`LaCG(DegL): the set of labelled connected graphs with `): print(`degree sequence DegL on {1,...,n} where n=nops(DegL)`): print(`for example try`): print(`LaCG([2,2,2]);`): elif nargs=1 and args[1]=MatToG then print(`MatToG(M): inputs a symmetric matrix and outputs`): print(` it as a list whose i-th term is`): print(`the set of neighbors of vertex i with the multiplicity`): print(`of the edge joining them in the form {[j,a]}`): print(`For example, try:`): print(`MatToG([[0,1],[1,0]]);`): elif nargs=1 and args[1]=NIVecs then print(`NIVecs(k,N,a): the set of weakly-decreading`): print(`vectors of non-neg. integers of length k <=a`): print(`that add-up to N `): print(` Try:`): print(`NIVecs(3,4,a);`): elif nargs=1 and args[1]=NIVecs1 then print(`NIVecs1(k,a,N): the set of weakly-decreading`): print(`vectors of non-neg. integers of length k`): print(`that add-up to N and that start with a`): print(` Try:`): print(`NIVecs1(3,2,4);`): elif nargs=1 and args[1]=SM then print(`SM(DegL): inputs a degree sequence DegL, a weakly-decreasing`): print(`vectors of non-negative integers and outputs all`): print(`symmetric matrices with that row sums. For example`): print(`try SM([2,2,2]);`): elif nargs=1 and args[1]=UCG then print(`UCG(DegL): the set of unlabelled connected graphs with `): print(`degree sequence DegL on {1,...,n} where n=nops(DegL)`): print(`for example try`): print(`UCG([2,2,2]);`): elif nargs=1 and args[1]=UCGcno then print(`UCGcno(DegL,m1,m2,m3): the set of unlabelled connected graphs with `): print(`degree sequence DegL on {1,...,n} where n=nops(DegL)`): print(`and m1 Carbon atoms, m2 Nitrogen atoms and m3 Oxygen atoms`): print(`for example try`): print(`UCGcno([2,2,2],1,1,1);`): elif nargs=1 and args[1]=Vecs then print(`Vecs(k,N): the set of vectors of non-neg. integers of length k`): print(`that add-up to N. Try`): print(`Vecs(3,4);`): else print(`There is no such thing as`, args): fi: end: with(combinat): #Vecs(k,N): the set of vectors of non-neg. integers of length k #that add-up to N. Try: #Vecs(3,4); Vecs:=proc(k,N) local gu,gu1,a,g1: option remember: if k=0 then if N=0 then RETURN({[]}): else RETURN({}): fi: fi: gu:={}: for a from 0 to N do gu1:=Vecs(k-1,N-a): gu:=gu union {seq([op(g1),a],g1 in gu1)}: od: gu: end: #NIVecs1(k,a,N): the set of weakly-decreading #vectors of non-neg. integers of length k #that add-up to N and that start with a # Try: #NIVecs1(3,2,4); NIVecs1:=proc(k,a,N) local gu,gu1,a1,g1: option remember: if a>N then RETURN({}): fi: if k<=0 then RETURN({}): fi: if k=1 then if N=a then RETURN({[a]}): else RETURN({}): fi: fi: gu:={}: for a1 from 0 to a do gu1:=NIVecs1(k-1,a1,N-a): gu:=gu union {seq([a,op(g1)],g1 in gu1)}: od: gu: end: #NIVecs(k,N,a): the set of weakly-decreading #vectors of non-neg. integers <=a of length k #that add-up to N # Try: #NIVecs(3,4,2); NIVecs:=proc(k,N,a) local a1: if k=0 then if N=0 then RETURN({[]}): else RETURN({}): fi: fi: {seq(op(NIVecs1(k,a1,N)), a1 =0..min(a,N))}: end: #DegSeqsCH(m,n): all the possible degree-sequences of #a compound with m Carbon and n Hydroden atoms. #For example, for Benzene, do: #DegSeqsCH(6,6); DegSeqsCH:=proc(m,n) local gu,g,i1: gu:=NIVecs(m,n,4): {seq([seq(4-g[nops(g)+1-i1],i1=1..nops(g))],g in gu)}: end: #FindPlaces(L,a): finds the list of places in the list L that equal a #For example, try: #FindPlaces([3,3,1,2,3]); FindPlaces:=proc(L,a) local gu,i: gu:=[]: for i from 1 to nops(L) do if L[i]=a then gu:=[op(gu),i]: fi: od: gu: end: #Sader(L): Sader:=proc(L) local L1,gu,L1a,rish,gu1,i: L1:=sort(L): L1:=[seq(L1[nops(L1)+1-i],i=1..nops(L1))]: L1a:=L1: gu:=[]: while L1a<>[] do rish:=L1a[1]: gu1:=FindPlaces(L,rish): gu:=[op(gu),op(gu1)]: L1a:=[op(nops(gu1)+1..nops(L1a),L1a)]: od: L1,gu: end: #SM(DegL): inputs a degree sequence DegL, a weakly-decreasing #vectors of non-negative integers and outputs all #symmetric matrices with that row sums. For example #try SM([2,2,2]); SM:=proc(DegL) local gu,d1,n,mu,mu1,DegL1,gu1,DegL1a,pi,gu2, i1,j1,i,g1: option remember: if nops(DegL)=1 then if DegL[1]=0 then RETURN({[[0]]}): else RETURN({}): fi: fi: n:=nops(DegL): d1:=DegL[1]: mu:=Vecs(n-1,d1): gu:={}: for mu1 in mu do DegL1:=[seq(DegL[i]-mu1[i-1],i=2..n)]: DegL1a:=Sader(DegL1): pi:=DegL1a[2]: DegL1a:=DegL1a[1]: gu1:=SM(DegL1a): for g1 in gu1 do for i1 from 1 to nops(g1) do for j1 from 1 to nops(g1) do gu2[pi[i1],pi[j1]]:=g1[i1][j1]: od: od: gu2:=[seq([seq(gu2[i1,j1],j1=1..nops(g1))],i1=1..nops(g1))]: gu2:=[[0,op(mu1)],seq([mu1[i1],op(gu2[i1])],i1=1..nops(gu2))]: gu:=gu union {gu2}: od: od: gu: end: #MatToG(M): inputs a symmetric matrix and outputs # it as a list whose i-th term is #the set of neighbors of vertex i with the multiplicity #of the edge joining them {[j,a]} MatToG:=proc(M) local n,gu,i,j,gu1: n:=nops(M): gu:=[]: for i from 1 to n do gu1:={}: for j from 1 to n do if M[i][j]>0 then gu1:=gu1 union {[j,M[i][j]]}: fi: od: gu:=[op(gu),gu1]: od: gu: end: #MatToG1(M): inputs a symmetric matrix and outputs # it as a list whose i-th term is #the set of neighbors of vertex i #of the edge (ignoring multiple edges) MatToG1:=proc(M) local n,gu,i,j,gu1: n:=nops(M): gu:=[]: for i from 1 to n do gu1:={}: for j from 1 to n do if M[i][j]>0 then gu1:=gu1 union {j}: fi: od: gu:=[op(gu),gu1]: od: gu: end: #IsConn(G): Given a simple graph G given in terms #of a list of neighbors of vertices {1,...n} #with nops(G)=n, decides whether it is connected. Try: #IsConn([{2},{1}]); IsConn:=proc(G) local mu,mu1,NewGuys,m1,n,i: n:=nops(G): mu:={1}: mu1:={1} union G[1]: while mu<>mu1 do NewGuys:= {seq(op(G[m1]),m1 in mu1)}: mu:=mu1: mu1:=mu1 union NewGuys: od: evalb(mu1={seq(i,i=1..n)}): end: #LaCG(DegL): the set of labelled connected graphs with #degree sequence DegL on {1,...,n} where n=nops(DegL) #for example try #LaCG([2,2,2]); LaCG:=proc(DegL) local gu,mu,m: option remember: mu:=SM(DegL): gu:={}: for m in mu do if IsConn(MatToG1(m)) then gu:=gu union {MatToG(m)}: fi: od: gu: end: #LaT(DegL): the set of labelled connected graphs with #degree sequence DegL on {1,...,n} where n=nops(DegL) #for example try #LaT([2,2,2]); LaT:=proc(DegL) local gu,mu,m,lu,n,i: option remember: n:=nops(DegL): mu:=SM(DegL): gu:={}: for m in mu do lu:=MatToG1(m): if IsConn(lu) and add(nops(lu[i]),i=1..nops(lu))=2*n-2 then gu:=gu union {MatToG(m)}: fi: od: gu: end: #Image1(G,pi): the images of the graph G #under the permutation pi #try Image1([{[2,2]},{[1,2]}],[2,1]); Image1:=proc(G,pi) local n, T,lu,lu1,i: n:=nops(G): for i from 1 to n do lu:=G[i]: T[pi[i]]:={seq([pi[lu1[1]],lu1[2]],lu1 in lu)}: od: [seq(T[i],i=1..n)]: end: #Images(G): all the images of the graph G, try: #Images([{[2,2]},{[1,2]}]); Images:=proc(G) local n,gu,pi: n:=nops(G): gu:=permute(n): {seq(Image1(G,pi), pi in gu)}: end: #UCG(DegL): the set of unlabelled connected graphs with #degree sequence DegL on {1,...,n} where n=nops(DegL) #for example try #UCG([2,2,2]); UCG:=proc(DegL) local gu,mu: option remember: mu:=LaCG(DegL): gu:={}: while mu<>{} do gu:=gu union {mu[1]}: mu:=mu minus Images(mu[1]): od: gu: end: #UCT(DegL): the set of unlabelled trees with #degree sequence DegL on {1,...,n} where n=nops(DegL) #for example try #UCT([2,2,2]); UCT:=proc(DegL) local gu,mu: option remember: mu:=LaT(DegL): gu:={}: while mu<>{} do gu:=gu union {mu[1]}: mu:=mu minus Images(mu[1]): od: gu: end: #CH(m,n): all the isomers of a molecule with m Carbons #and n Hydrogens. For example, for the set if isomers #of Benzene, type: #CH(6,6); CH:=proc(m,n) local gu,g: option remember: gu:=DegSeqsCH(m,n): {seq(op(UCG(g)), g in gu)}: end: #CHtrees(m,n): all the isomers of a molecule with m Carbons #and n Hydrogens that are trees. For example, for the set if isomers #of Benzene, type: #CH(6,6); CHtrees:=proc(m,n) local gu,g: option remember: gu:=DegSeqsCH(m,n): {seq(op(UCT(g)), g in gu)}: end: ####New stuff with generalized C,N,O and H #DegSeqsCNOH(m1,m2,m3,n): all the possible degree-sequences #for triples of degree-sequences #a compound with m1 Carbon, m2 Nitrogen, and m3 oxygen #and n Hydroden atoms. #For example, for C4NOH6 do: #DegSeqsCNOH(4,1,1,6); DegSeqsCNOH:=proc(m1,m2,m3,n) local mu,mu1,lu1a,lu2a,lu3a,gu1,lu1,lu2, lu3,i1,gu: gu:={}: mu:=Vecs(3,n): for mu1 in mu do lu1:=NIVecs(m1,mu1[1],4): lu2:=NIVecs(m2,mu1[2],3): lu3:=NIVecs(m3,mu1[3],2): for lu1a in lu1 do for lu2a in lu2 do for lu3a in lu3 do gu1:=[seq(4-lu1a[i1],i1=1..m1),seq(3-lu2a[i1],i1=1..m2), seq(2-lu3a[i1],i1=1..m3)]: gu:=gu union {gu1}: od: od: od: od: gu: end: #UCGcno(DegL,m1,m2,m3): the set of unlabelled connected graphs with #degree sequence DegL on {1,...,n} where n=nops(DegL) #and m1 Carbon atoms, m2 Nitrogen atoms and m3 Oxygen atoms #for example try #UCGcno([2,2,2],1,1,1); UCGcno:=proc(DegL,m1,m2,m3) local gu,mu: option remember: if nops(DegL)<>m1+m2+m3 then print(`bad input`): RETURN(FAIL): fi: mu:=LaCG(DegL): gu:={}: while mu<>{} do gu:=gu union {mu[1]}: mu:=mu minus ImagesCNO(mu[1],m1,m2,m3): od: gu: end: #ImagesCNO(G,m1,m2,m3): all the images of the graph G, #with m1 Carbon, m2 Nitrogen, and m3 Oxygen atoms. Try: #ImagesCNO([{[2,2]},{[1,2]}],1,1,0); ImagesCNO:=proc(G,m1,m2,m3) local lu1,lu2,lu3,i,pi,lu1a,lu2a,lu3a,gu: if nops(G)<>m1+m2+m3 then RETURN(FAIL): fi: lu1:=permute([seq(i,i=1..m1)]): lu2:=permute([seq(i,i=m1+1..m1+m2)]): lu3:=permute([seq(i,i=m1+m2+1..m1+m2+m3)]): gu:={}: for lu1a in lu1 do for lu2a in lu2 do for lu3a in lu3 do pi:=[op(lu1a),op(lu2a),op(lu3a)]: gu:=gu union {Image1(G,pi)}: od: od: od: gu: end: #CNOH(m1,m2,m3,n): all the isomers of a molecule with m1 Carbons #m2 Nitrogens and m3 Hydrogens and #and n Hydrogens. For example, for the set if isomers #of Benzene, type: #CNOH(6,0,0,6); CNOH:=proc(m1,m2,m3,n) local gu,g: option remember: gu:=DegSeqsCNOH(m1,m2,m3,n): {seq(op(UCGcno(g,m1,m2,m3)), g in gu)}: end: #CHNO(m1,n,m2,m3): all the isomers of a molecule with m1 Carbons #n Hydrogens, m2 Nitrogens and m3 Hydrogens and #and . For example, for the set if isomers #of Benzene, type: #CHNO(6,6,0,0); CHNO:=proc(m1,n,m2,m3): CNOH(m1,m2,m3,n): end: #UCTcno(DegL,m1,m2,m3): the set of unlabelled connected trees with #degree sequence DegL on {1,...,n} where n=nops(DegL) #and m1 Carbon atoms, m2 Nitrogen atoms and m3 Oxygen atoms #for example try #UCTcno([2,2,2],1,1,1); UCTcno:=proc(DegL,m1,m2,m3) local gu,mu: option remember: if nops(DegL)<>m1+m2+m3 then print(`bad input`): RETURN(FAIL): fi: mu:=LaT(DegL): gu:={}: while mu<>{} do gu:=gu union {mu[1]}: mu:=mu minus ImagesCNO(mu[1],m1,m2,m3): od: gu: end: #CNOHtrees(m1,m2,m3,n): all the isomers of a molecule with m1 Carbons #m2 Nitrogens and m3 Oxygens #and n Hydrogens that are trees. For example, for the set if #tree-isomers #of Benzene, type: #CNOHtrees(6,0,0,6); CNOHtrees:=proc(m1,m2,m3,n) local gu,g: option remember: gu:=DegSeqsCNOH(m1,m2,m3,n): {seq(op(UCTcno(g,m1,m2,m3)), g in gu)}: end: