fix: providing self.box from sim.py when .pdb does not have it

calvados/build.py

@@ -294,7 +294,7 @@ def build_xyzgrid(N,box):
     return xyz
 
 # FOLDED
-def geometry_from_pdb(pdb,use_com=False):
+def geometry_from_pdb(pdb, use_com=False, auto_box_margin=100., verbose=1):   # [A]
     """ positions in nm"""
     with catch_warnings():
         simplefilter("ignore")
@@ -310,7 +310,19 @@ def geometry_from_pdb(pdb,use_com=False):
     else:
         cas = u.select_atoms('name CA')
         pos = cas.positions / 10.
-    box = np.append(u.dimensions[:3]/10.,u.dimensions[3:])
+
+    if(u.dimensions is None):   # the case when .pdb does not provide CRYST cell
+        if(auto_box_margin >= 0):
+            if(verbose > 0):
+                print(f'WARNING: did not find dimensions in PDB: {pdb}, automatically setting the box to the X-range + margin={auto_box_margin}[A]')
+            pos_size = np.amax(ag.positions, axis=0) - np.amin(ag.positions, axis=0)
+            u.dimensions = [pos_size[0] + auto_box_margin, pos_size[1] + auto_box_margin, pos_size[2] + auto_box_margin, 90, 90, 90]
+
+    if(u.dimensions is None):
+        box = None
+    else:
+        box = np.append(u.dimensions[:3]/10., u.dimensions[3:])
+
     return pos, box
 
 def geometry_from_pdb_rna(pdb,use_com=False):

calvados/components.py

@@ -112,11 +112,17 @@ class Protein(Component):
     def __init__(self, name: str, comp_dict: dict, defaults: dict):
         super().__init__(name, comp_dict, defaults)
 
-    def calc_x_from_pdb(self):
+    def calc_x_from_pdb(self, verbose: bool = False):
         """ Calculate protein positions from pdb. """
 
         input_pdb = f'{self.pdb_folder}/{self.name}.pdb'
-        self.xinit, self.dimensions = build.geometry_from_pdb(input_pdb,use_com=self.use_com) # read from pdb
+        if(hasattr(self, 'dimensions')):
+            self.xinit, _ = build.geometry_from_pdb(input_pdb, use_com=self.use_com, auto_box_margin=-1, verbose=verbose) # read from pdb
+            if(self.to_check_box > 0):
+                backbone_Xrange = np.amax(self.xinit, axis=0) - np.amin(self.xinit, axis=0)
+                assert(np.all(backbone_Xrange < np.array(self.dimensions[:3]))), f'ERROR: X-range from the PDB "{input_pdb}" is ({backbone_Xrange[0], backbone_Xrange[1], backbone_Xrange[2]})[A], which does not fit into the provided box ({box[0], box[1], box[2]})[A]'
+        else:
+            self.xinit, self.dimensions = build.geometry_from_pdb(input_pdb, use_com=self.use_com, verbose=verbose)
 
     def calc_ssdomains(self):
         """ Get bounds for restraints (harmonic). """
@@ -167,7 +173,6 @@ def calc_go_scale(self, bscale_shift = 0.1, bscale_width = 80):#,
 
     def calc_properties(self, pH: float = 7.0, verbose: bool = False, comp_setup: str = 'spiral'):
         """ Protein properties. """
-
         super().calc_properties(pH=pH, verbose=verbose)
 
         # fix charge and mw of termini
@@ -178,7 +183,7 @@ def calc_properties(self, pH: float = 7.0, verbose: bool = False, comp_setup: st
 
         if self.restraint:
             # self.init_restraint_force() # Done via sim.py
-            self.calc_x_from_pdb()
+            self.calc_x_from_pdb(verbose=verbose)
             self.calc_dmap()
             if self.restraint_type == 'harmonic':
                 self.calc_ssdomains()

calvados/sim.py

@@ -35,6 +35,9 @@ def __init__(self,path,config,components):
         self.comp_defaults = components['defaults']
 
         self.box = np.array(self.box)
+        if('dimensions' not in self.comp_defaults):   # Protein has attr 'dimensions', not 'box'
+            self.comp_defaults['dimensions'] = list(self.box) + [90.] * 3
+            self.comp_defaults['to_check_box'] = 1
         self.eps_lj *= 4.184 # kcal to kJ/mol
 
         if self.restart == 'checkpoint' and os.path.isfile(f'{self.path}/{self.frestart}'):
@@ -141,7 +144,6 @@ def build_system(self):
         a, b, c = build.build_box(self.box[0],self.box[1],self.box[2])
         self.system.setDefaultPeriodicBoxVectors(a, b, c)
 
-
         # init interaction parameters (required before make components)
         self.eps_yu, self.k_yu = interactions.genParamsDH(self.temp,self.ionic)