You can't change which color is opposite which color on the cubes. So each cube boils down to the 3 pairs of colors that are opposite each other on that cube. If we label the stack of cubes to be built with X, Y, and Z axes, then you can position each cube such that any of the three pairs shows on the X axis in either order and either of the remaining pairs shows on the Y axis in either order. It doesn't do any good to change the order in which the cubes are stacked. For the first cube, the only positions that matter are 3 of them, each cooresponding to a different pair of faces being hidden on the Z axis (the other positions would just boil down to duplicates where the stack is rotated or all 4 cubes on the stack are flipped along the same axis).

@cubes contains the four cubes, each as an array of color abbreviations grouped in the array as pairs of colors that are on opposite faces of the cube.

So we let the first cubes assume 3 different positions. We let each subsequent cube assume all 24 possible positions. For simplicity and speed, I precalculate arrays that show which face goes on which side. @base holds three positions that matter for the first cube. @full contains the mappings for all 24 orientations for subsequent cubes.

@side contains the colors showing on the sides of the stack as it is built, the first two colors showing on opposite sides of the X axis, the last two on the Y axis. setSides() adds the next 4 faces, adding one color to each side. If a duplicate color is found, then 0 is returned immediately. Otherwise, 1 is returned after all of the sides have been updated.

#!/usr/bin/perl -w
use strict;

my @cubes= map [/./g], 'ppgygr', 'rrprgy', 'prpyyg', 'ryrgpy';
my @sides= ( '', '', '', '' );
my @base= map [/./g], '0123', '0145', '2345';
my @full;
# Populate @full:
for my $base (  @base  ) {
    my @perm= @$base;
    for( 0..1 ) {
        for( 0..1 ) {
            for( 0..1 ) {
                push @full, [ @perm ];
                @perm[0,1,2,3]= @perm[2,3,0,1]; # Rotate 90deg; swap X
+/Y
            }
            @perm[2,3]= @perm[3,2]; # Flip to reverse faces on Y axis
        }
        @perm[0,1]= @perm[1,0];     # Flip to reverse faces on X axis
    }
    die   if  "@perm" ne "@$base";
}

# Try to explain format of the solution output:
print "For cubes 1, 2, 3, and 4.  The colors showing:\n";
print "frnt back left right\n";
print "1234 1234 1234 1234\n";

# Build the stack one cube at-a-time:
for my $perm ( @base ) {
    next   if  ! setSides( 0, $perm );
    for my $perm ( @full ) {
        next   if  ! setSides( 1, $perm );
        for my $perm ( @full ) {
            next   if  ! setSides( 2, $perm );
            for my $perm ( @full ) {
                if(  setSides( 3, $perm )  ) {
                    print "@sides\n";
                }
            }
        }
    }
}

sub setSides {
    my( $pos, $perm )= @_;
    for my $cube (  $cubes[$pos]  ) {
        for my $idx (  0..3  ) {
            for my $side (  $sides[$idx]  ) {
                for my $face (  $cube->[ $perm->[$idx] ]  ) {
                    substr( $side, $pos )= $face;
                    return 0   if  index( $side, $face ) < $pos;
                }
            }
        }
    }
    return 1;
}
[download]

Since the number of possible configurations for each of the 4 cubes is small (24), I generated them all up front and limited the data to only the 4 faces that would actually show (not the top or bottom faces). I then used a brute force method to search all possible combinations of configurations, but I didn't permute the order of the blocks (which would add another factor of 24 to the total number of solutions found). I also used extensive short-circuiting to avoid checking combinations that were already known to be failures.

I printed all of the solutions to STDOUT as well as the total number found at the end. The format for each solution is as follows:

ggry
ypgr
rypp
pryg
[download]

As you noted, 8 solutions are reported for each unique solution.

use strict;
use warnings;

# cube strings are the colors for the left, front, right, back, top
# and bottom faces, respectively
my @cubes = ( 'pgpygr', 'rprrgy', 'ppryyg', 'rrygpy' );
my ( $num_solutions, @opts );

