polyconf.PDB has missing unit cell dimensions in some cases, possibly related to polymer.copy() #55
Description
I can't figure out what's going on here, but sometimes PDB.cleanup() is not correctly assigning to box size to the universe
Example code uses the input files in polyconf_examples
Here's a process that works
polymer1=Polymer(Monomer('PEI_start.pdb'))
adds=126
for i in range (0,adds):
polymer1.extend( # extend with one monomer, aligned along this step's linearization vector
Monomer('PEI_monomer.pdb'), # extend with this monomer
n=polymer1.maxresid(), # we will allways add onto the existing monomer with the highest resid
nn=polymer1.newresid(), # the incoming monomer needs a new resid
names=dict(P1='CX',P2='C1',Q1='N1',Q2='NX'), # C1_i+1 fit to CX_i, then rotate so NX_i+1 fit to N1_i
joins=[('N1','C1')],# new connection between N1_i and C1_i+1
)
polymer1.extend(
Monomer('PEI_end.pdb'),
n=polymer1.maxresid(),
nn=polymer1.newresid(),
names=dict(P1='CX',P2='C1',Q1='N1',Q2='NX'),
joins=[('N1','C1')],
)
Saver = PDB(polymer1)
Saver.cleanup()
Saver.save(dummyAtoms='CX NX',fname='polymer_01_vanilla-extend')
This saves normally
But this doesn't.
polymer4=polymer1.copy()
CN_dihedrals=polymer4.gen_pairlist(a1='C2',a2='N1',first_resid=1,last_resid=128,mult=3)
cutoff = 0.7
polymer4.dihedral_solver(CN_dihedrals,dummy='CX,NX',cutoff=cutoff)
Saver = PDB(polymer4)
# for debug:
print( (Saver.atoms.positions).max(axis=0), (Saver.atoms.positions).min(axis=0))
print( (Saver.atoms.positions).max(axis=0)- (Saver.atoms.positions).min(axis=0))
Saver.cleanup()
Saver.save(dummyAtoms='CX NX',fname='polymer_04_solved') # save, excluding dummy atoms
This returns an error
UserWarning: Unit cell dimensions not found. CRYST1 record set to unitary values.
And it sets the unit cell to 1,1,1,90,90,90
The start of the output file is
0.00 90.00 P 1 1
REMARK 285 UNITARY VALUES FOR THE UNIT CELL AUTOMATICALLY SET
REMARK 285 BY MDANALYSIS PDBWRITER BECAUSE UNIT CELL INFORMATION
REMARK 285 WAS MISSING.
REMARK 285 PROTEIN DATA BANK CONVENTIONS REQUIRE THAT
REMARK 285 CRYST1 RECORD IS INCLUDED, BUT THE VALUES ON
REMARK 285 THIS RECORD ARE MEANINGLESS.
The output files look correct, apart from the box dimensions.
I think what is happening is that the way we handle polymers means we don't unit cell define dimensions at all, unless they happen to be in the monomer input file. polymer.extend() assigns unit cell dimensions based on the size of the polymer after extension
from Polymer.py lines 252 -253 (inside polymer.extend() )
new.dimensions = list(new.atoms.positions.max(axis=0) + [0.5,0.5,0.5]) + [90]*3
self.polymer = new.copy()
but cleanup doesn't seem to do that correctly, even though there's a step in there to do it
from PDB.py lines 22 to 30
def cleanup(self):
"""
Adjust the box size and center the polymer in the box
"""
box = (self.atoms.positions).max(axis=0)- (self.atoms.positions).min(axis=0) + [10,10,10]
self.dimensions = list(box) + [90]*3 <--- this line should be resetting the dimensions
cog = self.atoms.center_of_geometry(wrap=False)
box_center = box / 2
self.atoms.translate(box_center - cog)
Luna, are you able to have a look at this and see if you can use your object oriented magic to figure out what the fuck is going wrong? I am guessing either polymer.copy() or PDB.cleanup() is not inheriting something from MDAnalysis, but I can't really tell.