From 225c92a717a6846c40144825de6fc358afa91754 Mon Sep 17 00:00:00 2001 From: Oliver Henrich Date: Tue, 1 Oct 2024 14:42:35 +0100 Subject: [PATCH 1/4] Restored PR#19 --- examples/PACKAGES/cgdna/util/lj2real.py | 79 ++++++++++++++++++------- 1 file changed, 59 insertions(+), 20 deletions(-) diff --git a/examples/PACKAGES/cgdna/util/lj2real.py b/examples/PACKAGES/cgdna/util/lj2real.py index 69076e9b010..24b3ff58fa3 100644 --- a/examples/PACKAGES/cgdna/util/lj2real.py +++ b/examples/PACKAGES/cgdna/util/lj2real.py @@ -143,10 +143,7 @@ def check_datafile_header(line: str, sections: Sections): Args: line (str): The line to check - masses_section (bool): If the current section is the masses section - atoms_section (bool): If the current section is the atoms section - velocities_section (bool): If the current section is the velocities section - ellipsoids_section (bool): If the current section is the ellipsoids section + sections (Sections): The Sections object to modify """ if any(header in line for header in ["xlo", "xhi", "ylo", "yhi", "zlo", "zhi"]): @@ -192,6 +189,7 @@ def modify_datafile(datafile_path: str, conversion_factors: ConversionFactors): Args: datafile_path (str): The path to the file to modify + conversion_factors (ConversionFactors): The conversion factors to use """ lines_changed = 0 current_section = Sections(False, False, False, False, False) @@ -290,6 +288,7 @@ def modify_inputfile(inputfile_path: str, conversion_factors: ConversionFactors) Args: inputfile_path (str): The path to the input file to modify + conversion_factors (ConversionFactors): The conversion factors to use """ lines_changed = 0 @@ -310,6 +309,10 @@ def modify_inputfile(inputfile_path: str, conversion_factors: ConversionFactors) "Warning: Both oxdna2 and oxrna2 found in input file. Output will likely be incorrect." ) + elements = line.split() + if not elements: + continue + if "variable T" in line: old_value = line.split()[3] @@ -327,7 +330,6 @@ def modify_inputfile(inputfile_path: str, conversion_factors: ConversionFactors) lines_changed += 1 elif "atom_modify" in line: - elements = line.split() elements[3] = str( round(float(elements[3]) * conversion_factors.length_conv_factor, 3) ) @@ -335,7 +337,6 @@ def modify_inputfile(inputfile_path: str, conversion_factors: ConversionFactors) lines_changed += 1 elif "neighbor" in line: - elements = line.split() elements[1] = str( round(float(elements[1]) * conversion_factors.length_conv_factor, 3) ) @@ -343,18 +344,15 @@ def modify_inputfile(inputfile_path: str, conversion_factors: ConversionFactors) lines_changed += 1 elif "read_data" in line: - elements = line.split() if conversion_factors.inverted: elements[1] = elements[1] + "_lj" else: - elements[1] = ( - elements[1] + "_real" - ) # naming convention of datafile after conversion + elements[1] = elements[1] + "_real" + # naming convention of datafile after conversion lines[i] = " ".join(elements) + "\n" lines_changed += 1 elif "mass" in line: - elements = line.split() elements[4] = str( round(float(elements[4]) * conversion_factors.mass_conv_factor, 4) ) @@ -362,16 +360,25 @@ def modify_inputfile(inputfile_path: str, conversion_factors: ConversionFactors) lines_changed += 1 elif "bond_coeff" in line or "pair_coeff" in line: - if ".lj" in line or ".real" in line: + if ".cgdna" in line: # potential files if conversion_factors.inverted: - line = line.replace(".real", ".lj") + line = line.replace("real.cgdna", "lj.cgdna") else: - line = line.replace(".lj", ".real") + line = line.replace("lj.cgdna", "real.cgdna") lines[i] = line lines_changed += 1 + # converting modifiable paramters outside of potential files if "stk" in line and "xstk" not in line and "coaxstk" not in line: elements = line.split() + if elements[5] != "${T}": + elements[5] = str( # convert T + round( + float(elements[5]) + * conversion_factors.temp_conv_factor, + 1, + ) + ) elements[6] = str( # convert xi round( float(elements[6]) * conversion_factors.energy_conv_factor, @@ -383,9 +390,19 @@ def modify_inputfile(inputfile_path: str, conversion_factors: ConversionFactors) ) lines[i] = " ".join(elements) + "\n" - else: - elements = line.split() + elif "dh" in line: + elements = line.split() + if elements[4] != "${T}": + elements[4] = str( # convert T + round( + float(elements[4]) + * conversion_factors.temp_conv_factor, + 1, + ) + ) + lines[i] = " ".join(elements) + "\n" + else: # non-potential files if "bond_coeff" in line: if oxdna2_flag: elements[2:] = conversion_factors.oxdna2_fene_string.split() @@ -427,6 +444,14 @@ def modify_inputfile(inputfile_path: str, conversion_factors: ConversionFactors) elements[4:] = conversion_factors.oxdna_xstk_string.split() else: # stk + if elements[5] != "${T}": + elements[5] = str( # convert T + round( + float(elements[5]) + * conversion_factors.temp_conv_factor, + 1, + ) + ) elements[6] = str( # convert xi round( float(elements[6]) @@ -472,11 +497,27 @@ def modify_inputfile(inputfile_path: str, conversion_factors: ConversionFactors) elements[5:] = ( conversion_factors.oxdna_hbond_1_4_2_3_string.split() ) + elif "dh" in line: + if elements[4] != "${T}": + elements[4] = str( + round( + float(elements[4]) + * conversion_factors.temp_conv_factor, + 1, + ) + ) lines[i] = " ".join(elements) + "\n" lines_changed += 1 - elif "langevin" in line: - elements = line.split() + elif "langevin" in line: # compatible with fix langevin + if elements[4] != "${T}": + elements[4] = str( + round(float(elements[4]) * conversion_factors.temp_conv_factor, 1) + ) + if elements[5] != "${T}": + elements[5] = str( + round(float(elements[5]) * conversion_factors.temp_conv_factor, 1) + ) elements[6] = str( round(float(elements[6]) * conversion_factors.time_conv_factor, 2) ) @@ -484,7 +525,6 @@ def modify_inputfile(inputfile_path: str, conversion_factors: ConversionFactors) lines_changed += 1 elif "timestep" in line: - elements = line.split() elements[1] = str( round(float(elements[1]) * conversion_factors.time_conv_factor, 5) ) @@ -492,7 +532,6 @@ def modify_inputfile(inputfile_path: str, conversion_factors: ConversionFactors) lines_changed += 1 elif "comm_modify" in line: - elements = line.split() elements[2] = str( round(float(elements[2]) * conversion_factors.length_conv_factor, 1) ) From 013910ce59e98ee174f7860bb4f94bc8a74e58be Mon Sep 17 00:00:00 2001 From: Oliver Henrich Date: Tue, 1 Oct 2024 14:43:42 +0100 Subject: [PATCH 2/4] Restored ring generation functionality, debugged unique base pairing --- examples/PACKAGES/cgdna/util/generate.py | 322 ++++++++++++++++++++--- 1 file changed, 285 insertions(+), 37 deletions(-) diff --git a/examples/PACKAGES/cgdna/util/generate.py b/examples/PACKAGES/cgdna/util/generate.py index 0d8630c88b1..2dfc4b6bc13 100644 --- a/examples/PACKAGES/cgdna/util/generate.py +++ b/examples/PACKAGES/cgdna/util/generate.py @@ -22,13 +22,17 @@ For dsDNA the sequence should be preceded by the 'DOUBLE' keyword. Usage: -$$ python generate.py box_offset box_length sequence_file +$$ python generate.py box_offset box_length sequence_file [topology (strand OR ring)] [twist ((p-p_0)/p_0) for strand OR Delta_Lk for ring] """ # for python2/3 compatibility from __future__ import print_function #!/usr/bin/env python +# number of ACGT base type groups +# default value 16 protects against asymmetric base pairing over 1.5 turns on the helix +N_BASE_TYPE_GROUPS = 16 + """ Import basic modules """ @@ -54,9 +58,21 @@ box_offset = float(sys.argv[1]) box_length = float(sys.argv[2]) infile = sys.argv[3] + # default is strand with no twist, i.e. pitch = 10.5 bp + if len(sys.argv) == 4: + topo = 'strand' + twist = 0 + # default is no twist, i.e. pitch = 10.5 bp + elif len(sys.argv) == 5: + topo = sys.argv[4] + twist = 0 + elif len(sys.argv) == 6: + topo = sys.argv[4] + twist = float(sys.argv[5]) + except: - print( "Usage: %s <%s> <%s> <%s>" % (sys.argv[0], \ - "box offset", "box length", "file with sequences"), file=sys.stderr) + print("Usage: %s <%s> <%s> <%s> <%s> <%s> " % (sys.argv[0], \ + "box offset", "box length", "file with sequences", "[topology (strand OR ring)]", "[twist ((p-p_0)/p_0) for strand OR Delta_Lk for ring]"), file=sys.stderr) sys.exit(1) box = np.array ([box_length, box_length, box_length]) @@ -253,11 +269,48 @@ def add_strands (mynewpositions, mynewa1s, mynewa3s): return True +""" +Calculate angle of rotation site to site +""" +def get_angle(bp): + #n, default number of bases per turn + n = 10.55 + np = n + nm = n + found = False + angle = 0 + while found == False: + + turnsp = bp/np + turnsm = bp/nm + + diffp = abs(turnsp - round(turnsp)) + diffm = abs(turnsm - round(turnsm)) + + if diffp < 0.03: + found = True + turnsp = round(turnsp)+int(twist) + angle = (360*turnsp)/bp + angle = round(angle,2) + break + + elif diffm < 0.03: + found = True + turnsm = round(turnsm)+int(twist) + angle = (360*turnsm)/bp + angle = round(angle,2) + break + + else: + np += 0.001 + nm -= 0.001 + + return angle """ Returns the rotation matrix defined by an axis and angle """ -def get_rotation_matrix(axis, anglest): +def get_rotation_matrix(axis, anglest, nbp=0): # The argument anglest can be either an angle in radiants # (accepted types are float, int or np.float64 or np.float64) # or a tuple [angle, units] where angle is a number and @@ -266,10 +319,19 @@ def get_rotation_matrix(axis, anglest): if not isinstance (anglest, (np.float64, np.float32, float, int)): if len(anglest) > 1: if anglest[1] in ["degrees", "deg", "o"]: - #angle = np.deg2rad (anglest[0]) angle = (np.pi / 180.) * (anglest[0]) elif anglest[1] in ["bp"]: - angle = int(anglest[0]) * (np.pi / 180.) * (35.9) + + # strand: angle is 360/pitch + # pitch is p = 10.5 * (1-twist) + if nbp == 0: + ang = 360./(10.5*(1-twist)) + + # ring: angle depends on Delta_Lk = Lk - Lk_0 + else: + ang = get_angle(nbp) + + angle = int(anglest[0]) * (np.pi / 180.) * (ang) else: angle = float(anglest[0]) else: @@ -323,6 +385,8 @@ def generate_strand(bp, sequence=None, start_pos=np.array([0, 0, 0]), \ # if not provided switch off random orientation if perp is None or perp is False: v1 = np.random.random_sample(3) + # comment in to suppress randomised base vector + v1 = [1,0,0] v1 -= dir * (np.dot(dir, v1)) v1 /= np.sqrt(sum(v1*v1)) else: @@ -347,10 +411,10 @@ def generate_strand(bp, sequence=None, start_pos=np.array([0, 0, 0]), \ a3 = dir for i in range(bp): # work out the position of the centre of mass of the nucleotide - rcdm = rb - CM_CENTER_DS * a1 + rcom = rb - CM_CENTER_DS * a1 # append to newpositions - mynewpositions.append(rcdm) + mynewpositions.append(rcom) mynewa1s.append(a1) mynewa3s.append(a3) @@ -367,8 +431,8 @@ def generate_strand(bp, sequence=None, start_pos=np.array([0, 0, 0]), \ a3 = -dir R = R.transpose() for i in range(bp): - rcdm = rb - CM_CENTER_DS * a1 - mynewpositions.append (rcdm) + rcom = rb - CM_CENTER_DS * a1 + mynewpositions.append (rcom) mynewa1s.append (a1) mynewa3s.append (a3) a1 = np.dot(R, a1) @@ -378,6 +442,125 @@ def generate_strand(bp, sequence=None, start_pos=np.array([0, 0, 0]), \ return [mynewpositions, mynewa1s, mynewa3s] +""" +Generates the position and orientation vectors of a +(single or double) ring from a sequence string +""" +def generate_ring(bp, sequence=None, start_pos=np.array([0, 0, 0]), \ + dir=np.array([0, 0, 1]), perp=False, double=True, rot=0.): + # generate empty arrays + mynewpositions, mynewa1s, mynewa3s = [], [], [] + + # cast the provided start_pos array into a numpy array + start_pos = np.array(start_pos, dtype=float) + + # overall direction of the helix + dir = np.array(dir, dtype=float) + if sequence == None: + sequence = np.random.randint(1, 5, bp) + + # the elseif here is most likely redundant + elif len(sequence) != bp: + n = bp - len(sequence) + sequence += np.random.randint(1, 5, n) + print("sequence is too short, adding %d random bases" % n, file=sys.stderr) + + # normalize direction + dir_norm = np.sqrt(np.dot(dir,dir)) + if dir_norm < 1e-10: + print("direction must be a valid vector,\ + defaulting to (0, 0, 1)", file=sys.stderr) + dir = np.array([0, 0, 1]) + else: dir /= dir_norm + + # find a vector orthogonal to dir to act as helix direction, + # if not provided switch off random orientation + if perp is None or perp is False: + v1 = np.random.random_sample(3) + # comment in to suppress randomised base vector + v1 = [1,0,0] + v1 -= dir * (np.dot(dir, v1)) + v1 /= np.sqrt(sum(v1*v1)) + else: + v1 = perp; + + + #work out v2 (lab frame) (Orthogonal vector to dir and v1) - probably need to translate by R + v2 = np.cross(dir,v1) + v2 /= np.sqrt(np.dot(v2,v2)) + + # generate rotational matrix representing the overall rotation of the helix + R0 = get_rotation_matrix(dir, rot) + + # rotation matrix corresponding to one step along the helix + R = get_rotation_matrix(dir, [1, "bp"],bp) + + #rotation matrix for ring creation, rotation around lab v2 by 2 pi/ len(bp) + + ang =( 2*PI)/(bp) + + Rv2 = get_rotation_matrix(v2,ang) + + # set the vector a1 (backbone to base) to v1 + a1 = v1 + + # apply the global rotation to a1 + a1 = np.dot(R0, a1) + + # set the position of the fist backbone site to start_pos + rb = np.array(start_pos) + + # set a3 to the direction of the helix + a3 = dir + for i in range(bp): + + R = get_rotation_matrix(a3,[1,"bp"],bp) + + rcom = rb - CM_CENTER_DS * a1 + + # append to newpositions + mynewpositions.append(rcom) + mynewa1s.append(a1) + mynewa3s.append(a3) + if i != bp -1 : + rb += a3*BASE_BASE + #calculate a3p, a1p by rotating a3, a1 around v2 + a1 = np.dot(R0,a1) + a3p = np.dot(Rv2,a3) + a1p = np.dot(Rv2,a1) + a1pp = np.dot(R,a1p) + a3 = a3p + a1 = a1pp + + if double == True: + cor=(bp+1)/bp + cor=1 + ang =(-2*PI)/(bp)*cor + + a1 = -a1 + a3 = -a3 + R = R.transpose() + Rv2 = get_rotation_matrix(v2,ang) + + for i in range(bp): + + rcom = rb - CM_CENTER_DS * a1 + + mynewpositions.append (rcom) + mynewa1s.append (a1) + mynewa3s.append(a3) + + a3p = np.dot(Rv2,a3) + a1p = np.dot(Rv2,a1) + a1pp = np.dot(R,a1p) + a3 = a3p + a1 = a1pp + rb += a3 * BASE_BASE + R = get_rotation_matrix(a3,[1,"bp"],bp) + + assert (len (mynewpositions) > 0) + + return [mynewpositions, mynewa1s, mynewa3s] """ Main function for this script. @@ -417,14 +600,21 @@ def read_strands(filename): print( "## Found duplex of %i base pairs" % length, file=sys.stdout) nnucl += 2*length nstrands += 2 - nbonds += (2*length-2) + if topo == 'ring': + nbonds+= 2*length + else: + nbonds += (2*length-2) else: line = line.split()[0] length = len(line) print( "## Found single strand of %i bases" % length, file=sys.stdout) nnucl += length nstrands += 1 - nbonds += length-1 + if topo == 'ring': + nbonds =+ length + else: + nbonds += length-1 + # rewind the sequence input file infile.seek(0) @@ -458,7 +648,17 @@ def read_strands(filename): length = len(line) seq = [(base_to_number[x]) for x in line] + # modify default type for unique base pairing + group = 0 + for s in range(len(seq)): + seq[s] += group*4 + group += 1 + # reset according to maximal number of base type groups + if group == N_BASE_TYPE_GROUPS: + group = 0 + myns += 1 + for b in range(length): basetype.append(seq[b]) strandnum.append(myns) @@ -466,14 +666,30 @@ def read_strands(filename): for b in range(length-1): bondpair = [noffset + b, noffset + b + 1] bonds.append(bondpair) - noffset += length - # create the sequence of the second strand as made of - # complementary bases - seq2 = [5-s for s in seq] - seq2.reverse() + if topo== 'ring': + bondpair = [noffset,noffset+length-1] + bonds.append(bondpair) + noffset += length + + # create the complementary sequence of the second strand + + seq2 = seq + for s in range(len(seq2)): + if seq2[s]%4 == 1: + seq2[s] += 3 + elif seq2[s]%4 == 2: + seq2[s] += 1 + elif seq2[s]%4 == 3: + seq2[s] -= 1 + elif seq2[s]%4 == 0: + seq2[s] -= 3 + + seq2 = seq2[::-1] + myns += 1 + for b in range(length): basetype.append(seq2[b]) strandnum.append(myns) @@ -481,31 +697,48 @@ def read_strands(filename): for b in range(length-1): bondpair = [noffset + b, noffset + b + 1] bonds.append(bondpair) + + if topo == 'ring': + bondpair = [noffset, noffset + length-1] + bonds.append(bondpair) + noffset += length print( "## Created duplex of %i bases" % (2*length), file=sys.stdout) # generate random position of the first nucleotide - cdm = box_offset + np.random.random_sample(3) * box + com = box_offset + np.random.random_sample(3) * box + # comment out to randomise + com = [0,0,0] # generate the random direction of the helix axis = np.random.random_sample(3) + # comment out to randomise + axis = [0,0,1] axis /= np.sqrt(np.dot(axis, axis)) - # use the generate function defined above to create + # use the generate function defined above to create # the position and orientation vector of the strand - newpositions, newa1s, newa3s = generate_strand(len(line), \ - sequence=seq, dir=axis, start_pos=cdm, double=True) + if topo == 'ring': + newpositions, newa1s, newa3s = generate_ring(len(line), \ + sequence=seq, dir=axis, start_pos=com, double=True) + else: + newpositions, newa1s, newa3s = generate_strand(len(line), \ + sequence=seq, dir=axis, start_pos=com, double=True) # generate a new position for the strand until it does not overlap # with anything already present start = timer() while not add_strands(newpositions, newa1s, newa3s): - cdm = box_offset + np.random.random_sample(3) * box + com = box_offset + np.random.random_sample(3) * box axis = np.random.random_sample(3) axis /= np.sqrt(np.dot(axis, axis)) - newpositions, newa1s, newa3s = generate_strand(len(line), \ - sequence=seq, dir=axis, start_pos=cdm, double=True) + if topo == 'ring': + newpositions, newa1s, newa3s = generate_ring(len(line), \ + sequence=seq, dir=axis, start_pos=com, double=True) + else: + newpositions, newa1s, newa3s = generate_strand(len(line), \ + sequence=seq, dir=axis, start_pos=com, double=True) print( "## Trying %i" % i, file=sys.stdout) end = timer() print( "## Added duplex of %i bases (line %i/%i) in %.2fs, now at %i/%i" % \ @@ -525,26 +758,43 @@ def read_strands(filename): for b in range(length-1): bondpair = [noffset + b, noffset + b + 1] bonds.append(bondpair) + + if topo == 'ring': + bondpair = [noffset, noffset + length-1] + bonds.append(bondpair) + noffset += length # generate random position of the first nucleotide - cdm = box_offset + np.random.random_sample(3) * box + com = box_offset + np.random.random_sample(3) * box # generate the random direction of the helix axis = np.random.random_sample(3) axis /= np.sqrt(np.dot(axis, axis)) print("## Created single strand of %i bases" % length, file=sys.stdout) - - newpositions, newa1s, newa3s = generate_strand(length, \ - sequence=seq, dir=axis, start_pos=cdm, double=False) + if topo == 'ring': + newpositions, newa1s, newa3s = generate_ring(length, \ + sequence=seq, dir=axis, start_pos=com, double=False) + else: + newpositions, newa1s, newa3s = generate_strand(length, \ + sequence=seq, dir=axis, start_pos=com, double=False) start = timer() while not add_strands(newpositions, newa1s, newa3s): - cdm = box_offset + np.random.random_sample(3) * box + com = box_offset + np.random.random_sample(3) * box + # comment out to randomise + com = [0,0,0] axis = np.random.random_sample(3) + # comment out to randomise + axis = [0,0,1] axis /= np.sqrt(np.dot(axis, axis)) - newpositions, newa1s, newa3s = generate_strand(length, \ - sequence=seq, dir=axis, start_pos=cdm, double=False) + if topo == 'ring': + newpositions, newa1s, newa3s = generate_ring(length, \ + sequence=seq, dir=axis, start_pos=com, double=False) + + else: + newpositions, newa1s, newa3s = generate_strand(length, \ + sequence=seq, dir=axis, start_pos=com, double=False) print >> sys.stdout, "## Trying %i" % (i) end = timer() print( "## Added single strand of %i bases (line %i/%i) in %.2fs, now at %i/%i" % \ @@ -562,7 +812,7 @@ def read_strands(filename): out.write('%d ellipsoids\n' % nnucl) out.write('%d bonds\n' % nbonds) out.write('\n') - out.write('4 atom types\n') + out.write('%d atom types\n' % (N_BASE_TYPE_GROUPS*4)) out.write('1 bond types\n') out.write('\n') out.write('# System size\n') @@ -573,10 +823,8 @@ def read_strands(filename): out.write('\n') out.write('Masses\n') out.write('\n') - out.write('1 3.1575\n') - out.write('2 3.1575\n') - out.write('3 3.1575\n') - out.write('4 3.1575\n') + for i in range(N_BASE_TYPE_GROUPS*4): + out.write('%d 3.1575\n' % (i+1)) # for each nucleotide print a line under the headers # Atoms, Velocities, Ellipsoids and Bonds @@ -591,7 +839,7 @@ def read_strands(filename): % (i+1, basetype[i], positions[i][0], positions[i][1], positions[i][2], strandnum[i])) out.write('\n') - out.write('# Atom-ID, translational, rotational velocity\n') + out.write('# Atom-ID, translational velocity, angular momentum\n') out.write('Velocities\n') out.write('\n') From e52ffcbae0c900ef8a08222e8461dd34a7a9a2ca Mon Sep 17 00:00:00 2001 From: Oliver Henrich Date: Tue, 1 Oct 2024 14:45:46 +0100 Subject: [PATCH 3/4] Set default number of base type groups to 1 --- examples/PACKAGES/cgdna/util/generate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/PACKAGES/cgdna/util/generate.py b/examples/PACKAGES/cgdna/util/generate.py index 2dfc4b6bc13..3431db41dbd 100644 --- a/examples/PACKAGES/cgdna/util/generate.py +++ b/examples/PACKAGES/cgdna/util/generate.py @@ -31,7 +31,7 @@ # number of ACGT base type groups # default value 16 protects against asymmetric base pairing over 1.5 turns on the helix -N_BASE_TYPE_GROUPS = 16 +N_BASE_TYPE_GROUPS = 1 """ Import basic modules From 48cb166f76b62c58164da9c5787d308299e961c1 Mon Sep 17 00:00:00 2001 From: Oliver Henrich Date: Fri, 1 Nov 2024 15:10:54 +0000 Subject: [PATCH 4/4] Restored ring functionality, corrected angle twist formula --- examples/PACKAGES/cgdna/util/generate.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/examples/PACKAGES/cgdna/util/generate.py b/examples/PACKAGES/cgdna/util/generate.py index 3431db41dbd..ba9b65becc1 100644 --- a/examples/PACKAGES/cgdna/util/generate.py +++ b/examples/PACKAGES/cgdna/util/generate.py @@ -22,7 +22,7 @@ For dsDNA the sequence should be preceded by the 'DOUBLE' keyword. Usage: -$$ python generate.py box_offset box_length sequence_file [topology (strand OR ring)] [twist ((p-p_0)/p_0) for strand OR Delta_Lk for ring] +$$ python generate.py box_offset box_length sequence_file [topology (strand OR ring)] [twist ((Lk-Lk_0)/Lk_0) for strand OR Delta_Lk=(Lk-Lk_0) for ring] """ # for python2/3 compatibility @@ -58,11 +58,11 @@ box_offset = float(sys.argv[1]) box_length = float(sys.argv[2]) infile = sys.argv[3] - # default is strand with no twist, i.e. pitch = 10.5 bp + # default is strand with no twist, i.e. pitch = 10.55 bp if len(sys.argv) == 4: topo = 'strand' twist = 0 - # default is no twist, i.e. pitch = 10.5 bp + # default is no twist, i.e. pitch = 10.55 bp elif len(sys.argv) == 5: topo = sys.argv[4] twist = 0 @@ -72,7 +72,7 @@ except: print("Usage: %s <%s> <%s> <%s> <%s> <%s> " % (sys.argv[0], \ - "box offset", "box length", "file with sequences", "[topology (strand OR ring)]", "[twist ((p-p_0)/p_0) for strand OR Delta_Lk for ring]"), file=sys.stderr) + "box offset", "box length", "file with sequences", "[topology (strand OR ring)]", "[twist ((Lk-Lk_0)/Lk_0) for strand OR Delta_Lk=(Lk-Lk_0) for ring]"), file=sys.stderr) sys.exit(1) box = np.array ([box_length, box_length, box_length]) @@ -322,10 +322,10 @@ def get_rotation_matrix(axis, anglest, nbp=0): angle = (np.pi / 180.) * (anglest[0]) elif anglest[1] in ["bp"]: - # strand: angle is 360/pitch - # pitch is p = 10.5 * (1-twist) + # strand: angle is 360 devided by default no. bpspitch + # positive/negative supercoiling increases/decreases the angle if nbp == 0: - ang = 360./(10.5*(1-twist)) + ang = 360./10.55*(1+twist) # ring: angle depends on Delta_Lk = Lk - Lk_0 else: