Added routines for applying grid potentials

mrdna/model/arbdmodel.py

@@ -843,12 +843,14 @@ diffusion {diffusivity}
 """.format(**particleParams))
                 if 'grid' in particleParams:
                     if not isinstance(pt.grid, list): pt.grid = [pt.grid]
-                    for g in pt.grid:
+                    for g,s in pt.grid:
+                        ## TODO, use Path.relative_to?
+                        try:
+                            fh.write("gridFile {}\n".format(g.relative_to(os.getcwd())))
+                        except:
                             fh.write("gridFile {}\n".format(g))
-                    ## TODO: make this prettier? multiple scaling factors?
-                    gridFileScale = 1.0
-                    if 'gridFileScale' in pt.__dict__: gridFileScale = pt.gridFileScale
-                    fh.write("gridFileScale {}\n".format(gridFileScale))
+
+                        fh.write("gridFileScale {}\n".format(s))
 
                 else:
                     fh.write("gridFile {}/null.dx\n".format(self.potential_directory))

mrdna/segmentmodel.py

 import pdb
+from pathlib import Path
 import numpy as np
 import random
 from .model.arbdmodel import PointParticle, ParticleType, Group, ArbdModel
@@ -1386,7 +1387,7 @@ class SegmentModel(ArbdModel):
 
         self._bonded_potential = dict() # cache for bonded potentials
         self._generate_strands()
-        self.potentials = []
+        self.grid_potentials = []
         self._generate_bead_model( max_basepairs_per_bead, max_nucleotides_per_bead, local_twist, escapable_twist)
 
         self.useNonbondedScheme( nbDnaScheme )
@@ -1908,6 +1909,7 @@ class SegmentModel(ArbdModel):
         for b,c in zip(beads,clusters):
             _assign_bead_type(b, cluster_size[c-1], decimals=order)
 
+        self._apply_grid_potentials_to_beads(beadtype_s)
 
         # for bead in [b for s in segments for b in s]:
         #     num_nt0 = bead.num_nt
@@ -2728,6 +2730,27 @@ proc calcforces {} {
             dx = dy = dz = max((dx,dy,dz))
         self.dimensions = [dx,dy,dz]
 
+    def add_grid_potential(self, grid_file, scale=1, per_nucleotide=True):
+        grid_file = Path(grid_file)
+        if not grid_file.is_file():
+            raise ValueError("Grid file {} does not exist".format(grid_file))
+        if not grid_file.is_absolute():
+            grid_file = Path.cwd() / grid_file
+        self.grid_potentials.append((grid_file,scale,per_nucleotide))
+
+    def _apply_grid_potentials_to_beads(self, bead_type_dict):
+        if len(self.grid_potentials) > 1:
+            raise NotImplementedError("Multiple grid potentials are not yet supported")
+
+        for grid_file, scale, per_nucleotide in self.grid_potentials:
+            for key,particle_type in bead_type_dict.items():
+                if particle_type.name[0] == "O": continue
+                s = scale*particle_type.nts if per_nucleotide else scale
+                try:
+                    particle_type.grid = particle_type.grid + (grid_file, s)
+                except:
+                    particle_type.grid = tuple((grid_file, s))
+
     def _generate_atomic_model(self, scale=1):
         ## TODO: deprecate
         self.generate_atomic_model(scale=scale)