beadModelTwist.py

            ## Write out bonds
            bonds = self.bonds
            fh.write("{:>8d} !NBOND\n".format(len(bonds)))
            counter = 0
            for n1,n2,pot in bonds:
                fh.write( "{:d} {:d} ".format(n1.idx+1,n2.idx+1) )
                counter += 1
                if counter == 3:
                    fh.write("\n")
                    counter = 0
            fh.write("\n")

        return

    def writeArbdFiles(self, prefix, numSteps=100000000, timestep=100e-6):
        ## TODO: save and reference directories and prefixes using member data
        d = "potentials"
        self._writeArbdCoordFile( prefix + ".coord.txt" )
        self._writeArbdBondFile(  prefix, directory = d )
        self._writeArbdAngleFile( prefix, directory = d )
        self._writeArbdDihedralFile( prefix, directory = d )
        self._writeArbdExclFile(  prefix + ".excludes.txt" )
        self._writeArbdPotentialFiles( prefix, directory = d )
        self._writeArbdConf( prefix, numSteps, timestep, "%s/%s-" % (d,prefix) )
        
    def _writeArbdCoordFile(self, filename):
        with open(filename,'w') as fh:
            for n,hid,zid in self.particles:
                fh.write("%f %f %f\n" % tuple(x for x in n.position))
        
    def _writeArbdConf(self, prefix, numSteps=100000000, timestep=100e-6, potentialPrefix='' ):
        ## TODO: raise exception if _writeArbdPotentialFiles has not been called
        filename = "%s.bd" % prefix
        with open(filename,'w') as fh:
            fh.write("""seed 1234
timestep %f
steps %d
numberFluct 0
interparticleForce 1
fullLongRange 0
temperature 291
electricField 0.0

outputPeriod 1000
outputEnergyPeriod 1000
outputFormat dcd

decompPeriod 100000
cutoff 40.0
pairlistDistance 80
""" % (timestep, numSteps))
            
            for x in self.getParticleTypesAndCounts():
                fh.write("\nparticle %s\nnum %d\n" % x)
                ## TODO: look up better values in dictionary for particle types
                fh.write("gridFile null.dx\ndiffusion 150\n")

            fh.write("\ninputCoordinates %s.coord.txt\n"  % prefix )
            if os.path.exists("test.0.restart"):
                fh.write("restartCoordinates test.0.restart\n" )
            fh.write("""\n## Interaction potentials
tabulatedPotential  1
## The i@j@file syntax means particle type i will have NB interactions with particle type j using the potential in file
""")
            for pair,f in zip(self._particleTypePairIter(), self._nbParamFiles):
                i,j,t1,t2 = pair
                fh.write("tabulatedFile %d@%d@%s\n" % (i,j,f))


            fh.write("\n")
            for f in list(set([b[-1].filename(potentialPrefix) for b in self.bonds])):
                fh.write("tabulatedBondFile %s\n" % f)

            fh.write("\n")
            for f in list(set([b[-1].filename(potentialPrefix) for b in self.angles])):
                fh.write("tabulatedAngleFile %s\n" % f)

            fh.write("\n")
            for f in list(set([b[-1].filename(potentialPrefix) for b in self.dihedrals])):
                fh.write("tabulatedDihedralFile %s\n" % f)

            fh.write("""\n## Files that specify connectivity of particles
inputBonds {prefix}.bonds.txt
inputAngles {prefix}.angles.txt
inputDihedrals {prefix}.dihedrals.txt
inputExcludes {prefix}.excludes.txt
""".format( prefix=prefix ))
        
        with open("null.dx",'w') as fh:
                fh.write("""object 1 class gridpositions counts  2 2 2
origin -4000.00000 -4000.00000 -4000.00000
delta  8000.00000 0.000000 0.000000
delta  0.000000 8000.00000 0.000000
delta  0.000000 0.000000 8000.00000
object 2 class gridconnections counts  2 2 2
object 3 class array type float rank 0 items 8 data follows
0.0	0.0	0.0	
0.0	0.0	0.0	
0.0	0.0	
attribute "dep" string "positions"
object "density" class field 
component "positions" value 1
component "connections" value 2
component "data" value 3
""")

    def getParticleTypesAndCounts(self):
        return sorted( self.particleTypeCounts.items(), key=lambda x: x[0] )

    def _particleTypePairIter(self):
        typesAndCounts = self.getParticleTypesAndCounts()
        for i in range(len(typesAndCounts)):
            t1 = typesAndCounts[i][0]
            for j in range(i,len(typesAndCounts)):
                t2 = typesAndCounts[j][0]
                yield( (i,j,t1,t2) )
    
    def _writeArbdPotentialFiles(self, prefix, directory = "potentials"):
        ## TODO: remove reduncant directory calls
        try: 
            os.makedirs(directory)
        except OSError:
            if not os.path.isdir(directory):
                raise

        pathPrefix = "%s/%s-" % (directory,prefix)
        self._writeNonbondedParameterFiles( pathPrefix + "nb" )
        self._writeBondParameterFiles( pathPrefix )
        self._writeAngleParameterFiles( pathPrefix )
        self._writeDihedralParameterFiles( pathPrefix )
                
    def _writeHarmonicPotentialFile(self, filename, k, x0, resolution=0.1, xmin=0, xmax=35, maxForce=None, periodicity=None):
        x = np.arange( xmin, xmax+resolution*2, resolution )
        if periodicity is None:
            dx = x-x0
        else:
            dx = np.mod( x-x0 + 0.5*periodicity, periodicity) - 0.5*periodicity 
        u = 0.5*k*dx**2
        if maxForce is not None:
            assert(maxForce > 0)
            f = np.diff(u)/np.diff(x)
            f[f>maxForce] = maxForce
            f[f<-maxForce] = -maxForce
            u[0] = 0
            u[1:] = np.cumsum(f*np.diff(x))
        
        np.savetxt( filename, np.array([x, u]).T, fmt="%f" )

    def _writeNonbondedParameterFiles(self, prefix):
        x = np.arange(0, 50, 0.1)
        for i,j,t1,t2 in self._particleTypePairIter():
            f = "%s.%s-%s.dat" % (prefix, t1, t2)

            if t1 == "O" or t2 == "O":
                y = np.zeros(np.shape(x))
            else:
                bps1,bps2 = [float( t[1:] )/10 for t in (t1,t2)]
                y = nbPot.nbPot(x, bps1, bps2)

            np.savetxt( f, np.array([x, y]).T )
            self._nbParamFiles.append(f)

    def _writeBondParameterFiles(self, prefix):
        for pot in list(set([item[-1] for item in self.bonds])):
            pot.write_file(prefix)

    def _writeAngleParameterFiles(self, prefix):
        for pot in list(set([item[-1] for item in self.angles])):
            pot.write_file(prefix)

    def _writeDihedralParameterFiles(self, prefix):
        for pot in list(set([item[-1] for item in self.dihedrals])):
            pot.write_file(prefix)
            
    def addBond(self, *args):
        self.bonds.add(args)
    def addAngle(self, *args):
        self.angles.add(args)
    def addDihedral(self, *args):
        self.dihedrals.add(args)

    def _buildBonds(self, prefix, directory="potentials"):
        self.bonds = set()
        
        ## Get intrahelical bonds
        for nodes,seps in self._getIntrahelicalBonds():
            n1,n2 = nodes
            sep,  = seps
            if n1.type[0] == "d" and n2.type[0] == "d": 
                k = 10.0/sqrt(sep) # TODO: determine from simulations
                d = 3.4*sep
            else:
                ## TODO: get correct numbers from ssDNA model
                k = 1.0/sqrt(sep)
                d = 5*sep   
            self.addBond(n1, n2, Bond(k, d))

        ## Get crossover bonds
        for nodes,fwds in self._getCrossoverBonds():
            n1,n2 = nodes
            self.addBond(n1, n2, Bond(4, 18.5))

        ## Get crossover bonds
        for nodes,fwds in self._getSsCrossoverBonds():
            n1,n2 = nodes
            self.addBond(n1, n2, Bond(1, 5))

        ## Get crossover bonds
        for nodes,seps in self._getOrientationBonds():
            n1,n2 = nodes
            self.addBond(n1, n2, Bond(30, 1))  # TODO: improve params

    def _buildAngles(self, prefix, directory="potentials"):
        kT = 0.58622522         # kcal/mol
        for nodes,seps in self._getIntrahelicalAngles():
            n1,n2,n3 = nodes
            sep1,sep2  = seps
            sep = sep1+sep2
            if n1.type[0] == "d" and n2.type[0] == "d" and n3.type[0] == "d":
                ## <cos(q)> = exp(-s/Lp) = integrate( x^4 exp(-A x^2) / 2, {x, 0, pi} ) / integrate( x^2 exp(-A x^2), {x, 0, pi} )
                ## <cos(q)> ~ 1 - 3/4A                                                                                            
                ## where A = k_spring / (2 kT)                                                                                    
                k = 1.5 * kT * (1.0 / (1-exp(-float(sep)/147))) * 0.00030461742; # kcal_mol/degree^2
                # k *= 5
            else:
                ## TODO: get correct number from ssDNA model
                k = 1.5 * kT * (1.0 / (1-exp(-float(sep)/3))) * 0.00030461742; # kcal_mol/degree^2                

            ## Intrahelical 180 deg orientation angles
            if None not in [n.orientationNode for n in nodes]:
                k *= 0.5    # halve spring constant because using 2 springs
                args = [n.orientationNode for n in nodes]
                args.append( Angle(k,180) )
                self.addAngle( *args )

            self.addAngle( n1,n2,n3,Angle(k,180) )
            
        a,d = self._getCrossoverAnglesAndDihedrals()
        for nodes,sep in a:
            n1,n2,n3 = nodes
            k = (1.0/2) * 1.5 * kT * (1.0 / (1-exp(-float(sep)/147))) * 0.00030461742; # kcal_mol/degree^2                
            self.addAngle( n1,n2,n3,Angle(k,90) )

        ## Intrahelical 90 deg orientation angles
        for nodes,seps in self._getOrientationAngles():
            n1,n2,n3 = nodes
            sep = np.sum(seps)
            k = (1.0/2) * 1.5 * kT * (1.0 / (1-exp(-float(sep)/147))) * 0.00030461742; # kcal_mol/degree^2                
            self.addAngle( n1,n2,n3,Angle(k,90) )

        ## Crossover orientation angles
        for nodes,fwds in self._getCrossoverBonds():
            n1,n2 = nodes
            f1,f2 = fwds
            o1,o2 = [n.orientationNode for n in nodes]

            k = (1.0/2) * 1.5 * kT * (1.0 / (1-exp(-float(1)/147))) * 0.00030461742; # kcal_mol/degree^2                

            if o1 is not None:
                t0 = 90 + 60
                if f1: t0 -= 120
                self.addAngle( o1,n1,n2,Angle(k,t0) )

            if o2 is not None:
                t0 = 90 + 60
                if f2: t0 -= 120
                self.addAngle( n1,n2,o2,Angle(k,t0) )

    def _buildDihedrals(self, prefix, directory="potentials"):
        kT = 0.58622522         # kcal/mol
        a,d = self._getCrossoverAnglesAndDihedrals()
        for nodes,sep,isFwd1,isFwd2 in d:
            n1,n2,n3,n4 = nodes
            ## <cos(q)> = exp(-s/Lp) = integrate( cos[x] exp(-A x^2), {x, 0, pi} ) / integrate( exp(-A x^2), {x, 0, pi} )
            ##   Assume A is small
            ## int[B_] :=  Normal[Integrate[ Series[Cos[x] Exp[-B x^2], {B, 0, 1}], {x, 0, \[Pi]}]/
            ##             Integrate[Series[Exp[-B x^2], {B, 0, 1}], {x, 0, \[Pi]}]]
            ## Actually, without assumptions I get fitFun below

            ## From http://www.annualreviews.org/doi/pdf/10.1146/annurev.bb.17.060188.001405
            ##   units "3e-19 erg cm/ 295 k K" "nm" =~ 73
            Lp = self.twistPersistenceLength/0.34  # set semi-arbitrarily as there is a large spread in literature

            fitFun = lambda x: np.real(erf( (4*np.pi*x + 1j)/(2*np.sqrt(x)) )) * np.exp(-1/(4*x)) / erf(2*np.sqrt(x)*np.pi) - exp(-sep/Lp)
            k = opt.leastsq( fitFun, x0=exp(-sep/Lp) )
            k = k[0][0] * 2*kT*0.00030461742

            # intrinsicDegrees=30
            # fitFun = lambda x: (1.0/(2*x) - 2*np.sqrt(np.pi)*np.exp(-4*np.pi**2*x) / (np.sqrt(x)*erf(2*np.pi*np.sqrt(x))) ) - \
            #          ( (intrinsicDegrees*np.pi/180)**2 + 2*(1-exp(-sep/Lp)) )
            # k = opt.leastsq( fitFun, x0=1/(1-exp(-sep/Lp)) )
            # k = k[0][0] * 2*kT*0.00030461742

            t0 = sep*(360.0/10.5)

            # pdb.set_trace()
            if isFwd1[0]: t0 -= 120
            if isFwd2[0]: t0 += 120

            t0 = t0 % 360

            # if n2.idx == 0:
            #     print( n1.idx,n2.idx,n3.idx,n4.idx,k,t0,sep )
            self.addDihedral( n1,n2,n3,n4,Dihedral(k,t0) )

        for nodes,seps in self._getOrientationDihedrals():
            n1,n2,n3,n4 = nodes
            sep = seps[1]
            t0 = sep*(360.0/10.5)

            Lp = self.twistPersistenceLength/0.34  # set semi-arbitrarily as there is a large spread in literature

            fitFun = lambda x: np.real(erf( (4*np.pi*x + 1j)/(2*np.sqrt(x)) )) * np.exp(-1/(4*x)) / erf(2*np.sqrt(x)*np.pi) - exp(-sep/Lp)
            k = opt.leastsq( fitFun, x0=exp(-sep/Lp) )
            k = k[0][0] * 2*kT*0.00030461742
            # k *= 0.1
            # k *= 0

            self.addDihedral( n1,n2,n3,n4,Dihedral(k,t0,maxPotential=1) )

        ## Crossover dihedral angles
        for nodes,fwds in self._getCrossoverBonds():
            n1,n2 = nodes
            f1,f2 = fwds
            o1,o2 = [n.orientationNode for n in nodes]
            a1,a2 = [n.nodeAbove for n in nodes]
            b1,b2 = [n.nodeBelow for n in nodes]

            k = (1.0/2) * 1.5 * kT * (1.0 / (1-exp(-float(1)/147))) * 0.00030461742; # kcal_mol/degree^2                
            if o1 is not None:
                t0 = 90
                # if f1: t0 = -90
                if a2 is not None:
                    self.addDihedral( o1,n1,n2,a2,Dihedral(k,t0) )

            if o2 is not None:
                t0 = 90
                # if f2: t0 = -90
                if a1 is not None:
                    self.addDihedral( o2,n2,n1,a1,Dihedral(k,t0) )

            if o1 is not None and o2 is not None:
                if a1 is not None and a2 is not None:
                    t0 = 0
                    self.addDihedral( a1,n1,n2,a2,Dihedral(k,t0) )
                elif b1 is not None and b2 is not None:
                    t0 = 0
                    self.addDihedral( b1,n1,n2,b2,Dihedral(k,t0) )
                elif b1 is not None and a2 is not None:
                    t0 = 180
                    self.addDihedral( b1,n1,n2,a2,Dihedral(k,t0) )
                elif a1 is not None and b2 is not None:
                    t0 = 180
                    self.addDihedral( a1,n1,n2,b2,Dihedral(k,t0) )


    def _writeArbdBondFile(self, prefix, directory="potentials"):
        filename = prefix + ".bonds.txt"
        prefix = "%s/%s-" % (directory,prefix)
        with open(filename,'w') as fh:        
            for n1,n2,pot in self.bonds:
                fh.write("BOND ADD %d %d %s\n" % (n1.idx, n2.idx, pot.filename(prefix)))

    def _writeArbdAngleFile(self, prefix, directory="potentials"):
        filename = prefix + ".angles.txt"
        prefix = "%s/%s-" % (directory,prefix)
        with open(filename,'w') as fh:        
            for n1,n2,n3,pot in self.angles:
                fh.write("ANGLE %d %d %d %s\n" % (n1.idx, n2.idx, n3.idx, pot.filename(prefix)))

    def _writeArbdDihedralFile(self, prefix, directory="potentials"):
        filename = prefix + ".dihedrals.txt"
        prefix = "%s/%s-" % (directory,prefix)
        with open(filename,'w') as fh:        
            for n1,n2,n3,n4,pot in self.dihedrals:
                fh.write("DIHEDRAL %d %d %d %d %s\n" % (n1.idx, n2.idx, n3.idx, n4.idx, pot.filename(prefix)))

    def _writeArbdExclFile(self, filename):
        ## Exclude all 1-4 intrahelical nodes
        # e = 4
        e = 8
        exclusions = { (nodes[i],nodes[j]) 
                       for nodes,seps in self._getIntrahelicalNodeSeries(e)
                       for i in range(e-1) 
                       for j in range(i,e) }
        
        ## TODO, make exclusions depend on distance
        
        ## Exclude ssDNA contacts
        for nodes,seps in self._getSsCrossoverBonds():
            n1,n2 = nodes       # recall that nodes is sorted by .idx
            exclusions.add( nodes )
            exclusions.update( [(n1,n) for n in (n2.nodeBelow,n2.nodeAbove) if n is not None] )
            exclusions.update( [(n,n2) for n in (n1.nodeBelow,n1.nodeAbove) if n is not None] )

        ## Exclude crossovers and nearby
        for nodes,fwds in self._getCrossoverBonds():
            n1,n2 = nodes       # recall that nodes is sorted by .idx
            exclusions.add( nodes )
            exclusions.update( [(n1,n) for n in (n2.nodeBelow,n2.nodeAbove) if n is not None] )
            exclusions.update( [(n,n2) for n in (n1.nodeBelow,n1.nodeAbove) if n is not None] )

        ## Write exclusions
        with open(filename,'w') as fh:        
            for n1,n2 in exclusions:
                fh.write( "EXCLUDE %d %d\n" % (n1.idx,n2.idx) )

    def _getNonbondedPotential(self,x,a,b):
        return a*(np.exp(-x/b))