foreach my $cube ( @cubes )
{
    push( @opts, get_all_opts( $cube ) );
}

foreach my $cube0str ( @{ $opts[0] } )
{
    foreach my $cube1str ( @{ $opts[1] } )
    {
        next if has_duplicate_faces( $cube0str, $cube1str );

        foreach my $cube2str ( @{ $opts[2] } )
        {
            next if has_duplicate_faces( $cube0str, $cube2str );
            next if has_duplicate_faces( $cube1str, $cube2str );

            foreach my $cube3str ( @{ $opts[3] } )
            {
                next if has_duplicate_faces( $cube0str, $cube3str );
                next if has_duplicate_faces( $cube1str, $cube3str );
                next if has_duplicate_faces( $cube2str, $cube3str );

                # we have a winner
                print "**********\n";
                print join( "\n", $cube0str, $cube1str,
                                  $cube2str, $cube3str ), "\n";

                $num_solutions++;
            }
        }
    }
}

printf "\n\nFound $num_solutions %s (%d unique %s),\n",
      $num_solutions == 1 ? 'solution' : 'solutions',
      $num_solutions / 8,
      $num_solutions / 8 == 1 ? 'solution' : 'solutions';
print "but by permuting the block order this total can be ",
      "increased by a factor of 24\n";

sub get_all_opts
{
    my ( $cubestr ) = @_;

    my @opts;

    # generate all 8 options for the rings around the
    # X, Y, and Z axes of the cube
    my @faces = split( //, $cubestr );
    push( @opts, permute_ring( @faces[ 0,1,2,3 ] ) );
    push( @opts, permute_ring( @faces[ 4,1,5,3 ] ) );
    push( @opts, permute_ring( @faces[ 0,4,2,5 ] ) );

    return( \@opts );
}

sub permute_ring
{
    my ( @faces ) = @_;

    # rotate and reverse the 4 faces
    my @opts;
    for ( 1 .. 4 )
    {
        push( @opts, join( '', @faces ) );
        push( @opts, scalar reverse join( '', @faces ) );
        push( @faces, shift( @faces ) );
    }

    return( @opts );
}

sub has_duplicate_faces
{
    my ( $cube1, $cube2 ) = @_;

    for( 0 .. 3 )
    {
        if( substr( $cube1, $_, 1 ) eq substr( $cube2, $_, 1 ) )
        {
            return 1;
        }
    }

    return 0;
}
[download]

#!/usr/bin/perl

use warnings;
use strict;

print ieronim_pure();

sub ieronim_pure {
    #my @cubes = qw/pgpygr rprrgy ppryyg rrygpy/;
    my @cubes = qw/rprrgy ppryyg rrygpy/;
    my %seen;
    my $str = join("", map {$_, scalar reverse $_} "pgpy","ggry","pygr
+")
                ."!".join "!", 
                map { join "", grep {!$seen{$_}++} map {join "", @$_} 
+@$_} 
                map { [ rotate([@{$_}[0,1,2,3]]),
                        rotate([@{$_}[4,1,5,3]]),
                        rotate([@{$_}[0,4,2,5]]) ] } map { [split //] 
+} @cubes;

    my @result;
    my $re = qr/^(?:[^!]{4})*
                 (?:
                    ([^!])
                    ([^!])
                    ([^!])
                    ([^!])
                 )
                 (?:[^!]{4})*!
                 (?:[^!]{4})*
                 (?:
                    (?!\1)([^!])
                    (?!\2)([^!])
                    (?!\3)([^!])
                    (?!\4)([^!])
                 )
                 (?:[^!]{4})*!
                 (?:[^!]{4})*
                 (?:
                    (?!\1|\5)([^!])
                    (?!\2|\6)([^!])
                    (?!\3|\7)([^!])
                    (?!\4|\8)([^!])
                 )
                 (?:[^!]{4})*!
                 (?:[^!]{4})*
                 (?:
                    (?!\1|\5|\9)([^!])
                    (?!\2|\6|\10)([^!])
                    (?!\3|\7|\11)([^!])
                    (?!\4|\8|\12)([^!])
                 )(?{
                    push @result, [join("", $1,$2,$3,$4), 
                                   join("", $5,$6,$7,$8), 
                                   join("", $9,$10,$11,$12),
                                   join("", $13,$14,$15,$16)]
                })
                (?!)
                /x;
    #print $str;

    $str =~ /$re/;
    return map { join "\n", @$_, "", ""} @result;
}


sub rotate {
    my @strips;
    $strips[0] = shift;
    $strips[1] = [reverse @{$strips[0]}];
    my @rotated;
    foreach (0..$#{$strips[0]}) {
        foreach my $strip (@strips) {
            push @$strip, shift @$strip;
            push @rotated, [@$strip];
        }
    }
    return @rotated;
}
[download]

#!/usr/bin/perl
use strict;
use warnings;

for (qw/pgpy ggyr gprp ypgp rygg prpg/) {    
    my @rotate = rotations($_);
    for my $cube1 (@rotate) {
        for (qw/rprr pgry gryr rrpr yrgp ryrg/) {
            my @rotate = rotations($_);
            for my $cube2 (@rotate) {
                for (qw/ppry pyyg yrgp yrpp gyyp pgry/) {
                    my @rotate = rotations($_);
                    for my $cube3 (@rotate) {
                        for (qw/rryg rpgy pyyr gyrr ygpr ryyp/) {
                            my @rotate = rotations($_);
                            for my $cube4 (@rotate) {
                                my $sol = solution($cube1, $cube2, $cu
+be3, $cube4);
                                print "$sol\n\n" if $sol;
                            }
                        }
                    }
                }
            }
        }
    }
}

sub rotations {
    my @rotate = shift @_;
    my @color  = split //, $rotate[0];
    (push @color, shift @color) && (push @rotate, join '', @color) for
+ 1 .. 3;
    return @rotate;
}

sub solution {
    my ($cube1, $cube2, $cube3, $cube4) = @_;
    for my $pos (0 .. 3) {
        my %uniq = map { substr($_, $pos, 1) => undef } $cube1, $cube2
+, $cube3, $cube4;
        return 0 if keys %uniq != 4;
    }
    return join "\n", $cube1, $cube2, $cube3, $cube4;
}
[download]

#!/usr/bin/perl
use strict;
use warnings;

my @cube = map rotations($_), qw/pgpygr rprrgy ppryyg rrygpy/;

for my $c1 (@{$cube[0]}) {
    for my $c2 (@{$cube[1]}) {
        for my $c3 (@{$cube[2]}) {
            for my $c4 (@{$cube[3]}) {
                my $sol = solution($c1, $c2, $c3, $c4);
                print "$sol\n\n" if $sol;
            }
        }
    }
}

sub rotations {
    my (%seen, @rot);
    my @cube = split //, shift @_;
    for ([0 .. 3], [1, 4, 3, 5], [4, 0, 5, 2]) {
        my @col = @cube[@$_];
        push @rot, map {push @col, shift @col; $seen{"@col"}++ ? () : 
+[@col]} 1..4;
        @col = reverse @{$rot[-1]};
        push @rot, map {push @col, shift @col; $seen{"@col"}++ ? () : 
+[@col]} 1..4;
    }
    return \@rot;
}

sub solution {
    my ($cube1, $cube2, $cube3, $cube4) = @_;
    for my $i (0 .. 3) {
        my %uniq = map { $_->[$i] => undef } $cube1, $cube2, $cube3, $
+cube4;
        return 0 if keys %uniq != 4;
    }
    return "@$cube1\n@$cube2\n@$cube3\n@$cube4";
}
[download]

#!/usr/bin/perl
use strict;
use warnings;

my @cube = map rotations($_), qw/pgpygr rprrgy ppryyg rrygpy/;
my @used;
for my $c1 (@{$cube[0]}) {
    @used = map {{$_ => 1}} @$c1;
    CUBE2:
    for my $c2 (@{$cube[1]}) {
        $used[$_]{$c2->[$_]} && next CUBE2 for 0 .. 3;
        $used[$_]{$c2->[$_]} = 1 for 0 .. 3;
        CUBE3:
        for my $c3 (@{$cube[2]}) {
            $used[$_]{$c3->[$_]} && next CUBE3 for 0 .. 3;
            $used[$_]{$c3->[$_]} = 1 for 0 .. 3;
            CUBE4:
            for my $c4 (@{$cube[3]}) {
                $used[$_]{$c4->[$_]} && next CUBE4 for 0 .. 3;
                print "@$c1\n@$c2\n@$c3\n@$c4\n\n";
            }
            $used[$_]{$c3->[$_]} = 0 for 0 .. 3;
        }
        $used[$_]{$c2->[$_]} = 0 for 0 .. 3;
    }
}

sub rotations {
    my (%seen, @rot);
    my @cube = split //, shift @_;
    for ([0 .. 3], [1, 4, 3, 5], [4, 0, 5, 2]) {
        my @col = @cube[@$_];
        push @rot, map {push @col, shift @col; $seen{"@col"}++ ? () : 
+[@col]} 1..4;
        @col = reverse @{$rot[-1]};
        push @rot, map {push @col, shift @col; $seen{"@col"}++ ? () : 
+[@col]} 1..4;
    }
    return \@rot;
}
[download]

#!/perl -w
use strict;

$|=1;

my @cubes = ('pgpygr', 'rprrgy', 'ppryyg', 'rrygpy');

foreach my $c1perm (permuate_colours($cubes[0])) {
    foreach my $c2perm (permuate_colours($cubes[1])) {
        foreach my $c3perm (permuate_colours($cubes[2])) {
            foreach my $c4perm (permuate_colours($cubes[3])) {
                if ( validate($c1perm, $c2perm, $c3perm, $c4perm) ) {
                    print "Solution: $c1perm, $c2perm, $c3perm, $c4per
+m\n";
                }
            }
        }
    }
}

sub validate {
    my (@colours) = @_;
    my %seen;
    foreach my $row (0..3) {
        foreach my $cube (0..3) {
            my $side_colour = substr($colours[$cube], $row, 1);
            if (exists $seen{$row}{$side_colour}) {
                return 0;
            }
            $seen{$row}{$side_colour}++;
        }
    }
    return 1;
}

sub permuate_colours {
    my ($colours) = @_;
    my @perms;
    for (1..4) {
        for (1..4) {
            push(@perms, $colours);
            $colours =~ s/(.)(.)(.)(.)(.)(.)/$2$3$4$1$5$6/; # spin
        }
        $colours =~ s/(.)(.)(.)(.)(.)(.)/$1$6$3$5$2$4/; # roll
    }

    $colours =~ s/(.)(.)(.)(.)(.)(.)/$2$3$4$1$5$6/; # spin
    $colours =~ s/(.)(.)(.)(.)(.)(.)/$1$6$3$5$2$4/; # roll
    for (1..4) {
        push(@perms, $colours);
        $colours =~ s/(.)(.)(.)(.)(.)(.)/$2$3$4$1$5$6/; # spin
    }
    $colours =~ s/(.)(.)(.)(.)(.)(.)/$1$6$3$5$2$4/; # roll
    $colours =~ s/(.)(.)(.)(.)(.)(.)/$1$6$3$5$2$4/; # roll
    for (1..4) {
        push(@perms, $colours);
        $colours =~ s/(.)(.)(.)(.)(.)(.)/$2$3$4$1$5$6/; # spin
    }
    return @perms;
}
[download]

NB. Very inefficent wrt memory and somewhat inefficent wrt speed,
NB. but no need to optimize on a modern machine.

NB. a cube is a list of 6 sides (each atoms): 
NB. left front right back top bottom.
NB. the first four of these sides are visible on the sides in a tower.
NB. a towers is a list of 4 cubes.

NB. colorings.
input =: >;: 'pgpygr rprrgy ppryyg rrygpy' 

NB. all 24 rotations of a cube.
gens =: (i.6),'NSKTFL'i.'FSLTKN',:'SKTNFL'
rots =: ([:~.[:/:~[:,/{"1/)^:_~ gens
   
NB. all 24 rotations of each 4 cube.
posn =: rots {"_ 1 input

NB. stack a list of 4 lists of 24 cubes to a tower of 4 cubes in each 
+24^4 way,
NB. order of 4 cubes doesn't matter.
stack =: [:>[:>[:,[:{[:<"1<"1
towers =: stack posn

NB. check if a cube is good or bad -- rank 2, gives bool atom.
good =: [: (*./"1) 4 ({."1) 4 = [: (#@:~."1) |:"2
solns =: (good # ]) towers

NB. each solution gives 8 trivially equivalent solutions because you c
+an 
NB. rotate all cubes simultanously around the vertical axis or twice a
+round
NB. the horizontal axis.  (this may be less than 8 for degenerate cube
+s.)
NB. rota y gives all 8 rotations of a tower y
rott =: (2{gens)&({"_ 1)
rotu =: ({~1{gens)&({"_ 1)
rota =: ([: ,/ [: rott^:(<4) rotu^:(<2))"2
NB. verify that each rotation of each solution is indeed among the sol
+utions.
assert =: 3 :'assert. y'
assert *./ solns e.~ ,/ rota solns
NB. now normalize the solutions by changing a tower to the 
NB. lexicographically first out of these eight variants,
NB. and weed out repetitions.
norml =: ([:{.[:/:~rota)"2
solns1 =: ~./:~ norml solns

NB. output the solutions. there's just one if you use the sample
NB. cube colorings in the post.
echo solns1
echo #solns1

NB. end
[download]

ggrypp
ypgrrr
ryppyg
prygyr
1
[download]

Most of the solutions just output the solution in a format similar to ours, except sometimes only the four visible sides are shown, so I visually compared one of their outputs to mine directly, and guessed the rest of the solutions they find would match as well. This table form is handy.

   <"2 solns
+------+------+------+------+------+------+------+------+
|ggyrpp|grygpp|yggrpp|yrggpp|ggrypp|gyrgpp|rggypp|ryggpp|
|pyrgrr|pgryrr|rypgrr|rgpyrr|ypgrrr|yrgprr|gpyrrr|gryprr|
|yrppgy|yppryg|prypyg|ppyrgy|ryppyg|rppygy|pyrpgy|ppryyg|
|rpgyry|rygpyr|gpryyr|gyrpry|prygyr|pgyrry|yrpgry|ygpryr|
+------+------+------+------+------+------+------+------+
[download]

To compare with the output of tye's script, we convert each solution to the format of his output. We can't use just the normalized solution because tye's script normalizes in a different way.

   ,/"2,"1&' ' |:"2 (1 3 0 2){"1 solns
gyrp rgpy gpyr yrpg 
rgpy gyrp gpyr yrpg 
gyrp rgpy yrpg gpyr 
rgpy gyrp yrpg gpyr 
gpyr yrpg gyrp rgpy 
yrpg gpyr gyrp rgpy 
gpyr yrpg rgpy gyrp 
yrpg gpyr rgpy gyrp
[download]

Indeed, output for index 4 is equal to tye's output.

I've checked the output of most solutions in this thread, and each seems to give at least one correct solution, or at most sneaked in so typos so smartly that I didn't notice them in my not too rigorous comparision. I don't yet understand the output of solutions in L~R's second order reply.