This example script relies on the example script Create mni-aligned grids in individual head_space. But for Neuromag data there are some differences. First make a mni template as is done in the above mentioned example script.

%%batch create_headmodel_neuromag
%==================================================================
%Loading mri of single subject and make a single shell head model
%==================================================================

filename           = 'mri.fif';
mri                = ft_read_mri(filename);
save mri mri

%The following makes the mri into ctf coordinates.. you have to select the
%fiducials again yourself
cfg                = [];
cfg.interactive    = 'yes';
mri                = ft_volumerealign(cfg,mri);
save mri_volreal mri

%segment the brain into gray, white and csf matter to later make a single
%shell model.
cfg                = [];
cfg.template       = '/data/apps/spm/spm2/templates/single_subject_T1.mnc';
cfg.coordinates    = 'ctf';
cfg.write          = 'no';
cfg.name           = 'temp';
[segmentedmri]     = ft_volumesegment(cfg, mri)
save segmentedmri segmentedmri

%check how it looks; does the segmented mri fit into the mri? Probably not
%because of neuromag coordinates (x and y are swapped) and a bug in volume
%segment.
cfg                = [];
test               = segmentedmri;
test.avg.pow       = test.gray+test.white+test.csf;
load mri_volreal
test.anatomy       = mri.anatomy;
cfg.funparameter   = avg.pow';
cfg.interactive    = no';
ft_sourceplot(cfg,test);

%------NEUROMAG SPECIFIC
% Swap the x and y coordinates back 
test.gray          = permute(segmentedmri.gray,[3 2 1]);
test.white         = permute(segmentedmri.white,[3 2 1]);
test.csf           = permute(segmentedmri.csf,[3 2 1]);
test.dim           = [256 256 192];
%------END
%and flip the dimensions back
for t = 1:3;
   test.gray       = flipdim(test.gray,t);
   test.csf        = flipdim(test.csf,t);
   test.white      = flipdim(test.white,t);
end

%check if it worked
test.avg.pow       = test.gray+test.white+test.csf;
figure;
cfg.interactive    = 'no';
ft_sourceplot(cfg,test);

test               = rmfield(test,'anatomy');
test               = rmfield(test,'avg');
segmentedmri       = test;
save flippedsegmentedmri segmentedmri

%make the single_shell headmodel
cfg                = [];
hdm                = ft_prepare_singleshell(cfg,segmentedmri);

%------NEUROMAG SPECIFIC
%Flip x and y coordinates to fit to the elekcta neuromag MEG helmet!
xcoords            = hdm.bnd.pnt(:,2);
ycoords            = hdm.bnd.pnt(:,1);
hdm.bnd.pnt(:,1)   = xcoords;
hdm.bnd.pnt(:,2)   = ycoords;
%------END

%make the leadfield normalised to the mni template
%normalise the mri first
cfg                = [];
cfg.template       = '/data/apps/spm/spm2/templates/T1.mnc'; %is in MNI coordinates, from templates/spm2
load mri_volreal
cfg.downsample     = 2;
cfg.coordinates    = 'ctf';
cfg.nonlinear      = 'no';
norm               = ft_volumenormalise(cfg,mri); % subjects m
%make the leadfield
%Use the final transform matrix of norm to make the grid.pos (and the template_grid)
load Template_brain/template_grid.mat
grid               = [];
grid.pos           = warp_apply(inv(norm.cfg.final), template_grid.pos, 'homogenous')/10; %in cm
%------NEUROMAG SPECIFIC
%Flip x and y coordinates to fit the elekcta neuromag MEG helmet!
x                  = grid.pos(:,2);
y                  = grid.pos(:,1);
grid.pos(:,1)      = x;
grid.pos(:,2)      = y;
%------END
grid.inside        = template_grid.inside;
grid.outside       = template_grid.outside;
grid.dim           = template_grid.dim;
clear template_grid

% also remember the normalization transformation matrix
grid.transform     = norm.cfg.final;

%to see if everything has worked:
%make a figure of the single subject headmodel, sensor positions and grid positions
cfg                = [];
cfg.vol            = hdm;
cfg.grid           = grid;
% cfg.hdmfile      = hdmfile{i};  % an alternative to cfg.vol, this uses the single sphere model
cfg.gradfile       = 'grad.mat';
figure;
ft_headmodelplot(cfg);

Figure 1 The results

Doing source-analysis with the created headmodel

load f
cfg                = [];
cfg.grid           = grid;
cfg.frequency      = 10;
cfg.vol            = hdm;
cfg.gradfile       = 'grad.mat';
cfg.projectnoise   = 'yes';
cfg.keeptrials     = 'no';
cfg.keepfilter     = 'yes';
cfg.keepcsd        = 'yes';
cfg.keepmom        = 'yes';
cfg.lambda         = 0.1 * mean(f.powspctrm(:,9),1);
cfg.method         = 'dics';
cfg.feedback       = 'textbar';
source             = ft_sourceanalysis(cfg,f);

save source source

load grid
cfg                = [];
source.dim         = grid.dim;
sd                 = ft_sourcedescriptives(cfg,source);
%-----MNI SPECIFIC
%Because sourceanalysis worked in NM coordinates and the interpolation goes
%in MNI coordinates, we have to replace the sd.pos by the
%template_grid.pos
load Template_brain/template_grid
sd.pos             = template_grid.pos;
%-----END

load norm
cfg                = [];
cfg.parameter      = 'avg.nai';
cfg.voxelcoord     = 'no';
cfg.interpmethod   = 'linear';
sdint              = ft_sourceinterpolate(cfg, sd, norm);

%plotting the results
cfg                = [];
cfg.surfdownsample = 2;
cfg.downsample     = 2;
cfg.funparameter   = 'avg.nai';
cfg.anaparameter   = 'anatomy';
cfg.funcolormap    = 'jet';
cfg.method         = 'ortho';

figure;
ft_sourceplot(cfg,sdint)