diff --git a/.gitignore b/.gitignore
index 6641a666..8c34210a 100755
--- a/.gitignore
+++ b/.gitignore
@@ -5,3 +5,6 @@
 .idea/
 venvs/
 _build/
+build/
+.vscode/
+
diff --git a/docs/heuristics.rst b/docs/heuristics.rst
index 874cc9b4..ec9f0e76 100644
--- a/docs/heuristics.rst
+++ b/docs/heuristics.rst
@@ -84,4 +84,66 @@ or::
     def grouping(files, dcmfilter, seqinfo):
         seqinfos = collections.OrderedDict()
         ...
-        return seqinfos  # ordered dict containing seqinfo objects: list of DICOMs
\ No newline at end of file
+        return seqinfos  # ordered dict containing seqinfo objects: list of DICOMs
+
+
+-------------------------------
+``POPULATE_INTENDED_FOR_OPTS``
+-------------------------------
+
+Dictionary to specify options to populate the ``'IntendedFor'`` field of the ``fmap``
+jsons.
+
+When a BIDS session has ``fmaps``, they can automatically be assigned to be used for
+susceptibility distortion correction of other non-``fmap`` images in the session
+(populating the ``'IntendedFor'`` field in the ``fmap`` json file).
+
+For this automated assignment, ``fmaps`` are taken as groups (``_phase`` and ``_phasediff``
+images and the corresponding ``_magnitude`` images; consecutive Spin-Echo images collected
+with opposite phase encoding polarity (``pepolar`` case); etc.).
+
+This is achieved by checking, for every non-``fmap`` image in the session, which ``fmap``
+groups are suitable candidates to correct for distortions in that image.  Then, if there is
+more than one candidate (e.g., if there was a ``fmap`` collected at the beginning of the
+session and another one at the end), the user can specify which one to use.
+
+The parameters that can be specified and the allowed options are defined in ``bids.py``:
+ - ``'matching_parameter'``: The imaging parameter that needs to match between the ``fmap``
+   and an image for the ``fmap`` to be considered as a suitable to correct that image.
+   Allowed options are:
+
+   * ``'Shims'``: ``heudiconv`` will check the ``ShimSetting`` in the ``.json`` files and
+     will only assign ``fmaps`` to images if the ``ShimSettings`` are identical for both.
+   * ``'ImagingVolume'``: both ``fmaps`` and images will need to have the same the imaging
+     volume (the header affine transformation: position, orientation and voxel size, as well
+     as number of voxels along each dimensions).
+   * ``'ModalityAcquisitionLabel'``: it checks for what modality (``anat``, ``func``, ``dwi``) each
+     ``fmap`` is intended by checking the ``_acq-`` label in the ``fmap`` filename and finding
+     corresponding modalities (e.g. ``_acq-fmri``, ``_acq-bold`` and ``_acq-func`` will be matched
+     with the ``func`` modality)
+   * ``'CustomAcquisitionLabel'``: it checks for what modality images each  ``fmap`` is intended
+     by checking the ``_acq-`` custom label (e.g. ``_acq-XYZ42``) in the ``fmap`` filename, and
+     matching it with:
+     - the corresponding modality image ``_acq-`` label for modalities other than ``func``
+     (e.g. ``_acq-XYZ42`` for ``dwi`` images)
+     - the corresponding image ``_task-`` label for the ``func`` modality (e.g. ``_task-XYZ42``)
+   * ``'Force'``: forces ``heudiconv`` to consider any ``fmaps`` in the session to be
+     suitable for any image, no matter what the imaging parameters are.
+
+
+ - ``'criterion'``: Criterion to decide which of the candidate ``fmaps`` will be assigned to
+   a given file, if there are more than one. Allowed values are:
+
+   * ``'First'``: The first matching ``fmap``.
+   * ``'Closest'``: The closest in time to the beginning of the image acquisition.
+
+.. note::
+  Example::
+
+    POPULATE_INTENDED_FOR_OPTS = {
+            'matching_parameters': ['ImagingVolume', 'Shims'],
+            'criterion': 'Closest'
+    }
+
+If ``POPULATE_INTENDED_FOR_OPTS`` is not present in the heuristic file, ``IntendedFor``
+will not be populated automatically.
diff --git a/heudiconv/bids.py b/heudiconv/bids.py
index 2043629c..40deb1e5 100644
--- a/heudiconv/bids.py
+++ b/heudiconv/bids.py
@@ -4,12 +4,14 @@
 import os
 import os.path as op
 import logging
+import numpy as np
 import re
 from collections import OrderedDict
 from datetime import datetime
 import csv
 from random import sample
 from glob import glob
+import errno
 
 from .external.pydicom import dcm
 
@@ -19,9 +21,12 @@
     save_json,
     create_file_if_missing,
     json_dumps,
+    update_json,
     set_readonly,
     is_readonly,
     get_datetime,
+    remove_suffix,
+    remove_prefix,
 )
 from . import __version__
 
@@ -51,6 +56,25 @@ class BIDSError(Exception):
 
 BIDS_VERSION = "1.4.1"
 
