@property
def bonds():
    ''''''

def run(fragObj, bonds='harm', angles='harm', dihedrals='cos', vdw='12-6', parameterise=None, unit='kcal'):
    ''''''

    from numpy     import argsort, array, dtype, full, savetxt, vstack, recarray, unique, fromiter, searchsorted, zeros, __version__
    # TD 18/06/2024 : Numpy v2.0.0 version dependence
    from numpy.lib import NumpyVersion
    if NumpyVersion(__version__) >= '2.0.0b1':
        from numpy._core import defchararray
    else:
        from numpy.core import defchararray
    from io        import StringIO
    from itertools import product
    from ..data    import atomicdata, forcefield

    strbuff = 'generated by ChemShell\n'
    strbuff += 'UNITS %s\n'%unit
    strbuff += 'molecular types 1\n'
    strbuff += 'molecule 1\n'
    strbuff += 'nummols 1\n'


    ### atoms ###

    strbuff += 'atoms %d\n'%fragObj.natoms
    dt = dtype([('names'  ,'8U'),
                ('masses' ,float),
                ('charges',float)])
    abuff = recarray(dtype=dt, shape=(fragObj.natoms,))
    abuff.names[:]   = fragObj.names.astype(str, copy=False)
    abuff.masses[:]  = fragObj.masses
    abuff.charges[:] = fragObj.charges
    with StringIO() as fp:
        savetxt(fp, abuff, fmt='%-8s %9.6f %13.8f  1  0  0')
        strbuff += fp.getvalue()

    ### bonds ###

    strbuff += 'bonds %d\n'%fragObj.nbonds
    iznums = fragObj.znums[fragObj.allbonds[:,0]]
    jznums = fragObj.znums[fragObj.allbonds[:,1]]
    dt = dtype([('i'      , int),
                ('j'      , int),
                ('contsts', float),
                ('bonds'  , float)])
    abuff = recarray(dtype=dt, shape=(fragObj.nbonds,))
    abuff.i[:]       = fragObj.allbonds[:,0]+1
    abuff.j[:]       = fragObj.allbonds[:,1]+1
    abuff.contsts[:] = forcefield.bonds[iznums,jznums,1]
    abuff.bonds[:]   = forcefield.bonds[iznums,jznums,0]
    with StringIO() as fp:
        savetxt(fp, abuff, fmt=bonds+' %6d %6d %13.6f %13.6f')
        strbuff += fp.getvalue()

    ### angles ###

    strbuff += 'angles %d\n'%fragObj.nangles
    iznums = fragObj.znums[fragObj.allangles[:,0]]
    jznums = fragObj.znums[fragObj.allangles[:,1]]
    kznums = fragObj.znums[fragObj.allangles[:,2]]
    dt = dtype([('i'      , int),
                ('j'      , int),
                ('k'      , int),
                ('contsts', float),
                ('angles' , float)])
    abuff = recarray(dtype=dt, shape=(fragObj.nangles,))
    abuff.i[:]       = fragObj.allangles[:,0]+1
    abuff.j[:]       = fragObj.allangles[:,1]+1
    abuff.k[:]       = fragObj.allangles[:,2]+1
    abuff.contsts[:] = forcefield.angles[iznums,jznums,kznums,1]
    abuff.angles[:]  = forcefield.angles[iznums,jznums,kznums,0]
    with StringIO() as fp:
        savetxt(fp, abuff, fmt=angles+' %6d %6d %6d %13.6f %13.6f')
        strbuff += fp.getvalue()

    ### dihedrals ###

    strbuff += 'dihedrals %d\n'%fragObj.ndihedrals
    iznums = fragObj.znums[fragObj.dihedrals[:,0]]
    jznums = fragObj.znums[fragObj.dihedrals[:,1]]
    kznums = fragObj.znums[fragObj.dihedrals[:,2]]
    lznums = fragObj.znums[fragObj.dihedrals[:,3]]
    dt = dtype([('i'        , int),
                ('j'        , int),
                ('k'        , int),
                ('l'        , int),
                ('contsts'  , float),
                ('dihedrals', float),
                ('m'        , int)])
    abuff = recarray(dtype=dt, shape=(fragObj.ndihedrals,))
    abuff.i[:]          = fragObj.dihedrals[:,0]+1
    abuff.j[:]          = fragObj.dihedrals[:,1]+1
    abuff.k[:]          = fragObj.dihedrals[:,2]+1
    abuff.l[:]          = fragObj.dihedrals[:,3]+1
    abuff.contsts[:]    = forcefield.dihedrals[iznums,jznums,kznums,lznums,1]
    abuff.dihedrals[:]  = forcefield.dihedrals[iznums,jznums,kznums,lznums,0]
    abuff.m[:]          = forcefield.dihedrals[iznums,jznums,kznums,lznums,2]
    with StringIO() as fp:
        savetxt(fp, abuff, fmt=dihedrals+' %6d %6d %6d %6d %13.9e %13.9e %6d 1.e+00 1.e+00')
        strbuff += fp.getvalue()

    strbuff += 'finish\n'

    ### van der Waals ###
