Super

Some code examples.

A protein backbone superimposer.

The following will superimpose the proteins contained in the pdb files passed as arguments to the program.

For a Quick&Dirty(tm), throw-away program, most errors can just be ignored:

    package main
    
    import (
    	"fmt"
    	"github.com/rmera/gochem"
    	"os"
    	"strings"
    )
    
//This will Superimpose the backbone of two molecules
func main() {
	backbone := []string{"C", "CA", "O", "N"} //The PDB name of the atoms in the backbone
	mol1, _ := chem.PdbRead(os.Args[1], true) //true means that we try to read the symbol from the PDB file
	mol2, _ := chem.PdbRead(os.Args[2], true)
	mols := []*chem.Molecule{mol1, mol2}
	superlist := make([][]int, 2, 2)
	//We collect the atoms that are part of the backbone.
	for molnumber, mol := range mols {
		for atomindex, atom := range mol.Atoms {
			if chem.IsIn(atom.Name, backbone) != -1 {
				superlist[molnumber] = append(superlist[molnumber], atomindex)
			}
		}
	}
	mol1.Coords[0], err = chem.Super(mol1.Coords[0], mol2.Coords[0], superlist[0], superlist[1])
	if err != nil {
		panic(err.Error())
	}
	newname := strings.Replace(os.Args[1], ".pdb", "_super.pdb", -1)
	chem.PdbWrite(mol1, newname)
}

Were we have left only the error-check that ensures the superposition itself was successful. It is possible that the best rotation matrix that superimposes the structures turns out to be a reflexion. This is the case when the structures are enantiomers of each other. In this case the function Super will return and error and no structure.

Note that the program can be made to handle non-protein molecules just by changing the criteria to include indexes in superlist, and perhaps using XyzRead/XyzWrite instead of their PDB counterparts.

If the function Super is called with empty atom lists it will attempt to use the whole corresponding molecule and return an error if the atom numbers of both structures don't match.

A more proper version of the previous program, with error checking, is slightly more complex, but still very readable and concise:

    package main
    
    import (
    	"fmt"
    	"github.com/rmera/gochem"
    	"os"
    	"strings"
    )
    
//This will Superimpose the backbone of two molecules
func main() {
	backbone := []string{"C", "CA", "O", "N"} //The PDB name of the atoms in the backbone
	if len(os.Args) < 3 || !strings.Contains(os.Args[1], ".pdb") || !strings.Contains(os.Args[2], ".pdb") {
		fmt.Printf("Usage: %s  mol1.pdb mol2.pdb\nNo caps in the pdb extension allowed.\n", os.Args[0])
		os.Exit(1)
	}
	mol1, err := chem.PdbRead(os.Args[1], true) //true means that we try to read the symbol from the PDB file
	mol2, err2 := chem.PdbRead(os.Args[2], true)
	if err != nil || err2 != nil {
		panic("Unable to open input files!")
	}
	mols := []*chem.Molecule{mol1, mol2}
	superlist := make([][]int, 2, 2)
	//We collect the atoms that are part of the backbone.
	for molnumber, mol := range mols {
		for atomindex, atom := range mol.Atoms {
			if chem.IsIn(atom.Name, backbone) != -1 {
				superlist[molnumber] = append(superlist[molnumber], atomindex)
			}
		}
	}
	mol1.Coords[0], err = chem.Super(mol1.Coords[0], mol2.Coords[0], superlist[0], superlist[1])
	if err != nil {
		panic(err.Error())
	}
	newname := strings.Replace(os.Args[1], ".pdb", "_super.pdb", -1)
	chem.PdbWrite(mol1, newname)
}