+# List defining allowed parameter matching for fmap assignment:
+SHIM_KEY = 'ShimSetting'
+AllowedFmapParameterMatching = [
+    'Shims',
+    'ImagingVolume',
+    'ModalityAcquisitionLabel',
+    'CustomAcquisitionLabel',
+    'Force',
+]
+# Key info returned by get_key_info_for_fmap_assignment when
+# matching_parameter = "Force"
+KeyInfoForForce = "Forced"
+# List defining allowed criteria to assign a given fmap to a non-fmap run
+# among the different fmaps with matching parameters:
+AllowedCriteriaForFmapAssignment = [
+    'First',
+    'Closest',
+]
+
 
 def maybe_na(val):
     """Return 'n/a' if non-None value represented as str is not empty
@@ -193,7 +217,7 @@ def populate_aggregated_jsons(path):
                 #   the next for loop iteration:
                 continue
 
-        events_file = fpath[:-len(suf)] + '_events.tsv'
+        events_file = remove_suffix(fpath, suf) + '_events.tsv'
         # do not touch any existing thing, it may be precious
         if not op.lexists(events_file):
             lgr.debug("Generating %s", events_file)
@@ -366,7 +390,7 @@ def save_scans_key(item, bids_files):
         subj_, ses_ = find_subj_ses(f_name)
         if not subj_:
             lgr.warning(
-                "Failed to detect fullfilled BIDS layout.  "
+                "Failed to detect fulfilled BIDS layout.  "
                 "No scans.tsv file(s) will be produced for %s",
                 ", ".join(bids_files)
             )
@@ -496,3 +520,507 @@ def convert_sid_bids(subject_id):
     lgr.warning('{0} contained nonalphanumeric character(s), subject '
                 'ID was cleaned to be {1}'.format(subject_id, sid))
     return sid, subject_id
+
+
+def get_shim_setting(json_file):
+    """
+    Gets the "ShimSetting" field from a json_file.
+    If no "ShimSetting" present, return error
+
+    Parameters:
+    ----------
+    json_file : str
+
+    Returns:
+    -------
+    str with "ShimSetting" value
+    """
+    data = load_json(json_file)
+    try:
+        shims = data[SHIM_KEY]
+    except KeyError as e:
+        lgr.error('File %s does not have "%s". '
+                  'Please use a different "matching_parameters" in your heuristic file',
+                   json_file, SHIM_KEY)
+        raise KeyError
+    return shims
+
+
+def find_fmap_groups(fmap_dir):
+    """
+    Finds the different fmap groups in a fmap directory.
+    By groups here we mean fmaps that are intended to go together
+    (with reversed PE polarity, magnitude/phase, etc.)
+
+    Parameters:
+    ----------
+    fmap_dir : str or os.path
+        path to the session folder (or to the subject folder, if there are no
+        sessions).
+
+    Returns:
+    -------
+    fmap_groups : dict
+        key: prefix common to the group (e.g. no "dir" entity, "_phase"/"_magnitude", ...)
+        value: list of all fmap paths in the group
+    """
+    if op.basename(fmap_dir) != 'fmap':
+        lgr.error('%s is not a fieldmap folder', fmap_dir)
+
+    # Get a list of all fmap json files in the session:
+    fmap_jsons = sorted(glob(op.join(fmap_dir, '*.json')))
+
+    # RegEx to remove fmap-specific substrings from fmap file names
+    # "_phase[1,2]", "_magnitude[1,2]", "_phasediff", "_dir-<label>", ...
+    fmap_regex = re.compile(
+        '(_dir-[0-9,a-z,A-Z]*)*'  # for pepolar case
+        '(_phase[12])*'  # for phase images
+        '(_phasediff)*'  # for phasediff images
+        '(_magnitude[12])*'  # for magnitude images
+        '(_fieldmap)*'  # for actual fieldmap images
+    )
+
+    # Find the unique prefixes ('splitext' removes the extension):
+    prefixes = sorted(
+        set(fmap_regex.sub('', remove_suffix(op.basename(fm), '.json')) for fm in fmap_jsons)
+    )
+    fmap_groups = OrderedDict()
+    for k in prefixes:
+        fmap_groups[k] = [
+            fm for fm in fmap_jsons if fmap_regex.sub('', remove_suffix(op.basename(fm), '.json')) == k
+        ]
+    return fmap_groups
+
+
+def get_key_info_for_fmap_assignment(json_file, matching_parameter):
+    """
+    Gets key information needed to assign fmaps to other modalities.
+    (Note: It is the responsibility of the calling function to make sure
+    the arguments are OK)
+
+    Parameters:
+    ----------
+    json_file : str or os.path
+        path to the json file
+    matching_parameter : str in AllowedFmapParameterMatching
+        matching_parameter that will be used to match runs
+
+    Returns:
+    -------
+    key_info : dict
+        part of the json file that will need to match between the fmap and
+        the other image
+    """
+    if not op.exists(json_file):
+        raise FileNotFoundError(
+            errno.ENOENT, os.strerror(errno.ENOENT), json_file
+        )
+
+    # loop through the possible criteria and extract the info needed
+    if matching_parameter == 'Shims':
+        key_info = [get_shim_setting(json_file)]
+    elif matching_parameter == 'ImagingVolume':
+        from nibabel import load as nb_load
+        nifti_file = glob(remove_suffix(json_file, '.json') + '.nii*')
+        assert len(nifti_file) == 1
+        nifti_file = nifti_file[0] 
+        nifti_header = nb_load(nifti_file).header
+        key_info = [nifti_header.get_best_affine(), nifti_header.get_data_shape()[:3]]
+    elif matching_parameter == 'ModalityAcquisitionLabel':
+        # Check the acq label for the fmap and the modality for others:
+        modality = op.basename(op.dirname(json_file))
+        if modality == 'fmap':
+            # extract the <acq> entity:
+            acq_label = BIDSFile.parse(op.basename(json_file))['acq']
+            if any(s in acq_label.lower() for s in ['fmri', 'bold', 'func']):
+                key_info = ['func']
+            elif any(s in acq_label.lower() for s in ['diff', 'dwi']):
+                key_info = ['dwi']
+            elif any(s in acq_label.lower() for s in ['anat', 'struct']):
+                key_info = ['anat']
+        else:
+            key_info = [modality]
+    elif matching_parameter == 'CustomAcquisitionLabel':
+        modality = op.basename(op.dirname(json_file))
+        if modality  == 'func':
+            # extract the <task> entity:
+            custom_label = BIDSFile.parse(op.basename(json_file))['task']
+        else:
+            # extract the <acq> entity:
+            custom_label = BIDSFile.parse(op.basename(json_file))['acq']
+        # Get the custom acquisition label, acq_label is None if no custom field found
+        key_info = [custom_label]
+    elif matching_parameter == 'Force':
+        # We want to force the matching, so just return some string
+        # regardless of the image
+        key_info = [KeyInfoForForce]
+    else:
+        # fallback:
+        key_info = []
+
+    return key_info
+
+def find_compatible_fmaps_for_run(json_file, fmap_groups, matching_parameters):
+    """
+    Finds compatible fmaps for a given run, for populate_intended_for.
+    (Note: It is the responsibility of the calling function to make sure
+    the arguments are OK)
+
+    Parameters:
+    ----------
+    json_file : str or os.path
+        path to the json file
+    fmap_groups : dict
+        key: prefix common to the group
+        value: list of all fmap paths in the group
+    matching_parameters : list of str from AllowedFmapParameterMatching
+        matching_parameters that will be used to match runs
+
+    Returns:
+    -------
+    compatible_fmap_groups : dict
+        Subset of the fmap_groups which match json_file, according
+        to the matching_parameters.
+        key: prefix common to the group
+        value: list of all fmap paths in the group
+    """
+    lgr.debug('Looking for fmaps for %s', json_file)
+    json_info = {}
+    for param in matching_parameters:
+        json_info[param] = get_key_info_for_fmap_assignment(json_file, param)
+
+    compatible_fmap_groups = {}
+    for fm_key, fm_group in fmap_groups.items():
+        # check the key_info (for all parameters) for one (the first) of
+        # the fmaps in the group:
+        compatible = False
+        for param in matching_parameters:
+            json_info_1st_item = json_info[param][0]
+            fm_info = get_key_info_for_fmap_assignment(fm_group[0], param)
+            # for the case in which key_info is a list of strings:
+            if isinstance(json_info_1st_item, str):
+                compatible = json_info[param] == fm_info
+            # for the case when no key info was found (e.g. "acq" field does not exist)
+            elif json_info_1st_item is None:
+                compatible = False
+            else:
+                # allow for tiny differences between the affines etc
+                compatible = all(np.allclose(x, y) for x, y in zip(json_info[param], fm_info))
+            if not compatible:
+                continue     # don't bother checking more params
+        if compatible:
+            compatible_fmap_groups[fm_key] = fm_group
+
+    return compatible_fmap_groups
+
+
+def find_compatible_fmaps_for_session(path_to_bids_session, matching_parameters):
+    """
+    Finds compatible fmaps for all non-fmap runs in a session.
+    (Note: It is the responsibility of the calling function to make sure
+    the arguments are OK)
+
+    Parameters:
+    ----------
+    path_to_bids_session : str or os.path
+        path to the session folder (or to the subject folder, if there are no
+        sessions).
+    matching_parameters : list of str from AllowedFmapParameterMatching
+        matching_parameters that will be used to match runs
+
+    Returns:
+    -------
+    compatible_fmap : dict
+        Dict of compatible_fmaps_groups (values) for each non-fmap run (keys)
+    """
+    lgr.debug('Looking for fmaps for session: %s', path_to_bids_session)
+
+    # Resolve path (eliminate '..')
+    path_to_bids_session = op.abspath(path_to_bids_session)
+
+    # find the different groups of fmaps:
+    fmap_dir = op.join(path_to_bids_session, 'fmap')
+    if not op.exists(fmap_dir):
+        lgr.warning('We cannot add the IntendedFor field: no fmap/ in %s', path_to_bids_session)
+        return
+    fmap_groups = find_fmap_groups(fmap_dir)
+
+    # Get a set with all non-fmap json files in the session (exclude SBRef files).
+    session_jsons = [
+        j for j in glob(op.join(path_to_bids_session, '*/*.json')) if not (
+            op.basename(op.dirname(j)) == 'fmap'
+            or remove_suffix(j, '.json').endswith('_sbref')
+        )
+    ]
+
+    # Loop through session_jsons and find the compatible fmap_groups for each
+    compatible_fmaps = {
+        j: find_compatible_fmaps_for_run(j, fmap_groups, matching_parameters)
+        for j in session_jsons
+    }
+    return compatible_fmaps
+
+
+def select_fmap_from_compatible_groups(json_file, compatible_fmap_groups, criterion):
+    """
+    Selects the fmap that will be used to correct for distortions in json_file
+    from the compatible fmap_groups list, based on the given criterion
+    (Note: It is the responsibility of the calling function to make sure
+    the arguments are OK)
+
+    Parameters:
+    ----------
+    json_file : str or os.path
+        path to the json file
+    compatible_fmap_groups : dict
+        fmap_groups that are compatible with the specific json_file
+    criterion : str in ['First', 'Closest']
+        matching_parameters that will be used to decide which fmap to use
+
+    Returns:
+    -------
+    selected_fmap_key : str or os.path
+        key from the compatible_fmap_groups for the selected fmap group
+    """
+    if len(compatible_fmap_groups) == 0:
+        return None
+    # if compatible_fmap_groups has only one entry, that's it:
+    elif len(compatible_fmap_groups) == 1:
+        return list(compatible_fmap_groups.keys())[0]
+
+    # get the modality folders, then session folder:
+    modality_folders = set(
+        op.dirname(fmap) for v in compatible_fmap_groups.values() for fmap in v
+    )      # there should be only one value, ending in 'fmap'
+    sess_folder = set(op.dirname(k) for k in modality_folders)
+    if len(sess_folder) > 1:
+        # for now, we only deal with single sessions:
+        raise RuntimeError
+    # if we made it here, we have only one session:
+    sess_folder = list(sess_folder)[0]
+
+    # get acquisition times from '_scans.tsv':
+    try:
+        scans_tsv = glob(op.join(sess_folder,'*_scans.tsv'))[0]
+    except IndexError:
+        raise FileNotFoundError("No '*_scans' file found for session %s" % sess_folder)
+    with open(scans_tsv) as f:
+        # read the contents, splitting by lines and by tab separators:
+        scans_tsv_content = [line.split('\t') for line in f.read().splitlines()]
+    # get column indices for filename and acq_time from the first line:
+    (fname_idx, time_idx) = (scans_tsv_content[0].index(k) for k in ['filename', 'acq_time'])
+    acq_times = {line[fname_idx]: line[time_idx] for line in scans_tsv_content[1:]}
+    # acq_times for the compatible fmaps:
+    acq_times_fmaps = {
+        k: acq_times[
+            # remove session folder and '.json', add '.nii.gz':
+            remove_suffix(remove_prefix(v[0], sess_folder + op.sep), '.json') + '.nii.gz'
+            ]
+        for k, v in compatible_fmap_groups.items()
+    }
+
+    if criterion == 'First':
+        # find the first acquired fmap_group from the compatible_fmap_groups:
+        first_acq_time = sorted(acq_times_fmaps.values())[0]
+        selected_fmap_key = [
+            k for k, v in acq_times_fmaps.items() if v == first_acq_time
+        ][0]
+    elif criterion == 'Closest':
+        json_acq_time = datetime.strptime(
+            acq_times[
+                # remove session folder and '.json', add '.nii.gz':
+                remove_suffix(remove_prefix(json_file, sess_folder + op.sep), '.json') + '.nii.gz'
+            ],
+            "%Y-%m-%dT%H:%M:%S.%f"
+        )
+        # differences in acquisition time (abs value):
+        diff_fmaps_acq_times = {
+            k: abs(datetime.strptime(v, "%Y-%m-%dT%H:%M:%S.%f")-json_acq_time)
+            for k, v in acq_times_fmaps.items()
+        }
+        min_diff_acq_times = sorted(diff_fmaps_acq_times.values())[0]
+        selected_fmap_key = [
+            k for k, v in diff_fmaps_acq_times.items() if v == min_diff_acq_times
+        ][0]
+
+    return selected_fmap_key
+
+
+def populate_intended_for(path_to_bids_session, matching_parameters, criterion):
+    """
+    Adds the 'IntendedFor' field to the fmap .json files in a session folder.
+    It goes through the session folders and for every json file, it finds
+    compatible_fmaps: fmaps that have the same matching_parameters as the json
+    file (e.g., same 'Shims').
+
+    If there are more than one compatible_fmaps, it will use the criterion
+    specified by the user (default: 'Closest' in time).
+
+    Because fmaps come in groups (with reversed PE polarity, or magnitude/
+    phase), we work with fmap_groups.
+
+    Parameters:
+    ----------
+    path_to_bids_session : str or os.path
+        path to the session folder (or to the subject folder, if there are no
+        sessions).
+    matching_parameters : list of str from AllowedFmapParameterMatching
+        matching_parameters that will be used to match runs
+    criterion : str in ['First', 'Closest']
+        matching_parameters that will be used to decide which of the matching
+        fmaps to use
+    """
+
+    if not isinstance(matching_parameters, list):
+        assert isinstance(matching_parameters, str), "matching_parameters must be a str or a list, got %s" % matching_parameters
+        matching_parameters = [matching_parameters]
+    for param in matching_parameters:
+        if param not in AllowedFmapParameterMatching:
+            raise ValueError(
+                "Fmap matching_parameter %s not allowed." % param
+            )
+    if criterion not in AllowedCriteriaForFmapAssignment:
+        raise ValueError(
+            "Fmap assignment criterion '%s' not allowed." % criterion
+        )
+
+    lgr.info('Adding "IntendedFor" to the fieldmaps in %s.', path_to_bids_session)
+
+    # Resolve path (eliminate '..')
+    path_to_bids_session = op.abspath(path_to_bids_session)
+
+    # Get the subject folder (if "path_to_bids_session" includes the session,
+    # remove it). "IntendedFor" paths will be relative to it.
+    if op.basename(path_to_bids_session).startswith('ses-'):
+        subj_folder = op.dirname(path_to_bids_session)
+    else:
+        subj_folder = path_to_bids_session
+
+    fmap_dir = op.join(path_to_bids_session, 'fmap')
+    if not op.exists(fmap_dir):
+        lgr.warning('We cannot add the IntendedFor field: no fmap/ in %s', path_to_bids_session)
+        return
+
+    compatible_fmaps = find_compatible_fmaps_for_session(
+        path_to_bids_session,
+        matching_parameters=matching_parameters
+    )
+    selected_fmaps = {}
+    for json_file, fmap_groups in compatible_fmaps.items():
+        if not op.dirname(json_file).endswith('fmap'):
+            selected_fmaps[json_file] = select_fmap_from_compatible_groups(
+                json_file,
+                fmap_groups,
+                criterion=criterion
+            )
+
+    # Loop through all the unique fmap_groups in compatible_fmaps:
+    unique_fmap_groups = {}
+    for cf in compatible_fmaps.values():
+        for (key, values) in cf.items():
+            if key not in unique_fmap_groups:
+                unique_fmap_groups[key] = values
+
+    for fmap_group in unique_fmap_groups:
+        intended_for = []
+        for json_file, selected_fmap_group in selected_fmaps.items():
+            if selected_fmap_group and (fmap_group in selected_fmap_group):
+                intended_for.append(
+                    op.relpath(remove_suffix(json_file, '.json') + '.nii.gz', start=subj_folder)
+                )
+        if intended_for:
+            intended_for = sorted(str(f) for f in intended_for)
+            # Add this intended_for to all fmap files in the fmap_group:
+            for fm_json in unique_fmap_groups[fmap_group]:
+                update_json(fm_json, {"IntendedFor": intended_for}, pretty=True)
+
+
+class BIDSFile(object):
+    """ as defined in https://bids-specification.readthedocs.io/en/stable/99-appendices/04-entity-table.html
+    which might soon become machine readable
+    order matters
+    """
+
+    _known_entities = ['sub', 'ses', 'task', 'acq', 'ce', 'rec', 'dir', 'run', 'mod',
+                       'echo', 'flip', 'inv', 'mt', 'part', 'recording',
+                       ]
+
+    def __init__(self, entities, suffix, extension):
+        self._entities = entities
+        self._suffix = suffix
+        self._extension = extension
+
+    def __eq__(self, other):
+        if not isinstance(other, self.__class__):
+            return False
+        if (
+            all([other[k] == v for k, v in self._entities.items()])
+            and self.extension == other.extension
+            and self.suffix == other.suffix
+        ):
+            return True
+        else:
+            return False
+
+    @classmethod
+    def parse(cls, filename):
+        """ Parse the filename for BIDS entities, suffix and extension """
+        # use re.findall to find all lower-case-letters + '-' + alphanumeric + '_' pairs:
+        entities_list = re.findall('([a-z]+)-([a-zA-Z0-9]+)[_]*', filename)
+        # keep only those in the _known_entities list:
+        entities = {k: v for k, v in entities_list if k in BIDSFile._known_entities}
+        # get whatever comes after the last key-value pair, and remove any '_' that
+        # might come in front:
+        ending = filename.split('-'.join(entities_list[-1]))[-1]
+        ending = remove_prefix(ending, '_')
+        # the first dot ('.') separates the suffix from the extension:
+        if '.' in ending:
+            suffix, extension = ending.split('.', 1)
+        else:
+            suffix, extension = ending, None
+        return BIDSFile(entities, suffix, extension)
+
+    def __str__(self):
+        """ reconstitute in a legit BIDS filename using the order from entity table """
+        if 'sub' not in self._entities:
+            raise ValueError('The \'sub-\' entity is mandatory')
+        # reconstitute the ending for the filename:
+        suffix = '_' + self.suffix if self.suffix else ''
+        extension = '.' + self.extension if self.extension else ''
+        return '_'.join(
+            ['-'.join([e, self._entities[e]]) for e in self._known_entities if e in self._entities]
+        ) + suffix + extension
+
+    def __getitem__(self, entity):
+        return self._entities[entity] if entity in self._entities else None
+
+    def __setitem__(self, entity, value):  # would puke with some exception if already known
+        return self.set(entity, value, overwrite=False)
+
+    def set(self, entity, value, overwrite=True):
+        if entity not in self._entities:
+            # just set it; no complains here
+            self._entities[entity] = value
+        elif overwrite:
+            lgr.warning("Overwriting the entity %s from %s to %s for file %s",
+                        str(entity),
+                        str(self[entity]),
+                        str(value),
+                        self.__str__()
+                        )
+            self._entities[entity] = value
+        else:
+            # if it already exists, and overwrite is false:
+            lgr.warning("Setting the entity %s to %s for file %s failed",
+                        str(entity),
+                        str(value),
+                        self.__str__()
+                        )
+
+    @property  # as needed make them RW
+    def suffix(self):
+        return self._suffix
+
+    @property
+    def extension(self):
+        return self._extension
diff --git a/heudiconv/cli/run.py b/heudiconv/cli/run.py
old mode 100755
new mode 100644
index 8a381829..42bc85d0
--- a/heudiconv/cli/run.py
+++ b/heudiconv/cli/run.py
@@ -136,6 +136,7 @@ def get_parser():
             'heuristics', 'heuristic-info',
             'ls', 'populate-templates',
             'sanitize-jsons', 'treat-jsons',
+            'populate-intended-for'
         ),
         help='Custom action to be performed on provided files instead of '
              'regular operation.')
diff --git a/heudiconv/convert.py b/heudiconv/convert.py
index 2cbdf055..bf9e6def 100644
--- a/heudiconv/convert.py
+++ b/heudiconv/convert.py
@@ -26,6 +26,7 @@
 from .bids import (
     convert_sid_bids,
     populate_bids_templates,
+    populate_intended_for,
     save_scans_key,
     tuneup_bids_json_files,
     add_participant_record,
@@ -207,6 +208,7 @@ def prep_conversion(sid, dicoms, outdir, heuristic, converter, anon_sid,
                 converter=converter,
                 scaninfo_suffix=getattr(heuristic, 'scaninfo_suffix', '.json'),
                 custom_callable=getattr(heuristic, 'custom_callable', None),
+                populate_intended_for_opts=getattr(heuristic, 'POPULATE_INTENDED_FOR_OPTS', None),
                 with_prov=with_prov,
                 bids_options=bids_options,
                 outdir=tdir,
@@ -448,7 +450,7 @@ def update_uncombined_name(metadata, filename, channel_names):
 
 def convert(items, converter, scaninfo_suffix, custom_callable, with_prov,
             bids_options, outdir, min_meta, overwrite, symlink=True, prov_file=None,
-            dcmconfig=None):
+            dcmconfig=None, populate_intended_for_opts={}):
     """Perform actual conversion (calls to converter etc) given info from
     heuristic's `infotodict`
 
