1#!/usr/bin/env python 2 3# Author: Andrew Jewett (jewett.aij at g mail) 4# License: MIT License (See LICENSE.md) 5# Copyright (c) 2017, California Institute of Technology 6# All rights reserved. 7 8""" 9 Generate a moltemplate (.lt) file containing a definition of a polymer 10 molecule whose monomers are located at the positions specified in 11 "coords.raw" (a 3-column text file). Monomers will be rotated so 12 that they point in the direction connecting neighbors (r[i+1]-r[i]) 13 The user can specify the subunits to use when building the polymer, 14 the atoms to to build bonds (and angles, and dihedrals) between monomers 15 and the helical pitch of the polymer. The output of this program is 16 a text file in moltemplate (.lt) format containing the sequence of 17 moltemplate commands needed to build this polymer molecule(s). (One must 18 then run moltemplate on this file to build the LAMMPS simulation files.) 19 Multiple Polymers: 20 To make it easier to create polymer melts, multiple polymers can be created 21 from coordinates in the same file by using the "-cuts" command line argument. 22 Encapsulation: 23 If the "-polymer-name PolyName" command line option is given, then these 24 moltemplate commands will be nested within the definition of a moltemplate 25 object (named "PolyName", in this example. Later in your moltemplate files, 26 you must remember to instantiate a copy of this moltemplate object using 27 a command like "polymer = new PolyName" Atoms within this object will 28 share the same molecule-ID number.) If multiple polymers are requested, then 29 each of them will have their own polymer object. 30 31""" 32 33 34g_usage_msg = """ 35Usage: 36 37 genpoly_lt.py \\ 38 [-polymer-name pname] \\ 39 [-monomer-name mname] \\ 40 [-sequence sequence.txt] \\ 41 [-bond a1 a2] \\ 42 [-bond btype a1 a2] \\ 43 [-angle atype a1 a2 a3 i1 i2 i3] \\ 44 [-dihedral dtype a1 a2 a3 a4 i1 i2 i3 i4] \\ 45 [-improper itype a1 a2 a3 a4 i1 i2 i3 i4] \\ 46 [-inherits ForceFieldObject] \\ 47 [-header "import \"monomer.lt\""] \\ 48 [-helix deltaphi] \\ 49 [-axis x,y,z] \\ 50 [-circular yes/no/connected] \\ 51 [-cuts cuts.txt] \\ 52 [-polymer-directions polarities.txt \\ 53 [-dir-indices ia ib] \\ 54 [-padding paddingX,paddingY,paddingZ] \\ 55 < coords.raw > polymer.lt 56 57""" 58 59 60import sys 61import random 62from math import * 63 64 65class InputError(Exception): 66 """ A generic exception object containing a string for error reporting. 67 (Raising this exception implies that the caller has provided 68 a faulty input file or argument.) 69 70 """ 71 72 def __init__(self, err_msg): 73 self.err_msg = err_msg 74 75 def __str__(self): 76 return self.err_msg 77 78 def __repr__(self): 79 return str(self) 80 81 82 83class GPSettings(object): 84 85 def __init__(self): 86 self.infile_name = '' 87 self.direction_orig = [1.0, 0.0, 0.0] 88 self.is_circular = False 89 self.connect_ends = False 90 self.delta_phi = 0.0 91 #self.header = 'import \"forcefield.lt\"' 92 self.header = '' 93 self.name_monomer = 'Monomer' 94 self.name_sequence_file = '' 95 self.name_sequence_multi = [] 96 self.name_polymer = '' 97 self.inherits = '' 98 self.dir_index_offsets = (-1,1) 99 self.cuts = [] 100 self.reverse_polymer_directions = [] 101 self.box_padding = None 102 self.bonds_name = [] 103 self.bonds_type = [] 104 self.bonds_atoms = [] 105 self.bonds_index_offsets = [] 106 self.bonds_name_notype = [] 107 self.bonds_atoms_notype = [] 108 self.bonds_index_offsets_notype = [] 109 self.angles_name = [] 110 self.angles_type = [] 111 self.angles_atoms = [] 112 self.angles_index_offsets = [] 113 self.dihedrals_name = [] 114 self.dihedrals_type = [] 115 self.dihedrals_atoms = [] 116 self.dihedrals_index_offsets = [] 117 self.impropers_name = [] 118 self.impropers_type = [] 119 self.impropers_atoms = [] 120 self.impropers_index_offsets = [] 121 self.helix_angles_file = '' 122 self.orientations_file = '' 123 self.orientations_use_quats = False 124 125 def ParseArgs(self, argv): 126 i = 1 127 while i < len(argv): 128 #sys.stderr.write('argv['+str(i)+'] = \"'+argv[i]+'\"\n') 129 130 if ((argv[i].lower() == '-in') or 131 (argv[i].lower() == '-i')): 132 if i + 1 >= len(argv): 133 raise InputError( 134 'Error: The ' + argv[i] + ' flag should be followed by a file name.\n') 135 self.infile_name = argv[i + 1] 136 del(argv[i:i + 2]) 137 elif argv[i].lower() == '-bond': 138 if i + 2 >= len(argv): 139 raise InputError( 140 'Error: The ' + argv[i] + ' flag should be followed by at 2 or 3 strings.\n') 141 has_third_arg = True 142 if ((i + 3 >= len(argv)) or (argv[i+3][0:1] == '-')): 143 has_third_arg = False 144 # Did the user specify the type of the bond? 145 if has_third_arg: 146 # then the type was specified (put in "Data Bonds") 147 self.bonds_type.append(argv[i + 1]) 148 self.bonds_atoms.append((argv[i + 2], 149 argv[i + 3])) 150 self.bonds_index_offsets.append((0, 1)) 151 del(argv[i:i + 4]) 152 else: 153 # then the type was not specified (put in "Data Bond List") 154 self.bonds_atoms_notype.append((argv[i + 1], 155 argv[i + 2])) 156 self.bonds_index_offsets_notype.append((0, 1)) 157 del(argv[i:i + 3]) 158 elif argv[i].lower() == '-angle': 159 if i + 7 >= len(argv): 160 raise InputError( 161 'Error: The ' + argv[i] + ' flag should be followed by 4 strings and 3 integers.\n') 162 # self.angles_name.append(argv[i+1]) 163 self.angles_type.append(argv[i + 1]) 164 self.angles_atoms.append((argv[i + 2], 165 argv[i + 3], 166 argv[i + 4])) 167 self.angles_index_offsets.append((int(argv[i + 5]), 168 int(argv[i + 6]), 169 int(argv[i + 7]))) 170 if ((self.angles_index_offsets[-1][0] < 0) or 171 (self.angles_index_offsets[-1][1] < 0) or 172 (self.angles_index_offsets[-1][2] < 0)): 173 raise InputError( 174 'Error: ' + argv[i] + ' indices (i1 i2 i3) must be >= 0\n') 175 del(argv[i:i + 8]) 176 elif argv[i].lower() == '-dihedral': 177 if i + 9 >= len(argv): 178 raise InputError( 179 'Error: The ' + argv[i] + ' flag should be followed by 5 strings and 4 integers.\n') 180 # self.dihedrals_name.append(argv[i+1]) 181 self.dihedrals_type.append(argv[i + 1]) 182 self.dihedrals_atoms.append((argv[i + 2], 183 argv[i + 3], 184 argv[i + 4], 185 argv[i + 5])) 186 self.dihedrals_index_offsets.append((int(argv[i + 6]), 187 int(argv[i + 7]), 188 int(argv[i + 8]), 189 int(argv[i + 9]))) 190 if ((self.dihedrals_index_offsets[-1][0] < 0) or 191 (self.dihedrals_index_offsets[-1][1] < 0) or 192 (self.dihedrals_index_offsets[-1][2] < 0) or 193 (self.dihedrals_index_offsets[-1][3] < 0)): 194 raise InputError( 195 'Error: ' + argv[i] + ' indices (i1 i2 i3 i4) must be >= 0\n') 196 del(argv[i:i + 10]) 197 elif argv[i].lower() == '-improper': 198 if i + 9 >= len(argv): 199 raise InputError( 200 'Error: The ' + argv[i] + ' flag should be followed by 5 strings and 4 integers.\n') 201 # self.impropers_name.append(argv[i+1]) 202 self.impropers_type.append(argv[i + 1]) 203 self.impropers_atoms.append((argv[i + 2], 204 argv[i + 3], 205 argv[i + 4], 206 argv[i + 5])) 207 self.impropers_index_offsets.append((int(argv[i + 6]), 208 int(argv[i + 7]), 209 int(argv[i + 8]), 210 int(argv[i + 9]))) 211 if ((self.impropers_index_offsets[-1][0] < 0) or 212 (self.impropers_index_offsets[-1][1] < 0) or 213 (self.impropers_index_offsets[-1][2] < 0) or 214 (self.impropers_index_offsets[-1][3] < 0)): 215 raise InputError( 216 'Error: ' + argv[i] + ' indices (i1 i2 i3 i4) must be >= 0\n') 217 del(argv[i:i + 10]) 218 elif (argv[i].lower() == '-monomer-name'): 219 if i + 1 >= len(argv): 220 raise InputError( 221 'Error: The ' + argv[i] + ' flag should be followed by a string\n') 222 self.name_monomer = argv[i + 1] 223 del(argv[i:i + 2]) 224 elif (argv[i].lower() == '-sequence'): 225 if i + 1 >= len(argv): 226 raise InputError( 227 'Error: The ' + argv[i] + ' flag should be followed by a file name\n') 228 self.name_sequence_file = argv[i + 1] 229 del(argv[i:i + 2]) 230 231 elif (argv[i].lower() == '-cuts'): 232 if i + 1 >= len(argv): 233 raise InputError( 234 'Error: The ' + argv[i] + ' flag should be followed by a file name\n') 235 try: 236 f = open(argv[i + 1], "r") 237 except IOError: 238 raise InputError( 239 'Error: File "' + argv[i + 1] + '" could not be opened for reading\n') 240 for line_orig in f: 241 line = line_orig.strip() 242 ic = line.find('#') 243 if ic != -1: 244 line = line[:ic] 245 else: 246 line = line.strip() 247 if len(line) > 0: 248 try: 249 self.cuts.append(int(line)) 250 except ValueError: 251 raise InputError( 252 'Error: file ' + argv[i + 1] + ' should contain only nonnegative integers.\n') 253 del(argv[i:i + 2]) 254 elif (argv[i].lower() == '-polymer-directions'): 255 if i + 1 >= len(argv): 256 raise InputError( 257 'Error: The ' + argv[i] + ' flag should be followed by a file name\n') 258 try: 259 f = open(argv[i + 1], "r") 260 except IOError: 261 raise InputError( 262 'Error: File "' + argv[i + 1] + '" could not be opened for reading\n') 263 for line_orig in f: 264 line = line_orig.strip() 265 ic = line.find('#') 266 if ic != -1: 267 line = line[:ic] 268 else: 269 line = line.strip() 270 if len(line) > 0: 271 try: 272 entry = line 273 if ((entry == '1') or 274 (entry == '+1') or 275 (entry == 'true') or 276 (entry == 'increasing')): 277 self.reverse_polymer_directions.append(False) 278 else: 279 self.reverse_polymer_directions.append(True) 280 except ValueError: 281 raise InputError( 282 'Error: file ' + argv[i + 1] + ' should contain only \"+1\" or \"-1\" on each line.\n') 283 del(argv[i:i + 2]) 284 elif (argv[i].lower() == '-polymer-name'): 285 if i + 1 >= len(argv): 286 raise InputError( 287 'Error: The ' + argv[i] + ' flag should be followed by a string\n') 288 self.name_polymer = argv[i + 1] 289 del(argv[i:i + 2]) 290 elif (argv[i].lower() == '-inherits'): 291 if i + 1 >= len(argv): 292 raise InputError( 293 'Error: The ' + argv[i] + ' flag should be followed by a string\n') 294 self.inherits = argv[i + 1] 295 if self.inherits.find('inherits ') == 0: 296 self.inherits = ' ' + self.inherits 297 else: 298 self.inherits = ' inherits ' + self.inherits 299 if self.name_polymer == '': 300 self.name_polymer = 'Polymer' # supply a default name 301 del(argv[i:i + 2]) 302 elif (argv[i].lower() == '-header'): 303 if i + 1 >= len(argv): 304 raise InputError( 305 'Error: The ' + argv[i] + ' flag should be followed by a string (usually in quotes)\n') 306 if self.header != '': 307 self.header += '\n' 308 self.header += argv[i + 1] 309 del(argv[i:i + 2]) 310 elif argv[i].lower() == '-axis': 311 if i + 1 >= len(argv): 312 raise InputError('Error: The ' + argv[i] + ' flag should be followed ' + 313 'by 3 numbers separated by commas (no spaces)\n') 314 self.direction_orig = list(map(float, argv[i + 1].split(','))) 315 del(argv[i:i + 2]) 316 elif argv[i].lower() == '-circular': 317 if i + 1 >= len(argv): 318 raise InputError('Error: The ' + argv[i] + ' flag should be followed by an argument\n' + 319 ' ("yes", "no", or "connected")\n') 320 if argv[i + 1].lower() == 'yes': 321 self.connect_ends = True 322 self.is_circular = True 323 elif argv[i + 1].lower() == 'connected': 324 self.connect_ends = True 325 self.is_circular = False 326 elif argv[i + 1].lower() == 'no': 327 self.connect_ends = False 328 self.is_circular = False 329 else: 330 raise InputError('Error: The ' + argv[i] + ' flag should be followed by an argument\n' + 331 ' ("yes", "no", or "connected")\n') 332 del(argv[i:i + 2]) 333 elif argv[i].lower() == '-helix': 334 if i + 1 >= len(argv): 335 raise InputError( 336 'Error: The ' + argv[i] + ' flag should be followed by a number (angle in degrees)\n') 337 self.delta_phi = float(argv[i + 1]) 338 del(argv[i:i + 2]) 339 elif (argv[i].lower() == '-dir-indices'): 340 if i + 2 >= len(argv): 341 raise InputError( 342 'Error: The ' + argv[i] + ' flag should be followed by two integers\n') 343 self.dir_index_offsets = (int(argv[i + 1]), int(argv[i + 2])) 344 if self.dir_index_offsets[0] == self.dir_index_offsets[1]: 345 raise InputError( 346 'Error: The two numbers following ' + argv[i] + ' must not be equal.\n') 347 del(argv[i:i + 3]) 348 elif ((argv[i].lower() == '-padding') or 349 (argv[i].lower() == '-box')): 350 if i + 1 >= len(argv): 351 raise InputError('Error: The ' + argv[i] + ' flag should be followed\n' + 352 'by 3 numbers separated by commas (no spaces)\n') 353 self.box_padding = list(map(float, argv[i + 1].split(','))) 354 if len(self.box_padding) == 1: 355 self.box_padding = [self.box_padding] * 3 356 del(argv[i:i + 2]) 357 elif (argv[i].lower() in ('-orientations','-orientations')): 358 self.orientations_file = argv[i+1] 359 self.orientations_use_quats = False 360 del(argv[i:i + 2]) 361 elif (argv[i].lower() in ('-quaternion','-quaternions')): 362 self.orientations_file = argv[i+1] 363 self.orientations_use_quats = True 364 del(argv[i:i + 2]) 365 elif (argv[i].lower() in ('-helix-angle','-helix-angles')): 366 self.helix_angles_file = argv[i+1] 367 del(argv[i:i + 2]) 368 369 # elif ((argv[i][0] == '-') and (__name__ == '__main__')): 370 # 371 # raise InputError('Error('+g_program_name+'):\n'+\ 372 # 'Unrecogized command line argument \"'+argv[i]+\ 373 # '\"\n\n'+\ 374 # __doc__) 375 else: 376 i += 1 377 378 379 380 381 382 383class WrapPeriodic(object): 384 """ Wrap() calculates the remainder of i % N. 385 It turns out to be convenient to do this multiple times and later 386 query whether i/N != 0 in any of them once (by checking bounds_err). 387 388 """ 389 bounds_err = False 390 391 @classmethod 392 def Wrap(obj, i, N): 393 if i // N != 0: 394 obj.bounds_err = True 395 return i % N 396 397 def WrapF(obj, x, L): 398 i = floor(x / L) 399 if i != 0: 400 obj.bounds_err = True 401 return x - i * L 402 403 404class GenPoly(object): 405 """ 406 Read coordinates from a file, and generate a list of \"new\" commands 407 in moltemplate format with the position of each monomer located 408 at these positions, oriented appropriately, with bonds (and angles, 409 dihedrals, etc...) connecting successive monomers together. 410 By default, only a single polymer is created. 411 However this class can create multiple polymers of different lengths. 412 The list of coordinates for each polymer are saved separately within 413 the "self.coords_multi" member. 414 415 """ 416 417 def __init__(self): 418 self.settings = GPSettings() 419 self.coords_multi = [] # a list-of-list-of-lists of numbers Nxnx3 420 self.name_sequence_multi = [] 421 self.direction_vects = [] 422 self.box_bounds_min = None 423 self.box_bounds_max = None 424 self.N = 0 425 self.orientations_multi = [] 426 self.helix_angles_multi = [] 427 428 429 def ParseArgs(self, argv): 430 """ 431 parse the argument list 432 """ 433 # The command above will remove arguments from argv which are 434 # understood by GPSettings.ParseArgs(argv). 435 # The remaining arguments will be handled below. 436 self.settings.ParseArgs(argv) 437 438 439 # ---- Parsing is finished. Now load files and cleanup. ---- 440 441 if self.settings.infile_name != '': 442 infile = open(self.settings.infile_name, 'r') 443 self.coords_multi = self.ReadCoords(infile) 444 infile.close() 445 446 if ((len(self.settings.reverse_polymer_directions) != 0) and 447 (len(self.settings.reverse_polymer_directions) != 448 len(self.cuts) + 1)): 449 raise InputError('Error: The number of entries in the file you supplied to "-polymer-directions"\n' 450 ' does not equal the number of polymers (which is either 1, or 1 + the\n' 451 ' number of entries in the file you supplied to "-cuts", if applicable)\n') 452 453 454 for b in range(0, len(self.settings.bonds_type)): 455 if len(self.settings.bonds_type) > 1: 456 self.settings.bonds_name.append('genp_bond' + str(b + 1) + '_') 457 else: 458 self.settings.bonds_name.append('genp_bond_') 459 for b in range(0, len(self.settings.bonds_atoms_notype)): 460 if len(self.settings.bonds_atoms_notype) > 1: 461 self.settings.bonds_name_notype.append('genp_bond_nt' + str(b + 1) + '_') 462 else: 463 self.settings.bonds_name_notype.append('genp_bond_nt') 464 for b in range(0, len(self.settings.angles_type)): 465 if len(self.settings.angles_type) > 1: 466 self.settings.angles_name.append('genp_angle' + str(b + 1) + '_') 467 else: 468 self.settings.angles_name.append('genp_angle_') 469 for b in range(0, len(self.settings.dihedrals_type)): 470 if len(self.settings.dihedrals_type) > 1: 471 self.settings.dihedrals_name.append('genp_dihedral' + str(b + 1) + '_') 472 else: 473 self.settings.dihedrals_name.append('genp_dihedral_') 474 for b in range(0, len(self.settings.impropers_type)): 475 if len(self.settings.impropers_type) > 1: 476 self.settings.impropers_name.append('genp_improper' + str(b + 1) + '_') 477 else: 478 self.settings.impropers_name.append('genp_improper_') 479 480 481 482 # Did the user ask us to read a custom sequence of monomer type names? 483 name_sequence_file = None 484 if isinstance(self.settings.name_sequence_file, str): 485 if self.settings.name_sequence_file != '': 486 name_sequence_file = open(self.settings.name_sequence_file,'r') 487 if name_sequence_file: 488 # Note: This will fill the contents of self.name_sequence_multi 489 self.ReadSequence(name_sequence_file) 490 491 # Did the user specify a file with a list of orientations? 492 orientations_file = None 493 if isinstance(self.settings.orientations_file, str): 494 if self.settings.orientations_file != '': 495 orientations_file = open(self.settings.orientations_file, 'r') 496 else: 497 orientations_file = self.settings.orientations_file 498 if orientations_file: 499 if self.settings.orientations_use_quats: 500 self.orientations_multi = self.ReadCoords(orientations_file, 4) 501 else: 502 self.orientations_multi = self.ReadCoords(orientations_file, 9) 503 504 # Did the user specify a file with a list of helix angles? 505 helix_angles_file = None 506 if isinstance(self.settings.helix_angles_file, str): 507 if self.settings.helix_angles_file != '': 508 helix_angles_file = open(self.settings.helix_angles_file, 'r') 509 else: 510 helix_angles_file = self.settings.helix_angles_file 511 if helix_angles_file: 512 self.helix_angles_multi = self.ReadCoords(self.settings.helix_angles_file, 1) 513 # Because I borrowed the ReadCoords() function to read the file, 514 # each "angle" is ended up being a list containing only one element 515 # (the angle). Convert this list into the number stored there. 516 for i in range(0, len(self.helix_angles_multi)): 517 for j in range(0, len(self.helix_angles_multi[i])): 518 assert(len(self.helix_angles_multi[i][j]) == 1) 519 self.helix_angles_multi[i][j] = self.helix_angles_multi[i][j][0] 520 521 522 523 524 def ReadCoords(self, infile, ncolumns=3): 525 """ 526 Read x y z coordinates from a multi-column ASCII text file. 527 Store the coordinates in the self.coords_multi[][][] list. 528 (Note: Information stored in self.settings.cuts and 529 self.settings.reverse_polymer_directions 530 is used to split the coordinate list into multiple lists 531 and to determine the order that coordinates are read.) 532 This function can also be used to read other (multi-column) ASCII text 533 files by overriding the default value for "ncolumns". (For example, 534 to read files with 4 numbers on each line, set ncolumns=4.) 535 """ 536 537 filename = '' 538 if isinstance(infile, str): 539 filename = infile 540 try: 541 infile = open(filename, 'r') 542 except IOError: 543 raise InputError( 544 'Error: file ' + filename + 545 ' could not be opened for reading\n') 546 547 coords_multi = [] # (do not confuse this with "self.coords_multi") 548 coords = [] 549 550 lines = infile.readlines() 551 for i in range(0, len(lines)): 552 tokens = lines[i].strip().split() 553 if (len(tokens) == ncolumns): 554 coords.append(list(map(float, tokens))) 555 556 self.N = len(coords) 557 if self.N < 2: 558 err_msg = 'Error: Coordinate file must have at least 2 positions' 559 raise InputError(err_msg+'.\n') 560 561 # Did the caller ask us to split the polymer into multiple polymers? 562 if len(self.settings.cuts) > 0: 563 self.settings.cuts.append(self.N) 564 self.settings.cuts.sort() 565 i = 0 566 for j in self.settings.cuts: 567 if ((j-i < 1) or (j > self.N)): 568 err_msg = 'Error in "-cuts" argument: You have bad numbers in the "-cuts" file. This\n' + \ 569 ' could cause one or more the polymers to have zero length. The numbers\n' + \ 570 ' in the "-cuts" file must be in increasing order and must be in\n' + \ 571 ' the range from 1 to N-1 (where "N" is the sum of the number of monomers\n' + \ 572 ' in all of the polymers, which is '+str(self.N)+' in this case).\n' + \ 573 ' No integer can be listed more than once.\n' 574 raise InputError(err_msg+'.\n') 575 coords_multi.append(coords[i:j]) 576 i = j 577 else: 578 coords_multi.append(coords) 579 580 # Did the caller ask us to reverse the direction of any polymers? 581 for i in range(0, len(self.settings.reverse_polymer_directions)): 582 assert(i < len(self.coords_multi)) 583 if self.settings.reverse_polymer_directions[i]: 584 self.coords_multi[i].reverse() 585 586 if filename != '': 587 # Then we opened a new file with that name. We should close it now. 588 infile.close() 589 590 return coords_multi 591 592 593 594 def ReadSequence(self, infile): 595 """ 596 Read a sequence of monomer type names from a file. 597 This function is similar to ReadCoords(). 598 """ 599 name_sequence = [] 600 601 filename = '' 602 if isinstance(infile, str): 603 filename = infile 604 try: 605 infile = open(filename, 'r') 606 except IOError: 607 raise InputError( 608 'Error: file ' + filename + 609 ' could not be opened for reading\n') 610 611 for line_orig in infile: 612 line = line_orig.strip() 613 ic = line.find('#') 614 if ic != -1: 615 line = line[:ic] 616 else: 617 line = line.strip() 618 if len(line) > 0: 619 name_sequence.append(line) 620 621 N = len(name_sequence) 622 623 # Did the caller ask us to split the polymer into multiple polymers? 624 if len(self.settings.cuts) > 0: 625 if (self.settings.cuts[-1] < N): 626 self.settings.cuts.append(N) 627 self.settings.cuts.sort() 628 i = 0 629 for j in self.settings.cuts: 630 self.name_sequence_multi.append(name_sequence[i:j]) 631 i = j 632 else: 633 self.name_sequence_multi.append(name_sequence) 634 635 # Did the caller ask us to reverse the direction of any polymers? 636 for i in range(0, len(self.settings.reverse_polymer_directions)): 637 if self.settings.reverse_polymer_directions[i]: 638 self.name_sequence_multi[i].reverse() 639 640 if filename != '': 641 # Then we opened a new file with that name. We should close it now. 642 infile.close() 643 644 645 646 def ChooseDirections(self, coords): 647 """ 648 Calculate the direction each monomer subunit should be pointing at: 649 650 """ 651 652 N = len(coords) 653 654 if N == 1: 655 self.direction_vects = [[1.0, 0.0, 0.0]] 656 return 657 658 self.direction_vects = [[0.0, 0.0, 0.0] for i in range(0, N + 1)] 659 660 if self.settings.is_circular: 661 for i in range(0, N): 662 # By default, the direction that monomer "i" is pointing is 663 # determined by the position of the monomers before and after it 664 # (at index i-1, and i+1). More generally, we allow the user 665 # to choose what these offsets are ("dir_index_offsets[") 666 ia = WrapPeriodic.Wrap(i + self.settings.dir_index_offsets[0], 667 N) 668 ib = WrapPeriodic.Wrap(i + self.settings.dir_index_offsets[1], 669 N) 670 for d in range(0, 3): 671 self.direction_vects[i][d] = coords[ib][d] - coords[ia][d] 672 673 else: 674 for i in range(1, N - 1): 675 for d in range(0, 3): 676 self.direction_vects[i][d] = (coords[i + self.settings.dir_index_offsets[1]][d] 677 - 678 coords[i + self.settings.dir_index_offsets[0]][d]) 679 680 for d in range(0, 3): 681 self.direction_vects[0][d] = coords[1][d] - coords[0][d] 682 self.direction_vects[N-1][d] = coords[N-1][d] - coords[N-2][d] 683 684 # Optional: normalize the direction vectors 685 686 for i in range(0, N): 687 direction_len = 0.0 688 for d in range(0, 3): 689 direction_len += (self.direction_vects[i][d])**2 690 direction_len = sqrt(direction_len) 691 for d in range(0, 3): 692 self.direction_vects[i][d] /= direction_len 693 694 # Special case: self.direction_vects[-1] is the direction that the original monomer 695 # in "monomer.lt" was pointing. (By default, 1,0,0 <--> the "x" 696 # direction) 697 698 self.direction_vects[-1] = self.settings.direction_orig 699 700 701 702 def WriteLTFile(self, outfile): 703 """ Write an moltemplate (.lt) file containing the definition of 704 this polymer object. (If multiple polymer objects were requested by 705 the user (using the -cuts argument), then their definitions will 706 appear nested within this object, and each of them will be 707 instantiated once when the parent object is instantiated.) 708 709 """ 710 711 # make sure len(self.orientations_multi) == len(self.coords_multi) 712 if len(self.orientations_multi) == 0: 713 self.orientations_multi = [[] for entry in self.coords_multi] 714 # make sure len(self.helix_angles_multi) == len(self.coords_multi) 715 if len(self.helix_angles_multi) == 0: 716 self.helix_angles_multi = [[] for entry in self.coords_multi] 717 718 outfile.write(self.settings.header + "\n\n\n") 719 720 if len(self.coords_multi) == 1: 721 self.WritePolymer(outfile, 722 self.settings.name_polymer + 723 self.settings.inherits, 724 self.coords_multi[0], 725 self.name_sequence_multi[0], 726 self.orientations_multi[0], 727 self.helix_angles_multi[0]) 728 else: 729 if self.settings.name_polymer != '': 730 outfile.write(self.settings.name_polymer + " {\n\n") 731 outfile.write('# Definitions of individual polymers to follow\n\n') 732 for i in range(0, len(self.coords_multi)): 733 # insert ".../" in front of the inherits string 734 ih_str = self.settings.inherits 735 ih = ih_str.find('inherits ') 736 if ih != -1: 737 ih += len('inherits ') 738 ih_str = ih_str[0:ih] + '.../' + ih_str[ih:] 739 self.WritePolymer(outfile, 740 self.settings.name_polymer + '_sub' + str(i + 1) + 741 ih_str, 742 self.coords_multi[i], 743 self.name_sequence_multi[i], 744 self.orientations_multi[i], 745 self.helix_angles_multi[i]) 746 outfile.write('\n\n') 747 outfile.write('# Now instantiate all the polymers (once each)\n\n') 748 749 for i in range(0, len(self.coords_multi)): 750 outfile.write('polymers[' + str(i) + '] = new ' + 751 self.settings.name_polymer + '_sub' + str(i + 1) + '\n') 752 753 if self.settings.name_polymer != '': 754 outfile.write('\n\n' 755 '} # ' + self.settings.name_polymer + '\n\n') 756 757 if self.settings.box_padding != None: 758 for i in range(0, len(self.coords_multi)): 759 # calculate the box big enough to collectively enclose 760 # all of the coordinates (even multiple coordinate sets) 761 self.CalcBoxBoundaries(self.coords_multi[i]) 762 self.WriteBoxBoundaries(outfile) 763 764 765 766 def WritePolymer(self, 767 outfile, 768 name_polymer, 769 coords, 770 names_monomers, 771 orientations=[], 772 helix_angles=[]): 773 """ Write a single polymer object to a file. 774 This function is invoked by WriteLTFile() 775 776 """ 777 N = len(coords) 778 779 self.ChooseDirections(coords) 780 781 if name_polymer != '': 782 outfile.write(name_polymer + ' {\n' 783 '\n\n\n' 784 ' create_var {$mol}\n' 785 ' # The line above forces all monomer subunits to share the same molecule-ID\n' 786 ' # (Note: The molecule-ID number is optional and is usually ignored by LAMMPS.)\n\n\n\n') 787 788 outfile.write(""" 789 # ----- list of monomers: ----- 790 # 791 # (Note: move(), rot(), and rotvv() commands control the position 792 # of each monomer. (See the moltemplate manual for an explanation 793 # of what they do.) Commands enclosed in push() are cumulative 794 # and remain in effect until removed by pop().) 795 796 797 798""" 799 ) 800 801 outfile.write(" push(move(0,0,0))\n") 802 803 phi = 0.0 804 805 for i in range(0, N): 806 #im1 = i-1 807 # if im1 < 0 or self.settings.connect_ends: 808 # if im1 < 0: 809 # im1 += N 810 outfile.write(" pop()\n") 811 if len(orientations) > 0: 812 assert(len(orientations) == N) 813 if self.settings.orientations_use_quats: 814 assert(len(orientations[i]) == 4) 815 outfile.write(" push(quat(" + 816 str(orientations[i][0]) + "," + 817 str(orientations[i][1]) + "," + 818 str(orientations[i][2]) + "," + 819 str(orientations[i][3]) + "))\n") 820 else: 821 assert(len(orientations[i]) == 9) 822 outfile.write(" push(matrix(" + 823 str(orientations[i][0]) + "," + 824 str(orientations[i][1]) + "," + 825 str(orientations[i][2]) + "," + 826 str(orientations[i][3]) + "," + 827 str(orientations[i][4]) + "," + 828 str(orientations[i][5]) + "," + 829 str(orientations[i][6]) + "," + 830 str(orientations[i][7]) + "," + 831 str(orientations[i][8]) + "))\n") 832 else: 833 # Otherwise, if no orientations were explicitly specified, then 834 # infer the orientation from the direction of the displacement. 835 outfile.write(" push(rotvv(" + 836 str(self.direction_vects[i - 1][0]) + "," + 837 str(self.direction_vects[i - 1][1]) + "," + 838 str(self.direction_vects[i - 1][2]) + "," + 839 str(self.direction_vects[i][0]) + "," + 840 str(self.direction_vects[i][1]) + "," + 841 str(self.direction_vects[i][2]) + "))\n") 842 outfile.write(" push(move(" + 843 str(coords[i][0]) + "," + 844 str(coords[i][1]) + "," + 845 str(coords[i][2]) + "))\n") 846 847 outfile.write(" mon[" + str(i) + "] = new " + 848 names_monomers[i]) 849 850 # If requested, apply additional rotations about the polymer axis 851 phi = 0.0 852 if len(helix_angles) > 0: 853 phi += helix_angles[i] 854 elif self.settings.delta_phi != 0.0: 855 phi = i * self.settings.delta_phi 856 # Recall that self.direction_vects[-1] = 857 # self.settings.direction_orig (usually 1,0,0) 858 outfile.write(".rot(" + str(phi) + 859 "," + str(self.settings.direction_orig[0]) + 860 "," + str(self.settings.direction_orig[1]) + 861 "," + str(self.settings.direction_orig[2]) + 862 ")\n") 863 if len(orientations) > 0: 864 outfile.write(" pop()\n") 865 866 867 outfile.write(" pop()\n") 868 if len(orientations) == 0: 869 outfile.write(" pop()\n"*N) 870 871 assert(len(self.settings.bonds_name) == 872 len(self.settings.bonds_type) == 873 len(self.settings.bonds_atoms) == 874 len(self.settings.bonds_index_offsets)) 875 876 if len(self.settings.bonds_type) > 0: 877 outfile.write("\n" 878 "\n" 879 " write(\"Data Bonds\") {\n") 880 881 WrapPeriodic.bounds_err = False 882 for i in range(0, N): 883 test = False 884 for b in range(0, len(self.settings.bonds_type)): 885 I = i + self.settings.bonds_index_offsets[b][0] 886 J = i + self.settings.bonds_index_offsets[b][1] 887 I = WrapPeriodic.Wrap(I, N) 888 J = WrapPeriodic.Wrap(J, N) 889 if WrapPeriodic.bounds_err: 890 WrapPeriodic.bounds_err = False 891 if not self.settings.connect_ends: 892 continue 893 outfile.write( 894 " $bond:" + self.settings.bonds_name[b] + str(i + 1)) 895 outfile.write(" @bond:" + self.settings.bonds_type[b] + " $atom:mon[" + str(I) + "]/" + self.settings.bonds_atoms[ 896 b][0] + " $atom:mon[" + str(J) + "]/" + self.settings.bonds_atoms[b][1] + "\n") 897 if len(self.settings.bonds_type) > 0: 898 outfile.write(" } # write(\"Data Bonds\") {...\n\n\n") 899 900 901 assert(len(self.settings.bonds_name_notype) == 902 len(self.settings.bonds_atoms_notype) == 903 len(self.settings.bonds_index_offsets_notype)) 904 905 if len(self.settings.bonds_atoms_notype) > 0: 906 outfile.write("\n" 907 "\n" 908 " write(\"Data Bond List\") {\n") 909 910 WrapPeriodic.bounds_err = False 911 for i in range(0, N): 912 test = False 913 for b in range(0, len(self.settings.bonds_atoms_notype)): 914 I = i + self.settings.bonds_index_offsets_notype[b][0] 915 J = i + self.settings.bonds_index_offsets_notype[b][1] 916 I = WrapPeriodic.Wrap(I, N) 917 J = WrapPeriodic.Wrap(J, N) 918 if WrapPeriodic.bounds_err: 919 WrapPeriodic.bounds_err = False 920 if not self.settings.connect_ends: 921 continue 922 outfile.write( 923 " $bond:" + self.settings.bonds_name_notype[b] + str(i + 1)) 924 outfile.write(" $atom:mon[" + str(I) + "]/" + self.settings.bonds_atoms_notype[ 925 b][0] + " $atom:mon[" + str(J) + "]/" + self.settings.bonds_atoms_notype[b][1] + "\n") 926 if len(self.settings.bonds_atoms_notype) > 0: 927 outfile.write(" } # write(\"Data Bond List\") {...\n\n\n") 928 929 930 assert(len(self.settings.angles_name) == 931 len(self.settings.angles_type) == 932 len(self.settings.angles_atoms) == 933 len(self.settings.angles_index_offsets)) 934 if len(self.settings.angles_type) > 0: 935 outfile.write("\n" 936 "\n" 937 " write(\"Data Angles\") {\n") 938 for i in range(0, N): 939 for b in range(0, len(self.settings.angles_type)): 940 I = i + self.settings.angles_index_offsets[b][0] 941 J = i + self.settings.angles_index_offsets[b][1] 942 K = i + self.settings.angles_index_offsets[b][2] 943 I = WrapPeriodic.Wrap(I, N) 944 J = WrapPeriodic.Wrap(J, N) 945 K = WrapPeriodic.Wrap(K, N) 946 if WrapPeriodic.bounds_err: 947 WrapPeriodic.bounds_err = False 948 if not self.settings.connect_ends: 949 continue 950 outfile.write( 951 " $angle:" + self.settings.angles_name[b] + str(i + 1)) 952 outfile.write(" @angle:" + self.settings.angles_type[b] + 953 " $atom:mon[" + str(I) + "]/" + self.settings.angles_atoms[b][0] + 954 " $atom:mon[" + str(J) + "]/" + self.settings.angles_atoms[b][1] + 955 " $atom:mon[" + str(K) + "]/" + self.settings.angles_atoms[b][2] + 956 "\n") 957 if len(self.settings.angles_type) > 0: 958 outfile.write(" } # write(\"Data Angles\") {...\n\n\n") 959 960 961 assert(len(self.settings.dihedrals_name) == 962 len(self.settings.dihedrals_type) == 963 len(self.settings.dihedrals_atoms) == 964 len(self.settings.dihedrals_index_offsets)) 965 if len(self.settings.dihedrals_type) > 0: 966 outfile.write("\n" 967 "\n" 968 " write(\"Data Dihedrals\") {\n") 969 for i in range(0, N): 970 for b in range(0, len(self.settings.dihedrals_type)): 971 I = i + self.settings.dihedrals_index_offsets[b][0] 972 J = i + self.settings.dihedrals_index_offsets[b][1] 973 K = i + self.settings.dihedrals_index_offsets[b][2] 974 L = i + self.settings.dihedrals_index_offsets[b][3] 975 I = WrapPeriodic.Wrap(I, N) 976 J = WrapPeriodic.Wrap(J, N) 977 K = WrapPeriodic.Wrap(K, N) 978 L = WrapPeriodic.Wrap(L, N) 979 if WrapPeriodic.bounds_err: 980 WrapPeriodic.bounds_err = False 981 if not self.settings.connect_ends: 982 continue 983 outfile.write(" $dihedral:" + 984 self.settings.dihedrals_name[b] + str(i + 1)) 985 outfile.write(" @dihedral:" + self.settings.dihedrals_type[b] + 986 " $atom:mon[" + str(I) + "]/" + self.settings.dihedrals_atoms[b][0] + 987 " $atom:mon[" + str(J) + "]/" + self.settings.dihedrals_atoms[b][1] + 988 " $atom:mon[" + str(K) + "]/" + self.settings.dihedrals_atoms[b][2] + 989 " $atom:mon[" + str(L) + "]/" + self.settings.dihedrals_atoms[b][3] + 990 "\n") 991 if len(self.settings.dihedrals_type) > 0: 992 outfile.write(" } # write(\"Data Dihedrals\") {...\n\n\n") 993 994 995 assert(len(self.settings.impropers_name) == 996 len(self.settings.impropers_type) == 997 len(self.settings.impropers_atoms) == 998 len(self.settings.impropers_index_offsets)) 999 if len(self.settings.impropers_type) > 0: 1000 outfile.write("\n" 1001 "\n" 1002 " write(\"Data Impropers\") {\n") 1003 for i in range(0, N): 1004 for b in range(0, len(self.settings.impropers_type)): 1005 I = i + self.settings.impropers_index_offsets[b][0] 1006 J = i + self.settings.impropers_index_offsets[b][1] 1007 K = i + self.settings.impropers_index_offsets[b][2] 1008 L = i + self.settings.impropers_index_offsets[b][3] 1009 I = WrapPeriodic.Wrap(I, N) 1010 J = WrapPeriodic.Wrap(J, N) 1011 K = WrapPeriodic.Wrap(K, N) 1012 L = WrapPeriodic.Wrap(L, N) 1013 if WrapPeriodic.bounds_err: 1014 WrapPeriodic.bounds_err = False 1015 if not self.settings.connect_ends: 1016 continue 1017 outfile.write(" $improper:" + 1018 self.settings.impropers_name[b] + str(i + 1)) 1019 outfile.write(" @improper:" + self.settings.impropers_type[b] + 1020 " $atom:mon[" + str(I) + "]/" + self.settings.impropers_atoms[b][0] + 1021 " $atom:mon[" + str(J) + "]/" + self.settings.impropers_atoms[b][1] + 1022 " $atom:mon[" + str(K) + "]/" + self.settings.impropers_atoms[b][2] + 1023 " $atom:mon[" + str(L) + "]/" + self.settings.impropers_atoms[b][3] + 1024 "\n") 1025 if len(self.settings.impropers_type) > 0: 1026 outfile.write(" } # write(\"Data Impropers\") {...\n\n\n") 1027 1028 if name_polymer != '': 1029 outfile.write("} # " + name_polymer + "\n\n\n\n") 1030 1031 def CalcBoxBoundaries(self, coords): 1032 N = len(coords) 1033 for i in range(0, N): 1034 for d in range(0, 3): 1035 if not self.box_bounds_min: 1036 assert(not self.box_bounds_max) 1037 self.box_bounds_min = [xd for xd in coords[i]] 1038 self.box_bounds_max = [xd for xd in coords[i]] 1039 else: 1040 if coords[i][d] > self.box_bounds_max[d]: 1041 self.box_bounds_max[d] = coords[i][d] 1042 if coords[i][d] < self.box_bounds_min[d]: 1043 self.box_bounds_min[d] = coords[i][d] 1044 1045 def WriteBoxBoundaries(self, outfile): 1046 for d in range(0, 3): 1047 self.box_bounds_min[d] -= self.settings.box_padding[d] 1048 self.box_bounds_max[d] += self.settings.box_padding[d] 1049 outfile.write("\n# ---------------- simulation box -----------------\n" 1050 1051 "# Now define a box big enough to hold a polymer with this (initial) shape\n" 1052 "# (The user can override this later on. This is the default box size.)" 1053 "\n\n" 1054 "write_once(\"Data Boundary\") {\n" 1055 + str(self.box_bounds_min[0]) + " " + 1056 str(self.box_bounds_max[0]) + " xlo xhi\n" 1057 + str(self.box_bounds_min[1]) + " " + 1058 str(self.box_bounds_max[1]) + " ylo yhi\n" 1059 + str(self.box_bounds_min[2]) + " " + 1060 str(self.box_bounds_max[2]) + " zlo zhi\n" 1061 "}\n\n\n") 1062 1063 1064def main(): 1065 try: 1066 g_program_name = __file__.split('/')[-1] 1067 g_version_str = '0.1.10' 1068 g_date_str = '2021-12-28' 1069 sys.stderr.write(g_program_name + ' v' + 1070 g_version_str + ' ' + g_date_str + '\n') 1071 argv = [arg for arg in sys.argv] 1072 genpoly = GenPoly() 1073 genpoly.ParseArgs(argv) 1074 # Any remain arguments? 1075 if len(argv) > 1: 1076 raise InputError('Error(' + g_program_name + '):\n' + 1077 'Unrecogized command line argument \"' + argv[1] + 1078 '\"\n\n' + 1079 g_usage_msg) 1080 1081 if genpoly.settings.infile_name == '': 1082 # If the user did not specify the name of the file where the 1083 # coordinates are stored, then the user wants us to read 1084 # them from the standard (terminal) input (sys.stdin) 1085 genpoly.coords_multi = genpoly.ReadCoords(sys.stdin) 1086 1087 outfile = sys.stdout 1088 1089 # Read the coordinates 1090 1091 1092 if genpoly.name_sequence_multi == []: 1093 # In that case it means we should assume that each monomer in 1094 # the polymer(s) is of the same type. So we fill the 1095 # genpoly.name_sequence_multi array with the same entry. Recall that 1096 # variables ending in "_multi" are two-dimensional lists take two 1097 # indices [i][j] 1098 # The first index [i] specifies the polymer 1099 # The second index [j] specifies the monomer within that polymer. 1100 # I will use the following ugly list comprehension to create a 1101 # two-dimensional list-of-lists whose shape and size is consistent 1102 # with the size of the genpoly.coords_multi array. 1103 # And I will fill it with "genpoly.settings.name_monomer" 1104 genpoly.name_sequence_multi = [[genpoly.settings.name_monomer 1105 for j in 1106 range(0, len(genpoly.coords_multi[i]))] 1107 for i in range(0, len(genpoly.coords_multi))] 1108 1109 # Now, check for polymer and sequence length inconsistencies: 1110 if (len(genpoly.coords_multi) != len(genpoly.name_sequence_multi)): 1111 raise InputError('Error(' + 1112 g_program_name + '):\n' + 1113 ' The coordinate file and sequence file have different lengths.\n') 1114 if ((len(genpoly.orientations_multi) > 0) and 1115 (len(genpoly.orientations_multi) != len(genpoly.coords_multi))): 1116 raise InputError('Error(' + 1117 g_program_name + '):\n' + 1118 ' The coordinate file and orientations/quats file have different lengths.\n') 1119 if ((len(genpoly.helix_angles_multi) > 0) and 1120 (len(genpoly.helix_angles_multi) != len(genpoly.coords_multi))): 1121 raise InputError('Error(' + 1122 g_program_name + '):\n' + 1123 ' The coordinate file and helix_angles file have different lengths.\n') 1124 for i in range(0, len(genpoly.coords_multi)): 1125 if len(genpoly.name_sequence_multi[i]) > 0: 1126 if (len(genpoly.coords_multi[i]) != 1127 len(genpoly.name_sequence_multi[i])): 1128 raise InputError('Error(' + 1129 g_program_name + '):\n' + 1130 ' The coordinate file and sequence file have different lengths.\n') 1131 if ((len(genpoly.orientations_multi) > 0) and 1132 (len(genpoly.orientations_multi[i]) > 0)): 1133 if (len(genpoly.coords_multi[i]) != 1134 len(genpoly.orientations_multi[i])): 1135 raise InputError('Error(' + 1136 g_program_name + '):\n' + 1137 ' The coordinate file and orientations/quats file have different lengths.\n') 1138 if ((len(genpoly.helix_angles_multi) > 0) and 1139 (len(genpoly.helix_angles_multi[i]) > 0)): 1140 if (len(genpoly.coords_multi[i]) != 1141 len(genpoly.helix_angles_multi[i])): 1142 raise InputError('Error(' + 1143 g_program_name + '):\n' + 1144 ' The coordinate file and helix_angles file have different lengths.\n') 1145 1146 1147 1148 # Convert all of this information to moltemplate (LT) format: 1149 genpoly.WriteLTFile(outfile) 1150 1151 1152 except (ValueError, InputError) as err: 1153 sys.stderr.write('\n' + str(err) + '\n') 1154 sys.exit(-1) 1155 1156 return 1157 1158if __name__ == '__main__': 1159 main() 1160