CJam, 189 187 bytes

This one's gonna be tough to explain... Time complexity is guaranteed to be O(scary). While writing the explanation I'm seeing a few golfs that I overlooked, but overall I doubt I'll make it much shorter.

qi:N_3>{,aN*]N({{:L;N,X)-e!{X)_@+L@@t}%{X2+<z{_fe=:(:+}%:+!},}%:+}fX{:G;N3m*{_~{G@==}:F~F\1m>~F\F=}%:*},:L,({LX=LX)>1$f{\_@a\a+Ne!\f{\:M;~{M\f=z}2*\Mff==}:|{;}|}\a+}fX]:~$e`{0=1=},,}{!!}? 

If you're brave enough, [try it online][1]. On my crappy laptop I can get up to 6 with the Java interpreter or 5 in the online interpreter.

Explanation

I don't have a big math background (just finished high school, starting CS at uni next week). So bear with me if I make mistakes, state the obvious, or do things in horribly inefficent ways.

My approach is a brute force, though I tried to make it a little more clever. The main steps are:

1. Generate all the possible operands ∗ for a group of order n (i.e., enumerate all groups of order n);
2. Generate all the possible bijections φ between two groups of order n;
3. Using the results from steps 1 and 2, determine all isomorphisms between two groups of order n;
4. Using the result from step 3, count the number of groups up to isomorphism.

Before looking at the how each step is done, let's get some trivial code out of the way:

qi:N_ e# Get input as integer, store in N, make a copy 3>{...} ? e# If N > 3, do... (see below) {!!} e# Else, push !!N (0 if N=0, 1 otherwise) 

The following algorithm doesn't work correctly with n < 4, the cases from 0 to 3 are handled with a double negation.

From now on, the elements of a group will be written as {1, a, b, c, ...}, where 1 is the identity element. In the CJam implementation, the corresponding elements are {0, 1, 2, 3, ...}, where 0 is the identity element.

Let's start from step 1. Writing all possible operators for a group of order n is equivalent to generating all the valid n×n [Cayley tables][2]. I'll use an order 4 group as example. The first row and column are trivial: they are both {1, a, b c} (left-to-right, up-to-down).

 e# N is on the stack (duplicated before the if) ,a e# Generate first row [0 1 2 3 ...] and wrap it in a list N* e# Repeat row N times (placeholders for next rows) ] e# Wrap everything in a list e# First column will be taken care of later 

Knowing that a Cayley table is also a reduced [Latin square][3] (due to the cancellation property) allows to generate the possible tables row-by-row. Starting from the second row (index 1), we generate all the unique permutations for that row, keeping the first column fixed to the value of the index.

N({ }fX e# For X in [0 ... N-2]: { }% e# For each table in the list: :L; e# Assign the table to L and pop it off the stack N, e# Push [0 ... N-1] X) e# Push X+1 - e# Remove X+1 from [0 ... N-1] e! e# Generate all the unique permutations of this list { }% e# For each permutation: X)_ e# Push two copies of X+1 @+ e# Prepend X+1 to the permutation L@@t e# Store the permutation at index X+1 in L {...}, e# Filter permutations (see below) :+ e# Concatenate the generated tables to the table list 

Not all those permutations are valid, of course: each row and column must contain all the elements exactly one time. A filter block is used for this purpose (a truthy value keeps the permutation, a falsy one removes it):

X2+ e# Push X+2 < e# Slice the permutations to the first X+2 rows z e# Transpose rows and columns { }% e# For each column: _fe= e# Count occurences of each element :( e# Subtract 1 from counts :+ e# Sum counts together :+ e# Sum counts from all columns together ! e# Negate count sum: e# if the sum is 0 (no duplicates) the permutation is kept e# if the sum is not zero the permutation is filtered away 

Note that I am filtering inside the generation loop: this makes the code a bit longer (compared to distinct generation and filtering), but greatly improves performance. The number of permutations of a set of size n is n!, so the shorter solution would require a lot of memory and time.

A list of valid Cayley tables is a great step towards enumerating the operators, but being a 2D structure, it can't check for associativity, which is a 3D property. So next step is filtering out non-associative functions.

{ }, e# For each table, keep table if result is true: :G; e# Assign table to G, pop it off the stack N3m* e# Generate triples [0 ... N-1]^3 { }% e# For each triple [a b c]: _~ e# Make a copy, unwrap top one { }:F e# Define function F(x,y): G@== e# x∗y (using table G) ~F e# Push a∗(b∗c) \1m> e# Rotate triple right by 1 ~ e# Unwrap rotated triple F\F e# Push (a∗b)∗c = e# Push 1 if a∗(b∗c) == (a∗b)∗c (associative), 0 otherwise :* e# Multiply all the results together e# 1 (true) only if F was associative for every [a b c] 

I'll be saving periodically as I work. To do: bijection generation, isomorphism testing, isomorphism counting.

[1]: http://cjam.aditsu.net/#code=qi%3AN_3%3E%7B%2CaN*%5DN(%7B%7B%3AL%3BN%2CX)-e!%7BX)%40%2BL%40%40t%7D%25%7BX2%2B%3Cz%7B_fe%3D%3A(%3A%2B%7D%25%3A%2B!%7D%2C%7D%25%3A%2B%7DfX%7B%3AG%3BN3m*%7B~%7BG%40%3D%3D%7D%3AF~F%5C1m%3E~F%5CF%3D%7D%25%3A*%7D%2C%3AL%2C(%7BLX%3DLX)%3E1%24f%7B%5C_%40a%5Ca%2BNe!%5Cf%7B%5C%3AM%3B~%7BM%5Cf%3Dz%7D2*%5CMff%3D%3D%7D%3A%7C%7B%3B%7D%7C%7D%5Ca%2B%7DfX%5D%3A~%24e%60%7B0%3D1%3D%7D%2C%2C%7D%7B!!%7D%3F&input=4 [2]: https://en.wikipedia.org/wiki/Cayley_table [3]: https://en.wikipedia.org/wiki/Latin_square