diff --git a/common/python/nomadcore/md_data_access/MDDataAccess.py b/common/python/nomadcore/md_data_access/MDDataAccess.py
index 715854a262b927492bb9e847a376b46736a71f47..d762bc805aff5c48fde1eca2b7d1f5282c8c902e 100644
--- a/common/python/nomadcore/md_data_access/MDDataAccess.py
+++ b/common/python/nomadcore/md_data_access/MDDataAccess.py
@@ -2269,194 +2269,3 @@ if __name__ == "__main__":
             atom = MDdata.iread()
     print(i)
 
-    if "mdanalysis" in MDdata.topocode: 
-
-        def getseg(seg):
-            regex = re.compile(r"^seg_[0-9]+_b")
-            if regex.findall(seg.segid):
-                segname = bytes(re.sub(r"^seg_[0-9]+_b","",seg.segid).replace("'",""), "utf-8").decode("utf-8")
-            else:
-                segname = seg.segid
-            return [seg.ix,segname]
-
-        def get_atomseg(atom):
-            return [atom.segment.ix,atom.segment.segname.decode('utf-8')]
-        
-        def getres(atom):
-            return [atom.resname.decode('utf-8'),atom.resid]
-
-        def getatom(atom):
-            atmid = atom.ix
-            atmname = atom.name.decode('utf-8')
-            atmtyp = atom.type.decode('utf-8')
-            atmres = atom.resname.decode('utf-8')
-            atmresid = atom.resid
-            atmsegid = atom.segid
-            return [atmid,atmname,atmtyp,atmres,atmresid,atmsegid]
-
-        def residuepairs(atom1,atom2):
-            return [getres(atom1),getres(atom2)]
-
-        def atompairs(atom1,atom2):
-            return [getatom(atom1),getatom(atom2)]
-
-        print(MDdata.topohandler.atoms)
-        attrlist = dir(MDdata.topohandler.atoms)
-        if "names" in attrlist:
-            print("Atom Names:",MDdata.topohandler.atoms.names)
-        if "types" in attrlist:
-            print("Atom Types:",MDdata.topohandler.atoms.types)
-        if "masses" in attrlist:
-            print("Atom Masses:",MDdata.topohandler.atoms.masses)
-        if "charges" in attrlist:
-            print("Atom Charges:",MDdata.topohandler.atoms.charges)
-        if "radiuses" in attrlist:
-            print("Atom Radiuses:",MDdata.topohandler.atoms.radiuses)
-        if "bfactors" in attrlist:
-            print("Atom bFactors:",MDdata.topohandler.atoms.bfactors)
-        print("Atom List:",[mda.topology.guessers.guess_atom_element(getatom(a)[1]) for a in MDdata.topohandler.atoms])
-        print("Atom Residues:",MDdata.topohandler.atoms.residues)
-        print("Atom Residue ids:",MDdata.topohandler.atoms.residues.ix)
-        print("Atom Residue Names:",[a.resname for a in MDdata.topohandler.atoms.residues])
-        print("Atom Segments:",MDdata.topohandler.atoms.segments)
-        print("Atom Bonds:",MDdata.topohandler.atoms.bonds)
-        print("Atom Angles:",MDdata.topohandler.atoms.angles)
-        print("Atom Dihedrals:",MDdata.topohandler.atoms.dihedrals)
-        print("Atom Impropers:",MDdata.topohandler.atoms.impropers)
-        print(MDdata.topohandler.bonds)
-        topd = MDdata.topohandler.bonds.topDict
-        topdk = topd.keys()
-        for key in topdk:
-            topt = topd[key]
-            for b in topt:
-                reslist=residuepairs(b.atoms[0],b.atoms[1])
-                atmlist=atompairs(b.atoms[0],b.atoms[1])
-                if str(reslist[0][0]) != str(reslist[1][0]):
-                    print("RES:",reslist,
-                          "ATM:",atmlist)
-
-
-        #print("Bonds:",[ [a[0].ix,a[1].ix] for a in MDdata.topohandler.bonds])
-        print(MDdata.topohandler.angles)
-        #print("Angles:",[ [a[0].ix,a[1].ix,a[2].ix] for a in MDdata.topohandler.angles])
-        print(MDdata.topohandler.residues)
-        #print("Residue bonds:",[a.bonds for a in MDdata.topohandler.residues])
-        #for r in MDdata.topohandler.residues:
-        #    print(r.resname,r.resid)
-        print([getseg(seg)[1] for seg in MDdata.topohandler.segments])
-        #print("Segment Names:",[a.segment for a in MDdata.topohandler.segments])
-        #print("Residue Names:",MDdata.topohandler.residues.names)
-        #print("Residue Types:",MDdata.topohandler.residues.names)
-        #print("Residue Bonds:",MDdata.topohandler.residues.bonds)
-        #print("Residue Angles:",MDdata.topohandler.residues.angles)
-
-
-        print(mda.topology.guessers.guess_atom_element("BERK"))
-
-        #angles = mdtraj.compute_angles(mytopology)
-        #dihedrals = mdtraj.compute_dihedrals(mytopology)
-        #print(dihedrals)
-        #print(' Chains: %s' % [res for res in mytopology.chains])
-        #print(' Residues: %s' % [res for res in mytopology.residues])
-        #print(' Atoms: %s' % [res for res in mytopology.atoms])
-        #print(' Bonds: %s' % [res for res in mytopology.bonds])
-        #print(mytopology.chain(0))
-        table, bonds = mytopology.to_dataframe()
-        print(bonds)
-        print(len(bonds))
-        array = table.get("resSeq")
-        #print([res for res in array])
-        mydict = table.to_dict(orient='list')
-        mylist = mydict["resSeq"]
-        print('mylist')
-        print(mylist)
-        print('mylist')
-        atomindex = np.arange(len(mylist))
-        passarray = np.zeros((len(mylist), 2), dtype=int)
-        passarray[:,0] = atomindex
-        passarray[:,1] = mylist
-        print(table)
-        print(passarray)
-        #print(bonds)
-        #for bond in bonds:
-        #    print(bond)
-        mol_to_mol_bond = []
-        mol_to_mol_bond_num = []
-        atom_in_mol_bond = []
-        atom_in_mol_bond_num = []
-        index=0
-        print(len([res for res in mytopology.bonds]))
-        for res in mytopology.bonds:
-            molname1, molatom1 = str(res[0]).split('-')
-            molname2, molatom2 = str(res[1]).split('-')
-            if molname1 in molname2: 
-                atom_in_mol_bond.append(res)
-                atom_in_mol_bond_num.append(bonds[index])
-            else:
-                mol_to_mol_bond.append(res)
-                mol_to_mol_bond_num.append(bonds[index])
-            index += 1
-        print(mol_to_mol_bond)
-        print(np.array(mol_to_mol_bond_num))
-
-        atom_list = list(set(mydict["name"]))
-        atom_type_dict = {}
-        print(atom_list)
-        print(len(atom_list))
-        masses = {}
-        for ielm in range(len(atom_list)):
-            elm = atom_list[ielm]
-            atom_type_dict.update({elm : ielm})
-            for atom in mytopology.atoms:
-                if elm == atom.name:
-                    masses.update({atom.name : atom.element})
-        print([atom_list, masses])
-        interDict = {}
-        interTypeDict = {}
-        for elm in atom_list:
-            bondList = []
-            typeList = []
-            bondid = 0
-            for bond in mytopology.bonds:
-                molname1, molatom1 = str(bond[0]).split('-')
-                molname2, molatom2 = str(bond[1]).split('-')
-                if (elm == str(molatom1) or elm == str(molatom2)):
-                    bondList.append(list(bonds[bondid]))
-                    #typeList.append(list([molatom1,molatom2]))
-                    typeList.append(list([
-                        atom_type_dict[str(molatom1)],
-                        atom_type_dict[str(molatom2)]
-                        ]))
-                    interDict.update({elm : bondList})
-                    interTypeDict.update({elm : typeList})
-                bondid += 1
-
-    #    print(interDict)
-        print(interDict["CH3"])
-        print(interTypeDict["CH3"])
-        print(len(interDict["CH3"][0]))
-
-    #    print([[atom.element.radius, atom.index, atom.name, atom.residue.name, atom.element] for atom in mytopology.atoms])
-
-    # Atom labels:
-        print(table)
-    #    for atom in mytopology.atoms:
-    #        print(atom)
-        labellist = mydict["element"]
-        #print(labellist)
-    #    atomlabels = np.arange(len(labellist))
-    #    labelarray = np.zeros((len(labellist), 2), dtype=int)
-    #    labelarray[:,0] = atomlabels
-    #    labelarray[:,1] = labellist
-    #    print(labelarray)
-    #    print(readdata.unitcell_lengths)
-        latticevectors = readdata.iread().unitcell_vectors
-        box = readdata.iread().unitcell_lengths
-        step = 0
-        for latvec in latticevectors:
-            print(step, latvec)
-            step += 1
-        for pbc in box:
-            print(step, box)
-            step += 1
-