@@ -575,6 +577,30 @@ def convert(items, converter, scaninfo_suffix, custom_callable, with_prov,
         if custom_callable is not None:
             custom_callable(*item)
 
+    # Populate "IntendedFor" for fmap files if requested in heuristic
+    if populate_intended_for_opts is not None:
+        # Because fmap files can only be used to correct for distortions in images
+        # collected within the same scanning session, find unique subject/session
+        # combinations from the outname in each item:
+        outnames = [item[0] for item in items]
+        # - grab "sub-<sID>[/ses-<ses>]", and keep only unique ones:
+        try:
+            sessions = set(
+                re.search(
+                    'sub-(?P<subj>[a-zA-Z0-9]*)([{0}_]ses-(?P<ses>[a-zA-Z0-9]*))?'.format(op.sep),
+                    oname
+                ).group(0)
+                for oname in outnames
+            )
+        except AttributeError:
+            # "sub-<sID>[/ses-<ses>]" is not present, so this is not BIDS compliant
+            # and it doesn't make sense to add "IntendedFor":
+            sessions = set()
+
+        for ses in sessions:
+            session_path = op.join(outdir, ses)
+            populate_intended_for(session_path, **populate_intended_for_opts)
+
 
 def convert_dicom(item_dicoms, bids_options, prefix,
                   outdir, tempdirs, symlink, overwrite):
diff --git a/heudiconv/heuristics/example.py b/heudiconv/heuristics/example.py
index 5c2bc5ca..dd152caf 100644
--- a/heudiconv/heuristics/example.py
+++ b/heudiconv/heuristics/example.py
@@ -1,6 +1,15 @@
 import os
 
 
+# Dictionary to specify options for the `populate_intended_for`.
+# Valid options are defined in 'bids.py' (for 'matching_parameters':
+# ['Shims', 'ImagingVolume',]; for 'criterion': ['First', 'Closest']
+POPULATE_INTENDED_FOR_OPTS = {
+    'matching_parameters': 'ImagingVolume',
+    'criterion': 'Closest'
+}
+
+
 def create_key(template, outtype=('nii.gz',), annotation_classes=None):
     if template is None or not template:
         raise ValueError('Template must be a valid format string')
diff --git a/heudiconv/main.py b/heudiconv/main.py
index 8b65f45e..5c6f37d5 100644
--- a/heudiconv/main.py
+++ b/heudiconv/main.py
@@ -1,9 +1,10 @@
 import logging
 import os.path as op
 import sys
+from glob import glob
 
 from . import __version__, __packagename__
-from .bids import populate_bids_templates, tuneup_bids_json_files
+from .bids import populate_bids_templates, tuneup_bids_json_files, populate_intended_for
 from .convert import prep_conversion
 from .due import due, Doi
 from .parser import get_study_sessions
@@ -47,13 +48,13 @@ def process_extra_commands(outdir, command, files, dicom_dir_template,
                            heuristic, session, subjs, grouping):
     """
     Perform custom command instead of regular operations. Supported commands:
-    ['treat-json', 'ls', 'populate-templates']
+    ['treat-json', 'ls', 'populate-templates', 'populate-intended-for']
 
     Parameters
     ----------
     outdir : str
         Output directory
-    command : {'treat-json', 'ls', 'populate-templates'}
+    command : {'treat-json', 'ls', 'populate-templates', 'populate-intended-for'}
         Heudiconv command to run
     files : list of str
         List of files
@@ -108,6 +109,47 @@ def process_extra_commands(outdir, command, files, dicom_dir_template,
         ensure_heuristic_arg(heuristic)
         from .utils import get_heuristic_description
         print(get_heuristic_description(heuristic, full=True))
+    elif command == 'populate-intended-for':
+        kwargs = {}
+        if heuristic:
+            heuristic = load_heuristic(heuristic)
+            kwargs = getattr(heuristic, 'POPULATE_INTENDED_FOR_OPTS', {})
+        if not subjs:
+            subjs = [
+                # search outdir for 'sub-*'; if it is a directory (not a regular file), remove
+                # the initial 'sub-':
+                op.basename(s)[len('sub-'):] for s in glob(op.join(outdir, 'sub-*')) if op.isdir(s)
+            ]
+            # read the subjects from the participants.tsv file to compare:
+            participants_tsv = op.join(outdir, 'participants.tsv')
+            if op.lexists(participants_tsv):
+                with open(participants_tsv, 'r') as f:
+                    # read header line and find index for 'participant_id':
+                    participant_id_index = f.readline().split('\t').index('participant_id')
+                    # read all participants, removing the initial 'sub-':
+                    known_subjects = [
+                        l.split('\t')[participant_id_index][len('sub-'):] for l in f.readlines()
+                    ]
+                if not set(subjs) == set(known_subjects):
+                    # issue a warning, but continue with the 'subjs' list (the subjects for
+                    # which there is data):
+                    lgr.warning(
+                        "'participants.tsv' contents are not identical to subjects found "
+                        "in the BIDS dataset %s", outdir
+                    )
+
+        for subj in subjs:
+            subject_path = op.join(outdir, 'sub-' + subj)
+            if session:
+                session_paths = [op.join(subject_path, 'ses-' + session)]
+            else:
+                # check to see if the data for this subject is organized by sessions; if not
+                # just use the subject_path
+                session_paths = [
+                    s for s in glob(op.join(subject_path, 'ses-*')) if op.isdir(s)
+                ] or [subject_path]
+            for session_path in session_paths:
+                populate_intended_for(session_path, **kwargs)
     else:
         raise ValueError("Unknown command %s" % command)
     return
@@ -203,7 +245,7 @@ def workflow(*, dicom_dir_template=None, files=None, subjs=None,
     debug : bool, optional
         Do not catch exceptions and show exception traceback. Default is False.
     command : {'heuristics', 'heuristic-info', 'ls', 'populate-templates',
-               'sanitize-jsons', 'treat-jsons', None}, optional
+               'sanitize-jsons', 'treat-jsons', 'populate-intended-for', None}, optional
         Custom action to be performed on provided files instead of regular
         operation. Default is None.
     grouping : {'studyUID', 'accession_number', 'all', 'custom'}, optional
diff --git a/heudiconv/tests/data/sample_nifti.nii.gz b/heudiconv/tests/data/sample_nifti.nii.gz
new file mode 100644
index 00000000..53b90ba3
Binary files /dev/null and b/heudiconv/tests/data/sample_nifti.nii.gz differ
diff --git a/heudiconv/tests/data/sample_nifti_params.txt b/heudiconv/tests/data/sample_nifti_params.txt
new file mode 100644
index 00000000..251c704d
--- /dev/null
+++ b/heudiconv/tests/data/sample_nifti_params.txt
@@ -0,0 +1,5 @@
+-1.99887669 0.0670080110 0.00144872535 109.236588
+0.0650372952 1.94962955 -0.441264838 -71.7261658
+0.0161963794 0.440969884 1.95071352 -64.1231308
+ 0 0 0 1
+ 5 5 5 1
diff --git a/heudiconv/tests/test_bids.py b/heudiconv/tests/test_bids.py
index a980e359..8716bb59 100644
--- a/heudiconv/tests/test_bids.py
+++ b/heudiconv/tests/test_bids.py
@@ -1,11 +1,63 @@
 """Test functions in heudiconv.bids module.
 """
 
+import re
+import os
+import os.path as op
+from pathlib import Path
+from random import (random,
+                    shuffle,
+                    choice,
+                    seed
+                    )
+from datetime import (datetime,
+                      timedelta,
+                      )
+from collections import (OrderedDict,
+                         )
+from glob import glob
+
+import nibabel
+import string
+from numpy import testing as np_testing
+
+from heudiconv.utils import (
+    load_json,
+    save_json,
+    create_tree,
+    remove_suffix,
+)
+
 from heudiconv.bids import (
     maybe_na,
     treat_age,
+    find_fmap_groups,
+    populate_intended_for,
+    get_shim_setting,
+    get_key_info_for_fmap_assignment,
+    find_compatible_fmaps_for_run,
+    find_compatible_fmaps_for_session,
+    select_fmap_from_compatible_groups,
+    SHIM_KEY,
+    AllowedCriteriaForFmapAssignment,
+    KeyInfoForForce,
+    BIDSFile,
 )
 
+from .utils import (
+    TESTS_DATA_PATH,
+)
+
+import pytest
+
+def gen_rand_label(label_size, label_seed, seed_stdout=True):
+    seed(label_seed)
+    rand_char = ''.join(choice(string.ascii_letters) for _ in range(label_size-1))
+    seed(label_seed)
+    rand_num = choice(string.digits)
+    if seed_stdout:
+        print(f'Seed used to generate custom label: {label_seed}')
+    return rand_char + rand_num
 
 def test_maybe_na():
     for na in '', ' ', None, 'n/a', 'N/A', 'NA':
@@ -23,4 +75,1044 @@ def test_treat_age():
     assert treat_age('1M') == '0.08'
     assert treat_age('12M') == '1'
     assert treat_age('0000.1') == '0.1'
-    assert treat_age(0000.1) == '0.1'
\ No newline at end of file
+    assert treat_age(0000.1) == '0.1'
+
+
+SHIM_LENGTH = 6
+TODAY = datetime.today()
+LABEL_SEED = int.from_bytes(os.urandom(8), byteorder="big")
+
+A_SHIM = [random() for i in range(SHIM_LENGTH)]
+def test_get_shim_setting(tmpdir):
+    """ Tests for get_shim_setting """
+    json_dir = op.join(str(tmpdir), 'foo')
+    if not op.exists(json_dir):
+        os.makedirs(json_dir)
+    json_name = op.join(json_dir, 'sub-foo.json')
+    # 1) file with no "ShimSetting", should return None
+    save_json(json_name, {})
+    with pytest.raises(KeyError):
+        assert get_shim_setting(json_name)
+
+    # -file with "ShimSetting" field
+    save_json(json_name, {SHIM_KEY: A_SHIM})
+    assert get_shim_setting(json_name) == A_SHIM
+
+
+def test_get_key_info_for_fmap_assignment(tmpdir, label_size=4, label_seed=LABEL_SEED):
+    """
+    Test get_key_info_for_fmap_assignment.
+    
+    label_size and label_seed are used for the "CustomAcquisitionLabel" matching
+    parameter. label_size is the size of the random label while label_seed is 
+    the seed for the random label creation.
+    """
+
+    nifti_file = op.join(TESTS_DATA_PATH, 'sample_nifti.nii.gz')
+    # Get the expected parameters from the NIfTI header:
+    MY_HEADER = nibabel.ni1.np.loadtxt(
+        op.join(TESTS_DATA_PATH, remove_suffix(nifti_file, '.nii.gz') + '_params.txt')
+    )
+    json_name = op.join(TESTS_DATA_PATH, remove_suffix(nifti_file, '.nii.gz') + '.json')
+
+    # 1) Call for a non-existing file should give an error:
+    with pytest.raises(FileNotFoundError):
+        assert get_key_info_for_fmap_assignment('foo.json', 'ImagingVolume')
+
+    # 2) matching_parameters = 'Shims'
+    json_name = op.join(TESTS_DATA_PATH, remove_suffix(nifti_file, '.nii.gz') + '.json')
+    save_json(json_name, {SHIM_KEY: A_SHIM})      # otherwise get_key_info_for_fmap_assignment will give an error
+    key_info = get_key_info_for_fmap_assignment(
+        json_name, matching_parameter='Shims'
+    )
+    assert key_info == [A_SHIM]
+
+    # 3) matching_parameters = 'ImagingVolume'
+    key_info = get_key_info_for_fmap_assignment(
+        json_name, matching_parameter='ImagingVolume'
+    )
+    np_testing.assert_almost_equal(key_info[0], MY_HEADER[:4], decimal=6)
+    np_testing.assert_almost_equal(key_info[1], MY_HEADER[4][:3], decimal=6)
+
+    # 4) matching_parameters = 'Force'
+    key_info = get_key_info_for_fmap_assignment(
+        json_name, matching_parameter='Force'
+    )
+    assert key_info == [KeyInfoForForce]
+
+    # 5) matching_parameter = 'ModalityAcquisitionLabel'
+    for d in ['fmap', 'func', 'dwi', 'anat']:
+        Path(op.join(str(tmpdir), d)).mkdir(parents=True, exist_ok=True)
+    for (dirname, fname, expected_key_info) in [
+        ('fmap', 'sub-foo_acq-fmri_epi.json', 'func'),
+        ('fmap', 'sub-foo_acq-bold_epi.json', 'func'),
+        ('fmap', 'sub-foo_acq-func_epi.json', 'func'),
+        ('fmap', 'sub-foo_acq-diff_epi.json', 'dwi'),
+        ('fmap', 'sub-foo_acq-anat_epi.json', 'anat'),
+        ('fmap', 'sub-foo_acq-struct_epi.json', 'anat'),
+        ('func', 'sub-foo_bold.json', 'func'),
+        ('dwi', 'sub-foo_dwi.json', 'dwi'),
+        ('anat', 'sub-foo_T1w.json', 'anat'),
+    ]:
+        json_name = op.join(str(tmpdir), dirname, fname)
+        save_json(json_name, {SHIM_KEY: A_SHIM})
+        assert [expected_key_info] == get_key_info_for_fmap_assignment(
+            json_name, matching_parameter='ModalityAcquisitionLabel'
+        )
+
+    # 6) matching_parameter = 'CustomAcquisitionLabel'
+    A_LABEL = gen_rand_label(label_size, label_seed)
+    for d in ['fmap', 'func', 'dwi', 'anat']:
+        Path(op.join(str(tmpdir), d)).mkdir(parents=True, exist_ok=True)
+        
+    for (dirname, fname, expected_key_info) in [
+        ('fmap', f'sub-foo_acq-{A_LABEL}_epi.json', A_LABEL),
+        ('func', f'sub-foo_task-{A_LABEL}_acq-foo_bold.json', A_LABEL),
+        ('dwi', f'sub-foo_acq-{A_LABEL}_dwi.json', A_LABEL),
+        ('anat', f'sub-foo_acq-{A_LABEL}_T1w.json', A_LABEL),
+    ]:
+        json_name = op.join(str(tmpdir), dirname, fname)
+        save_json(json_name, {SHIM_KEY: A_SHIM})
+        assert [expected_key_info] == get_key_info_for_fmap_assignment(
+            json_name, matching_parameter='CustomAcquisitionLabel'
+        )
+
+    # Finally: invalid matching_parameters:
+    assert get_key_info_for_fmap_assignment(
+            json_name, matching_parameter='Invalid'
+        ) == []
+
+
+def generate_scans_tsv(session_struct):
+    """
+    Generates the contents of the "_scans.tsv" file, given a session structure.
+    Currently, it will have the columns "filename" and "acq_time".
+    The acq_time will increase by one minute from run to run.
+
+    Parameters:
+    ----------
+    session_struct : dict
+        structure for the session, as a dict with modality: files
+
+    Returns:
+    -------
+    scans_file_content : str
+        multi-line string with the content of the file
+    """
+    # for each modality in session_struct (k), get the filenames:
+    scans_fnames = [
+        op.join(k, vk)
+            for k, v in session_struct.items()
+            for vk in sorted(v.keys())
+            if vk.endswith('.nii.gz')
+    ]
+    # for each file, increment the acq_time by one minute:
+    scans_file_content = ['filename\tacq_time'] + [
+        '%s\t%s' % (fn, (TODAY + timedelta(minutes=i)).isoformat()) for fn, i in
+        zip(scans_fnames, range(len(scans_fnames)))
+    ]
+    # convert to multiline string:
+    return "\n".join(scans_file_content)
+
+
+def create_dummy_pepolar_bids_session(session_path):
+    """
+    Creates a dummy BIDS session, with slim json files and empty nii.gz
+    The fmap files are pepolar
+    The json files have ShimSettings
+
+    Parameters:
+    ----------
+    session_path : str or os.path
+        path to the session (or subject) level folder
+
+    Returns:
+    -------
+    session_struct : dict
+        Structure of the directory that was created
+    expected_result : dict
+        dictionary with fmap names as keys and the expected "IntendedFor" as
+        values.
+    expected_fmap_groups : dict
+        dictionary with the expected fmap groups
+    expected_compatible_fmaps : dict
+        dictionary with the expected fmap groups for each non-fmap run in the
+        session
+    """
+    session_parent, session_basename = op.split(session_path.rstrip(op.sep))
+    if session_basename.startswith('ses-'):
+        prefix = op.split(session_parent)[1] + '_' + session_basename
+    else:
+        prefix = session_basename
+
+    # 1) Simulate the file structure for a session:
+
+    # Generate some random ShimSettings:
+    anat_shims = [random() for i in range(SHIM_LENGTH)]
+    dwi_shims = [random() for i in range(SHIM_LENGTH)]
+    func_shims_A = [random() for i in range(SHIM_LENGTH)]
+    func_shims_B = [random() for i in range(SHIM_LENGTH)]
+
+    # Dict with the file structure for the session:
+    # -anat:
+    anat_struct = {
+        '{p}_{m}.{e}'.format(p=prefix, m=mod, e=ext): dummy_content
+        for ext, dummy_content in zip(['nii.gz', 'json'], ['', {'ShimSetting': anat_shims}])
+        for mod in ['T1w', 'T2w']
+    }
+    # -dwi:
+    dwi_struct = {
+        '{p}_acq-A_run-{r}_dwi.nii.gz'.format(p=prefix, r=runNo): '' for runNo in [1, 2]
+    }
+    dwi_struct.update({
+        '{p}_acq-A_run-{r}_dwi.json'.format(p=prefix, r=runNo): {'ShimSetting': dwi_shims} for runNo in [1, 2]
+    })
+    # -func:
+    func_struct = {
+        '{p}_acq-{a}_bold.nii.gz'.format(p=prefix, a=acq): '' for acq in ['A', 'B', 'unmatched']
+    }
+    func_struct.update({
+        '{p}_acq-A_bold.json'.format(p=prefix): {'ShimSetting': func_shims_A},
+        '{p}_acq-B_bold.json'.format(p=prefix): {'ShimSetting': func_shims_B},
+        '{p}_acq-unmatched_bold.json'.format(p=prefix): {
+            'ShimSetting': [random() for i in range(SHIM_LENGTH)]
+        },
+    })
+    # -fmap:
+    #  * NIfTI files:
+    fmap_struct = {
+        '{p}_acq-{a}_dir-{d}_run-{r}_epi.nii.gz'.format(p=prefix, a=acq, d=d, r=r): ''
+        for acq in ['dwi', 'fMRI']
+        for d in ['AP', 'PA']
+        for r in [1, 2]
+    }
+    #  * dwi shims:
+    expected_fmap_groups = {
+        '{p}_acq-dwi_run-{r}_epi'.format(p=prefix, r=r): [
+            '{p}_acq-dwi_dir-{d}_run-{r}_epi.json'.format(
+                p=op.join(session_path, 'fmap', prefix), d=d, r=r
+            ) for d in ['AP', 'PA']
+        ]
+        for r in [1, 2]
+    }
+    fmap_struct.update({
+        '{p}_acq-dwi_dir-{d}_run-{r}_epi.json'.format(p=prefix, d=d, r=r): {'ShimSetting': dwi_shims}
+        for d in ['AP', 'PA']
+        for r in [1, 2]
+    })
+    #  * func_shims (_A and _B):
+    expected_fmap_groups.update({
+        '{p}_acq-fMRI_run-{r}_epi'.format(p=prefix, r=r): [
+            '{p}_acq-fMRI_dir-{d}_run-{r}_epi.json'.format(
+                p=op.join(session_path, 'fmap', prefix), d=d, r=r
+            ) for d in ['AP', 'PA']
+        ]
+        for r in [1, 2]
+    })
+    fmap_struct.update({
+        '{p}_acq-fMRI_dir-{d}_run-{r}_epi.json'.format(p=prefix, d=d, r=r): {'ShimSetting': shims}
+        for r, shims in {'1': func_shims_A, '2': func_shims_B}.items()
+        for d in ['AP', 'PA']
+    })
+    # structure for the full session (init the OrderedDict as a list to preserve order):
+    session_struct = OrderedDict([
+        ('fmap', fmap_struct),
+        ('anat', anat_struct),
+        ('dwi', dwi_struct),
+        ('func', func_struct),
+    ])
+    # add "_scans.tsv" file to the session_struct
+    scans_file_content = generate_scans_tsv(session_struct)
+    session_struct.update({'{p}_scans.tsv'.format(p=prefix): scans_file_content})
+
+    create_tree(session_path, session_struct)
+
+    # 2) Now, let's create a dict with the fmap groups compatible for each run
+    # -anat: empty
+    expected_compatible_fmaps = {
+        '{p}_{m}.json'.format(p=op.join(session_path, 'anat', prefix), m=mod): {}
+        for mod in ['T1w', 'T2w']
+    }
+    # -dwi: each of the runs (1, 2) is compatible with both of the dwi fmaps (1, 2):
+    expected_compatible_fmaps.update({
+        '{p}_acq-A_run-{r}_dwi.json'.format(p=op.join(session_path, 'dwi', prefix), r=runNo): {
+            key: val for key, val in expected_fmap_groups.items() if key in [
+                '{p}_acq-dwi_run-{r}_epi'.format(p=prefix, r=r) for r in [1, 2]
+            ]
+        }
+        for runNo in [1, 2]
+    })
+    # -func: acq-A is compatible w/ fmap fMRI run 1; acq-2 w/ fmap fMRI run 2
+    expected_compatible_fmaps.update({
+        '{p}_acq-{a}_bold.json'.format(p=op.join(session_path, 'func', prefix), a=acq): {
+            key: val for key, val in expected_fmap_groups.items() if key in [
+                '{p}_acq-fMRI_run-{r}_epi'.format(p=prefix, r=runNo)
+            ]
+        }
+        for runNo, acq in {'1': 'A', '2': 'B'}.items()
+    })
+    # -func (cont): acq-unmatched is empty
+    expected_compatible_fmaps.update({
+        '{p}_acq-unmatched_bold.json'.format(p=op.join(session_path, 'func', prefix)): {}
+    })
+
+    # 3) Then, let's create a dict with what we expect for the "IntendedFor":
+
+    sub_match = re.findall('(sub-([a-zA-Z0-9]*))', session_path)
+    sub_str = sub_match[0][0]
+    expected_prefix = session_path.split(sub_str)[-1].split(op.sep)[-1]
+
+    # dict, with fmap names as keys and the expected "IntendedFor" as values.
+    expected_result = {
+        '{p}_acq-dwi_dir-{d}_run-{r}_epi.json'.format(p=prefix, d=d, r=runNo):
+        intended_for
+        # (runNo=1 goes with the long list, runNo=2 goes with None):
+        for runNo, intended_for in zip(
+            [1, 2],
+            [[op.join(expected_prefix, 'dwi', '{p}_acq-A_run-{r}_dwi.nii.gz'.format(p=prefix, r=r)) for r in [1,2]],
+             None]
+        )
+        for d in ['AP', 'PA']
+    }
+    expected_result.update(
+        {
+            '{p}_acq-fMRI_dir-{d}_run-{r}_epi.json'.format(p=prefix, d=d, r=runNo):
+            [
+                op.join(expected_prefix,
+                        'func',
+                        '{p}_acq-{a}_bold.nii.gz'.format(p=prefix, a=acq))
+            ]
+            # runNo=1 goes with acq='A'; runNo=2 goes with acq='B'
+            for runNo, acq in zip([1, 2], ['A', 'B'])
+            for d in ['AP', 'PA']
+        }
+    )
+    
+    return session_struct, expected_result, expected_fmap_groups, expected_compatible_fmaps
+
+
+def create_dummy_no_shim_settings_bids_session(session_path):
+    """
+    Creates a dummy BIDS session, with slim json files and empty nii.gz
+    The fmap files are pepolar
+    The json files don't have ShimSettings
+
+    Parameters:
+    ----------
+    session_path : str or os.path
+        path to the session (or subject) level folder
+
+    Returns:
+    -------
+    session_struct : dict
+        Structure of the directory that was created
+    expected_result : dict
+        dictionary with fmap names as keys and the expected "IntendedFor" as
+        values.
+    None
+        it returns a third argument (None) to have the same signature as
+        create_dummy_pepolar_bids_session
+    """
+    session_parent, session_basename = op.split(session_path.rstrip(op.sep))
+    if session_basename.startswith('ses-'):
+        prefix = op.split(session_parent)[1] + '_' + session_basename
+    else:
+        prefix = session_basename
+
+    # 1) Simulate the file structure for a session:
+
+    # Dict with the file structure for the session.
+    # All json files will be empty.
+    # -anat:
+    anat_struct = {
+        '{p}_{m}.{e}'.format(p=prefix, m=mod, e=ext): dummy_content
+        for ext, dummy_content in zip(['nii.gz', 'json'], ['', {}])
+        for mod in ['T1w', 'T2w']
+    }
+    # -dwi:
+    dwi_struct = {
+        '{p}_acq-A_run-{r}_dwi.{e}'.format(p=prefix, r=runNo, e=ext): dummy_content
+        for ext, dummy_content in zip(['nii.gz', 'json'], ['', {}])
+        for runNo in [1, 2]
+    }
+    # -func:
+    func_struct = {
+        '{p}_acq-{a}_bold.{e}'.format(p=prefix, a=acq, e=ext): dummy_content
+        for ext, dummy_content in zip(['nii.gz', 'json'], ['', {}])
+        for acq in ['A', 'B']
+    }
+    # -fmap:
+    fmap_struct = {
+        '{p}_acq-{a}_dir-{d}_run-{r}_epi.{e}'.format(p=prefix, a=acq, d=d, r=r, e=ext): dummy_content
+        for ext, dummy_content in zip(['nii.gz', 'json'], ['', {}])
+        for acq in ['dwi', 'fMRI']
+        for d in ['AP', 'PA']
+        for r in [1, 2]
+    }
+    expected_fmap_groups = {
+        '{p}_acq-{a}_run-{r}_epi'.format(p=prefix, a=acq, r=r): [
+            '{p}_acq-{a}_dir-{d}_run-{r}_epi.json'.format(
+                p=op.join(session_path, 'fmap', prefix), a=acq, d=d, r=r
+            ) for d in ['AP', 'PA']
+        ]
+        for acq in ['dwi', 'fMRI']
+        for r in [1, 2]
+    }
+
+    # structure for the full session (init the OrderedDict as a list to preserve order):
+    session_struct = OrderedDict([
+        ('fmap', fmap_struct),
+        ('anat', anat_struct),
+        ('dwi', dwi_struct),
+        ('func', func_struct),
+    ])
+    # add "_scans.tsv" file to the session_struct
+    scans_file_content = generate_scans_tsv(session_struct)
+    session_struct.update({'{p}_scans.tsv'.format(p=prefix): scans_file_content})
+
+    create_tree(session_path, session_struct)
+
+    # 2) Now, let's create a dict with the fmap groups compatible for each run
+    # -anat: empty
+    expected_compatible_fmaps = {
+        '{p}_{m}.json'.format(p=op.join(session_path, 'anat', prefix), m=mod): {}
+        for mod in ['T1w', 'T2w']
+    }
+    # -dwi: each of the runs (1, 2) is compatible with both of the dwi fmaps (1, 2):
+    expected_compatible_fmaps.update({
+        '{p}_acq-A_run-{r}_dwi.json'.format(p=op.join(session_path, 'dwi', prefix), r=runNo): {
+            key: val for key, val in expected_fmap_groups.items() if key in [
+                '{p}_acq-dwi_run-{r}_epi'.format(p=prefix, r=r) for r in [1, 2]
+            ]
+        }
+        for runNo in [1, 2]
+    })
+    # -func: each of the acq (A, B) is compatible w/ both fmap fMRI runs (1, 2)
+    expected_compatible_fmaps.update({
+        '{p}_acq-{a}_bold.json'.format(p=op.join(session_path, 'func', prefix), a=acq): {
+            key: val for key, val in expected_fmap_groups.items() if key in [
+                '{p}_acq-fMRI_run-{r}_epi'.format(p=prefix, r=r) for r in [1, 2]
+           ]
+        }
+        for acq in ['A', 'B']
+    })
+
+    # 3) Now, let's create a dict with what we expect for the "IntendedFor":
+    # NOTE: The "expected_prefix" (the beginning of the path to the
+    # "IntendedFor") should be relative to the subject level (see:
+    # https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/01-magnetic-resonance-imaging-data.html#fieldmap-data)
+
+    sub_match = re.findall('(sub-([a-zA-Z0-9]*))', session_path)
+    sub_str = sub_match[0][0]
+    expected_prefix = session_path.split(sub_str)[-1].split(op.sep)[-1]
+
+    # dict, with fmap names as keys and the expected "IntendedFor" as values.
+    expected_result = {
+        # (runNo=1 goes with the long list, runNo=2 goes with None):
+        '{p}_acq-dwi_dir-{d}_run-{r}_epi.json'.format(p=prefix, d=d, r=runNo): intended_for
+        for runNo, intended_for in zip(
+            [1, 2],
+            [[op.join(expected_prefix, 'dwi', '{p}_acq-A_run-{r}_dwi.nii.gz'.format(p=prefix, r=r)) for r in [1,2]],
+             None]
+        )
+        for d in ['AP', 'PA']
+    }
+    expected_result.update(
+        {
+            # The first "fMRI" run gets all files in the "func" folder;
+            # the second shouldn't get any.
+            '{p}_acq-fMRI_dir-{d}_run-{r}_epi.json'.format(p=prefix, d=d, r=runNo): intended_for
+            for runNo, intended_for in zip(
+                [1, 2],
+                [[op.join(expected_prefix, 'func', '{p}_acq-{a}_bold.nii.gz'.format(p=prefix, a=acq))
+                  for acq in ['A', 'B']],
+                  None]
+            )
+            for d in ['AP', 'PA']
+        }
+    )
+
+    return session_struct, expected_result, expected_fmap_groups, expected_compatible_fmaps
+
+def create_dummy_no_shim_settings_custom_label_bids_session(session_path, label_size=4, label_seed=LABEL_SEED):
+    """
+    Creates a dummy BIDS session, with slim json files and empty nii.gz
+    The fmap files are pepolar
+    The json files don't have ShimSettings
+    The fmap files have a custom ACQ label matching:
+        - TASK label for <func> modality
+        - ACQ label for any other modality (e.g. <dwi>)
+
+    Parameters:
+    ----------
+    session_path : str or os.path
+        path to the session (or subject) level folder
+    label_size : int, optional
+        size of the random label
+    label_seed : int, optional
+        seed for the random label creation
+
+    Returns:
+    -------
+    session_struct : dict
+        Structure of the directory that was created
+    expected_result : dict
+        dictionary with fmap names as keys and the expected "IntendedFor" as
+        values.
+    None
+        it returns a third argument (None) to have the same signature as
+        create_dummy_pepolar_bids_session
+    """
+    session_parent, session_basename = op.split(session_path.rstrip(op.sep))
+    if session_basename.startswith('ses-'):
+        prefix = op.split(session_parent)[1] + '_' + session_basename
+    else:
+        prefix = session_basename
+
+    # 1) Simulate the file structure for a session:
+
+    # Dict with the file structure for the session.
+    # All json files will be empty.
+    # -anat:
+    anat_struct = {
+        f'{prefix}_{mod}.{ext}': dummy_content
+        for ext, dummy_content in zip(['nii.gz', 'json'], ['', {}])
+        for mod in ['T1w', 'T2w']
+    }
+    # -dwi:
+    label_seed += 1
+    DWI_LABEL = gen_rand_label(label_size, label_seed)
+    dwi_struct = {
+        f'{prefix}_acq-{DWI_LABEL}_run-{runNo}_dwi.{ext}': dummy_content
+        for ext, dummy_content in zip(['nii.gz', 'json'], ['', {}])
+        for runNo in [1, 2]
+    }
+    # -func:
+    label_seed += 1
+    FUNC_LABEL = gen_rand_label(label_size, label_seed)
+    func_struct = {
+        f'{prefix}_task-{FUNC_LABEL}_acq-{acq}_bold.{ext}': dummy_content
+        for ext, dummy_content in zip(['nii.gz', 'json'], ['', {}])
+        for acq in ['A', 'B']
+    }
+    # -fmap:
+    fmap_struct = {
+        f'{prefix}_acq-{acq}_dir-{d}_run-{r}_epi.{ext}': dummy_content
+        for ext, dummy_content in zip(['nii.gz', 'json'], ['', {}])
+        for acq in [DWI_LABEL, FUNC_LABEL]
+        for d in ['AP', 'PA']
+        for r in [1, 2]
+    }
+    expected_fmap_groups = {
+        f'{prefix}_acq-{acq}_run-{r}_epi': [
+            f'{op.join(session_path, "fmap", prefix)}_acq-{acq}_dir-{d}_run-{r}_epi.json'
+            for d in ['AP', 'PA']
+        ]
+        for acq in [DWI_LABEL, FUNC_LABEL]
+        for r in [1, 2]
+    }
+
+    # structure for the full session (init the OrderedDict as a list to preserve order):
+    session_struct = OrderedDict([
+        ('fmap', fmap_struct),
+        ('anat', anat_struct),
+        ('dwi', dwi_struct),
+        ('func', func_struct),
+    ])
+    # add "_scans.tsv" file to the session_struct
+    scans_file_content = generate_scans_tsv(session_struct)
+    session_struct.update({'{p}_scans.tsv'.format(p=prefix): scans_file_content})
+
+    create_tree(session_path, session_struct)
+
+    # 2) Now, let's create a dict with the fmap groups compatible for each run
+    # -anat: empty
+    expected_compatible_fmaps = {
+        f'{op.join(session_path, "anat", prefix)}_{mod}.json': {}
+        for mod in ['T1w', 'T2w']
+    }
+    # -dwi: each of the runs (1, 2) is compatible with both of the dwi fmaps (1, 2):
+    expected_compatible_fmaps.update({
+        f'{op.join(session_path, "dwi", prefix)}_acq-{DWI_LABEL}_run-{runNo}_dwi.json': {
+            key: val for key, val in expected_fmap_groups.items() if key in [
+                f'{prefix}_acq-{DWI_LABEL}_run-{r}_epi' for r in [1, 2]
+            ]
+        }
+        for runNo in [1, 2]
+    })
+    # -func: each of the acq (A, B) is compatible w/ both fmap fMRI runs (1, 2)
+    expected_compatible_fmaps.update({
+        f'{op.join(session_path, "func", prefix)}_task-{FUNC_LABEL}_acq-{acq}_bold.json': {
+            key: val for key, val in expected_fmap_groups.items() if key in [
+                f'{prefix}_acq-{FUNC_LABEL}_run-{r}_epi' for r in [1, 2]
+           ]
+        }
+        for acq in ['A', 'B']
+    })
+
+    # 3) Now, let's create a dict with what we expect for the "IntendedFor":
+    # NOTE: The "expected_prefix" (the beginning of the path to the
+    # "IntendedFor") should be relative to the subject level (see:
+    # https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/01-magnetic-resonance-imaging-data.html#fieldmap-data)
+
+    sub_match = re.findall('(sub-([a-zA-Z0-9]*))', session_path)
+    sub_str = sub_match[0][0]
+    expected_prefix = session_path.split(sub_str)[-1].split(op.sep)[-1]
+
+    # dict, with fmap names as keys and the expected "IntendedFor" as values.
+    expected_result = {
+        # (runNo=1 goes with the long list, runNo=2 goes with None):
+        f'{prefix}_acq-{DWI_LABEL}_dir-{d}_run-{runNo}_epi.json': intended_for
+        for runNo, intended_for in zip(
+            [1, 2],
+            [[op.join(expected_prefix, 'dwi', f'{prefix}_acq-{DWI_LABEL}_run-{r}_dwi.nii.gz') for r in [1,2]],
+             None]
+        )
+        for d in ['AP', 'PA']
+    }
+    expected_result.update(
+        {
+            # The first "fMRI" run gets all files in the "func" folder;
+            # the second shouldn't get any.
+            f'{prefix}_acq-{FUNC_LABEL}_dir-{d}_run-{runNo}_epi.json': intended_for
+            for runNo, intended_for in zip(
+                [1, 2],
+                [[op.join(expected_prefix, 'func', f'{prefix}_task-{FUNC_LABEL}_acq-{acq}_bold.nii.gz')
+                  for acq in ['A', 'B']],
+                  None]
+            )
+            for d in ['AP', 'PA']
+        }
+    )
+
+    return session_struct, expected_result, expected_fmap_groups, expected_compatible_fmaps
+
+def create_dummy_magnitude_phase_bids_session(session_path):
+    """
+    Creates a dummy BIDS session, with slim json files and empty nii.gz
+    The fmap files are a magnitude/phase pair
+    The json files have ShimSettings
+    We just need to test a very simple case to make sure the mag/phase have
+    the same "IntendedFor" field:
+
+    Parameters:
+    ----------
+    session_path : str or os.path
+        path to the session (or subject) level folder
+
+    Returns:
+    -------
+    session_struct : dict
+        Structure of the directory that was created
+    expected_result : dict
+        dictionary with fmap names as keys and the expected "IntendedFor" as
+        values.
+    expected_fmap_groups : dict
+        dictionary with the expected fmap groups
+    """
+    session_parent, session_basename = op.split(session_path.rstrip(op.sep))
+    if session_basename.startswith('ses-'):
+        prefix = op.split(session_parent)[1] + '_' + session_basename
+    else:
+        prefix = session_basename
+
+    # 1) Simulate the file structure for a session:
+
+    # Generate some random ShimSettings:
+    dwi_shims = [random() for i in range(SHIM_LENGTH)]
+    func_shims_A = [random() for i in range(SHIM_LENGTH)]
+    func_shims_B = [random() for i in range(SHIM_LENGTH)]
+
+    # Dict with the file structure for the session:
+    # -dwi:
+    dwi_struct = {
+        '{p}_acq-A_run-{r}_dwi.nii.gz'.format(p=prefix, r=runNo): '' for runNo in [1, 2]
+    }
+    dwi_struct.update({
+        '{p}_acq-A_run-{r}_dwi.json'.format(p=prefix, r=runNo): {'ShimSetting': dwi_shims} for runNo in [1, 2]
+    })
+    # -func:
+    func_struct = {
+        '{p}_acq-{a}_bold.nii.gz'.format(p=prefix, a=acq): '' for acq in ['A', 'B', 'unmatched']
+    }
+    func_struct.update({
+        '{p}_acq-A_bold.json'.format(p=prefix): {'ShimSetting': func_shims_A},
+        '{p}_acq-B_bold.json'.format(p=prefix): {'ShimSetting': func_shims_B},
+        '{p}_acq-unmatched_bold.json'.format(p=prefix): {
+            'ShimSetting': [random() for i in range(SHIM_LENGTH)]
+        },
+    })
+    # -fmap:
+    #    * Case 1 in https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/01-magnetic-resonance-imaging-data.html#fieldmap-data
+    fmap_struct = {
+        '{p}_acq-case1_{s}.nii.gz'.format(p=prefix, s=suffix): ''
+        for suffix in ['phasediff', 'magnitude1', 'magnitude2']
+    }
+    expected_fmap_groups = {
+        '{p}_acq-case1'.format(p=prefix): [
+            '{p}_acq-case1_phasediff.json'.format(p=op.join(session_path, 'fmap', prefix))
+        ]
+    }
+    fmap_struct.update({
+        '{p}_acq-case1_phasediff.json'.format(p=prefix): {'ShimSetting': dwi_shims}
+    })
+    #    * Case 2:
+    fmap_struct.update({
+        '{p}_acq-case2_{s}.nii.gz'.format(p=prefix, s=suffix): ''
+        for suffix in ['magnitude1', 'magnitude2', 'phase1', 'phase2']
+    })
+    expected_fmap_groups.update({
+        '{p}_acq-case2'.format(p=prefix): [
+            '{p}_acq-case2_phase{n}.json'.format(
+                p=op.join(session_path, 'fmap', prefix), n=n
+            ) for n in [1, 2]
+        ]
+    })
+    fmap_struct.update({
+        '{p}_acq-case2_phase{n}.json'.format(p=prefix, n=n): {'ShimSetting': func_shims_A}
+        for n in [1, 2]
+    })
+    #    * Case 3:
+    fmap_struct.update({
+        '{p}_acq-case3_{s}.nii.gz'.format(p=prefix, s=suffix): ''
+        for suffix in ['magnitude', 'fieldmap']
+    })
+    expected_fmap_groups.update({
+        '{p}_acq-case3'.format(p=prefix): [
+            '{p}_acq-case3_fieldmap.json'.format(p=op.join(session_path, 'fmap', prefix))
+        ]
+    })
+    fmap_struct.update({
+        '{p}_acq-case3_fieldmap.json'.format(p=prefix): {'ShimSetting': func_shims_B}
+    })
+    # structure for the full session (init the OrderedDict as a list to preserve order):
+    session_struct = OrderedDict([
+        ('fmap', fmap_struct),
+        ('dwi', dwi_struct),
+        ('func', func_struct),
+    ])
+    # add "_scans.tsv" file to the session_struct
+    scans_file_content = generate_scans_tsv(session_struct)
+    session_struct.update({'{p}_scans.tsv'.format(p=prefix): scans_file_content})
+
+    create_tree(session_path, session_struct)
+
+    # 2) Now, let's create a dict with the fmap groups compatible for each run
+    # -dwi: each of the runs (1, 2) is compatible with case1 fmap:
+    expected_compatible_fmaps = {
+        '{p}_acq-A_run-{r}_dwi.json'.format(p=op.join(session_path, 'dwi', prefix), r=runNo): {
+            key: val for key, val in expected_fmap_groups.items() if key in [
+                '{p}_acq-case1'.format(p=prefix)
+            ]
+        }
+        for runNo in [1, 2]
+    }
+    # -func: acq-A is compatible w/ fmap case2; acq-B w/ fmap case3
+    expected_compatible_fmaps.update({
+        '{p}_acq-{a}_bold.json'.format(p=op.join(session_path, 'func', prefix), a=acq): {
+            key: val for key, val in expected_fmap_groups.items() if key in [
+                '{p}_acq-case{c}'.format(p=prefix, c=caseNo)
+            ]
+        }
+        for caseNo, acq in {'2': 'A', '3': 'B'}.items()
+    })
+    # -func (cont): acq-unmatched is empty
+    expected_compatible_fmaps.update({
+        '{p}_acq-unmatched_bold.json'.format(p=op.join(session_path, 'func', prefix)): {}
+    })
+
+    # 3) Now, let's create a dict with what we expect for the "IntendedFor":
+
+    sub_match = re.findall('(sub-([a-zA-Z0-9]*))', session_path)
+    sub_str = sub_match[0][0]
+    expected_prefix = session_path.split(sub_str)[-1].split(op.sep)[-1]
+
+    # dict, with fmap names as keys and the expected "IntendedFor" as values.
+    expected_result = {
+        '{p}_acq-case1_{s}.json'.format(p=prefix, s='phasediff'):
+            [op.join(expected_prefix, 'dwi', '{p}_acq-A_run-{r}_dwi.nii.gz'.format(p=prefix, r=r)) for r in [1, 2]]
+    }
+    expected_result.update({
+        '{p}_acq-case2_phase{n}.json'.format(p=prefix, n=n):
+            # populate_intended_for writes lists:
+            [op.join(expected_prefix, 'func', '{p}_acq-A_bold.nii.gz'.format(p=prefix))]
+        for n in [1, 2]
+    })
+    expected_result.update({
+        '{p}_acq-case3_fieldmap.json'.format(p=prefix):
+            # populate_intended_for writes lists:
+            [op.join(expected_prefix, 'func', '{p}_acq-B_bold.nii.gz'.format(p=prefix))]
+    })
+
+    return session_struct, expected_result, expected_fmap_groups, expected_compatible_fmaps
+
+
+# Test cases:
+# A) pepolar fmaps with ShimSetting in json files
+# B) same, with no ShimSetting
+# C) magnitude/phase, with ShimSetting
+@pytest.mark.parametrize(
+    "simulation_function", [create_dummy_pepolar_bids_session,
+                            create_dummy_no_shim_settings_bids_session,
+                            create_dummy_magnitude_phase_bids_session]
+)
+def test_find_fmap_groups(tmpdir, simulation_function):
+    """ Test for find_fmap_groups """
+    folder = op.join(str(tmpdir), 'sub-foo')
+    _, _, expected_fmap_groups, _ = simulation_function(folder)
+    fmap_groups = find_fmap_groups(op.join(folder, 'fmap'))
+    assert fmap_groups == expected_fmap_groups
+
+
+# Test cases:
+# A) pepolar fmaps with ShimSetting in json files
+# B) same, with no ShimSetting
+# C) magnitude/phase, with ShimSetting
+@pytest.mark.parametrize(
+    "simulation_function, match_param", [
+        (create_dummy_pepolar_bids_session, 'Shims'),
+        (create_dummy_no_shim_settings_bids_session, 'ModalityAcquisitionLabel'),
+        (create_dummy_no_shim_settings_custom_label_bids_session, 'CustomAcquisitionLabel'),
+        (create_dummy_magnitude_phase_bids_session, 'Shims')
+    ]
+)
+def test_find_compatible_fmaps_for_run(tmpdir, simulation_function, match_param):
+    """
+    Test find_compatible_fmaps_for_run.
+
+    Parameters:
+    ----------
+    tmpdir
+    simulation_function : function
+        function to create the directory tree and expected results
+    match_param : str
+        matching_parameter for assigning fmaps
+    """
+    folder = op.join(str(tmpdir), 'sub-foo')
+    _, _, expected_fmap_groups, expected_compatible_fmaps = simulation_function(folder)
+    for modality in ['anat', 'dwi', 'func']:
+        for json_file in glob(op.join(folder, modality, '*.json')):
+            compatible_fmaps = find_compatible_fmaps_for_run(
+                json_file,
+                expected_fmap_groups,
+                matching_parameters=[match_param]
+            )
+            assert compatible_fmaps == expected_compatible_fmaps[json_file]
+
+
+# Test two scenarios for each case:
+# -study without sessions
+# -study with sessions
+# Cases:
+# A) pepolar fmaps with ShimSetting in json files
+# B) same, with no ShimSetting
+# C) magnitude/phase, with ShimSetting
+@pytest.mark.parametrize(
+    "folder, expected_prefix, simulation_function, match_param", [
+        (folder, expected_prefix, sim_func, mp)
+        for folder, expected_prefix in zip(['no_sessions/sub-1', 'sessions/sub-1/ses-pre'], ['', 'ses-pre'])
+        for sim_func, mp in [
+            (create_dummy_pepolar_bids_session, 'Shims'),
+            (create_dummy_no_shim_settings_bids_session, 'ModalityAcquisitionLabel'),
+            (create_dummy_no_shim_settings_custom_label_bids_session, 'CustomAcquisitionLabel'),
+            (create_dummy_magnitude_phase_bids_session, 'Shims')
+        ]
+    ]
+)
+def test_find_compatible_fmaps_for_session(
+        tmpdir,
+        folder,
+        expected_prefix,
+        simulation_function,
+        match_param
+):
+    """
+    Test find_compatible_fmaps_for_session.
+
+    Parameters:
+    ----------
+    tmpdir
+    folder : str or os.path
+        path to BIDS study to be simulated, relative to tmpdir
+    expected_prefix : str
+        expected start of the "IntendedFor" elements
+    simulation_function : function
+        function to create the directory tree and expected results
+    match_param : str
+        matching_parameter for assigning fmaps
+    """
+    session_folder = op.join(str(tmpdir), folder)
+    _, _, _, expected_compatible_fmaps = simulation_function(session_folder)
+
+    compatible_fmaps = find_compatible_fmaps_for_session(session_folder, matching_parameters=[match_param])
+
+    assert compatible_fmaps == expected_compatible_fmaps
+
+
+# Test two scenarios for each case:
+# -study without sessions
+# -study with sessions
+# Cases:
+# A) pepolar fmaps with ShimSetting in json files
+# B) same, with no ShimSetting
+# C) magnitude/phase, with ShimSetting
+@pytest.mark.parametrize(
+    "folder, expected_prefix, simulation_function", [
+        (folder, expected_prefix, sim_func)
+        for folder, expected_prefix in zip(['no_sessions/sub-1', 'sessions/sub-1/ses-pre'], ['', 'ses-pre'])
+        for sim_func in [create_dummy_pepolar_bids_session,
+                         create_dummy_no_shim_settings_bids_session,
+                         create_dummy_magnitude_phase_bids_session]
+    ]
+)
+def test_select_fmap_from_compatible_groups(tmpdir, folder, expected_prefix, simulation_function):
+    """Test select_fmap_from_compatible_groups"""
+    session_folder = op.join(str(tmpdir), folder)
+    _, _, _, expected_compatible_fmaps = simulation_function(session_folder)
+
+    for json_file, fmap_groups in expected_compatible_fmaps.items():
+        for criterion in AllowedCriteriaForFmapAssignment:
+            if not op.dirname(json_file).endswith('fmap'):
+                selected_fmap = select_fmap_from_compatible_groups(
+                    json_file,
+                    fmap_groups,
+                    criterion=criterion
+                )
+            # when the criterion is 'First', you should get the first of
+            # the compatible_fmaps (for that json_file), if it is 'Closest',
+            # it should be the last one (the fmaps are "run" at the
+            # beginning of the session)
+            if selected_fmap:
+                if criterion == 'First':
+                    assert selected_fmap == sorted(expected_compatible_fmaps[json_file])[0]
+                elif criterion == 'Closest':
+                    assert selected_fmap == sorted(expected_compatible_fmaps[json_file])[-1]
+            else:
+                assert not expected_compatible_fmaps[json_file]
+
+
+# Test two scenarios for each case:
+# -study without sessions
+# -study with sessions
+# Cases:
+# A) pepolar fmaps with ShimSetting in json files
+# B) same, with no ShimSetting
+# C) magnitude/phase, with ShimSetting
+@pytest.mark.parametrize(
+    "folder, expected_prefix, simulation_function, match_param", [
+        (folder, expected_prefix, sim_func, mp)
+        for folder, expected_prefix in zip(['no_sessions/sub-1', 'sessions/sub-1/ses-pre'], ['', 'ses-pre'])
+        for sim_func, mp in [
+            (create_dummy_pepolar_bids_session, 'Shims'),
+            (create_dummy_no_shim_settings_bids_session, 'ModalityAcquisitionLabel'),
+            (create_dummy_no_shim_settings_custom_label_bids_session, 'CustomAcquisitionLabel'),
+            (create_dummy_magnitude_phase_bids_session, 'Shims')
+        ]
+    ]
+)
+def test_populate_intended_for(
+    tmpdir,
+    folder,
+    expected_prefix,
+    simulation_function,
+    match_param
+):
+    """
+    Test populate_intended_for.
+    Parameters:
+    ----------
+    tmpdir
+    folder : str or os.path
+        path to BIDS study to be simulated, relative to tmpdir
+    expected_prefix : str
+        expected start of the "IntendedFor" elements
+    simulation_function : function
+        function to create the directory tree and expected results
+    match_param : str
+        matching_parameter for assigning fmaps
+    """
+
+    session_folder = op.join(str(tmpdir), folder)
+    session_struct, expected_result, _, _ = simulation_function(session_folder)
+    populate_intended_for(session_folder, matching_parameters=match_param, criterion='First')
+
+    # Now, loop through the jsons in the fmap folder and make sure it matches
+    # the expected result:
+    fmap_folder = op.join(session_folder, 'fmap')
+    for j in session_struct['fmap'].keys():
+        if j.endswith('.json'):
+            assert j in expected_result.keys()
+            data = load_json(op.join(fmap_folder, j))
+            if expected_result[j]:
+                assert data['IntendedFor'] == expected_result[j]
+                # Also, make sure the run with random shims is not here:
+                # (It is assured by the assert above, but let's make it
+                # explicit)
+                run_prefix = j.split('_acq')[0]
+                assert '{p}_acq-unmatched_bold.nii.gz'.format(p=run_prefix) not in data['IntendedFor']
+            else:
+                assert 'IntendedFor' not in data.keys()
+
+
+def test_BIDSFile():
+    """ Tests for the BIDSFile class """
+
+    # define entities in the correct order:
+    sorted_entities = [
+        ('sub', 'Jason'),
+        ('acq', 'Treadstone'),
+        ('run', '2'),
+        ('echo', '1'),
+    ]
+    # 'sub-Jason_acq-Treadstone_run-2_echo-1':
+    expected_sorted_str = '_'.join(['-'.join(e) for e in sorted_entities])
+    # expected BIDSFile:
+    suffix = 'T1w'
+    extension = 'nii.gz'
+    expected_bids_file = BIDSFile(OrderedDict(sorted_entities), suffix, extension)
+
+    # entities in random order:
+    idcs = list(range(len(sorted_entities)))
+    shuffle(idcs)
+    shuffled_entities = [sorted_entities[i] for i in idcs]
+    shuffled_str = '_'.join(['-'.join(e) for e in shuffled_entities])
+
+    # Test __eq__ method.
+    # It should consider equal BIDSFiles with the same entities even in different order:
+    assert BIDSFile(OrderedDict(shuffled_entities), suffix, extension) == expected_bids_file
+
+    # Test __getitem__:
+    assert all([expected_bids_file[k] == v for k, v in shuffled_entities])
+
+    # Test filename parser and  __str__ method:
+    # Note: the __str__ method should return entities in the correct order
+    ending = '_T1w.nii.gz'            # suffix + extension
+    my_bids_file = BIDSFile.parse(shuffled_str + ending)
+    assert my_bids_file == expected_bids_file
+    assert str(my_bids_file) == expected_sorted_str + ending
+
+    ending = '.json'                  # just extension
+    my_bids_file = BIDSFile.parse(shuffled_str + ending)
+    assert my_bids_file.suffix == ''
+    assert str(my_bids_file) == expected_sorted_str + ending
+
+    ending = '_T1w'                   # just suffix
+    my_bids_file = BIDSFile.parse(shuffled_str + ending)
+    assert my_bids_file.extension is None
+    assert str(my_bids_file) == expected_sorted_str + ending
+
+    # Complain if entity 'sub' is not set:
+    with pytest.raises(ValueError) as e_info:
+        assert str(BIDSFile.parse('dir-reversed.json'))
+        assert 'sub-' in e_info.value
+
+    # Test set method:
+    # -for a new entity, just set it without a complaint:
+    my_bids_file['dir'] = 'AP'
+    assert my_bids_file['dir'] == 'AP'
+    # -for an existing entity, don't change it by default:
+    my_bids_file['echo'] = '2'
+    assert my_bids_file['echo'] == expected_bids_file['echo']  # still the original value
+    # -for an existing entity, you can overwrite it with "set":
+    my_bids_file.set('echo', '2')
+    assert my_bids_file['echo'] == '2'
diff --git a/heudiconv/tests/test_convert.py b/heudiconv/tests/test_convert.py
index f8316ee5..c9674c22 100644
--- a/heudiconv/tests/test_convert.py
+++ b/heudiconv/tests/test_convert.py
@@ -4,7 +4,9 @@
 from glob import glob
 
 import pytest
+import heudiconv.convert
 from heudiconv.bids import BIDSError
+from heudiconv.utils import load_heuristic
 from heudiconv.cli.run import main as runner
 from heudiconv.convert import (
     DW_IMAGE_IN_FMAP_FOLDER_WARNING,
@@ -163,3 +165,80 @@ def test_b0dwi_for_fmap(tmpdir, caplog):
         assert op.isdir(op.join(tmppath, 'sub-%s', mod) % (subID))
         for ext in ['bval', 'bvec']:
             assert glob(op.join(tmppath, 'sub-%s', mod, 'sub-%s_*.%s') % (subID, subID, ext))
+
+
+# Test two scenarios for each case:
+# -study without sessions
+# -study with sessions
+@pytest.mark.parametrize(
+    "subjects, sesID, expected_session_folder", [
+        (['Jason', 'Bourne'], None, 'sub-{sID}'),
+        (['Bourne'], 'Treadstone', op.join('sub-{{sID}}', 'ses-{{ses}}')),
+    ]
+)
+# Two possibilities: with or without heuristics:
+@pytest.mark.parametrize(
+    "heuristic", ['example', 'reproin', None]       # heuristics/example.py, heuristics/reproin.py
+)
+def test_populate_intended_for(tmpdir, monkeypatch, capfd,
+                 subjects, sesID, expected_session_folder,
+                 heuristic):
+    """
+    Test convert
+
+    For now, I'm just going to test that the call to populate_intended_for is
+    done with the correct argument.
+    More tests can be added here.
+    """
+
+    def mock_populate_intended_for(session, matching_parameters='Shims', criterion='Closest'):
+        """
+        Pretend we run populate_intended_for, but just print out the arguments.
+        """
+        print('session: {}'.format(session))
+        print('matching_parameters: {}'.format(matching_parameters))
+        print('criterion: {}'.format(criterion))
+        return
+    # mock the "populate_intended_for":
+    monkeypatch.setattr(
+        heudiconv.convert, "populate_intended_for", mock_populate_intended_for
+    )
+
+    outdir = op.join(str(tmpdir), 'foo')
+    outfolder = op.join(outdir, 'sub-{sID}', 'ses-{ses}') if sesID else op.join(outdir,'sub-{sID}')
+    sub_ses = 'sub-{sID}' + ('_ses-{ses}' if sesID else '')
+
+    # items are a list of tuples, with each tuple having three elements:
+    #   prefix, outtypes, item_dicoms
+    items = [
+        (op.join(outfolder, 'anat', sub_ses + '_T1w').format(sID=s, ses=sesID), ('',), [])
+        for s in subjects
+    ]
+
+    heuristic = load_heuristic(heuristic) if heuristic else None
+    heudiconv.convert.convert(items,
+            converter='',
+            scaninfo_suffix='.json',
+            custom_callable=None,
+            populate_intended_for_opts=getattr(heuristic, 'POPULATE_INTENDED_FOR_OPTS', None),
+            with_prov=None,
+            bids_options=[],
+            outdir=outdir,
+            min_meta=True,
+            overwrite=False)
+    output = capfd.readouterr()
+    # if the heuristic module has a 'POPULATE_INTENDED_FOR_OPTS' field, we expect
+    # to get the output of the mock_populate_intended_for, otherwise, no output:
+    if getattr(heuristic, 'POPULATE_INTENDED_FOR_OPTS', None):
+        assert all([
+            "\n".join([
+                "session: " + outfolder.format(sID=s, ses=sesID),
+                # "ImagingVolume" is defined in heuristic file; "Shims" is the default
+                "matching_parameters: " + "ImagingVolume",
+                "criterion: Closest"
+            ]) in output.out
+            for s in subjects
+        ])
+    else:
+        # If there was no heuristic, make sure populate_intended_for was not called
+        assert not output.out
diff --git a/heudiconv/tests/test_main.py b/heudiconv/tests/test_main.py
index 83998469..17ff08ae 100644
--- a/heudiconv/tests/test_main.py
+++ b/heudiconv/tests/test_main.py
@@ -1,7 +1,8 @@
 # TODO: break this up by modules
 
 from heudiconv.cli.run import main as runner
-from heudiconv.main import workflow
+from heudiconv.main import (workflow,
+                            process_extra_commands)
 from heudiconv import __version__
 from heudiconv.utils import (create_file_if_missing,
                              load_json,
@@ -303,3 +304,56 @@ def test_no_etelemetry():
     with patch.dict('sys.modules', {'etelemetry': None}):
         workflow(outdir='/dev/null', command='ls',
                  heuristic='reproin', files=[])
+
+
+# Test two scenarios:
+# -study without sessions
+# -study with sessions
+# The "expected_folder" is the session folder without the tmpdir
+@pytest.mark.parametrize(
+    "session, expected_folder", [
+        ('', 'foo/sub-{sID}'),
+        ('pre', 'foo/sub-{sID}/ses-pre')
+    ]
+)
+def test_populate_intended_for(tmpdir, session, expected_folder, capfd):
+    """
+    Tests for "process_extra_commands" when the command is
+    'populate-intended-for'
+    """
+    # Because the function .utils.populate_intended_for already has its own
+    # tests, here we just test that "process_extra_commands", when 'command'
+    # is 'populate_intended_for' does what we expect (loop through the list of
+    # subjects and calls 'populate_intended_for' using the 'POPULATE_INTENDED_FOR_OPTS'
+    # defined in the heuristic file 'example.py'). We call it using folders
+    # that don't exist, and check that the output is the expected.
+    bids_folder = expected_folder.split('sub-')[0]
+    subjects = ['1', '2']
+    process_extra_commands(bids_folder, 'populate-intended-for', [], '',
+                           'example', session, subjects, None)
+    captured_output = capfd.readouterr().err
+    for s in subjects:
+        expected_info = 'Adding "IntendedFor" to the fieldmaps in ' + expected_folder.format(sID=s)
+        assert expected_info in captured_output
+
+    # try the same, but without specifying the subjects or the session.
+    # the code in main should find any subject in the output folder and call
+    # populate_intended_for on each of them (or for each of the sessions, if
+    # the data for that subject is organized in sessions):
+    # TODO: Add a 'participants.tsv' file with one of the subjects missing;
+    #  the 'process_extra_commands' call should print out a warning
+    outdir = opj(str(tmpdir), bids_folder)
+    for subj in subjects:
+        subj_dir = opj(outdir, 'sub-' + subj)
+        print('Creating output dir: %s', subj_dir)
+        os.makedirs(subj_dir)
+        if session:
+            os.makedirs(opj(subj_dir, 'ses-' + session))
+    process_extra_commands(outdir, 'populate-intended-for', [], '',
+                           'example', [], [], None)
+    captured_output = capfd.readouterr().err
+    for s in subjects:
+        expected_info = 'Adding "IntendedFor" to the fieldmaps in ' + opj(str(tmpdir), expected_folder.format(sID=s))
+        assert expected_info in captured_output
+
+
diff --git a/heudiconv/tests/test_utils.py b/heudiconv/tests/test_utils.py
index f43a09b0..48448c25 100644
--- a/heudiconv/tests/test_utils.py
+++ b/heudiconv/tests/test_utils.py
@@ -12,7 +12,10 @@
     load_json,
     create_tree,
     save_json,
+    update_json,
     get_datetime,
+    remove_suffix,
+    remove_prefix,
     JSONDecodeError)
 
 import pytest
@@ -93,7 +96,7 @@ def test_load_json(tmpdir, caplog):
     vfname = "valid.json"
     valid_json_file = str(tmpdir / vfname)
     save_json(valid_json_file, vcontent)
-    
+
     assert load_json(valid_json_file) == vcontent
 
     calls = [0]
@@ -114,6 +117,32 @@ def json_load_patched(fp):
 
 
 
+def test_update_json(tmpdir):
+    """
+    Test utils.update_json()
+    """
+    dummy_json_file = str(tmpdir / 'dummy.json')
+    some_content = {"name": "Jason", "age": 30, "city": "New York"}
+    save_json(dummy_json_file, some_content, pretty=True)
+
+    added_content = {"LastName": "Bourne",
+                     "Movies": [
+                         "The Bourne Identity",
+                         "The Bourne Supremacy",
+                         "The Bourne Ultimatum",
+                         "The Bourne Legacy",
+                         "Jason Bourne"
+                     ]
+                     }
+    update_json(dummy_json_file, added_content)
+
+    # check that it was added:
+    with open(dummy_json_file) as f:
+        data = json.load(f)
+    some_content.update(added_content)
+    assert data == some_content
+
+
 def test_get_datetime():
     """
     Test utils.get_datetime()
@@ -121,3 +150,23 @@ def test_get_datetime():
     assert get_datetime('20200512', '162130') == '2020-05-12T16:21:30'
     assert get_datetime('20200512', '162130.5') == '2020-05-12T16:21:30.500000'
     assert get_datetime('20200512', '162130.5', microseconds=False) == '2020-05-12T16:21:30'
+
+
+def test_remove_suffix():
+    """
+    Test utils.remove_suffix()
+    """
+    s = 'jason.bourne'
+    assert remove_suffix(s, '') == s
+    assert remove_suffix(s, 'foo') == s
+    assert remove_suffix(s, '.bourne') == 'jason'
+
+
+def test_remove_prefix():
+    """
+    Test utils.remove_prefix()
+    """
+    s = 'jason.bourne'
+    assert remove_prefix(s, '') == s
+    assert remove_prefix(s, 'foo') == s
+    assert remove_prefix(s, 'jason') == '.bourne'
diff --git a/heudiconv/utils.py b/heudiconv/utils.py
index 35dae3e1..e3b74dfa 100644
--- a/heudiconv/utils.py
+++ b/heudiconv/utils.py
@@ -286,6 +286,34 @@ def json_dumps_pretty(j, indent=2, sort_keys=True):
     return js_
 
 
+def update_json(json_file, new_data, pretty=False):
+    """
+    Adds a given field (and its value) to a json file
+
+    Parameters:
+    -----------
+    json_file : str or Path
+        path for the corresponding json file
+    new_data : dict
+        pair of "key": "value" to add to the json file
+    pretty : bool
+        argument to be passed to save_json
+    """
+    for key, value in new_data.items():
+        lgr.debug(
+            'File "{f}": Setting {k} to {v}'.format(
+                f=json_file,
+                k=key,
+                v=value,
+            )
+        )
+
+    with open(json_file) as f:
+        data = json.load(f)
+    data.update(new_data)
+    save_json(json_file, data, pretty=pretty)
+
+
 def treat_infofile(filename):
     """Tune up generated .json file (slim down, pretty-print for humans).
     """
@@ -516,7 +544,11 @@ def create_tree(path, tree, archives_leading_dir=True):
             executable = False
             name = file_
         full_name = op.join(path, name)
-        if isinstance(load, (tuple, list, dict)):
+        if name.endswith('.json') and isinstance(load, dict):
+            # (For a json file, we expect the content to be a dictionary, so
+            #  don't continue creating a tree, but just write dict to file)
+            save_json(full_name, load)
+        elif isinstance(load, (tuple, list, dict)):
             # if name.endswith('.tar.gz') or name.endswith('.tar') or name.endswith('.zip'):
             #     create_tree_archive(path, name, load, archives_leading_dir=archives_leading_dir)
             # else:
@@ -575,3 +607,41 @@ def get_datetime(date, time, *, microseconds=True):
     if not microseconds:
         datetime_str = datetime_str.split('.', 1)[0]
     return datetime_str
+
+
+def remove_suffix(s, suf):
+    """
+    Remove suffix from the end of the string
+
+    Parameters:
+    ----------
+    s : str
+    suf : str
+
+    Returns:
+    -------
+    s : str
+        string with "suf" removed from the end (if present)
+    """
+    if suf and s.endswith(suf):
+        return s[:-len(suf)]
+    return s
+
+
+def remove_prefix(s, pre):
+    """
+    Remove prefix from the beginning of the string
+
+    Parameters:
+    ----------
+    s : str
+    pre : str
+
+    Returns:
+    -------
+    s : str
+        string with "pre" removed from the beginning (if present)
+    """
+    if pre and s.startswith(pre):
+        return s[len(pre):]
+    return s