Source code for arbdmodel.sali_polymer_model

# -*- coding: utf-8 -*-
## Test with `python -m arbdmodel.hps_polymer_model`

import numpy as np
import sys


## Local imports
from . import ParticleType, PointParticle, Group
from .logger import logger, get_resource_path    
from .model import ArbdModel
from .polymer import PolymerSection, PolymerGroup
from .interactions import AbstractPotential, HarmonicBond, HarmonicBondedPotential
from .coords import quaternion_to_matrix


"""Define particle types"""
type_ = ParticleType("AA",
                     # units "297.15 k K / (6 pi 0.92 mPa s 6 AA)" "AA**2/ns"
                     diffusivity = 39.429314
)

## Bonded potentials
[docs] class LinearBond(HarmonicBond): def __init__(self, k, r0, rRange=(0,50), resolution=0.1, maxForce=None, max_potential=None, prefix="potentials/"): HarmonicBondedPotential.__init__(self, k, r0, rRange, resolution, maxForce, max_potential, prefix) self.type_ = "linearbond" self.kscale_ = 1.0
[docs] def potential(self, dr): return self.k*np.abs(dr)
[docs] class SaliNonbonded(AbstractPotential): def __init__(self, resolution=0.1): AbstractPotential.__init__(self, resolution=resolution)
[docs] def potential(self, r): """ Constant force excluded volume """ force = 10 # kcal_mol/AA radius = 6 u = np.zeros(r.shape) s = r < 2*radius u[s] = (2*radius - r[s]) * force return u
[docs] class SaliBeadsFromPolymer(Group): # p = PointParticle(_P, (0,0,0), "P") # b = PointParticle(_B, (3,0,1), "B") # nt = Group( name = "nt", children = [p,b]) # nt.add_bond( i=p, j=b, bond = get_resource_path('two_bead_model/BPB.dat') ) def __init__(self, polymer, sequence=None, **kwargs): if sequence is None: raise NotImplementedError # ... set random sequence self.polymer = polymer self.sequence = sequence for prop in ('segname','chain'): if prop not in kwargs: # import pdb # pdb.set_trace() try: self.__dict__[prop] = polymer.__dict__[prop] except: pass if len(sequence) != polymer.num_monomers: raise ValueError("Length of sequence does not match length of polymer") Group.__init__(self, **kwargs) self.num_beads = int(np.ceil(polymer.num_monomers / 20)) def _clear_beads(self): ... def _generate_beads(self): # beads = self.children nb = self.num_beads for i in range(nb): # c = self.polymer.monomer_index_to_contour(i) c = float(i)/(nb-1) if nb > 1 else 0.5 r = self.polymer.contour_to_position(c) bead = PointParticle(type_, r, resid = i+1) self.add(bead) ## Two consecutive nts for i in range(len(self.children)-1): b1 = self.children[i] b2 = self.children[i+1] """ units "10 kJ/N_A" kcal_mol """ bond = LinearBond(k = 1, r0 = 18, rRange = (0,500), resolution = 0.01, maxForce = 30) self.add_bond( i=b1, j=b2, bond = bond, exclude=True )
[docs] class SaliModel(ArbdModel): def __init__(self, polymers, sequences = None, debye_length = 10, damping_coefficient = 10, DEBUG=False, **kwargs): """ [debye_length]: angstroms [damping_coefficient]: ns """ print("WARNING: diffusion coefficient arbitrarily set to 100 AA**2/ns in SaliModel") kwargs['timestep'] = 1047e-6 kwargs['cutoff'] = 18 if 'decompPeriod' not in kwargs: kwargs['decompPeriod'] = 1000 """ Assign sequences """ if sequences is None: raise NotImplementedError("HpsModel must be provided a sequences argument") self.polymer_group = PolymerGroup(polymers) self.sequences = sequences ArbdModel.__init__(self, [], **kwargs) # """ Update type diffusion coefficients """ # self.set_damping_coefficient( damping_coefficient ) """ Set up nonbonded interactions """ nonbonded = SaliNonbonded() self.useNonbondedScheme( nonbonded, typeA=type_, typeB=type_ ) """ Generate beads """ self.generate_beads()
[docs] def update_splines(self, coords): i = 0 for p,c in zip(self.polymer_group.polymers,self.children): n = len(c.children) p.set_splines(np.linspace(0,1,n), coords[i:i+n]) i += n self.clear_all() self.generate_beads()
## TODO Apply restraints, etc
[docs] def generate_beads(self): self.peptides = [SaliBeadsFromPolymer(p,s) for p,s in zip(self.polymer_group.polymers,self.sequences)] for s in self.peptides: self.add(s) s._generate_beads()
# def set_damping_coefficient(self, damping_coefficient): # for t in self.types: # t.damping_coefficient = damping_coefficient # # t.diffusivity = 831447.2 * temperature / (t.mass * damping_coefficient) if __name__ == "__main__": pass