arbdmodel.py 39.5 KB
Newer Older
 cmaffeo2 committed Feb 28, 2018 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 ``````# -*- coding: utf-8 -*- import abc import numpy as np from copy import copy, deepcopy from inspect import ismethod # from datetime import datetime # from math import pi,sqrt,exp,floor # from scipy.special import erf # import scipy.optimize as opt import os, sys, subprocess ## Abstract classes class Transformable(): def __init__(self, position, orientation=None): self.position = np.array(position) if orientation is not None: orientation = np.array(orientation) self.orientation = orientation def transform(self, R = ((1,0,0),(0,1,0),(0,0,1)), center = (0,0,0), offset = (0,0,0)): R,center,offset = [np.array(x) for x in (R,center,offset)] self.position = R.dot(self.position-center)+center+offset if self.orientation is not None: ## TODO: what if self.orientation is taken from parent?! self.orientation = self.orientation.dot(R) ... def collapsedPosition(self): # print("collapsedPosition called", type(self), self.name) if isinstance(self, Child): # print(self.parent, isinstance(self.parent,Transformable)) if isinstance(self.parent, Transformable): return self.applyOrientation(self.position) + self.parent.collapsedPosition() # if self.parent.orientation is not None: # return self.parent.collapsedOrientation().dot(self.position) + self.parent.collapsedPosition() return np.array(self.position) # return a copy def applyOrientation(self,obj): # print("applyOrientation called", self.name, obj) if isinstance(self, Child): # print( self.orientation, self.orientation is not None, None is not None ) # if self.orientation is not None: # # print("applyOrientation applying", self, self.name, self.orientation) # obj = self.orientation.dot(obj) if isinstance(self.parent, Transformable): if self.parent.orientation is not None: obj = self.parent.orientation.dot(obj) obj = self.parent.applyOrientation(obj) # print("applyOrientation returning", self.name, obj) return obj class Parent(): `````` cmaffeo2 committed Jul 03, 2018 60 `````` def __init__(self, children=None, remove_duplicate_bonded_terms=False): `````` cmaffeo2 committed Feb 28, 2018 61 `````` self.children = [] `````` cmaffeo2 committed Jul 03, 2018 62 63 64 `````` if children is not None: for x in children: self.add(x) `````` cmaffeo2 committed Jun 14, 2018 65 66 `````` self.remove_duplicate_bonded_terms = remove_duplicate_bonded_terms `````` cmaffeo2 committed Feb 28, 2018 67 68 69 `````` self.bonds = [] self.angles = [] self.dihedrals = [] `````` cmaffeo2 committed Mar 23, 2018 70 `````` self.impropers = [] `````` cmaffeo2 committed Feb 28, 2018 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 `````` self.exclusions = [] ## TODO: self.cacheInvalid = True # What will be in the cache? def add(self,x): ## TODO: check the parent-child tree to make sure there are no cycles if not isinstance(x,Child): raise Exception('Attempted to add an object to a group that does not inherit from the "Child" type') if x.parent is not None and x.parent is not self: raise Exception("Child {} already belongs to some group".format(x)) x.parent = self self.children.append(x) `````` cmaffeo2 committed Mar 20, 2018 86 87 88 89 90 91 92 93 `````` def clear_all(self, keep_children=False): if keep_children == False: for x in self.children: x.parent = None self.children = [] self.bonds = [] self.angles = [] self.dihedrals = [] `````` cmaffeo2 committed Mar 23, 2018 94 `````` self.impropers = [] `````` cmaffeo2 committed Mar 20, 2018 95 96 `````` self.exclusions = [] `````` cmaffeo2 committed Feb 28, 2018 97 98 99 100 101 102 103 `````` def remove(self,x): if x in self.children: self.children.remove(x) x.parent = None def add_bond(self, i,j, bond, exclude=False): ## TODO: how to handle duplicating and cloning bonds `````` cmaffeo2 committed Mar 30, 2018 104 105 `````` # beads = [b for b in self] # for b in (i,j): assert(b in beads) `````` cmaffeo2 committed Feb 28, 2018 106 107 108 `````` self.bonds.append( (i,j, bond, exclude) ) def add_angle(self, i,j,k, angle): `````` cmaffeo2 committed Mar 30, 2018 109 110 `````` # beads = [b for b in self] # for b in (i,j,k): assert(b in beads) `````` cmaffeo2 committed Feb 28, 2018 111 112 113 `````` self.angles.append( (i,j,k, angle) ) def add_dihedral(self, i,j,k,l, dihedral): `````` cmaffeo2 committed Mar 30, 2018 114 115 `````` # beads = [b for b in self] # for b in (i,j,k,l): assert(b in beads) `````` cmaffeo2 committed Feb 28, 2018 116 117 `````` self.dihedrals.append( (i,j,k,l, dihedral) ) `````` cmaffeo2 committed Mar 23, 2018 118 `````` def add_improper(self, i,j,k,l, dihedral): `````` cmaffeo2 committed Mar 30, 2018 119 120 `````` # beads = [b for b in self] # for b in (i,j,k,l): assert(b in beads) `````` cmaffeo2 committed Mar 23, 2018 121 122 `````` self.impropers.append( (i,j,k,l, dihedral) ) `````` cmaffeo2 committed Feb 28, 2018 123 124 `````` def add_exclusion(self, i,j): ## TODO: how to handle duplicating and cloning bonds `````` cmaffeo2 committed Mar 02, 2018 125 126 127 `````` ## TODO: perform following check elsewhere # beads = [b for b in self] # for b in (i,j): assert(b in beads) `````` cmaffeo2 committed Feb 28, 2018 128 129 130 131 132 133 `````` self.exclusions.append( (i,j) ) def get_bonds(self): ret = self.bonds for c in self.children: if isinstance(c,Parent): ret.extend( c.get_bonds() ) `````` cmaffeo2 committed Jun 14, 2018 134 135 136 137 138 `````` if self.remove_duplicate_bonded_terms: return list(set(ret)) else: return ret `````` cmaffeo2 committed Feb 28, 2018 139 140 141 142 143 `````` def get_angles(self): ret = self.angles for c in self.children: if isinstance(c,Parent): ret.extend( c.get_angles() ) `````` cmaffeo2 committed Jun 14, 2018 144 145 146 147 `````` if self.remove_duplicate_bonded_terms: return list(set(ret)) else: return ret `````` cmaffeo2 committed Feb 28, 2018 148 149 150 151 152 `````` def get_dihedrals(self): ret = self.dihedrals for c in self.children: if isinstance(c,Parent): ret.extend( c.get_dihedrals() ) `````` cmaffeo2 committed Jun 14, 2018 153 154 155 156 `````` if self.remove_duplicate_bonded_terms: return list(set(ret)) else: return ret `````` cmaffeo2 committed Feb 28, 2018 157 `````` `````` cmaffeo2 committed Mar 23, 2018 158 159 160 161 `````` def get_impropers(self): ret = self.impropers for c in self.children: if isinstance(c,Parent): ret.extend( c.get_impropers() ) `````` cmaffeo2 committed Jun 14, 2018 162 163 164 165 `````` if self.remove_duplicate_bonded_terms: return list(set(ret)) else: return ret `````` cmaffeo2 committed Mar 23, 2018 166 `````` `````` cmaffeo2 committed Feb 28, 2018 167 168 169 170 `````` def get_exclusions(self): ret = self.exclusions for c in self.children: if isinstance(c,Parent): ret.extend( c.get_exclusions() ) `````` cmaffeo2 committed Jun 14, 2018 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 `````` if self.remove_duplicate_bonded_terms: return list(set(ret)) else: return ret ## Removed because prohibitively slow # def remove_duplicate_terms(self): # for key in "bonds angles dihedrals impropers exclusions".split(): # self.remove_duplicate_item(key) # def remove_duplicate_item(self, dict_key, existing=None): # if existing is None: existing = [] # ret = [i for i in list(set(self.__dict__[dict_key])) if i not in existing] # self.__dict__[dict_key] = ret # existing.extend(ret) # for c in self.children: # if isinstance(c,Parent): # ret = ret + c.remove_duplicate_item(dict_key, existing) # return ret `````` cmaffeo2 committed Feb 28, 2018 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 `````` def __iter__(self): ## TODO: decide if this is the nicest way to do it! """Depth-first iteration through tree""" for x in self.children: if isinstance(x,Parent): if isinstance(x,Clone) and not isinstance(x.get_original_recursively(),Parent): yield x else: for y in x: yield y else: yield x def __len__(self): l = 0 for x in self.children: if isinstance(x,Parent): l += len(x) else: l += 1 return l def __getitem__(self, i): return self.children[i] def __setitem__(self, i, val): x = self.children[i] x.parent = None val.parent = self self.children[i] = val class Child(): def __init__(self, parent=None): self.parent = parent `````` cmaffeo2 committed Mar 22, 2018 226 227 228 `````` if parent is not None: assert( isinstance(parent, Parent) ) parent.children.append(self) `````` cmaffeo2 committed Feb 28, 2018 229 230 231 232 233 234 `````` def __getattr__(self, name): """ Try to get attribute from the parent """ # if self.parent is not None: `````` cmaffeo2 committed Mar 22, 2018 235 `````` if "parent" not in self.__dict__ or self.__dict__["parent"] is None or name is "children": `````` cmaffeo2 committed Feb 28, 2018 236 `````` raise AttributeError("'{}' object has no attribute '{}'".format(type(self).__name__, name)) `````` cmaffeo2 committed Mar 23, 2018 237 238 239 240 241 242 243 244 245 246 `````` ## TODO: determine if there is a way to avoid __getattr__ if a method is being looked up try: ret = getattr(self.parent,name) except: raise AttributeError("'{}' object has no attribute '{}'".format(type(self).__name__, name)) if ismethod(ret): raise AttributeError("'{}' object has no method '{}'".format(type(self).__name__, name)) return ret `````` cmaffeo2 committed Feb 28, 2018 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 `````` # def __getstate__(self): # print("Child getstate called", self) # print(self.__dict__) # return (self.__dict__,) # def __setstate__(self, state): # self.__dict__, = state class Clone(Transformable, Parent, Child): def __init__(self, original, parent=None, position = None, orientation = None): if position is None and original.position is not None: position = np.array( original.position ) if orientation is None and original.orientation is not None: orientation = np.array( original.orientation ) if parent is None: parent = original.parent self.original = original Child.__init__(self, parent) Transformable.__init__(self, position, orientation) ## TODO: keep own bond_list, etc, update when needed original changes if "children" in original.__dict__ and len(original.children) > 0: self.children = [Clone(c, parent = self) for c in original.children] else: self.children = [] def get_original_recursively(self): if isinstance(self.original, Clone): return self.original.get_original_recursively() else: return self.original def __getattr__(self, name): """ Try to get attribute from the original without descending the tree heirarchy, then look up parent TODO: handle PointParticle lookups into ParticleType """ # print("Clone getattr",name) if name in self.original.__dict__: return self.original.__dict__[name] else: if "parent" not in self.__dict__ or self.__dict__["parent"] is None: raise AttributeError("'{}' object has no attribute '{}'".format(type(self).__name__, name)) return getattr(self.parent, name) ## Particle classes class ParticleType(): """Class that hold common attributes that particles can point to""" excludedAttributes = ("idx","type_", "position", `````` cmaffeo2 committed Mar 22, 2018 304 `````` "children", `````` cmaffeo2 committed Feb 28, 2018 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 `````` "parent", "excludedAttributes", ) def __init__(self, name, charge=0, **kargs): # parent = None # if "parent" in kargs: parent = kargs["parent"] # Child.__init__(self, parent) self.name = name self.charge = charge for key in ParticleType.excludedAttributes: assert( key not in kargs ) for key,val in kargs.items(): self.__dict__[key] = val def __hash_key(self): l = [self.name,self.charge] for keyval in sorted(self.__dict__.items()): l.extend(keyval) return tuple(l) def __hash__(self): return hash(self.__hash_key()) def _equal_check(a,b): if a.name == b.name: if a.__hash_key() != b.__hash_key(): `````` cmaffeo2 committed Mar 02, 2018 333 `````` raise Exception("Two different ParticleTypes have same 'name' attribute") `````` cmaffeo2 committed Feb 28, 2018 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 `````` def __eq__(a,b): a._equal_check(b) return a.name == b.name def __lt__(a,b): a._equal_check(b) return a.name < b.name def __le__(a,b): a._equal_check(b) return a.name <= b.name def __gt__(a,b): a._equal_check(b) return a.name > b.name def __ge__(a,b): a._equal_check(b) return a.name >= b.name class PointParticle(Transformable, Child): `````` cmaffeo2 committed Mar 27, 2018 353 `````` def __init__(self, type_, position, name="A", **kwargs): `````` cmaffeo2 committed Mar 22, 2018 354 355 356 357 `````` parent = None if 'parent' in kwargs: parent = kwargs['parent'] Child.__init__(self, parent=parent) `````` cmaffeo2 committed Feb 28, 2018 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 `````` Transformable.__init__(self,position) self.type_ = type_ self.idx = None self.name = name self.counter = 0 for key,val in kwargs.items(): self.__dict__[key] = val def __getattr__(self, name): """ First try to get attribute from the parent, then type_ Note that this data structure seems to be fragile, can result in stack overflow """ # return Child.__getattr__(self,name) try: return Child.__getattr__(self,name) except Exception as e: if 'type_' in self.__dict__: return getattr(self.type_, name) else: raise AttributeError("'{}' object has no attribute '{}'".format(type(self).__name__, name)) `````` cmaffeo2 committed Sep 15, 2018 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 `````` def _get_psfpdb_dictionary(self): p = self try: segname = p.segname except: segname = "A" try: resname = p.resname except: resname = p.name[:3] try: resid = p.resid except: resid = p.idx+1 try: occ = p.occupancy except: occ = 0 try: beta = p.beta except: beta = 0 data = dict(segname = segname, resname = resname, name = str(p.name)[:4], chain = "A", resid = int(resid), idx = p.idx+1, type = p.type_.name[:4], charge = p.charge, mass = p.mass, occupancy = occ, beta = beta ) return data `````` cmaffeo2 committed Feb 28, 2018 423 424 `````` class Group(Transformable, Parent, Child): `````` cmaffeo2 committed Sep 15, 2018 425 `````` `````` cmaffeo2 committed Jul 03, 2018 426 `````` def __init__(self, name=None, children = None, parent=None, `````` cmaffeo2 committed Feb 28, 2018 427 428 `````` position = np.array((0,0,0)), orientation = np.array(((1,0,0),(0,1,0),(0,0,1))), `````` cmaffeo2 committed Jun 14, 2018 429 `````` remove_duplicate_bonded_terms = False, `````` cmaffeo2 committed Sep 15, 2018 430 431 `````` **kwargs): `````` cmaffeo2 committed Feb 28, 2018 432 433 `````` Transformable.__init__(self, position, orientation) Child.__init__(self, parent) # Initialize Child first `````` cmaffeo2 committed Jun 14, 2018 434 `````` Parent.__init__(self, children, remove_duplicate_bonded_terms) `````` cmaffeo2 committed Feb 28, 2018 435 436 437 `````` self.name = name self.isClone = False `````` cmaffeo2 committed Sep 15, 2018 438 439 440 441 `````` for key,val in kwargs.items(): self.__dict__[key] = val `````` cmaffeo2 committed Feb 28, 2018 442 443 444 445 446 447 448 449 450 451 452 453 454 `````` def clone(self): return Clone(self) g = copy(self) g.isClone = True # TODO: use? g.children = [copy(c) for c in g.children] for c in g.children: c.parent = g return g g = Group(position = self.position, orientation = self.orientation) g.children = self.children # lists point to the same object def duplicate(self): `````` cmaffeo2 committed Mar 23, 2018 455 456 457 458 `````` new = deepcopy(self) for c in new.children: c.parent = new return new `````` cmaffeo2 committed Feb 28, 2018 459 460 461 462 `````` # Group(position = self.position, # orientation = self.orientation) # g.children = deepcopy self.children.deepcopy() # lists are the same object `````` cmaffeo2 committed Apr 04, 2018 463 `````` ## TODO override deepcopy so parent can be excluded from copying? `````` cmaffeo2 committed Feb 28, 2018 464 465 466 467 468 469 470 471 472 473 `````` # def __getstate__(self): # return (self.children, self.parent, self.position, self.orientation) # def __setstate__(self, state): # self.children, self.parent, self.position, self.orientation = state class PdbModel(Transformable, Parent): `````` cmaffeo2 committed Jul 03, 2018 474 `````` def __init__(self, children=None, dimensions=None, remove_duplicate_bonded_terms=False): `````` cmaffeo2 committed Feb 28, 2018 475 `````` Transformable.__init__(self,(0,0,0)) `````` cmaffeo2 committed Jun 14, 2018 476 `````` Parent.__init__(self, children, remove_duplicate_bonded_terms) `````` cmaffeo2 committed Feb 28, 2018 477 478 479 480 481 482 483 484 485 486 487 488 `````` self.dimensions = dimensions self.particles = [p for p in self] self.cacheInvalid = True def _updateParticleOrder(self): pass def writePdb(self, filename): if self.cacheInvalid: self._updateParticleOrder() with open(filename,'w') as fh: ## Write header `````` 489 `````` fh.write("CRYST1{:>9.3f}{:>9.3f}{:>9.3f} 90.00 90.00 90.00 P 1 1\n".format( *self.dimensions )) `````` cmaffeo2 committed Feb 28, 2018 490 491 `````` ## Write coordinates `````` cmaffeo2 committed Sep 15, 2018 492 `````` formatString = "ATOM {idx:>6.6s} {name:^4.4s} {resname:3.3s} {chain:1.1s}{resid:>5.5s} {x:8.8s}{y:8.8s}{z:8.8s}{occupancy:6.2f}{beta:6.2f} {charge:2d}{segname:>6s}\n" `````` cmaffeo2 committed Feb 28, 2018 493 494 `````` for p in self.particles: ## http://www.wwpdb.org/documentation/file-format-content/format33/sect9.html#ATOM `````` cmaffeo2 committed Sep 15, 2018 495 496 497 `````` data = p._get_psfpdb_dictionary() idx = data['idx'] `````` cmaffeo2 committed Apr 04, 2018 498 `````` if np.log10(idx) >= 5: `````` cmaffeo2 committed Aug 24, 2018 499 `````` idx = " *****" `````` cmaffeo2 committed Apr 04, 2018 500 501 `````` else: idx = "{:>6d}".format(idx) `````` cmaffeo2 committed Sep 15, 2018 502 `````` data['idx'] = idx `````` cmaffeo2 committed Apr 04, 2018 503 504 `````` pos = p.collapsedPosition() `````` New Tbgl User committed Apr 05, 2018 505 `````` dig = [max(int(np.log10(np.abs(x)+1e-6)//1),0)+1 for x in pos] `````` cmaffeo2 committed Apr 04, 2018 506 507 508 509 `````` for d in dig: assert( d <= 7 ) # assert( np.all(dig <= 7) ) fs = ["{: %d.%df}" % (8,7-d) for d in dig] x,y,z = [f.format(x) for f,x in zip(fs,pos)] `````` cmaffeo2 committed Sep 15, 2018 510 511 512 513 514 515 516 `````` data['x'] = x data['y'] = y data['z'] = z assert(data['resid'] < 1e5) data['charge'] = int(data['charge']) data['resid'] = "{:<4d}".format(data['resid']) fh.write( formatString.format(**data) ) `````` cmaffeo2 committed Apr 04, 2018 517 `````` `````` cmaffeo2 committed Feb 28, 2018 518 519 520 521 522 523 524 525 `````` return def writePsf(self, filename): if self.cacheUpToDate == False: self._updateParticleOrder() with open(filename,'w') as fh: ## Write header fh.write("PSF NAMD\n\n") # create NAMD formatted psf `````` cmaffeo2 committed Mar 23, 2018 526 `````` fh.write("{:>8d} !NTITLE\n\n".format(0)) `````` cmaffeo2 committed Feb 28, 2018 527 528 529 530 531 532 `````` ## ATOMS section fh.write("{:>8d} !NATOM\n".format(len(self.particles))) ## From vmd/plugins/molfile_plugin/src/psfplugin.c ## "%d %7s %10s %7s %7s %7s %f %f" `````` cmaffeo2 committed Apr 01, 2018 533 534 `````` formatString = "{idx:>8d} {segname:7.7s} {resid:<10.10s} {resname:7.7s}" + \ " {name:7.7s} {type:7.7s} {charge:f} {mass:f}\n" `````` cmaffeo2 committed Feb 28, 2018 535 `````` for p in self.particles: `````` cmaffeo2 committed Sep 15, 2018 536 537 538 `````` data = p._get_psfpdb_dictionary() data['resid'] = "%d%c%c" % (data['resid']," "," ") # TODO: work with large indices fh.write( formatString.format(**data) ) `````` cmaffeo2 committed Feb 28, 2018 539 540 541 542 543 544 545 `````` fh.write("\n") ## Write out bonds bonds = self.get_bonds() fh.write("{:>8d} !NBOND\n".format(len(bonds))) counter = 0 for p1,p2,b,ex in bonds: `````` New Tbgl User committed Apr 05, 2018 546 `````` fh.write( "{:>8d}{:>8d}".format(p1.idx+1,p2.idx+1) ) `````` cmaffeo2 committed Feb 28, 2018 547 `````` counter += 1 `````` New Tbgl User committed Apr 05, 2018 548 `````` if counter == 4: `````` cmaffeo2 committed Feb 28, 2018 549 550 `````` fh.write("\n") counter = 0 `````` cmaffeo2 committed Jun 13, 2018 551 552 `````` else: fh.write(" ") `````` 553 `````` fh.write("\n" if counter == 0 else "\n\n") `````` cmaffeo2 committed Feb 28, 2018 554 `````` `````` cmaffeo2 committed Mar 23, 2018 555 556 557 558 559 `````` ## Write out angles angles = self.get_angles() fh.write("{:>8d} !NTHETA\n".format(len(angles))) counter = 0 for p1,p2,p3,a in angles: `````` New Tbgl User committed Apr 05, 2018 560 `````` fh.write( "{:>8d}{:>8d}{:>8d}".format(p1.idx+1,p2.idx+1,p3.idx+1) ) `````` cmaffeo2 committed Mar 23, 2018 561 `````` counter += 1 `````` New Tbgl User committed Apr 05, 2018 562 `````` if counter == 3: `````` cmaffeo2 committed Mar 23, 2018 563 564 `````` fh.write("\n") counter = 0 `````` cmaffeo2 committed Jun 13, 2018 565 566 `````` else: fh.write(" ") `````` 567 `````` fh.write("\n" if counter == 0 else "\n\n") `````` cmaffeo2 committed Mar 23, 2018 568 569 570 571 572 573 `````` ## Write out dihedrals dihedrals = self.get_dihedrals() fh.write("{:>8d} !NPHI\n".format(len(dihedrals))) counter = 0 for p1,p2,p3,p4,a in dihedrals: `````` New Tbgl User committed Apr 05, 2018 574 `````` fh.write( "{:>8d}{:>8d}{:>8d}{:>8d}".format(p1.idx+1,p2.idx+1,p3.idx+1,p4.idx+1) ) `````` cmaffeo2 committed Mar 23, 2018 575 576 577 578 `````` counter += 1 if counter == 2: fh.write("\n") counter = 0 `````` cmaffeo2 committed Jun 13, 2018 579 `````` else: `````` 580 581 `````` fh.write(" ") fh.write("\n" if counter == 0 else "\n\n") `````` cmaffeo2 committed Mar 23, 2018 582 583 584 585 586 587 `````` ## Write out impropers impropers = self.get_impropers() fh.write("{:>8d} !NIMPHI\n".format(len(impropers))) counter = 0 for p1,p2,p3,p4,a in impropers: `````` New Tbgl User committed Apr 05, 2018 588 `````` fh.write( "{:>8d}{:>8d}{:>8d}{:>8d}".format(p1.idx+1,p2.idx+1,p3.idx+1,p4.idx+1) ) `````` cmaffeo2 committed Mar 23, 2018 589 590 591 592 `````` counter += 1 if counter == 2: fh.write("\n") counter = 0 `````` cmaffeo2 committed Jun 13, 2018 593 594 `````` else: fh.write(" ") `````` 595 596 597 598 599 600 601 602 603 604 605 `````` fh.write("\n" if counter == 0 else "\n\n") fh.write("\n{:>8d} !NDON: donors\n\n\n".format(0)) fh.write("\n{:>8d} !NACC: acceptors\n\n\n".format(0)) fh.write("\n 0 !NNB\n\n") natoms = len(self.particles) for i in range(natoms//8): fh.write(" 0 0 0 0 0 0 0 0\n") for i in range(natoms-8*(natoms//8)): fh.write(" 0") fh.write("\n\n 1 0 !NGRP\n\n") `````` cmaffeo2 committed Feb 28, 2018 606 607 608 `````` class ArbdModel(PdbModel): `````` cmaffeo2 committed Jun 14, 2018 609 610 `````` def __init__(self, children, dimensions=(1000,1000,1000), temperature=291, timestep=50e-6, cutoff=50, decompPeriod=10000, pairlistDistance=None, nonbondedResolution=0.1, remove_duplicate_bonded_terms=True): PdbModel.__init__(self, children, dimensions, remove_duplicate_bonded_terms) `````` cmaffeo2 committed Feb 28, 2018 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 `````` self.temperature = temperature self.timestep = timestep self.cutoff = cutoff if pairlistDistance == None: pairlistDistance = cutoff+2 self.decompPeriod = decompPeriod self.pairlistDistance = pairlistDistance self.numParticles = 0 self.particles = [] self.type_counts = None self.nbSchemes = [] self._nbParamFiles = [] # This could be made more robust self.nbResolution = 0.1 self._written_bond_files = dict() self.cacheUpToDate = False `````` cmaffeo2 committed Mar 20, 2018 634 635 636 637 638 639 640 `````` def clear_all(self, keep_children=False): Parent.clear_all(self, keep_children=keep_children) self.particles = [] self.numParticles = 0 self.type_counts = None self._nbParamFiles = [] self._written_bond_files = dict() `````` cmaffeo2 committed Feb 28, 2018 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 `````` def _getNbScheme(self, typeA, typeB): scheme = None for s,A,B in self.nbSchemes: if A is None or B is None: if A is None and B is None: return s elif A is None and typeB == B: return s elif B is None and typeA == A: return s elif typeA == A and typeB == B: return s raise Exception("No nonbonded scheme found for %s and %s" % (typeA.name, typeB.name)) def _countParticleTypes(self): ## TODO: check for modifications to particle that require ## automatic generation of new particle type type_counts = dict() for p in self: t = p.type_ if t in type_counts: type_counts[t] += 1 else: type_counts[t] = 1 self.type_counts = type_counts def _updateParticleOrder(self): ## Create ordered list `````` cmaffeo2 committed Mar 20, 2018 670 671 `````` self.particles = [p for p in self] # self.particles = sorted(particles, key=lambda p: (p.type_, p.idx)) `````` cmaffeo2 committed Feb 28, 2018 672 673 674 675 676 677 678 679 680 681 682 683 `````` ## Update particle indices for p,i in zip(self.particles,range(len(self.particles))): p.idx = i # self.initialCoords = np.array([p.initialPosition for p in self.particles]) def useNonbondedScheme(self, nbScheme, typeA=None, typeB=None): self.nbSchemes.append( (nbScheme, typeA, typeB) ) if typeA != typeB: self.nbSchemes.append( (nbScheme, typeB, typeA) ) `````` cmaffeo2 committed Sep 10, 2018 684 `````` def simulate(self, output_name, output_directory='output', num_steps=100000000, timestep=None, gpu=0, output_period=1e4, arbd=None, directory='.'): `````` cmaffeo2 committed Feb 28, 2018 685 `````` assert(type(gpu) is int) `````` cmaffeo2 committed Sep 10, 2018 686 `````` num_steps = int(num_steps) `````` cmaffeo2 committed Feb 28, 2018 687 `````` `````` cmaffeo2 committed Sep 10, 2018 688 689 690 691 692 `````` d_orig = os.getcwd() try: if not os.path.exists(directory): os.makedirs(directory) os.chdir(directory) `````` New Tbgl User committed Apr 05, 2018 693 `````` `````` cmaffeo2 committed Sep 10, 2018 694 695 `````` if timestep is not None: self.timestep = timestep `````` cmaffeo2 committed Feb 28, 2018 696 `````` `````` cmaffeo2 committed Sep 10, 2018 697 698 699 700 `````` if self.cacheUpToDate == False: # TODO: remove cache? self._countParticleTypes() self._updateParticleOrder() `````` cmaffeo2 committed Sep 10, 2018 701 `````` if output_directory == '': output_directory='.' `````` cmaffeo2 committed Sep 10, 2018 702 703 704 705 706 707 708 709 710 711 `````` if arbd is None: for path in os.environ["PATH"].split(os.pathsep): path = path.strip('"') fname = os.path.join(path, "arbd") if os.path.isfile(fname) and os.access(fname, os.X_OK): arbd = fname break if arbd is None: raise Exception("ARBD was not found") `````` cmaffeo2 committed Feb 28, 2018 712 `````` `````` cmaffeo2 committed Sep 10, 2018 713 714 715 716 717 718 `````` if not os.path.exists(arbd): raise Exception("ARBD was not found") if not os.path.isfile(arbd): raise Exception("ARBD was not found") if not os.access(arbd, os.X_OK): raise Exception("ARBD is not executable") `````` cmaffeo2 committed Feb 28, 2018 719 `````` `````` cmaffeo2 committed Sep 10, 2018 720 721 722 723 `````` if not os.path.exists(output_directory): os.makedirs(output_directory) elif not os.path.isdir(output_directory): raise Exception("output_directory '%s' is not a directory!" % output_directory) `````` cmaffeo2 committed Sep 10, 2018 724 725 `````` `````` cmaffeo2 committed Sep 10, 2018 726 727 728 `````` self.writePdb( output_name + ".pdb" ) self.writePsf( output_name + ".psf" ) self.writeArbdFiles( output_name, numSteps=num_steps, outputPeriod=output_period ) `````` cmaffeo2 committed Sep 10, 2018 729 730 731 `````` os.sync() ## http://stackoverflow.com/questions/18421757/live-output-from-subprocess-command `````` cmaffeo2 committed Sep 10, 2018 732 `````` cmd = (arbd, '-g', "%d" % gpu, "%s.bd" % output_name, "%s/%s" % (output_directory, output_name)) `````` cmaffeo2 committed Sep 10, 2018 733 734 735 736 737 738 739 740 741 742 743 744 `````` cmd = tuple(str(x) for x in cmd) print("Running ARBD with: %s" % " ".join(cmd)) process = subprocess.Popen(cmd, stdout=subprocess.PIPE, universal_newlines=True) for line in process.stdout: sys.stdout.write(line) sys.stdout.flush() except: raise finally: os.chdir(d_orig) `````` cmaffeo2 committed Feb 28, 2018 745 746 747 748 749 750 751 752 `````` # -------------------------- # # Methods for printing model # # -------------------------- # def writeArbdFiles(self, prefix, numSteps=100000000, outputPeriod=10000): ## TODO: save and reference directories and prefixes using member data `````` cmaffeo2 committed Sep 06, 2018 753 `````` d = self.potential_directory = "potentials" `````` cmaffeo2 committed Mar 20, 2018 754 755 `````` if not os.path.exists(d): os.makedirs(d) `````` cmaffeo2 committed Feb 28, 2018 756 757 758 759 760 `````` self._bond_filename = "%s/%s.bonds.txt" % (d, prefix) self._angle_filename = "%s/%s.angles.txt" % (d, prefix) self._dihedral_filename = "%s/%s.dihedrals.txt" % (d, prefix) self._exclusion_filename = "%s/%s.exculsions.txt" % (d, prefix) `````` cmaffeo2 committed Mar 20, 2018 761 762 `````` # self._writeArbdCoordFile( prefix + ".coord.txt" ) self._writeArbdParticleFile( prefix + ".particles.txt" ) `````` cmaffeo2 committed Feb 28, 2018 763 764 765 766 767 768 769 `````` self._writeArbdBondFile() self._writeArbdAngleFile() self._writeArbdDihedralFile() self._writeArbdExclusionFile() self._writeArbdPotentialFiles( prefix, directory = d ) self._writeArbdConf( prefix, numSteps=numSteps, outputPeriod=outputPeriod ) `````` cmaffeo2 committed Mar 20, 2018 770 771 772 773 774 775 `````` # def _writeArbdCoordFile(self, filename): # with open(filename,'w') as fh: # for p in self.particles: # fh.write("%f %f %f\n" % tuple(x for x in p.collapsedPosition())) def _writeArbdParticleFile(self, filename): `````` cmaffeo2 committed Feb 28, 2018 776 777 `````` with open(filename,'w') as fh: for p in self.particles: `````` cmaffeo2 committed Jun 14, 2018 778 `````` data = tuple([p.idx,p.type_.name] + [x for x in p.collapsedPosition()]) `````` cmaffeo2 committed Mar 20, 2018 779 780 `````` fh.write("ATOM %d %s %f %f %f\n" % data) `````` cmaffeo2 committed Feb 28, 2018 781 782 783 784 785 786 787 788 789 790 791 792 793 794 `````` def _writeArbdConf(self, prefix, randomSeed=None, numSteps=100000000, outputPeriod=10000, restartCoordinateFile=None): ## TODO: raise exception if _writeArbdPotentialFiles has not been called filename = "%s.bd" % prefix ## Prepare a dictionary to fill in placeholders in the configuration file params = self.__dict__.copy() # get parameters from System object if randomSeed is None: params['randomSeed'] = "" else: params['randomSeed'] = "seed %s" % randomSeed params['numSteps'] = int(numSteps) `````` cmaffeo2 committed Mar 20, 2018 795 796 `````` # params['coordinateFile'] = "%s.coord.txt" % prefix params['particleFile'] = "%s.particles.txt" % prefix `````` cmaffeo2 committed Feb 28, 2018 797 798 799 800 801 802 803 804 805 `````` if restartCoordinateFile is None: params['restartCoordinates'] = "" else: params['restartCoordinates'] = "restartCoordinates %s" % restartCoordinateFile params['outputPeriod'] = outputPeriod for k,v in zip('XYZ', self.dimensions): params['origin'+k] = -v*0.5 params['dim'+k] = v `````` cmaffeo2 committed Apr 01, 2018 806 807 `````` params['pairlistDistance'] -= params['cutoff'] `````` cmaffeo2 committed Feb 28, 2018 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 `````` ## Actually write the file with open(filename,'w') as fh: fh.write("""{randomSeed} timestep {timestep} steps {numSteps} numberFluct 0 # deprecated interparticleForce 1 # other values deprecated fullLongRange 0 # deprecated temperature {temperature} outputPeriod {outputPeriod} ## Energy doesn't actually get printed! outputEnergyPeriod {outputPeriod} outputFormat dcd ## Infrequent domain decomposition because this kernel is still very slow decompPeriod {decompPeriod} cutoff {cutoff} pairlistDistance {pairlistDistance} origin {originX} {originY} {originZ} systemSize {dimX} {dimY} {dimZ} \n""".format(**params)) ## Write entries for each type of particle for pt,num in self.getParticleTypesAndCounts(): `````` cmaffeo2 committed Apr 01, 2018 836 `````` ## TODO create new particle types if existing has grid `````` cmaffeo2 committed Feb 28, 2018 837 838 839 840 841 842 843 844 `````` particleParams = pt.__dict__.copy() particleParams['num'] = num fh.write(""" particle {name} num {num} diffusion {diffusivity} """.format(**particleParams)) if 'grid' in particleParams: `````` cmaffeo2 committed Apr 01, 2018 845 `````` if not isinstance(pt.grid, list): pt.grid = [pt.grid] `````` cmaffeo2 committed Sep 17, 2018 846 847 848 849 850 851 852 853 `````` 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)) fh.write("gridFileScale {}\n".format(s)) `````` cmaffeo2 committed Apr 01, 2018 854 `````` `````` cmaffeo2 committed Feb 28, 2018 855 `````` else: `````` cmaffeo2 committed Sep 06, 2018 856 `````` fh.write("gridFile {}/null.dx\n".format(self.potential_directory)) `````` cmaffeo2 committed Feb 28, 2018 857 858 859 860 `````` ## Write coordinates and interactions fh.write(""" ## Input coordinates `````` cmaffeo2 committed Mar 20, 2018 861 ``````inputParticles {particleFile} `````` cmaffeo2 committed Feb 28, 2018 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 ``````{restartCoordinates} ## 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 """.format(**params)) for pair,f in zip(self._particleTypePairIter(), self._nbParamFiles): i,j,t1,t2 = pair fh.write("tabulatedFile %d@%d@%s\n" % (i,j,f)) ## Bonded interactions bonds = self.get_bonds() angles = self.get_angles() dihedrals = self.get_dihedrals() exclusions = self.get_exclusions() if len(bonds) > 0: for b in list(set([b for i,j,b,ex in bonds])): fh.write("tabulatedBondFile %s\n" % b) if len(angles) > 0: for b in list(set([b for i,j,k,b in angles])): fh.write("tabulatedAngleFile %s\n" % b) if len(dihedrals) > 0: for b in list(set([b for i,j,k,l,b in dihedrals])): fh.write("tabulatedDihedralFile %s\n" % b) fh.write("inputBonds %s\n" % self._bond_filename) fh.write("inputAngles %s\n" % self._angle_filename) fh.write("inputDihedrals %s\n" % self._dihedral_filename) fh.write("inputExcludes %s\n" % self._exclusion_filename) write_null_dx = False for pt,num in self.getParticleTypesAndCounts(): if "grid" not in pt.__dict__: `````` cmaffeo2 committed Sep 06, 2018 898 899 `````` gridfile = "{}/null.dx".format(self.potential_directory) with open(gridfile, 'w') as fh: `````` cmaffeo2 committed Feb 28, 2018 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 `````` fh.write("""object 1 class gridpositions counts 2 2 2 origin {originX} {originY} {originZ} delta {dimX} 0.000000 0.000000 delta 0.000000 {dimY} 0.000000 delta 0.000000 0.000000 {dimZ} 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 """.format(**params)) break def getParticleTypesAndCounts(self): ## TODO: remove(?) return sorted( self.type_counts.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"): 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 _writeNonbondedParameterFiles(self, prefix): x = np.arange(0, self.cutoff, self.nbResolution) for i,j,t1,t2 in self._particleTypePairIter(): f = "%s.%s-%s.dat" % (prefix, t1.name, t2.name) scheme = self._getNbScheme(t1,t2) scheme.write_file(f, t1, t2, rMax = self.cutoff) self._nbParamFiles.append(f) def _getNonbondedPotential(self,x,a,b): return a*(np.exp(-x/b)) def _writeArbdBondFile( self ): for b in list( set( [b for i,j,b,ex in self.get_bonds()] ) ): if type(b) is not str: b.write_file() with open(self._bond_filename,'w') as fh: for i,j,b,ex in self.get_bonds(): item = (i.idx, j.idx, str(b)) if ex: fh.write("BOND REPLACE %d %d %s\n" % item) else: fh.write("BOND ADD %d %d %s\n" % item) def _writeArbdAngleFile( self ): for b in list( set( [b for i,j,k,b in self.get_angles()] ) ): if type(b) is not str: b.write_file() with open(self._angle_filename,'w') as fh: for b in self.get_angles(): item = tuple([p.idx for p in b[:-1]] + [str(b[-1])]) fh.write("ANGLE %d %d %d %s\n" % item) def _writeArbdDihedralFile( self ): for b in list( set( [b for i,j,k,l,b in self.get_dihedrals()] ) ): if type(b) is not str: b.write_file() with open(self._dihedral_filename,'w') as fh: for b in self.get_dihedrals(): item = tuple([p.idx for p in b[:-1]] + [str(b[-1])]) fh.write("DIHEDRAL %d %d %d %d %s\n" % item) def _writeArbdExclusionFile( self ): with open(self._exclusion_filename,'w') as fh: for ex in self.get_exclusions(): item = tuple(int(p.idx) for p in ex) fh.write("EXCLUDE %d %d\n" % item) `````` cmaffeo2 committed Mar 23, 2018 992 `````` `````` 993 994 995 996 `````` def set_dimensions_from_structure( self, padding_factor=1.5, isotropic=False ): ## TODO: cache coordinates using numpy arrays for quick min/max raise(NotImplementedError) `````` cmaffeo2 committed Sep 10, 2018 997 `````` def write_namd_configuration( self, output_name, num_steps = 1e6, `````` 998 999 1000 `````` output_directory = 'output', update_dimensions=True, extrabonds=True ): ``````
For faster browsing, not all history is shown.