def make_gromacs(simulation, directory, clean=False):
    """Create gromacs directory structure"""
    if clean is False and os.path.exists(directory):
        raise ValueError(
            'Cannot override {}, use option clean=True'.format(directory))
    else:
        shutil.rmtree(directory, ignore_errors=True)
        os.mkdir(directory)

    # Check custom simulation potential
    
    if simulation.potential.intermolecular.type == 'custom':
        
        for pair in simulation.potential.intermolecular.special_pairs:
            table = to_table(simulation.potential.intermolecular.pair_interaction(*pair), 
                             simulation.cutoff)
            fname1 = os.path.join(directory, 
                                  'table_{}_{}.xvg'.format(pair[0], pair[1]))
            fname2 = os.path.join(directory, 
                                  'table_{}_{}.xvg'.format(pair[1], pair[0]))
            
            with open(fname1, 'w') as fd:
                fd.write(table)

            with open(fname2, 'w') as fd:
                fd.write(table)
            

        ndx = {'System' : np.arange(simulation.system.n_atoms, dtype='int')}
        for particle in simulation.potential.intermolecular.particles:
            idx = simulation.system.where(atom_name=particle)['atom'].nonzero()[0]
            ndx[particle] = idx
        
        with open(os.path.join(directory, 'index.ndx'), 'w') as fd:
            fd.write(to_ndx(ndx))
            
    # Parameter file
    mdpfile = to_mdp(simulation)
    with open(os.path.join(directory, 'grompp.mdp'), 'w') as fd:
        fd.write(mdpfile)

    # Topology file
    topfile = to_top(simulation.system, simulation.potential)
    with open(os.path.join(directory, 'topol.top'), 'w') as fd:
        fd.write(topfile)

    # Simulation file
    datafile(os.path.join(directory, 'conf.gro'),
             'w').write('system', simulation.system)

    return directory
Create gromacs directory structure