# TODO compare with Tcl-ChemShell

    strbuff += 'vdw %d\n'%fragObj.nvdwtypes
    # what is the dtype if we use fromiter()?
    types2znums = unique(vstack((fragObj.types.astype(str, copy=False),
                                 fragObj.znums)).T,
                         axis=0)

    # make a mapping between types and znums
    types, znums = types2znums[:,0], types2znums[:,1].astype(int, copy=False)
    
    # handle of fragObj.types and the sorted indices
    alltypes = fragObj.types.astype(str, copy=False)
    sorted = argsort(alltypes)

    # znums of the associated vdW pair types
    iznums = znums[searchsorted(types, fragObj.vdwtypes[:,0])]
    jznums = znums[searchsorted(types, fragObj.vdwtypes[:,1])]

    dt = dtype([('itypes' , '8U'),
                ('jtypes' , '8U'),
                ('a'      , float),
                ('b'      , float)])

    # array buff to hold the parameters
    abuff = recarray(dtype=dt, shape=(fragObj.nvdwtypes,))
    abuff.itypes = fragObj.vdwtypes[:,0]
    abuff.jtypes = fragObj.vdwtypes[:,1]

    ## interactions between non-H and non-H ##

    select = (iznums>1)*(jznums>1)
    abuff.a[select] = forcefield.vdw[vdw]['x-x'][iznums[select],jznums[select],0]
    abuff.b[select] = forcefield.vdw[vdw]['x-x'][iznums[select],jznums[select],1]

    # make a mapping between types and bonded hydrogens' znums (for example, it is 6 if a hydrogen is bonded with carbon)
    hbonded       = full((fragObj.natoms,), -1, dtype=int)
    znums_hbonded = zeros((fragObj.natoms,), dtype=int)

    # only interested in H
    hydrogens = fragObj.znums==1
    # a list of indices of hydrogens' bonded atoms
    hbonded[hydrogens] = fragObj.conn[hydrogens][:,0]
    # and the znums
    znums_hbonded[hbonded>=0] = fragObj.znums[hbonded[hbonded>=0]]
    # interactions between H and non-H
    select = (iznums==1)*(jznums>1)
    # znums of atoms bonded to iatoms (hydrogens) and znums of jatoms
    iznums_hx = znums_hbonded[sorted[searchsorted(alltypes, fragObj.vdwtypes[:,0][select], sorter=sorted)]]
    jznums_hx = fragObj.znums[sorted[searchsorted(alltypes, fragObj.vdwtypes[:,1][select], sorter=sorted)]]
    abuff.a[select] = forcefield.vdw[vdw]['h-x'][iznums_hx,jznums_hx,0]
    abuff.b[select] = forcefield.vdw[vdw]['h-x'][iznums_hx,jznums_hx,1]

    ## interactions between non-H and H ##

    select = (iznums>1)*(jznums==1)
    # znums of atoms bonded to iatoms (hydrogens) and znums of jatoms
    iznums_xh = fragObj.znums[sorted[searchsorted(alltypes, fragObj.vdwtypes[:,0][select], sorter=sorted)]]
    jznums_xh = znums_hbonded[sorted[searchsorted(alltypes, fragObj.vdwtypes[:,1][select], sorter=sorted)]]
    # NB: i- and jznums are swapped when looping up in the forcefield.vdw table because of the latter is fixed-formatted as H and non-H
    abuff.a[select] = forcefield.vdw[vdw]['h-x'][jznums_xh,iznums_xh,0]
    abuff.b[select] = forcefield.vdw[vdw]['h-x'][jznums_xh,iznums_xh,1]
    # znums of iatoms and znums of atoms bonded to jatoms (hydrogens)

    ## interactions between H and H ##

    select = (iznums==1)*(jznums==1)
    # znums of atoms bonded to iatoms (hydrogens) and znums of atoms bonded to jatoms (other hydrogens)
    iznums_hh = znums_hbonded[sorted[searchsorted(alltypes, fragObj.vdwtypes[:,0][select], sorter=sorted)]]
    jznums_hh = znums_hbonded[sorted[searchsorted(alltypes, fragObj.vdwtypes[:,1][select], sorter=sorted)]]
    abuff.a[select] = forcefield.vdw[vdw]['h-h'][iznums_hh,jznums_hh,0]
    abuff.b[select] = forcefield.vdw[vdw]['h-h'][iznums_hh,jznums_hh,1]
#    abuff.a[:] = forcefield.vdw['12-6']['h-h'][iznums,jznums]
    with StringIO() as fp:
        savetxt(fp, abuff, fmt='%6s %6s  '+vdw+' %15.6e %15.6e')
        strbuff += fp.getvalue()

    strbuff += 'close\n'

    return strbuff
