FSL's flirt operates only on single sub-bricks (i.e., a single volume or a single gradient) at one time, at least by default-- therefore, you would have to write an explicit loop in your processing script to use it, and stick those together. Actually, the FSL function called 'eddy_correct' is purely a wrapper to go through and do this linear registration process-- it does the same linear affine registration without changing the number of sub-bricks, nothing else is different.

The AFNI alignment tools, such as 3dAllineate, which I believe might be the default in align_epi_anat.py, will will go through sub-brick by sub-brick, preserving the total number of volumes.

Something to be careful of is that different registration programs will use different cost functions to say 'this is the rule I use to match these images'. That can have big consequences for the outputs, even though both programs may be 12 degree of freedom linear affine functions (which I believe both flirt and 3dAllineate are, by definition). It would probably be best to consider using the same function, particularly for DWI data whose contrasts change with each diffusion gradient-- that is, different parts of the brain will show up light and dark merely because of the applied direction of measuring diffusion.

Some usage of align_epi_anat.py (or other AFNI functions) has been addressed previously here:
[afni.nimh.nih.gov]
[blog.cogneurostats.com]
And I note in the latter that the 'mutualinfo' cost function is explicitly chosen in aligning DWIs.

Finally, one of the links above makes note of TORTOISE for processing, and that should be given very strong consideration. It does quite a comprehensive set of processing and correction even beyond motion distortion, such as for eddy currents and EPI distortion. It is slow to operate, as they themselves have noted in their acronym; and if you don't have an IDL license, it can't be batch processed by instead only operated one-at-a-time using a GUI, but it is definitely worth considering.

--pt