Index: gorithms.py
===================================================================
--- algorithms.py	(revision 5fef0d738f6c69046859b9b41da3ebeefc2efa61)
+++ 	(revision )
@@ -1,91 +1,0 @@
-# algorithms.py
-#
-# Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
-#                      Jan H. Meinke, Sandipan Mohanty
-#
-import copy
-import random
-import smmp_p as smmp
-from math import *
-from universe import *
-from protein import *
-
-class CanonicalMonteCarlo(object):
-    """An implementation of a canonical Monte Carlo run using Metropolis weights.
-    After a Monte-Carlo sweep over all internal variables, this implementation 
-    calculates the interactions between all Protein s in the Universe
-    """
-    
-    def __init__(self, myUniverse, nequi, sweeps, weight =  lambda x : x * smmp.bet.beta):
-        self.__defSweeps = sweeps
-        self.__defEqui = nequi
-        self.__univ = myUniverse
-        self.__acc = 0
-        self.__performedSweeps = 0
-        self.__minEnergy = myUniverse.energy()
-        self.__l = 0.1
-        self.__minCounter = 0
-        self.__weight = weight
-        
-    ## Fraction of accepted Monte Carlo moves of the total number of attempted moves.
-    #  @pre At least one move must have been performed before this function can be
-    #       called
-    #  @return fraction of accepted moves as a number between 0 and 1.
-    def acceptanceRate(self):
-        return float(self.__acc) / (self.__performedSweeps * smmp.mol_par.nvr)
-    ## The lowest energy ever seen at the end of a Monte Carlo move.
-    #  @return lowest energy seen at the end of a sweep, sweep
-    def minimumEnergy(self):
-        return self.__minEnergy, self.__minCounter
-    ## Number of sweeps performed so far
-    def performedSweeps(self):
-        return self.__performedSweeps
-        
-    def run(self):
-        """A complete Monte Carlo run including initial equilibration."""
-        self.sweeps(self.__defSweeps)
-        
-    def equilibrate(self):
-        """Do nequi sweeps to eqilibrate the system."""
-        self.sweeps(self.__defEqui)
-        
-    def sweeps(self, n):
-        """Perform n sweeps."""
-        for i in range(0, n):
-            self.sweep()
-        
-    def sweep(self):
-        """Performs a single Monte Carlo sweep."""
-        eol = self.__univ.energy()
-        #eol, self.__acc = smmp.metropolis(eol, self.__acc, self.__weight)
-        eol = self.metropolis()
-        self.__performedSweeps += 1
-        if eol < self.__minEnergy:
-            self.__minEnergy = eol
-            self.__minCounter = self.__performedSweeps
-            smmp.outpdb(0, "best.pdb")
-
-    def metropolis(self):
-        """Implementation of the Metropolis algorithm."""
-        eol = self.__univ.energy()
-        # Update the internal degrees of freedom first
-        for i in range(0, smmp.mol_par.nvr):
-            jv = smmp.var_i.idvr[i] - 1
-            vrol = smmp.var_r.vlvr[jv]
-            dv=smmp.var_r.axvr[jv] * (random.random() - 0.5)
-            smmp.var_r.vlvr[jv] = (vrol + dv + pi) % (2 * pi) - pi
-            etrial = self.__univ.energy()
-            dW = self.__weight(etrial) - self.__weight(eol)
-            if dW < 0:
-                eol = etrial
-                self.__acc += 1
-            else:
-                rd =random.random()
-                dW = min(dW, 100)
-                if  exp(-dW) > rd:
-                    eol = etrial
-                    self.__acc +=1
-                else:
-                    smmp.var_r.vlvr[jv] = vrol
-        # Update the global degrees of freedom
-        return eol
Index: in.py
===================================================================
--- main.py	(revision 5fef0d738f6c69046859b9b41da3ebeefc2efa61)
+++ 	(revision )
@@ -1,114 +1,0 @@
-#!/usr/bin/env python
-# main.py
-#
-# Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
-#                      Jan H. Meinke, Sandipan Mohanty
-
-## \mainpage pySMMP: Python bindings to SMMP v. 1.1.
-#    SMMP is a FORTRAN library for simulating proteins. It provides 
-#    minimization routines and various Monte Carlo methods. SMMP uses an all-atom
-#    representation of the protein. It stores the configuration within the 
-#    standard geometry model. In the standard geometry model the internal degrees
-#    of freedom, e.g., the dihedral angles, the bond lengths, etc. are stored and
-#    the cartesian coordinates are calculated on the fly. For more details on 
-#    the model SMMP uses and the functions it provides see the manual (manual.ps).
-#   
-#    pySMMP builds on SMMP. It provides python bindings to the functions and 
-#    global variables, but it uses classes to represent the Proteins and 
-#    algorithms. 
-#    
-#    \section sec_tutorial A Tutorial: A canonical Monte Carlo simulation of protein A
-#
-#    The basic setup is the same for most simulations. This tutorial takes you
-#    step by step to a running simulation. First we need to import the library
-#    \code 
-#         import smmp
-#     \endcode
-#    This makes the smmp routines available to the Python interpreter. Next we
-#    create the universe by initializing a universe.Universe object by first 
-#    importing the universe package and then setting the object:
-#    \code
-#        import universe
-#        myUniverse = universe.Universe(T=300)  
-#    \endcode
-#    This initializes the smmp library as well. We also set the initial 
-#    temperature of the universe to 300K. The Universe object initializes the 
-#    SMMP library and keeps track of global parameters, e.g, the temperature T.
-#    After you created a Universe, you can initialize a Protein. Proteins can
-#    be read from a sequence and a variable file or from a PDB file.
-#    \code
-#        import protein
-#        
-#        protA = protein.Protein("1bdd.seq", "1bdd.var")
-#        myUniverse.add(protA)
-#    \endcode
-#    Adding the protein to myUniverse makes the universe aware of it.
-#    Now you are ready to run a simulation. You can also querry the protein for 
-#    various properties, e.g., 
-#    \code
-#        print protA.rgyr()
-#    \endcode
-#    prints the radius of gyration of the protein.
-#    A canonical Monte Carlo algorithm is available in the algorithms package.
-#    \code
-#        import algorithms
-#        myMC = algorithms.CanonicalMonteCarlo(myUniverse, 1000, 10000)
-#    \endcode
-#    This sets up a simulation with 1000 steps for equilibration and a default 
-#    length of 10000 sweeps.
-#    The quickest way to run the simulation is by calling the run method of the 
-#    algorithm object myMC to do the entire Monte Carlo run.
-#    \code
-#        myMC.run()
-#    \endcode
-
-import sys, time
-from optparse import OptionParser
-from protein import *
-from universe import *
-from algorithms import CanonicalMonteCarlo
-
-if __name__ == "__main__":
-# Set up command line options    
-    usage = "usage: %prog [options] seq-filename"
-    parser = OptionParser(usage = usage, version="%prog 20050805")
-    parser.add_option("-v", "--variables", dest="varFileName",
-                  help="read amino-acid conformation from FILE", metavar="FILE")
-    parser.add_option("--equi", type="int", dest="nequi", default=100,
-                  help="number of sweeps for equilibration", metavar="n")
-    parser.add_option("--sweeps", type="int", dest="sweeps", default=10000,
-                  help="number of sweeps for annealing", metavar="n")
-    parser.add_option("--mes", type="int", dest="mes", default=10,
-                  help="number of sweeps between measurements", metavar="n")
-    parser.add_option("--Tmax", type="float", dest="tmax", default=1000,
-                  help="maximum temperature for annealing", metavar="T")
-    parser.add_option("--Tmin",  type="float", dest="tmin", default=100, 
-                  help="maximum temperature for annealing",metavar="T")
-
-    (options, args) = parser.parse_args()
-    if len(args) < 1:
-        parser.print_help()
-        sys.exit()
-    seqFileName = args[0]
-    
-# Start simulation setup.    
-    myUniv = Universe(T=300, st=0)   
-    p1 = Protein(seqFileName, options.varFileName)
-    p1.setOrigin((0, 0, 0), (0, 0, 0))
-    myUniv.add(p1)
-    T = 300.0
-    smmp.outpdb(1, "start.pdb")
-    print "E of starting configuration.", myUniv.energy()
-    ts = time.time()
-    smmp.minim(1, options.sweeps, 1.0E-12)
-
-    smmp.outpdb(1, "equi.pdb")
-    
-    t1 = time.clock()
-    print myUniv.energy(), myUniv.rgyr(), myUniv.helix(), p1.hbond(), (t1-ts)
-    p1.saveVariables("final.var")
-    smmp.outpdb(0, "final.pdb")
-    te = time.time()
-    print "The calculation took %ss." % (te - ts)
- 
-    
Index: otein.py
===================================================================
--- protein.py	(revision 5fef0d738f6c69046859b9b41da3ebeefc2efa61)
+++ 	(revision )
@@ -1,290 +1,0 @@
-#!/usr/bin/env python
-# protein.py
-#
-# Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
-#                      Jan H. Meinke, Sandipan Mohanty
-#
-import sys
-import numpy as np
-import smmp_p as smmp
-from math import *
-
-crd = 57.2957795130823
-
-def nursvr(ivr):
-    res = 0
-    for i in range(smmp.mol_i.ntlml - 1, -1, -1):
-        ifirs = smmp.mol_i.irsml1[i] - 1
-        if ivr + 1 > smmp.res_i.ivrrs1[ifirs] and smmp.mol_i.nvrml[i] > 0:
-            for j in range(smmp.mol_i.irsml2[i] - 1, ifirs - 1, -1):
-                if ivr + 1 > smmp.res_i.ivrrs1[j]  and smmp.res_i.nvrrs[j] > 0:
-                    res = j + 1
-                    return res
-    return res
-
-## Provides read-only access to SMMP's atoms. Usually, you won't initialize
-#  an Atom directly since you need to know the internal array index of the
-#  Atom in SMMP. You can get a list of all the Atoms in a Protein by using 
-#  protein.Protein.getAtoms().
-#  Various properties are accessible. See the functions' documentations for
-#  details.
-class Atom:
-    """Atoms are the basic chemical units."""
-    #TODO: insert namespace
-    
-    ## Initialize the atom with its ID. The ID is the array index of the atom
-    #  in SMMP. 
-    #       @param id index of atom in SMMP
-    def __init__(self, id):
-        self.__id = id
-    
-    def position(self):
-        """Cartesian coordinates of this atom."""
-        res = (smmp.atm_r.xat[self.__id], smmp.atm_r.yat[self.__id], smmp.atm_r.zat[self.__id])
-        return res
-    
-    def charge(self):
-        """Partial charge of this atom."""
-        return smmp.atm_r.cgat[self.__id]
-    
-    def radius(self):
-        """Van-der-Waals Radius of this atom."""
-        return smmp.solvent.rvdw[self.__id]
-        
-    def name(self):
-        return ''.join(smmp.atm_c.nmat[:,0].reshape([-1,4])[self.__id]).strip()
-        
-    def __str__(self):
-        return "%s at (%s, %s, %s)" % (self.name(), self.position()[0],  self.position()[1], self.position()[2])
-
-## Amino acids bound within a protein are often called residues. Strictly
-#  speaking the residue is the part that's connected to the Calpha and not part
-#  of the backbone, but we'll include the backbone atoms as well. Just like Atoms
-# Residues reference the corresponding data in SMMP and are not independent. 
-class Residue:
-    """An amino acid within a Protein."""
-    
-    ## Initialize the residue with the corresponding index in SMMP.
-    #  @param id index of residue in SMMP
-    def __init__(self, id):
-        self.__id = id
-
-    def atoms(self):
-        """Returns the list of atoms belonging to this residue."""
-        first = smmp.res_i.iatrs1[self.__id]
-        last = smmp.res_i.iatrs2[self.__id]
-        res = [Atom(i) for i in range(first - 1, last - 1)]
-        return res
-    
-        
-## Proteins are constructed from a sequence file and a configuration file or 
-#  read from a PDB file.
-class Protein:
-    """A protein is a polypeptide made of the 20 naturally occuring amino 
-    acids.
-    """
-    
-    def __init__(self, seq, var=' '):
-        """The structure of the protein is read in from a file. 
-        Pass the name of a textfile with a sequence string and--optionally--the name of a
-        variable file that describes the conformation. Or give the name of pdb
-        file---it must end with 'pdb'---that contains both information."""
-        
-        grpn = 'nh2'
-        grpc = 'cooh'
-        prevID = smmp.mol_i.ntlml - 1
-        if seq[-3:].lower() == 'pdb':
-            if smmp.mol_i.ntlml:
-                print "Sorry, but I can only read the first protein from a pdb file." 
-                raise ValueError
-            smmp.init_molecule(0, grpn, grpc, seq)
-        else:
-            smmp.init_molecule(1, grpn, grpc, seq, var)
-        self.__id = smmp.mol_i.ntlml - 1
-        if self.__id - prevID > 1:
-            print "Sorry, I don't know how to deal with multiple proteins in one sequence file, yet."
-            raise NotImplementedError
-    
-    def setOrigin(self, x, a):
-        """Sets the position of the N terminal of the protein in the global 
-        coordinate system and the orientation of the protein. Expects two 
-        tuples. The first contains the coordinates of the first atom (x, y,
-        and z) and the second contains three angles.
-        """
-        #TODO: Implement type checking
-        smmp.mol_r.gbpr[0][self.__id] = x[0]
-        smmp.mol_r.gbpr[1][self.__id] = x[1]
-        smmp.mol_r.gbpr[2][self.__id] = x[2]
-        smmp.mol_r.gbpr[3][self.__id] = a[0]
-        smmp.mol_r.gbpr[4][self.__id] = a[1]
-        smmp.mol_r.gbpr[5][self.__id] = a[2]
-
-    def origin(self):
-        r = [smmp.mol_r.gbpr[i][self.__id] for i in range(0, 3)]
-        a = [smmp.mol_r.gbpr[i][self.__id] for i in range(3, 6)]
-        return r, a
-    
-    def getDihedral(self, angle, residue):
-        try:
-            res = int(residue)
-            idx = res
-        except ValueError:
-            res = residue.split()
-            res[0] = res[0].lower()
-            res[1] = int(res[1])
-            idx = self.findResidue(res[0], res[1])
-        varIdx = ''
-        for i in range(smmp.res_i.ivrrs1[idx], smmp.res_i.ivrrs1[idx] + smmp.res_i.nvrrs[idx]):
-            varName = ''.join(list([str(smmp.var_c.nmvr[j][i]) for j in range(0, 3)]))
-            if varName == angle:
-                varIdx = i
-                break
-        if varIdx:
-            return smmp.var_r.vlvr[varIdx] * crd
-        raise KeyError, '%s not found for %s.' % (angle, residue)
-    
-    ## Residues are given by their 3-letter abbreviation and their index. 
-    #  For example, findResidue('gly', 5) returns the index of the fifth occurence
-    #  of glycine in the protein. If the residue is not found an exception is raised.
-    #
-    #  @param name 3-letter abbreviation of the residue
-    #  @param idx index of occurence of the amino acid given in name.
-    
-    def findResidue(self, name, idx = 1):
-        """Find the idxth residue of type name and return the index."""
-        res = ''.join(self.residues()).lower()
-        name = name.lower()
-        j = 0
-        try:
-            for i in range(0, idx):
-                k = res.index(name, j)
-                j = j + k
-        except:
-            raise KeyError, "%s%s not found." % (name, idx)
-        if j >=len(res):
-            raise KeyError, "%s%s not found." % (name, idx)
-        return j / 3
-        
-    def atoms(self):
-        """Returns the list of atoms belonging to this protein."""
-        first = smmp.res_i.iatrs1[smmp.mol_i.irsml1[self.__id] - 1]
-        last = smmp.res_i.iatrs2[smmp.mol_i.irsml2[self.__id] - 1]
-        res = [Atom(i) for i in range(first - 1, last - 1)]
-        return res
-    
-    def residues(self):
-        """Returns a list of the amino acids in this protein."""
-        seq = [''.join(l) for l in np.array([str(c[0]) for c in smmp.res_c.seq]).reshape([-1, 4])]
-        first = smmp.mol_i.irsml1[self.__id]- 1
-        last = smmp.mol_i.irsml2[self.__id] - 1
-        return seq[first:last]
-        
-    def id(self):
-        return self.__id
-        
-    def strand(self):
-        """Returns the number of residues in a beta-strand-like configuration."""
-        maxDev = 30.0
-        phiSheet = -150.0
-        psiSheet = 150.0
-        strandCount = 0
-        ifivr = smmp.mol_i.ivrml1[self.__id]
-        i = 0
-        iv = smmp.var_i.idvr[i]
-        while iv < ifivr:
-            i += 1
-            iv = smmp.var_i.idvr[i]
-        for j in range(i, ifivr + smmp.mol_i.nvrml[self.__id]):
-            iv = smmp.var_i.idvr[j] - 1
-            name =''. join([str(smmp.var_c.nmvr[iv * 3 + k][0]) for k in range(0,3)])
-            if name == 'phi':
-                phi = smmp.var_r.vlvr[iv] * crd
-                psi = smmp.var_r.vlvr[smmp.var_i.idvr[j+1] - 1] * crd
-                if abs(phi - phiSheet) < maxDev and abs(psi - psiSheet) < maxDev:
-                    strandCount += 1
-        return strandCount
-    
-    def directionVector(self):
-        """Calculates the normalized end-to-end vector for this molecule. 
-        The first and last c_alpha carbon is used as reference."""
-        atoms = self.atoms()
-        beginning = []
-        for a in atoms:
-            if a.name().strip() == 'ca':
-                if beginning:
-                    end = a.position()
-                else:
-                    beginning = a.position()
-        res = [end[i] - beginning[i] for i in range(0,3)] 
-        sum = sqrt(reduce(lambda x, y: x + y ** 2, res))
-        res = map(lambda x: x / sum, res)
-        return res
-        
-    def contacts(self):
-        """Returns a triplet containing the total number of contacts, the number
-        of native contacts, and the hemming distance. The latter two quantities 
-        are only meaningful if a reference structure has been defined.
-        """
-        return smmp.contacts()
-
-    def hbond(self):
-        return smmp.hbond(self.__id + 1, 0)
-        
-    def energy(self):
-        return smmp.energy()
-    
-    def rgyr(self):
-        return smmp.rgyr(self.__id + 1) 
-    
-    def rmsdfun(self):
-        return smmp.rmsdfun()
-        
-    def savePDB(self, file=""):
-       """Write the current configuration to file using PDB formatting."""
-       smmp.outpdb(self.__id, file)
-       
-    def saveVariables(self, file=""):
-        """Output the current conformation in the variable format. The output 
-        can be read in again and used to restart/continue a simulation.
-        
-        @param file if empty write to standard out
-        """
-        # Same functionality as smmp.outvar(self.__id + 1, 0)
-        if not file:
-            out = sys.stdout
-        else:
-            out = open(file, 'w')
-        for i in range(smmp.mol_i.ivrml1[self.__id] - 1 , smmp.mol_i.ivrml1[self.__id] + smmp.mol_i.nvrml[self.__id] - 1):
-            if not smmp.var_l.fxvr[i]:
-                name = ''.join([str(smmp.var_c.nmvr[i * 3 + j][0]) for j in range(0,3)])
-                print >>out, nursvr(i), ':', name, ':', smmp.var_r.vlvr[i] * crd
-            else:
-                it=smmp.var_i.ityvr[i]
-                if it == 3:
-                    print >>out, nursvr(i), ':', name, ':', smmp.var_r.vlvr[i] * crd, ' &'
-        
-    def zimmer(self):
-        first = smmp.mol_i.irsml1[self.__id] - 1
-        last = smmp.mol_i.irsml2[self.__id]
-        nrs = smmp.mol_i.irsml2[smmp.mol_i.ntlml - 1] - 1
-        smmp.zimmer(last)
-        return smmp.zimme.zimm[first:last]
-
-    def __len__(self):
-        return smmp.mol_i.irsml2[self.__id]  - smmp.mol_i.irsml1[self.__id] + 1
-
-
-if __name__ == "__main__":
-    libdir = 'SMMP/'
-
-    smmp.epar_l.flex = 0
-    smmp.epar_l.sh2 = 0
-    smmp.epar_l.epsd = 0
-    smmp.epar_l.tesgrd = 0
-    smmp.isolty.itysol = 0
-    smmp.init_energy(libdir)
-    
-    p1 = Protein('EXAMPLES/1bdd.pdb')
-    for i in p1.atoms():
-        print i
-    print p1.getDihedral('phi', 11)
Index: pysmmp/ParallelTempering.py
===================================================================
--- pysmmp/ParallelTempering.py	(revision 225475cd59f91e8a64cbbf82f4e34e4484191092)
+++ pysmmp/ParallelTempering.py	(revision 225475cd59f91e8a64cbbf82f4e34e4484191092)
@@ -0,0 +1,117 @@
+#!/usr/bin/env python
+
+#import smmp
+import sys, copy,  random,  math,  numpy as np
+
+from mpi4py import MPI 
+#from algorithms import CanonicalMonteCarlo
+
+class ParallelTempering(object):
+    """Parallel tempering is a generalized Monte Carlo algorithm. A set of
+    replica is simulated at different temperatures. Each replica performs a 
+    canoncial Monte Carlo simulation. After a given number of sweeps, an exchange
+    is attempted between neighboring replicas.
+    """
+    
+    def __init__(self, T,  myUniverse):
+        """Initialize ParallelTempering with a set of tempertures T.
+        """
+        self.kB = 1.98773e-3
+        
+        self._T = copy.deepcopy(T)
+        self._numberOfReplica = len(T)
+
+        # Calculate the number of processors available per replica
+        numberOfProcessors = MPI.COMM_WORLD.Get_size()
+        numberOfProcessorsPerReplica = numberOfProcessors // self._numberOfReplica
+        numberOfProcessors = numberOfProcessorsPerReplica * self._numberOfReplica
+        self._myReplica = MPI.COMM_WORLD.Get_rank() // numberOfProcessorsPerReplica
+
+        # Create the replica internal communicators for the energy function
+        groupWorld = MPI.COMM_WORLD.Get_group()
+        group = [MPI.Group.Range_incl(groupWorld, [[i * numberOfProcessorsPerReplica, (i + 1) * numberOfProcessorsPerReplica - 1, 1], ]) for i in xrange(self._numberOfReplica)]
+        comm = [MPI.COMM_WORLD.Create(g) for g in group]
+        
+        # Create the communicator needed to exchange information about replica
+        ptGroup = MPI.Group.Incl(groupWorld, range(0, numberOfProcessors, numberOfProcessorsPerReplica))
+        self._ptComm = MPI.COMM_WORLD.Create(ptGroup)
+        
+        self._localMaster = self._ptComm != MPI.COMM_NULL
+        if self._localMaster:
+            self._myPTRank = self._ptComm.Get_rank()
+            
+        self._oddCycle = False
+        
+        
+    def exchange(self):
+        """Collect temperatures and energies from all replicas, determine neighbor and do replica exchange move.
+        
+        The exchange step is the heart of parallel tempering. The probability of exchange between two 
+        temperatures i and j is given by P_{ij} = \Delta \beta \Delta E, where \beta = 1/k_B T.
+        """
+        if self._localMaster:
+#            energies = ptComm.allgather(self._myUniverse.energy())
+            myE = [250,  255,  75,  125, 120, 185][self._myReplica]
+            myT = [420,  370,  250,  275,  300,  330][self._myReplica]
+            energies = self._ptComm.allgather(myE)
+#            temperatures = ptComm.allgather(self._myUniverse.temperature())
+            temperatures = self._ptComm.allgather(myT)
+            
+            temperatureRank = zip(temperatures,  np.arange(self._numberOfReplica))
+            temperatureRank.sort()
+            
+            if (self._myPTRank % 2):
+                if self._oddCycle:
+                    rank0 = self._myPTRank
+                else:
+                    rank0 = self._myPTRank - 1
+                
+                rank1 = (rank0 + 1) % self._ptComm.Get_size()
+                
+                self._oddCycle = not self._oddCycle
+                
+                E0 = energies[temperatureRank[rank0][1]]
+                E1 = energies[temperatureRank[rank1][1]]
+                T0 = self._T[rank0]
+                T1 = self._T[rank1]
+                beta0 = 1.0 / (T0 * self.kB)
+                beta1 = 1.0 / (T1 * self.kB)
+                
+                P = math.exp((beta0 - beta1) * (E0 - E1))
+                print self._myPTRank,  P
+                
+                if P >= random.random(): # accept exchange
+                    rank = self._ptComm.allgather(np.array([temperatureRank[rank1][1],  temperatureRank[rank0][1]]))
+                else:
+                    rank = self._ptComm.allgather(np.array([temperatureRank[rank0][1],  temperatureRank[rank1][1]]))
+            else:
+                random.random() # Do one random number to stay in sync with the other processes for the same replica.
+                rank = self._ptComm.allgather(np.array([-1,  -1]))
+
+            rank = np.array(rank[1::2]).flatten()
+            print(rank)
+            rankTemperature = zip(rank,  self._T)
+            rankTemperature.sort()
+            print self._myPTRank,  rankTemperature[self._myReplica][1]
+#            self._myUniverse.setTemperature(rankTemperature[self._myReplica][1])
+        else:
+            random.random() # Do one random number to stay in sync with the master process for the same replica.
+            
+
+            
+ 
+if __name__ == "__main__":
+    
+    n = int(sys.argv[1])
+    T = [250,  275,  300,  330,  370,  420]
+    
+    myPT = ParallelTempering(T,  [])
+    myPT.exchange()
+    myPT.exchange()
+
+
+        
+        
+        
+
+         
Index: pysmmp/algorithms.py
===================================================================
--- pysmmp/algorithms.py	(revision 225475cd59f91e8a64cbbf82f4e34e4484191092)
+++ pysmmp/algorithms.py	(revision 225475cd59f91e8a64cbbf82f4e34e4484191092)
@@ -0,0 +1,91 @@
+# algorithms.py
+#
+# Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
+#                      Jan H. Meinke, Sandipan Mohanty
+#
+import copy
+import random
+import smmp_p as smmp
+from math import *
+from universe import *
+from protein import *
+
+class CanonicalMonteCarlo(object):
+    """An implementation of a canonical Monte Carlo run using Metropolis weights.
+    After a Monte-Carlo sweep over all internal variables, this implementation 
+    calculates the interactions between all Protein s in the Universe
+    """
+    
+    def __init__(self, myUniverse, nequi, sweeps, weight =  lambda x : x * smmp.bet.beta):
+        self.__defSweeps = sweeps
+        self.__defEqui = nequi
+        self.__univ = myUniverse
+        self.__acc = 0
+        self.__performedSweeps = 0
+        self.__minEnergy = myUniverse.energy()
+        self.__l = 0.1
+        self.__minCounter = 0
+        self.__weight = weight
+        
+    ## Fraction of accepted Monte Carlo moves of the total number of attempted moves.
+    #  @pre At least one move must have been performed before this function can be
+    #       called
+    #  @return fraction of accepted moves as a number between 0 and 1.
+    def acceptanceRate(self):
+        return float(self.__acc) / (self.__performedSweeps * smmp.mol_par.nvr)
+    ## The lowest energy ever seen at the end of a Monte Carlo move.
+    #  @return lowest energy seen at the end of a sweep, sweep
+    def minimumEnergy(self):
+        return self.__minEnergy, self.__minCounter
+    ## Number of sweeps performed so far
+    def performedSweeps(self):
+        return self.__performedSweeps
+        
+    def run(self):
+        """A complete Monte Carlo run including initial equilibration."""
+        self.sweeps(self.__defSweeps)
+        
+    def equilibrate(self):
+        """Do nequi sweeps to eqilibrate the system."""
+        self.sweeps(self.__defEqui)
+        
+    def sweeps(self, n):
+        """Perform n sweeps."""
+        for i in range(0, n):
+            self.sweep()
+        
+    def sweep(self):
+        """Performs a single Monte Carlo sweep."""
+        eol = self.__univ.energy()
+        #eol, self.__acc = smmp.metropolis(eol, self.__acc, self.__weight)
+        eol = self.metropolis()
+        self.__performedSweeps += 1
+        if eol < self.__minEnergy:
+            self.__minEnergy = eol
+            self.__minCounter = self.__performedSweeps
+            smmp.outpdb(0, "best.pdb")
+
+    def metropolis(self):
+        """Implementation of the Metropolis algorithm."""
+        eol = self.__univ.energy()
+        # Update the internal degrees of freedom first
+        for i in range(0, smmp.mol_par.nvr):
+            jv = smmp.var_i.idvr[i] - 1
+            vrol = smmp.var_r.vlvr[jv]
+            dv=smmp.var_r.axvr[jv] * (random.random() - 0.5)
+            smmp.var_r.vlvr[jv] = (vrol + dv + pi) % (2 * pi) - pi
+            etrial = self.__univ.energy()
+            dW = self.__weight(etrial) - self.__weight(eol)
+            if dW < 0:
+                eol = etrial
+                self.__acc += 1
+            else:
+                rd =random.random()
+                dW = min(dW, 100)
+                if  exp(-dW) > rd:
+                    eol = etrial
+                    self.__acc +=1
+                else:
+                    smmp.var_r.vlvr[jv] = vrol
+        # Update the global degrees of freedom
+        return eol
Index: pysmmp/interactionscan.py
===================================================================
--- pysmmp/interactionscan.py	(revision 225475cd59f91e8a64cbbf82f4e34e4484191092)
+++ pysmmp/interactionscan.py	(revision 225475cd59f91e8a64cbbf82f4e34e4484191092)
@@ -0,0 +1,100 @@
+# -*- coding: utf-8 -*-
+'''
+Created on Aug 31, 2010
+
+@author: Jan H. Meinke
+'''
+
+import smmp_p as smmp
+import numpy as np
+from math import sqrt, cos
+
+def scan():
+    """Iterate through all interaction pairs and return lists of 
+    interaction partners for each type of interaction.
+    
+    This implementation is almost equivalent to the Fortran 
+    implementation in enyshe.f.
+    
+    """
+    cutoff = 32
+    longDistancePartners = {}
+    oneFourPartners = {}
+    
+    EvdW = 0
+    Ec = 0
+    Ehb = 0
+    Etor = 0
+    ifivr = smmp.mol_i.ivrml1[0]
+    ntlvr = smmp.mol_par.nvr
+    i1s = smmp.mol_i.imsml1[0] + smmp.mol_i.nmsml[0]
+    # Loop over variables in reverse order.
+    for iv in smmp.var_i.iorvr[ifivr + ntlvr - 2:ifivr:-1] - 1:
+        ia = smmp.var_i.iatvr[iv] - 1
+        ic = smmp.var_i.iclvr[iv] - 1
+        if smmp.var_i.ityvr[iv] == 3:
+            Etor += smmp.epar_r.e0to[ic] * (1.0 + smmp.epar_r.sgto[ic] *
+                cos(smmp.atm_r.toat[ia] * smmp.epar_r.rnto[ic]))
+            
+        i2s = i1s - 1
+        i1s = smmp.var_i.imsvr1[iv - 1]
+        for i1, i2 in zip(smmp.var_i.latms1[i1s:i2s] - 1, 
+                          smmp.var_i.latms2[i1s:i2s]):
+            for i, ity, cqi, xi, yi, zi, ivw1, ivw2, i141, i142 in zip(
+                    np.arange(i1, i2), smmp.atm_i.ityat[i1:i2], 
+                    smmp.epar_r.conv * smmp.atm_r.cgat[i1:i2],
+                    smmp.atm_r.xat[i1:i2], 
+                    smmp.atm_r.yat[i1:i2], smmp.atm_r.zat[i1:i2], 
+                    smmp.eny_i.ivwat1[i1:i2] - 1, smmp.eny_i.ivwat2[i1:i2],
+                    smmp.eny_i.i14at1[i1:i2] - 1,smmp.eny_i.i14at1[i1:i2]):
+                longDistancePartners[i] = []
+                oneFourPartners[i] = []
+                for j, jty, xij, yij, zij, cqj in zip(np.arange(ivw1, ivw2), 
+                        smmp.atm_i.ityat[ivw1:ivw2],
+                        smmp.atm_r.xat[ivw1:ivw2] - xi, 
+                        smmp.atm_r.yat[ivw1:ivw2] - yi,
+                        smmp.atm_r.zat[ivw1:ivw2] - zi,
+                        smmp.atm_r.cgat[ivw1:ivw2]):
+                    if abs(i - j) < cutoff:
+                        longDistancePartners[i].append(j)
+                    rij2 = xij * xij + yij * yij + zij * zij
+                    rij4 = rij2 * rij2
+                    rij6 = rij4 * rij2
+                    rij = sqrt(rij2)
+                    Ec += cqi * cqj / rij
+                    if smmp.epar_i.ihbty[ity, jty]:
+                        Ehb += smmp.epar_r.ahb[ity, jty] / (rij6 * rij6) -\
+                            smmp.epar_r.chb[ity, jty] / (rij6 * rij4)
+                    else:
+                        EvdW += smmp.epar_r.aij[ity, jty] / (rij6 * rij6) -\
+                            smmp.epar_r.cij[ity, jty] / rij
+                
+                for j in smmp.eny_i.l14at[i141:i142]:
+                    jty = smmp.atm_i.ityat[j]
+                    xij = smmp.atm_r.xat[j] - xi
+                    yij = smmp.atm_r.xat[j] - yi
+                    zij = smmp.atm_r.xat[j] - zi
+                    rij2 = xij * xij + yij * yij + zij * zij
+                    rij4 = rij2 * rij2
+                    rij6 = rij4 * rij2
+                    rij = sqrt(rij2)
+                    Ec += cqi * cqj / rij
+                    if smmp.epar_i.ihbty[ity, jty]:
+                        Ehb += smmp.epar_r.ahb[ity, jty] / (rij6 * rij6) -\
+                            smmp.epar_r.chb[ity, jty] / (rij6 * rij4)
+                        longDistancePartners[i].append(j)
+                    else:
+                        oneFourPartners[i].append(j)
+                        EvdW += smmp.epar_r.a14[ity, jty] / (rij6 * rij6) -\
+                           smmp.epar_r.cij[ity, jty] / rij
+                           
+                if not longDistancePartners[i]:
+                    del(longDistancePartners[i])
+                if not oneFourPartners[i]:
+                    del(oneFourPartners[i])
+
+                    
+    print EvdW + Ec + Ehb + Etor, EvdW, Ec, Ehb, Etor
+    
+    return longDistancePartners, oneFourPartners
+                 
Index: pysmmp/main.py
===================================================================
--- pysmmp/main.py	(revision 225475cd59f91e8a64cbbf82f4e34e4484191092)
+++ pysmmp/main.py	(revision 225475cd59f91e8a64cbbf82f4e34e4484191092)
@@ -0,0 +1,114 @@
+#!/usr/bin/env python
+# main.py
+#
+# Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
+#                      Jan H. Meinke, Sandipan Mohanty
+
+## \mainpage pySMMP: Python bindings to SMMP v. 1.1.
+#    SMMP is a FORTRAN library for simulating proteins. It provides 
+#    minimization routines and various Monte Carlo methods. SMMP uses an all-atom
+#    representation of the protein. It stores the configuration within the 
+#    standard geometry model. In the standard geometry model the internal degrees
+#    of freedom, e.g., the dihedral angles, the bond lengths, etc. are stored and
+#    the cartesian coordinates are calculated on the fly. For more details on 
+#    the model SMMP uses and the functions it provides see the manual (manual.ps).
+#   
+#    pySMMP builds on SMMP. It provides python bindings to the functions and 
+#    global variables, but it uses classes to represent the Proteins and 
+#    algorithms. 
+#    
+#    \section sec_tutorial A Tutorial: A canonical Monte Carlo simulation of protein A
+#
+#    The basic setup is the same for most simulations. This tutorial takes you
+#    step by step to a running simulation. First we need to import the library
+#    \code 
+#         import smmp
+#     \endcode
+#    This makes the smmp routines available to the Python interpreter. Next we
+#    create the universe by initializing a universe.Universe object by first 
+#    importing the universe package and then setting the object:
+#    \code
+#        import universe
+#        myUniverse = universe.Universe(T=300)  
+#    \endcode
+#    This initializes the smmp library as well. We also set the initial 
+#    temperature of the universe to 300K. The Universe object initializes the 
+#    SMMP library and keeps track of global parameters, e.g, the temperature T.
+#    After you created a Universe, you can initialize a Protein. Proteins can
+#    be read from a sequence and a variable file or from a PDB file.
+#    \code
+#        import protein
+#        
+#        protA = protein.Protein("1bdd.seq", "1bdd.var")
+#        myUniverse.add(protA)
+#    \endcode
+#    Adding the protein to myUniverse makes the universe aware of it.
+#    Now you are ready to run a simulation. You can also querry the protein for 
+#    various properties, e.g., 
+#    \code
+#        print protA.rgyr()
+#    \endcode
+#    prints the radius of gyration of the protein.
+#    A canonical Monte Carlo algorithm is available in the algorithms package.
+#    \code
+#        import algorithms
+#        myMC = algorithms.CanonicalMonteCarlo(myUniverse, 1000, 10000)
+#    \endcode
+#    This sets up a simulation with 1000 steps for equilibration and a default 
+#    length of 10000 sweeps.
+#    The quickest way to run the simulation is by calling the run method of the 
+#    algorithm object myMC to do the entire Monte Carlo run.
+#    \code
+#        myMC.run()
+#    \endcode
+
+import sys, time
+from optparse import OptionParser
+from protein import *
+from universe import *
+from algorithms import CanonicalMonteCarlo
+
+if __name__ == "__main__":
+# Set up command line options    
+    usage = "usage: %prog [options] seq-filename"
+    parser = OptionParser(usage = usage, version="%prog 20050805")
+    parser.add_option("-v", "--variables", dest="varFileName",
+                  help="read amino-acid conformation from FILE", metavar="FILE")
+    parser.add_option("--equi", type="int", dest="nequi", default=100,
+                  help="number of sweeps for equilibration", metavar="n")
+    parser.add_option("--sweeps", type="int", dest="sweeps", default=10000,
+                  help="number of sweeps for annealing", metavar="n")
+    parser.add_option("--mes", type="int", dest="mes", default=10,
+                  help="number of sweeps between measurements", metavar="n")
+    parser.add_option("--Tmax", type="float", dest="tmax", default=1000,
+                  help="maximum temperature for annealing", metavar="T")
+    parser.add_option("--Tmin",  type="float", dest="tmin", default=100, 
+                  help="maximum temperature for annealing",metavar="T")
+
+    (options, args) = parser.parse_args()
+    if len(args) < 1:
+        parser.print_help()
+        sys.exit()
+    seqFileName = args[0]
+    
+# Start simulation setup.    
+    myUniv = Universe(T=300, st=0)   
+    p1 = Protein(seqFileName, options.varFileName)
+    p1.setOrigin((0, 0, 0), (0, 0, 0))
+    myUniv.add(p1)
+    T = 300.0
+    smmp.outpdb(1, "start.pdb")
+    print "E of starting configuration.", myUniv.energy()
+    ts = time.time()
+    smmp.minim(1, options.sweeps, 1.0E-12)
+
+    smmp.outpdb(1, "equi.pdb")
+    
+    t1 = time.clock()
+    print myUniv.energy(), myUniv.rgyr(), myUniv.helix(), p1.hbond(), (t1-ts)
+    p1.saveVariables("final.var")
+    smmp.outpdb(0, "final.pdb")
+    te = time.time()
+    print "The calculation took %ss." % (te - ts)
+ 
+    
Index: pysmmp/pPySMMP.py
===================================================================
--- pysmmp/pPySMMP.py	(revision 225475cd59f91e8a64cbbf82f4e34e4484191092)
+++ pysmmp/pPySMMP.py	(revision 225475cd59f91e8a64cbbf82f4e34e4484191092)
@@ -0,0 +1,40 @@
+# -*- coding: utf-8 -*-
+import smmp_p as smmp
+import mpi4py
+import mpi4py.MPI
+import time, random
+from universe import *
+from protein import *
+from algorithms import *
+
+smmp.paral.my_mpi_comm = mpi4py.MPI.COMM_WORLD.py2f()
+smmp.paral.myrank = mpi4py.MPI.COMM_WORLD.Get_rank()
+numberOfTasks = mpi4py.MPI.COMM_WORLD.Get_size()
+smmp.paral.no = numberOfTasks
+#smmp.init_energy("SMMP/")
+#smmp.init_molecule(1,'nh3+','coo-',seqfile='EXAMPLES/1bdd.seq',varfile='EXAMPLES/1bdd.var')
+myU = Universe(st=1)
+myP = Protein('EXAMPLES/1bdd.seq', 'EXAMPLES/1bdd.var')
+myU.add(myP)
+smmp.distributeworkload(numberOfTasks,1)
+E = myU.energy()
+if smmp.paral.myrank == 0:
+    print "The starting energy is %s." % E
+
+random.seed(82364)
+#myMC = CanonicalMonteCarlo(myU, 20, 20)
+
+t0 = time.time()
+for i in xrange(1000):
+    myU.energy()
+
+#myMC.run()
+t1 = time.time()
+E = myU.energy()
+if smmp.paral.myrank == 0:
+#    Emin = myMC.minimumEnergy()
+    print "The final energy is: %f. Calculated in %fs" % (E, (t1 -t0))
+#    print "The minimum energy found was %.2f at sweep %s." % Emin
+    
+print random.random()
+
Index: pysmmp/protein.py
===================================================================
--- pysmmp/protein.py	(revision 225475cd59f91e8a64cbbf82f4e34e4484191092)
+++ pysmmp/protein.py	(revision 225475cd59f91e8a64cbbf82f4e34e4484191092)
@@ -0,0 +1,290 @@
+#!/usr/bin/env python
+# protein.py
+#
+# Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
+#                      Jan H. Meinke, Sandipan Mohanty
+#
+import sys
+import numpy as np
+import smmp_p as smmp
+from math import *
+
+crd = 57.2957795130823
+
+def nursvr(ivr):
+    res = 0
+    for i in range(smmp.mol_i.ntlml - 1, -1, -1):
+        ifirs = smmp.mol_i.irsml1[i] - 1
+        if ivr + 1 > smmp.res_i.ivrrs1[ifirs] and smmp.mol_i.nvrml[i] > 0:
+            for j in range(smmp.mol_i.irsml2[i] - 1, ifirs - 1, -1):
+                if ivr + 1 > smmp.res_i.ivrrs1[j]  and smmp.res_i.nvrrs[j] > 0:
+                    res = j + 1
+                    return res
+    return res
+
+## Provides read-only access to SMMP's atoms. Usually, you won't initialize
+#  an Atom directly since you need to know the internal array index of the
+#  Atom in SMMP. You can get a list of all the Atoms in a Protein by using 
+#  protein.Protein.getAtoms().
+#  Various properties are accessible. See the functions' documentations for
+#  details.
+class Atom:
+    """Atoms are the basic chemical units."""
+    #TODO: insert namespace
+    
+    ## Initialize the atom with its ID. The ID is the array index of the atom
+    #  in SMMP. 
+    #       @param id index of atom in SMMP
+    def __init__(self, id):
+        self.__id = id
+    
+    def position(self):
+        """Cartesian coordinates of this atom."""
+        res = (smmp.atm_r.xat[self.__id], smmp.atm_r.yat[self.__id], smmp.atm_r.zat[self.__id])
+        return res
+    
+    def charge(self):
+        """Partial charge of this atom."""
+        return smmp.atm_r.cgat[self.__id]
+    
+    def radius(self):
+        """Van-der-Waals Radius of this atom."""
+        return smmp.solvent.rvdw[self.__id]
+        
+    def name(self):
+        return ''.join(smmp.atm_c.nmat[:,0].reshape([-1,4])[self.__id]).strip()
+        
+    def __str__(self):
+        return "%s at (%s, %s, %s)" % (self.name(), self.position()[0],  self.position()[1], self.position()[2])
+
+## Amino acids bound within a protein are often called residues. Strictly
+#  speaking the residue is the part that's connected to the Calpha and not part
+#  of the backbone, but we'll include the backbone atoms as well. Just like Atoms
+# Residues reference the corresponding data in SMMP and are not independent. 
+class Residue:
+    """An amino acid within a Protein."""
+    
+    ## Initialize the residue with the corresponding index in SMMP.
+    #  @param id index of residue in SMMP
+    def __init__(self, id):
+        self.__id = id
+
+    def atoms(self):
+        """Returns the list of atoms belonging to this residue."""
+        first = smmp.res_i.iatrs1[self.__id]
+        last = smmp.res_i.iatrs2[self.__id]
+        res = [Atom(i) for i in range(first - 1, last - 1)]
+        return res
+    
+        
+## Proteins are constructed from a sequence file and a configuration file or 
+#  read from a PDB file.
+class Protein:
+    """A protein is a polypeptide made of the 20 naturally occuring amino 
+    acids.
+    """
+    
+    def __init__(self, seq, var=' '):
+        """The structure of the protein is read in from a file. 
+        Pass the name of a textfile with a sequence string and--optionally--the name of a
+        variable file that describes the conformation. Or give the name of pdb
+        file---it must end with 'pdb'---that contains both information."""
+        
+        grpn = 'nh2'
+        grpc = 'cooh'
+        prevID = smmp.mol_i.ntlml - 1
+        if seq[-3:].lower() == 'pdb':
+            if smmp.mol_i.ntlml:
+                print "Sorry, but I can only read the first protein from a pdb file." 
+                raise ValueError
+            smmp.init_molecule(0, grpn, grpc, seq)
+        else:
+            smmp.init_molecule(1, grpn, grpc, seq, var)
+        self.__id = smmp.mol_i.ntlml - 1
+        if self.__id - prevID > 1:
+            print "Sorry, I don't know how to deal with multiple proteins in one sequence file, yet."
+            raise NotImplementedError
+    
+    def setOrigin(self, x, a):
+        """Sets the position of the N terminal of the protein in the global 
+        coordinate system and the orientation of the protein. Expects two 
+        tuples. The first contains the coordinates of the first atom (x, y,
+        and z) and the second contains three angles.
+        """
+        #TODO: Implement type checking
+        smmp.mol_r.gbpr[0][self.__id] = x[0]
+        smmp.mol_r.gbpr[1][self.__id] = x[1]
+        smmp.mol_r.gbpr[2][self.__id] = x[2]
+        smmp.mol_r.gbpr[3][self.__id] = a[0]
+        smmp.mol_r.gbpr[4][self.__id] = a[1]
+        smmp.mol_r.gbpr[5][self.__id] = a[2]
+
+    def origin(self):
+        r = [smmp.mol_r.gbpr[i][self.__id] for i in range(0, 3)]
+        a = [smmp.mol_r.gbpr[i][self.__id] for i in range(3, 6)]
+        return r, a
+    
+    def getDihedral(self, angle, residue):
+        try:
+            res = int(residue)
+            idx = res
+        except ValueError:
+            res = residue.split()
+            res[0] = res[0].lower()
+            res[1] = int(res[1])
+            idx = self.findResidue(res[0], res[1])
+        varIdx = ''
+        for i in range(smmp.res_i.ivrrs1[idx], smmp.res_i.ivrrs1[idx] + smmp.res_i.nvrrs[idx]):
+            varName = ''.join(list([str(smmp.var_c.nmvr[j][i]) for j in range(0, 3)]))
+            if varName == angle:
+                varIdx = i
+                break
+        if varIdx:
+            return smmp.var_r.vlvr[varIdx] * crd
+        raise KeyError, '%s not found for %s.' % (angle, residue)
+    
+    ## Residues are given by their 3-letter abbreviation and their index. 
+    #  For example, findResidue('gly', 5) returns the index of the fifth occurence
+    #  of glycine in the protein. If the residue is not found an exception is raised.
+    #
+    #  @param name 3-letter abbreviation of the residue
+    #  @param idx index of occurence of the amino acid given in name.
+    
+    def findResidue(self, name, idx = 1):
+        """Find the idxth residue of type name and return the index."""
+        res = ''.join(self.residues()).lower()
+        name = name.lower()
+        j = 0
+        try:
+            for i in range(0, idx):
+                k = res.index(name, j)
+                j = j + k
+        except:
+            raise KeyError, "%s%s not found." % (name, idx)
+        if j >=len(res):
+            raise KeyError, "%s%s not found." % (name, idx)
+        return j / 3
+        
+    def atoms(self):
+        """Returns the list of atoms belonging to this protein."""
+        first = smmp.res_i.iatrs1[smmp.mol_i.irsml1[self.__id] - 1]
+        last = smmp.res_i.iatrs2[smmp.mol_i.irsml2[self.__id] - 1]
+        res = [Atom(i) for i in range(first - 1, last - 1)]
+        return res
+    
+    def residues(self):
+        """Returns a list of the amino acids in this protein."""
+        seq = [''.join(l) for l in np.array([str(c[0]) for c in smmp.res_c.seq]).reshape([-1, 4])]
+        first = smmp.mol_i.irsml1[self.__id]- 1
+        last = smmp.mol_i.irsml2[self.__id] - 1
+        return seq[first:last]
+        
+    def id(self):
+        return self.__id
+        
+    def strand(self):
+        """Returns the number of residues in a beta-strand-like configuration."""
+        maxDev = 30.0
+        phiSheet = -150.0
+        psiSheet = 150.0
+        strandCount = 0
+        ifivr = smmp.mol_i.ivrml1[self.__id]
+        i = 0
+        iv = smmp.var_i.idvr[i]
+        while iv < ifivr:
+            i += 1
+            iv = smmp.var_i.idvr[i]
+        for j in range(i, ifivr + smmp.mol_i.nvrml[self.__id]):
+            iv = smmp.var_i.idvr[j] - 1
+            name =''. join([str(smmp.var_c.nmvr[iv * 3 + k][0]) for k in range(0,3)])
+            if name == 'phi':
+                phi = smmp.var_r.vlvr[iv] * crd
+                psi = smmp.var_r.vlvr[smmp.var_i.idvr[j+1] - 1] * crd
+                if abs(phi - phiSheet) < maxDev and abs(psi - psiSheet) < maxDev:
+                    strandCount += 1
+        return strandCount
+    
+    def directionVector(self):
+        """Calculates the normalized end-to-end vector for this molecule. 
+        The first and last c_alpha carbon is used as reference."""
+        atoms = self.atoms()
+        beginning = []
+        for a in atoms:
+            if a.name().strip() == 'ca':
+                if beginning:
+                    end = a.position()
+                else:
+                    beginning = a.position()
+        res = [end[i] - beginning[i] for i in range(0,3)] 
+        sum = sqrt(reduce(lambda x, y: x + y ** 2, res))
+        res = map(lambda x: x / sum, res)
+        return res
+        
+    def contacts(self):
+        """Returns a triplet containing the total number of contacts, the number
+        of native contacts, and the hemming distance. The latter two quantities 
+        are only meaningful if a reference structure has been defined.
+        """
+        return smmp.contacts()
+
+    def hbond(self):
+        return smmp.hbond(self.__id + 1, 0)
+        
+    def energy(self):
+        return smmp.energy()
+    
+    def rgyr(self):
+        return smmp.rgyr(self.__id + 1) 
+    
+    def rmsdfun(self):
+        return smmp.rmsdfun()
+        
+    def savePDB(self, file=""):
+       """Write the current configuration to file using PDB formatting."""
+       smmp.outpdb(self.__id, file)
+       
+    def saveVariables(self, file=""):
+        """Output the current conformation in the variable format. The output 
+        can be read in again and used to restart/continue a simulation.
+        
+        @param file if empty write to standard out
+        """
+        # Same functionality as smmp.outvar(self.__id + 1, 0)
+        if not file:
+            out = sys.stdout
+        else:
+            out = open(file, 'w')
+        for i in range(smmp.mol_i.ivrml1[self.__id] - 1 , smmp.mol_i.ivrml1[self.__id] + smmp.mol_i.nvrml[self.__id] - 1):
+            if not smmp.var_l.fxvr[i]:
+                name = ''.join([str(smmp.var_c.nmvr[i * 3 + j][0]) for j in range(0,3)])
+                print >>out, nursvr(i), ':', name, ':', smmp.var_r.vlvr[i] * crd
+            else:
+                it=smmp.var_i.ityvr[i]
+                if it == 3:
+                    print >>out, nursvr(i), ':', name, ':', smmp.var_r.vlvr[i] * crd, ' &'
+        
+    def zimmer(self):
+        first = smmp.mol_i.irsml1[self.__id] - 1
+        last = smmp.mol_i.irsml2[self.__id]
+        nrs = smmp.mol_i.irsml2[smmp.mol_i.ntlml - 1] - 1
+        smmp.zimmer(last)
+        return smmp.zimme.zimm[first:last]
+
+    def __len__(self):
+        return smmp.mol_i.irsml2[self.__id]  - smmp.mol_i.irsml1[self.__id] + 1
+
+
+if __name__ == "__main__":
+    libdir = 'SMMP/'
+
+    smmp.epar_l.flex = 0
+    smmp.epar_l.sh2 = 0
+    smmp.epar_l.epsd = 0
+    smmp.epar_l.tesgrd = 0
+    smmp.isolty.itysol = 0
+    smmp.init_energy(libdir)
+    
+    p1 = Protein('EXAMPLES/1bdd.pdb')
+    for i in p1.atoms():
+        print i
+    print p1.getDihedral('phi', 11)
Index: pysmmp/universe.py
===================================================================
--- pysmmp/universe.py	(revision 225475cd59f91e8a64cbbf82f4e34e4484191092)
+++ pysmmp/universe.py	(revision 225475cd59f91e8a64cbbf82f4e34e4484191092)
@@ -0,0 +1,116 @@
+# universe.py
+#
+# Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
+#                      Jan H. Meinke, Sandipan Mohanty
+#
+import smmp_p as smmp
+from math import *
+
+class Universe:
+    """Describes the environment of the simulation. Including the temperature,
+    the solvent, and the force field. This class keeps track of all the 
+    molecules.
+    """
+    
+    def __init__(self, T = 300, ff = 'ecepp3', st = 1, libdir = 'SMMP/'):
+        self.setTemperature(T)
+        if ff == 'ecepp3':
+            smmp.epar_l.flex = 0
+            smmp.epar_l.sh2 = 0
+            smmp.epar_l.epsd = 0
+            smmp.epar_l.ientyp=0
+        elif ff == 'ecepp2':
+            smmp.epar_l.flex = 0
+            smmp.epar_l.sh2 = 1
+            smmp.epar_l.epsd = 0
+            smmp.epar_l.ientyp=0
+        elif ff == 'flex':
+            smmp.epar_l.flex = 1
+            smmp.epar_l.sh2 = 0
+            smmp.epar_l.epsd = 0
+            smmp.epar_l.ientyp=0
+        smmp.epar_l.tesgrd = 0
+        smmp.isolty.itysol = st
+        smmp.init_energy(libdir)
+        self.__objects = []
+        self.__boxSize = 10
+        smmp.updchois.upchswitch = 0
+        smmp.updchois.rndord = 0
+        smmp.updchois.bgsprob = 0
+        
+    def add(self, m):
+        self.__objects.append(m)
+        
+    def energy(self):
+        """Return the total energy of this universe."""
+        return smmp.energy()
+    
+    def rgyr(self):
+        """Radius of gyration of the entire system."""
+        return smmp.rgyr(0)[0]
+        
+    def endToEnd(self):
+        return smmp.rgyr(0)[1]
+        
+    def hbond(self):
+        """Returns the total number of hydrogen bonds in the system"""
+        mhb = 0
+        for i in self.__objects:
+            mhb = mhb + i.hbond()
+        mhb = mhb + smmp.interhbond()
+        return mhb
+    
+    def helix(self):
+        return smmp.helix()
+    
+    def sheet(self):
+        """Checks if we have any parallel or antiparallel beta-sheets in the
+        system.
+        """
+        minStrandContent = 2
+        minHydrogenBonds = 3
+        res = []
+        for i in range(0, len(self.__objects)):
+            if self.__objects[i].strand() >= minStrandContent:
+                for j in range(i, len(self.__objects)):
+                    if smmp.h_bond.mmhb[self.__objects[i].id()][self.__objects[j].id()] >= minHydrogenBonds:
+                        if self.__objects[j].strand() >= minStrandContent:
+                            print "Found a candidate for a sheet between %s and %s." % (i, j)
+                            dir1 = self.__objects[i].directionVector()
+                            dir2 = self.__objects[j].directionVector()
+                            sum = 0
+                            for k in range(0, 3):
+                                sum += dir1[k] * dir2[k]
+                            if sum > 0.7:
+                                res.append([i, j, 1])
+                            elif sum < -0.7:
+                                res.append([i, j, -1])
+        return res 
+    
+    def temperature(self):
+        """Returns the current temperature of this Universe."""
+        return 1.0 / ( smmp.bet.beta * 1.98773e-3 )
+    ## @brief Sets the temperature value of this Universe.
+    #  Be aware that SMMP is not thread save. Changing the temperature may
+    #  lead to strange behavior in the current sweep.
+    #  @pre T > 0
+    #
+    #  @param T temperature in Kelvin.
+    def setTemperature(self, T):
+        if T> 0:
+            smmp.bet.beta = 1.0 / ( T * 1.98773e-3 )
+        
+    def boxSize(self):
+        return self.__boxSize
+
+    def objects(self):
+        """Returns the list of interacting things in the universe."""
+        res = self.__objects
+        return res
+        
+    def save(self):
+        """Saves the state of this Universe."""
+        pass
+
+if __name__ == '__main__':
+    u=Universe()
Index: iverse.py
===================================================================
--- universe.py	(revision 5fef0d738f6c69046859b9b41da3ebeefc2efa61)
+++ 	(revision )
@@ -1,116 +1,0 @@
-# universe.py
-#
-# Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
-#                      Jan H. Meinke, Sandipan Mohanty
-#
-import smmp_p as smmp
-from math import *
-
-class Universe:
-    """Describes the environment of the simulation. Including the temperature,
-    the solvent, and the force field. This class keeps track of all the 
-    molecules.
-    """
-    
-    def __init__(self, T = 300, ff = 'ecepp3', st = 1, libdir = 'SMMP/'):
-        self.setTemperature(T)
-        if ff == 'ecepp3':
-            smmp.epar_l.flex = 0
-            smmp.epar_l.sh2 = 0
-            smmp.epar_l.epsd = 0
-            smmp.epar_l.ientyp=0
-        elif ff == 'ecepp2':
-            smmp.epar_l.flex = 0
-            smmp.epar_l.sh2 = 1
-            smmp.epar_l.epsd = 0
-            smmp.epar_l.ientyp=0
-        elif ff == 'flex':
-            smmp.epar_l.flex = 1
-            smmp.epar_l.sh2 = 0
-            smmp.epar_l.epsd = 0
-            smmp.epar_l.ientyp=0
-        smmp.epar_l.tesgrd = 0
-        smmp.isolty.itysol = st
-        smmp.init_energy(libdir)
-        self.__objects = []
-        self.__boxSize = 10
-        smmp.updchois.upchswitch = 0
-        smmp.updchois.rndord = 0
-        smmp.updchois.bgsprob = 0
-        
-    def add(self, m):
-        self.__objects.append(m)
-        
-    def energy(self):
-        """Return the total energy of this universe."""
-        return smmp.energy()
-    
-    def rgyr(self):
-        """Radius of gyration of the entire system."""
-        return smmp.rgyr(0)[0]
-        
-    def endToEnd(self):
-        return smmp.rgyr(0)[1]
-        
-    def hbond(self):
-        """Returns the total number of hydrogen bonds in the system"""
-        mhb = 0
-        for i in self.__objects:
-            mhb = mhb + i.hbond()
-        mhb = mhb + smmp.interhbond()
-        return mhb
-    
-    def helix(self):
-        return smmp.helix()
-    
-    def sheet(self):
-        """Checks if we have any parallel or antiparallel beta-sheets in the
-        system.
-        """
-        minStrandContent = 2
-        minHydrogenBonds = 3
-        res = []
-        for i in range(0, len(self.__objects)):
-            if self.__objects[i].strand() >= minStrandContent:
-                for j in range(i, len(self.__objects)):
-                    if smmp.h_bond.mmhb[self.__objects[i].id()][self.__objects[j].id()] >= minHydrogenBonds:
-                        if self.__objects[j].strand() >= minStrandContent:
-                            print "Found a candidate for a sheet between %s and %s." % (i, j)
-                            dir1 = self.__objects[i].directionVector()
-                            dir2 = self.__objects[j].directionVector()
-                            sum = 0
-                            for k in range(0, 3):
-                                sum += dir1[k] * dir2[k]
-                            if sum > 0.7:
-                                res.append([i, j, 1])
-                            elif sum < -0.7:
-                                res.append([i, j, -1])
-        return res 
-    
-    def temperature(self):
-        """Returns the current temperature of this Universe."""
-        return 1.0 / ( smmp.bet.beta * 1.98773e-3 )
-    ## @brief Sets the temperature value of this Universe.
-    #  Be aware that SMMP is not thread save. Changing the temperature may
-    #  lead to strange behavior in the current sweep.
-    #  @pre T > 0
-    #
-    #  @param T temperature in Kelvin.
-    def setTemperature(self, T):
-        if T> 0:
-            smmp.bet.beta = 1.0 / ( T * 1.98773e-3 )
-        
-    def boxSize(self):
-        return self.__boxSize
-
-    def objects(self):
-        """Returns the list of interacting things in the universe."""
-        res = self.__objects
-        return res
-        
-    def save(self):
-        """Saves the state of this Universe."""
-        pass
-
-if __name__ == '__main__':
-    u=Universe()
