From 7f9e6fb077b7dae289a46ab4a5b014a32f125f95 Mon Sep 17 00:00:00 2001 From: Max Imakaev Date: Fri, 6 Jan 2023 21:45:10 -0500 Subject: [PATCH] black with 120 character limit; pyproject.toml --- examples/customIntegrators/activeBD.py | 20 +-- examples/customIntegrators/corr_noise.py | 8 +- .../directionalModel/flagshipNormLifetime.py | 12 +- .../directionalModel/showContactmap.py | 12 +- examples/loopExtrusion/newSimulatorCode.py | 58 ++------ polychrom/cli/traj_convert.py | 33 ++--- polychrom/contactmaps.py | 56 ++------ polychrom/contrib/integrators.py | 12 +- polychrom/forcekits.py | 18 +-- polychrom/forces.py | 100 ++++---------- polychrom/hdf5_format.py | 31 ++--- polychrom/legacy/contactmaps.py | 15 +- polychrom/legacy/forces.py | 112 ++++----------- polychrom/lib/extrusion.py | 31 +---- polychrom/polymer_analyses.py | 13 +- polychrom/polymerutils.py | 16 +-- polychrom/pymol_show.py | 42 ++---- polychrom/simulation.py | 130 ++++-------------- polychrom/starting_conformations.py | 5 +- pyproject.toml | 5 + tests/test_configurations.py | 26 +--- tests/test_contactmaps.py | 4 +- tests/test_forces.py | 4 +- tests/test_polymer_analyses.py | 20 +-- utilities/quickMovie.py | 5 +- utilities/showChainWithRasmol.py | 15 +- 26 files changed, 196 insertions(+), 607 deletions(-) create mode 100644 pyproject.toml diff --git a/examples/customIntegrators/activeBD.py b/examples/customIntegrators/activeBD.py index ad0dec7..bff2412 100644 --- a/examples/customIntegrators/activeBD.py +++ b/examples/customIntegrators/activeBD.py @@ -27,9 +27,7 @@ ids[int(N / 2) :] = 0 -def run_monomer_diffusion( - gpuid, N, ids, activity_ratio, timestep=170, nblocks=10, blocksize=100 -): +def run_monomer_diffusion(gpuid, N, ids, activity_ratio, timestep=170, nblocks=10, blocksize=100): """Run a single simulation on a GPU of a hetero-polymer with A monomers and B monomers. A monomers have a larger diffusion coefficient than B monomers, with an activity ratio of D_A / D_B. @@ -53,12 +51,8 @@ def run_monomer_diffusion( """ if len(ids) != N: - raise ValueError( - "The array of monomer identities must have length equal to the total number of monomers." - ) - D = np.ones( - (N, 3) - ) # Dx, Dy, Dz --> we assume the diffusion coefficient in each spatial dimension is the same + raise ValueError("The array of monomer identities must have length equal to the total number of monomers.") + D = np.ones((N, 3)) # Dx, Dy, Dz --> we assume the diffusion coefficient in each spatial dimension is the same # by default set the average diffusion coefficient to be 1 kT/friction # let D_A = 1 + Ddiff and D_B = 1 - Ddiff such that D_A / D_B is the given activity_ratio Ddiff = (activity_ratio - 1) / (activity_ratio + 1) @@ -79,9 +73,7 @@ def run_monomer_diffusion( particleD = unit.Quantity(D, kT / friction) integrator = ActiveBrownianIntegrator(timestep, collision_rate, particleD) gpuid = f"{gpuid}" - reporter = HDF5Reporter( - folder="active_inactive", max_data_length=100, overwrite=True - ) + reporter = HDF5Reporter(folder="active_inactive", max_data_length=100, overwrite=True) sim = simulation.Simulation( platform="CUDA", # for custom integrators, feed a tuple with the integrator class reference and a string specifying type, @@ -99,9 +91,7 @@ def run_monomer_diffusion( polymer = starting_conformations.grow_cubic(N, int(np.ceil(r))) sim.set_data(polymer, center=True) # loads a polymer, puts a center of mass at zero - sim.set_velocities( - v=np.zeros((N, 3)) - ) # initializes velocities of all monomers to zero (no inertia) + sim.set_velocities(v=np.zeros((N, 3))) # initializes velocities of all monomers to zero (no inertia) sim.add_force(forces.spherical_confinement(sim, density=density, k=5.0)) sim.add_force( forcekits.polymer_chains( diff --git a/examples/customIntegrators/corr_noise.py b/examples/customIntegrators/corr_noise.py index 2dfa420..3a2203f 100644 --- a/examples/customIntegrators/corr_noise.py +++ b/examples/customIntegrators/corr_noise.py @@ -88,16 +88,12 @@ def run_correlated_diffusion(gpuid, N, rhos, timestep=170, nblocks=10, blocksize """ if rhos.shape[1] != N: - raise ValueError( - "The array of monomer identities must have length equal to the total number of monomers." - ) + raise ValueError("The array of monomer identities must have length equal to the total number of monomers.") # monomer density in confinement in units of monomers/volume density = 0.224 r = (3 * N / (4 * 3.141592 * density)) ** (1 / 3) print(f"Radius of confinement: {r}") - D = np.ones( - (N, 3) - ) # Diffusion coefficients of N monomers in x,y,z spatial dimensions + D = np.ones((N, 3)) # Diffusion coefficients of N monomers in x,y,z spatial dimensions timestep = timestep # the monomer diffusion coefficient should be in units of kT / friction, where friction = mass*collision_rate collision_rate = 2.0 diff --git a/examples/loopExtrusion/directionalModel/flagshipNormLifetime.py b/examples/loopExtrusion/directionalModel/flagshipNormLifetime.py index 6c790ad..cb753c4 100644 --- a/examples/loopExtrusion/directionalModel/flagshipNormLifetime.py +++ b/examples/loopExtrusion/directionalModel/flagshipNormLifetime.py @@ -70,9 +70,7 @@ def averageContacts(contactFunction, inValues, N, **kwargs): filenameChunks = [filenames[i::nproc] for i in range(nproc)] - with closing( - mp.Pool(processes=nproc, initializer=init, initargs=sharedARrays_) - ) as p: + with closing(mp.Pool(processes=nproc, initializer=init, initargs=sharedARrays_)) as p: p.map(worker, filenameChunks) @@ -201,9 +199,7 @@ def doSim(i): hicdata = np.clip(hicdata, 0, np.percentile(hicdata, 99.99)) hicdata /= np.mean(np.sum(hicdata, axis=1)) - fmap( - doSim, range(30), n=20 - ) # number of threads to use. On a 20-core machine I use 20. + fmap(doSim, range(30), n=20) # number of threads to use. On a 20-core machine I use 20. arr = coarsegrain(arr, 20) arr = np.clip(arr, 0, np.percentile(arr, 99.9)) @@ -322,7 +318,5 @@ def doPolymerSimulation(steps, dens, stiff, folder): steps=5000, stiff=2, dens=0.2, - folder="flagshipLifetime{1}Mu{2}_try={0}".format( - sys.argv[1], sys.argv[2], sys.argv[3] - ), + folder="flagshipLifetime{1}Mu{2}_try={0}".format(sys.argv[1], sys.argv[2], sys.argv[3]), ) diff --git a/examples/loopExtrusion/directionalModel/showContactmap.py b/examples/loopExtrusion/directionalModel/showContactmap.py index 12e55c8..d8529ea 100644 --- a/examples/loopExtrusion/directionalModel/showContactmap.py +++ b/examples/loopExtrusion/directionalModel/showContactmap.py @@ -365,16 +365,8 @@ def showCmapNew(): if end > 75000: continue plt.xlim([st, end]) - plt.savefig( - os.path.join( - "heatmaps", "{0}_st={1}_end={2}_r=2.png".format(fname, st, end) - ) - ) - plt.savefig( - os.path.join( - "heatmaps", "{0}_st={1}_end={2}_r=2.pdf".format(fname, st, end) - ) - ) + plt.savefig(os.path.join("heatmaps", "{0}_st={1}_end={2}_r=2.png".format(fname, st, end))) + plt.savefig(os.path.join("heatmaps", "{0}_st={1}_end={2}_r=2.pdf".format(fname, st, end))) plt.clf() plt.show() diff --git a/examples/loopExtrusion/newSimulatorCode.py b/examples/loopExtrusion/newSimulatorCode.py index d9172a5..ee3809c 100644 --- a/examples/loopExtrusion/newSimulatorCode.py +++ b/examples/loopExtrusion/newSimulatorCode.py @@ -31,9 +31,7 @@ # new parameters because some things changed saveEveryBlocks = 10 # save every 10 blocks (saving every block is now too much almost) -skipSavedBlocksBeginning = ( - 20 # how many blocks (saved) to skip after you restart LEF positions -) +skipSavedBlocksBeginning = 20 # how many blocks (saved) to skip after you restart LEF positions totalSavedBlocks = 40 # how many blocks to save (number of blocks done is totalSavedBlocks * saveEveryBlocks) restartMilkerEveryBlocks = 100 @@ -49,16 +47,8 @@ savesPerMilker = restartMilkerEveryBlocks // saveEveryBlocks milkerInitsSkip = saveEveryBlocks * skipSavedBlocksBeginning // restartMilkerEveryBlocks -milkerInitsTotal = ( - (totalSavedBlocks + skipSavedBlocksBeginning) - * saveEveryBlocks - // restartMilkerEveryBlocks -) -print( - "Milker will be initialized {0} times, first {1} will be skipped".format( - milkerInitsTotal, milkerInitsSkip - ) -) +milkerInitsTotal = (totalSavedBlocks + skipSavedBlocksBeginning) * saveEveryBlocks // restartMilkerEveryBlocks +print("Milker will be initialized {0} times, first {1} will be skipped".format(milkerInitsTotal, milkerInitsSkip)) class smcTranslocatorMilker(object): @@ -93,9 +83,7 @@ def setup(self, bondForce, blocks=100, smcStepsPerBlock=1): """ if len(self.allBonds) != 0: - raise ValueError( - "Not all bonds were used; {0} sets left".format(len(self.allBonds)) - ) + raise ValueError("Not all bonds were used; {0} sets left".format(len(self.allBonds))) self.bondForce = bondForce @@ -115,11 +103,7 @@ def setup(self, bondForce, blocks=100, smcStepsPerBlock=1): self.curBonds = allBonds.pop(0) for bond in self.uniqueBonds: - paramset = ( - self.activeParamDict - if (bond in self.curBonds) - else self.inactiveParamDict - ) + paramset = self.activeParamDict if (bond in self.curBonds) else self.inactiveParamDict ind = bondForce.addBond(bond[0], bond[1], **paramset) self.bondInds.append(ind) self.bondToInd = {i: j for i, j in zip(self.uniqueBonds, self.bondInds)} @@ -133,9 +117,7 @@ def step(self, context, verbose=False): :return: (current bonds, previous step bonds); just for reference """ if len(self.allBonds) == 0: - raise ValueError( - "No bonds left to run; you should restart simulation and run setup again" - ) + raise ValueError("No bonds left to run; you should restart simulation and run setup again") pastBonds = self.curBonds self.curBonds = self.allBonds.pop(0) # getting current bonds @@ -153,12 +135,8 @@ def step(self, context, verbose=False): for bond, isAdd in zip(bondsToChange, bondsIsAdd): ind = self.bondToInd[bond] paramset = self.activeParamDict if isAdd else self.inactiveParamDict - self.bondForce.setBondParameters( - ind, bond[0], bond[1], **paramset - ) # actually updating bonds - self.bondForce.updateParametersInContext( - context - ) # now run this to update things in the context + self.bondForce.setBondParameters(ind, bond[0], bond[1], **paramset) # actually updating bonds + self.bondForce.updateParametersInContext(context) # now run this to update things in the context return self.curBonds, pastBonds @@ -195,9 +173,7 @@ def initModel(): # now polymer simulation code starts # ------------feed smcTran to the milker--- -SMCTran.steps( - 1000000 -) # first steps to "equilibrate" SMC dynamics. If desired of course. +SMCTran.steps(1000000) # first steps to "equilibrate" SMC dynamics. If desired of course. milker = smcTranslocatorMilker(SMCTran) # now feed this thing to milker (do it once!) # --------- end new code ------------ @@ -206,9 +182,7 @@ def initModel(): # simulation parameters are defined below a = Simulation(timestep=80, thermostat=0.01) - a.setup( - platform="CUDA", PBC=True, PBCbox=[box, box, box], GPU=0, precision="mixed" - ) # set up GPU here + a.setup(platform="CUDA", PBC=True, PBCbox=[box, box, box], GPU=0, precision="mixed") # set up GPU here a.saveFolder(folder) a.load(data) a.addHarmonicPolymerBonds(wiggleDist=0.1) @@ -234,13 +208,9 @@ def initModel(): ) # now only one step of SMC per step print("Restarting milker") - a.doBlock( - steps=steps, increment=False - ) # do block for the first time with first set of bonds in + a.doBlock(steps=steps, increment=False) # do block for the first time with first set of bonds in for i in range(restartMilkerEveryBlocks - 1): - curBonds, pastBonds = milker.step( - a.context - ) # this updates bonds. You can do something with bonds here + curBonds, pastBonds = milker.step(a.context) # this updates bonds. You can do something with bonds here if i % saveEveryBlocks == (saveEveryBlocks - 2): a.doBlock(steps=steps, increment=doSave) if doSave: @@ -250,9 +220,7 @@ def initModel(): open(os.path.join(a.folder, "SMC{0}.dat".format(a.step)), "wb"), ) else: - a.integrator.step( - steps - ) # do steps without getting the positions from the GPU (faster) + a.integrator.step(steps) # do steps without getting the positions from the GPU (faster) data = a.getData() # save data and step, and delete the simulation block = a.step diff --git a/polychrom/cli/traj_convert.py b/polychrom/cli/traj_convert.py index af884ce..53d6a57 100644 --- a/polychrom/cli/traj_convert.py +++ b/polychrom/cli/traj_convert.py @@ -103,9 +103,7 @@ def _find_matches(pat, filenames): a = re.search(pat, filename) if a is not None: if len(a.groups()) != 1: - raise ValueError( - "You should have one group in regex denoting the number of the file" - ) + raise ValueError("You should have one group in regex denoting the number of the file") assert len(a.groups()) == 1 gr = int(a.groups()[0]) result[filename] = gr @@ -176,12 +174,8 @@ def _find_matches(pat, filenames): help="allow blocks to be non-consecutive (1,2,3...)", ) @click.option("--verbose", is_flag=True) -@click.option( - "--round-to", default=2, show_default=True, help="round to this number of digits" -) -@click.option( - "--skip-files", default=1, show_default=True, help="save only every Nth file" -) +@click.option("--round-to", default=2, show_default=True, help="round to this number of digits") +@click.option("--skip-files", default=1, show_default=True, help="save only every Nth file") @click.option( "--HDF5-blocks-per-file", default=100, @@ -251,10 +245,7 @@ def trajcopy( if kwargs["input_style"] == "old": blocks = _find_matches(block_pattern, all_files) elif kwargs["input_style"] == "new": - blocks = { - i: j - for j, i in list_URIs(in_dir, empty_error=True, return_dict=True).items() - } + blocks = {i: j for j, i in list_URIs(in_dir, empty_error=True, return_dict=True).items()} else: raise ValueError("input-style should be 'old' or 'new'") @@ -300,15 +291,11 @@ def trajcopy( assert len(extra_loader) == len(extra_require) # matching patterns for extra files, populating the dataframe - for val_pat, val_name, require in zip( - extra_pattern, extra_pattern_name, extra_require - ): + for val_pat, val_name, require in zip(extra_pattern, extra_pattern_name, extra_require): datas = _find_matches(val_pat, all_files) if require: if len(datas) != len(blocks): - raise ValueError( - f"files missing for {val_name}: need {len(blocks)} found {len(datas)}" - ) + raise ValueError(f"files missing for {val_name}: need {len(blocks)} found {len(datas)}") if len(datas) > 0: datas = pd.Series(data=list(datas.keys()), index=list(datas.values())) datas.name = val_name @@ -327,9 +314,7 @@ def trajcopy( print(other[:: len(other) // 20 + 1]) print("Verify that none of these should be converted using extra_pattern") print("If not, increase max_unmatched_files") - raise ValueError( - "Limit exceeded: {0} files did not match anything".format(len(other)) - ) + raise ValueError("Limit exceeded: {0} files did not match anything".format(len(other))) # creating the reporter if (len(blocks) > 0) and (not kwargs["dry_run"]): @@ -361,9 +346,7 @@ def trajcopy( cur["block"] = i elif kwargs["input_style"] == "new": cur = load_URI(data) - cur["pos"] = np.round( - np.asarray(cur["pos"], dtype=np.float32), kwargs["round_to"] - ) + cur["pos"] = np.round(np.asarray(cur["pos"], dtype=np.float32), kwargs["round_to"]) # adding "extra" data in the dict to save for name, ldr in zip(extra_pattern_name, extra_loader): diff --git a/polychrom/contactmaps.py b/polychrom/contactmaps.py index 933aa0f..c8adb20 100644 --- a/polychrom/contactmaps.py +++ b/polychrom/contactmaps.py @@ -193,9 +193,7 @@ def averageContactsSimple(contactIterator, inValues, N, **kwargs): uniqueContacts = kwargs.get("uniqueContacts", False) contactProcessing = kwargs.get("contactProcessing", lambda x: x) finalSize = N * (N + 1) // 2 - sharedArrays = [ - np.zeros(finalSize, dtype=arrayDtype) - ] # just an array, not a shared array here bc 1 core + sharedArrays = [np.zeros(finalSize, dtype=arrayDtype)] # just an array, not a shared array here bc 1 core argset = list(sharedArrays) + [ contactProcessing, classInitArgs, @@ -228,9 +226,7 @@ def worker(x): sharedNumpy = list(map(tonumpyarray, sharedArrays__)) # shared numpy arrays allContacts = [] contactSum = 0 - myIterator = contactIterator__( - x, *classInitArgs__, **classInitKwargs__ - ) # acquiring and initializing iterator + myIterator = contactIterator__(x, *classInitArgs__, **classInitKwargs__) # acquiring and initializing iterator stopped = False while True: # main event loop try: @@ -246,9 +242,7 @@ def worker(x): except StopIteration: stopped = True - if ( - contactSum > contactBlock__ - ) or stopped: # we aggregated enough contacts. ready to dump them. + if (contactSum > contactBlock__) or stopped: # we aggregated enough contacts. ready to dump them. if len(allContacts) == 0: return # no contacts found at all - exiting (we must be stopped) contactSum = 0 @@ -258,30 +252,20 @@ def worker(x): if stopped: # contactProcessing killed all contacts? are we done? return # if yes, exiting continue # if not, going to the next bucket - ctrue = indexing( - contacts[:, 0], contacts[:, 1], N__ - ) # converting to 1D contacts + ctrue = indexing(contacts[:, 0], contacts[:, 1], N__) # converting to 1D contacts position, counts = np.unique(ctrue, return_counts=True) # unique contacts assert position[0] >= 0 assert position[-1] < N__ * (N__ + 1) // 2 # boundary check for contacts - chunks = np.array( - np.r_[0, np.cumsum(list(map(len, sharedArrays__)))], dtype=int - ) - inds = np.searchsorted( - position, chunks - ) # assinging contacts to chunks here + chunks = np.array(np.r_[0, np.cumsum(list(map(len, sharedArrays__)))], dtype=int) + inds = np.searchsorted(position, chunks) # assinging contacts to chunks here if inds[-1] != len(position): raise ValueError # last chunks boundary must be after all contacts) indszip = list(zip(inds[:-1], inds[1:])) # extents of contact buckets indarray = list(range(len(sharedArrays__))) - random.shuffle( - indarray - ) # shuffled array of contactmap bucket indices we are going to work with + random.shuffle(indarray) # shuffled array of contactmap bucket indices we are going to work with for j, (st, end) in enumerate(indszip): - position[st:end] -= chunks[ - j - ] # pre-subtracting offsets now - not to do it when the lock is being held + position[st:end] -= chunks[j] # pre-subtracting offsets now - not to do it when the lock is being held while len(indarray) > 0: # continue until all contacts are put in buckets for i in range(len(indarray)): # going over all buckets @@ -290,20 +274,12 @@ def worker(x): if i == len(indarray): # is this the last bucket? lock.acquire() # wait for it to be free, and work with it else: - if not lock.acquire( - 0 - ): # not the last bucket? Try to acquire the lock + if not lock.acquire(0): # not the last bucket? Try to acquire the lock continue # if failed, move to the next bucket - st, end = indszip[ - ind - ] # succeeded acquiring the lock? Then do things with our bucket - sharedNumpy[ind][position[st:end]] += counts[ - st:end - ] # add to the current bucket + st, end = indszip[ind] # succeeded acquiring the lock? Then do things with our bucket + sharedNumpy[ind][position[st:end]] += counts[st:end] # add to the current bucket lock.release() - indarray.pop( - i - ) # remove the index of the bucket because it's finished + indarray.pop(i) # remove the index of the bucket because it's finished break # back to the main loop allContacts = [] if stopped: @@ -396,9 +372,7 @@ def next(self): # actual realization of the self.next method N, ] - if ( - not useFmap - ): # for mp.map we need initializer because shared memory cannot be pickled + if not useFmap: # for mp.map we need initializer because shared memory cannot be pickled # # or passed as an argument in inValues with closing(mp.Pool(processes=nproc, initializer=init, initargs=argset)) as p: p.map(worker, inValues) @@ -519,9 +493,7 @@ def binnedContactMap( Nbase = len(bins) - 1 if Nbase > 25000: - warnings.warn( - UserWarning("very large contact map" " may be difficult to visualize") - ) + warnings.warn(UserWarning("very large contact map" " may be difficult to visualize")) chromosomeStarts = np.cumsum(chainBinNums) chromosomeStarts = np.hstack((0, chromosomeStarts)) diff --git a/polychrom/contrib/integrators.py b/polychrom/contrib/integrators.py index 0c1c740..0ea12d3 100644 --- a/polychrom/contrib/integrators.py +++ b/polychrom/contrib/integrators.py @@ -37,9 +37,7 @@ from openmmtools.integrators import PrettyPrintableIntegrator -class ActiveBrownianIntegrator( - utils.RestorableOpenMMObject, PrettyPrintableIntegrator, mm.CustomIntegrator -): +class ActiveBrownianIntegrator(utils.RestorableOpenMMObject, PrettyPrintableIntegrator, mm.CustomIntegrator): """Brownian dynamics integrator with monomer Diffusion coefficient that varies along the chain. Parameters @@ -67,9 +65,7 @@ def __init__(self, timestep, collision_rate, particleD): self.addConstrainPositions() -class CorrelatedNoiseIntegrator( - utils.RestorableOpenMMObject, PrettyPrintableIntegrator, mm.CustomIntegrator -): +class CorrelatedNoiseIntegrator(utils.RestorableOpenMMObject, PrettyPrintableIntegrator, mm.CustomIntegrator): """Brownian dynamics integrator with correlated active noise. To define the correlations, we define a set of k attributes that specify the monomer's identity, such @@ -103,9 +99,7 @@ def __init__(self, timestep, collision_rate, particleD, rhos): self.addPerDofVariable("noise", 0) # noise term in Langevin equation self.addPerDofVariable("D", 0) # monomer diffusion coefficient self.addPerDofVariable("sigma", 0) # standard deviation of noise - self.addGlobalVariable( - "ghostx", 0 - ) # the ghost random variable used to correlate noise + self.addGlobalVariable("ghostx", 0) # the ghost random variable used to correlate noise self.addGlobalVariable("ghosty", 0) self.addGlobalVariable("ghostz", 0) self.addPerDofVariable("maskx", 0) # mask y and z directions diff --git a/polychrom/forcekits.py b/polychrom/forcekits.py index 5e3fa21..6aa9532 100644 --- a/polychrom/forcekits.py +++ b/polychrom/forcekits.py @@ -63,11 +63,7 @@ def polymer_chains( force_list = [] - bonds = ( - [] - if ((extra_bonds is None) or len(extra_bonds) == 0) - else [tuple(b) for b in extra_bonds] - ) + bonds = [] if ((extra_bonds is None) or len(extra_bonds) == 0) else [tuple(b) for b in extra_bonds] triplets = extra_triplets if extra_triplets else [] newchains = [] @@ -116,9 +112,7 @@ def polymer_chains( if except_bonds: exc = list(set([tuple(i) for i in np.sort(np.array(bonds), axis=1)])) if hasattr(nb_force, "addException"): - print( - "Exclude neighbouring chain particles from {}".format(nb_force.name) - ) + print("Exclude neighbouring chain particles from {}".format(nb_force.name)) for pair in exc: nb_force.addException(int(pair[0]), int(pair[1]), 0, 0, 0, True) @@ -126,12 +120,8 @@ def polymer_chains( # The built-in LJ nonbonded force uses "exclusions" instead of "exceptions" elif hasattr(nb_force, "addExclusion"): - print( - "Exclude neighbouring chain particles from {}".format(nb_force.name) - ) - nb_force.createExclusionsFromBonds( - [(int(b[0]), int(b[1])) for b in bonds], int(except_bonds) - ) + print("Exclude neighbouring chain particles from {}".format(nb_force.name)) + nb_force.createExclusionsFromBonds([(int(b[0]), int(b[1])) for b in bonds], int(except_bonds)) # for pair in exc: # nb_force.addExclusion(int(pair[0]), int(pair[1])) num_exc = nb_force.getNumExclusions() diff --git a/polychrom/forces.py b/polychrom/forces.py index dad2f1b..18d3670 100644 --- a/polychrom/forces.py +++ b/polychrom/forces.py @@ -95,9 +95,7 @@ def _check_bonds(bonds, N): for bond in set(bonds): bonds.remove(bond) - raise ValueError( - f"Bonds {bonds} are repeated. Set override_checks=True to override this check." - ) + raise ValueError(f"Bonds {bonds} are repeated. Set override_checks=True to override this check.") # check that all monomers make at least one bond monomer_not_in_bond = ~np.zeros(N).astype(bool) @@ -122,9 +120,7 @@ def _check_angle_bonds(triplets): for triplet in set(triplets): triplets.remove(triplet) - raise ValueError( - f"Triplets {triplets} are repeated. Set override_checks=True to override this check." - ) + raise ValueError(f"Triplets {triplets} are repeated. Set override_checks=True to override this check.") # check that no triplet of the form (i, i, j) exists # check that no bonds of the form (i, i) exist @@ -204,9 +200,7 @@ def harmonic_bonds( force.name = name bondLength = _to_array_1d(bondLength, len(bonds)) * sim_object.length_scale - bondWiggleDistance = ( - _to_array_1d(bondWiggleDistance, len(bonds)) * sim_object.length_scale - ) + bondWiggleDistance = _to_array_1d(bondWiggleDistance, len(bonds)) * sim_object.length_scale # using kbondScalingFactor because force accepts parameters with units kbond = sim_object.kbondScalingFactor / (bondWiggleDistance**2) @@ -215,13 +209,10 @@ def harmonic_bonds( for bond_idx, (i, j) in enumerate(bonds): if (i >= sim_object.N) or (j >= sim_object.N): raise ValueError( - "\nCannot add bond with monomers %d,%d that" - "are beyound the polymer length %d" % (i, j, sim_object.N) + "\nCannot add bond with monomers %d,%d that" "are beyound the polymer length %d" % (i, j, sim_object.N) ) - force.addBond( - int(i), int(j), float(bondLength[bond_idx]), float(kbond[bond_idx]) - ) + force.addBond(int(i), int(j), float(bondLength[bond_idx]), float(kbond[bond_idx])) return force @@ -267,11 +258,7 @@ def constant_force_bonds( if not override_checks: _check_bonds(bonds, sim_object.N) - energy = ( - f"(1. / wiggle) * univK * " - f"(sqrt((r-r0 * conlen)* " - f" (r - r0 * conlen) + a * a) - a)" - ) + energy = f"(1. / wiggle) * univK * " f"(sqrt((r-r0 * conlen)* " f" (r - r0 * conlen) + a * a) - a)" force = openmm.CustomBondForce(energy) force.name = name @@ -282,15 +269,12 @@ def constant_force_bonds( force.addGlobalParameter("conlen", sim_object.conlen) bondLength = _to_array_1d(bondLength, len(bonds)) * sim_object.length_scale - bondWiggleDistance = ( - _to_array_1d(bondWiggleDistance, len(bonds)) * sim_object.length_scale - ) + bondWiggleDistance = _to_array_1d(bondWiggleDistance, len(bonds)) * sim_object.length_scale for bond_idx, (i, j) in enumerate(bonds): if (i >= sim_object.N) or (j >= sim_object.N): raise ValueError( - "\nCannot add bond with monomers %d,%d that" - "are beyound the polymer length %d" % (i, j, sim_object.N) + "\nCannot add bond with monomers %d,%d that" "are beyound the polymer length %d" % (i, j, sim_object.N) ) force.addBond( @@ -302,9 +286,7 @@ def constant_force_bonds( return force -def angle_force( - sim_object, triplets, k=1.5, theta_0=np.pi, name="angle", override_checks=False -): +def angle_force(sim_object, triplets, k=1.5, theta_0=np.pi, name="angle", override_checks=False): """Adds harmonic angle bonds. k specifies energy in kT at one radian If k is an array, it has to be of the length N. Xth value then specifies stiffness of the angle centered at @@ -353,9 +335,7 @@ def angle_force( return force -def polynomial_repulsive( - sim_object, trunc=3.0, radiusMult=1.0, name="polynomial_repulsive" -): +def polynomial_repulsive(sim_object, trunc=3.0, radiusMult=1.0, name="polynomial_repulsive"): """This is a simple polynomial repulsive potential. It has the value of `trunc` at zero, stays flat until 0.6-0.7 and then drops to zero together with its first derivative at r=1.0. @@ -461,9 +441,7 @@ def smooth_square_well( force.addGlobalParameter("REPsigma", repulsionRadius * sim_object.conlen) force.addGlobalParameter("ATTRe", attractionEnergy * sim_object.kT) - force.addGlobalParameter( - "ATTRdelta", sim_object.conlen * (attractionRadius - repulsionRadius) / 2.0 - ) + force.addGlobalParameter("ATTRdelta", sim_object.conlen * (attractionRadius - repulsionRadius) / 2.0) # Coefficients for the minimum of x^12*(x*x-1) force.addGlobalParameter("emin12", 46656.0 / 823543.0) force.addGlobalParameter("rmin12", np.sqrt(6.0 / 7.0)) @@ -574,9 +552,7 @@ def selective_SSW( force.addGlobalParameter("ATTRe", attractionEnergy * sim_object.kT) force.addGlobalParameter("ATTReAdd", selectiveAttractionEnergy * sim_object.kT) - force.addGlobalParameter( - "ATTRdelta", sim_object.conlen * (attractionRadius - repulsionRadius) / 2.0 - ) + force.addGlobalParameter("ATTRdelta", sim_object.conlen * (attractionRadius - repulsionRadius) / 2.0) # Coefficients for x^12*(x*x-1) force.addGlobalParameter("emin12", 46656.0 / 823543.0) @@ -697,13 +673,11 @@ def heteropolymer_SSW( "ATTReTot = ATTRe" ) if len(indexpairs) > 0: - energy += ( - " + ATTReAdd*(delta(type1-{0:d})*delta(type2-{1:d})" "*INT_{0:d}_{1:d}" - ).format(indexpairs[0][0], indexpairs[0][1]) + energy += (" + ATTReAdd*(delta(type1-{0:d})*delta(type2-{1:d})" "*INT_{0:d}_{1:d}").format( + indexpairs[0][0], indexpairs[0][1] + ) for i, j in indexpairs[1:]: - energy += "+delta(type1-{0:d})*delta(type2-{1:d})*INT_{0:d}_{1:d}".format( - i, j - ) + energy += "+delta(type1-{0:d})*delta(type2-{1:d})*INT_{0:d}_{1:d}".format(i, j) energy += ")" energy += ( ";" @@ -729,18 +703,14 @@ def heteropolymer_SSW( force.addGlobalParameter("ATTRe", attractionEnergy * sim_object.kT) force.addGlobalParameter("ATTReAdd", selectiveAttractionEnergy * sim_object.kT) - force.addGlobalParameter( - "ATTRdelta", sim_object.conlen * (attractionRadius - repulsionRadius) / 2.0 - ) + force.addGlobalParameter("ATTRdelta", sim_object.conlen * (attractionRadius - repulsionRadius) / 2.0) # Coefficients for x^12*(x*x-1) force.addGlobalParameter("emin12", 46656.0 / 823543.0) force.addGlobalParameter("rmin12", np.sqrt(6.0 / 7.0)) for i, j in indexpairs: - force.addGlobalParameter( - "INT_{0:d}_{1:d}".format(i, j), interactionMatrix[i, j] - ) + force.addGlobalParameter("INT_{0:d}_{1:d}".format(i, j), interactionMatrix[i, j]) force.addPerParticleParameter("type") force.addPerParticleParameter("ExtraHard") @@ -751,9 +721,7 @@ def heteropolymer_SSW( return force -def cylindrical_confinement( - sim_object, r, bottom=None, k=0.1, top=9999, name="cylindrical_confinement" -): +def cylindrical_confinement(sim_object, r, bottom=None, k=0.1, top=9999, name="cylindrical_confinement"): """As it says.""" if bottom == True: @@ -771,8 +739,7 @@ def cylindrical_confinement( ) else: force = openmm.CustomExternalForce( - "kt * k * step(dr) * (sqrt(dr*dr + t*t) - t);" - "dr = sqrt(x^2 + y^2 + tt^2) - r + 10*t" + "kt * k * step(dr) * (sqrt(dr*dr + t*t) - t);" "dr = sqrt(x^2 + y^2 + tt^2) - r + 10*t" ) force.name = name @@ -857,9 +824,7 @@ def spherical_confinement( return force -def spherical_well( - sim_object, particles, r, center=[0, 0, 0], width=1, depth=1, name="spherical_well" -): +def spherical_well(sim_object, particles, r, center=[0, 0, 0], width=1, depth=1, name="spherical_well"): """ A spherical potential well, suited for example to simulate attraction to a lamina. @@ -904,9 +869,7 @@ def spherical_well( return force -def tether_particles( - sim_object, particles, *, pbc=False, k=30, positions="current", name="Tethers" -): +def tether_particles(sim_object, particles, *, pbc=False, k=30, positions="current", name="Tethers"): """tethers particles in the 'particles' array. Increase k to tether them stronger, but watch the system! @@ -984,18 +947,14 @@ def pull_force(sim_object, particles, force_vecs, name="Pull"): force.addPerParticleParameter("fy") force.addPerParticleParameter("fz") - for num, force_vec in itertools.zip_longest( - particles, force_vecs, fillvalue=force_vecs[-1] - ): + for num, force_vec in itertools.zip_longest(particles, force_vecs, fillvalue=force_vecs[-1]): force_vec = [float(f) * (sim_object.kT / sim_object.conlen) for f in force_vec] force.addParticle(int(num), force_vec) return force -def grosberg_polymer_bonds( - sim_object, bonds, k=30, name="grosberg_polymer", override_checks=False -): +def grosberg_polymer_bonds(sim_object, bonds, k=30, name="grosberg_polymer", override_checks=False): """Adds FENE bonds according to Halverson-Grosberg paper. (Halverson, Jonathan D., et al. "Molecular dynamics simulation study of nonconcatenated ring polymers in a melt. I. Statics." @@ -1023,16 +982,13 @@ def grosberg_polymer_bonds( force = openmm.CustomBondForce(equation) force.name = name - force.addGlobalParameter( - "k", k * sim_object.kT / (sim_object.conlen * sim_object.conlen) - ) + force.addGlobalParameter("k", k * sim_object.kT / (sim_object.conlen * sim_object.conlen)) force.addGlobalParameter("r0", sim_object.conlen * 1.5) for bond_idx, (i, j) in enumerate(bonds): if (i >= sim_object.N) or (j >= sim_object.N): raise ValueError( - "\nCannot add bond with monomers %d,%d that" - "are beyound the polymer length %d" % (i, j, sim_object.N) + "\nCannot add bond with monomers %d,%d that" "are beyound the polymer length %d" % (i, j, sim_object.N) ) force.addBond(int(i), int(j)) @@ -1040,9 +996,7 @@ def grosberg_polymer_bonds( return force -def grosberg_angle( - sim_object, triplets, k=1.5, name="grosberg_angle", override_checks=False -): +def grosberg_angle(sim_object, triplets, k=1.5, name="grosberg_angle", override_checks=False): """ Adds stiffness according to the Grosberg paper. (Halverson, Jonathan D., et al. "Molecular dynamics simulation study of diff --git a/polychrom/hdf5_format.py b/polychrom/hdf5_format.py index 0f06510..19fee68 100644 --- a/polychrom/hdf5_format.py +++ b/polychrom/hdf5_format.py @@ -186,9 +186,7 @@ def list_URIs(folder, empty_error=True, read_error=True, return_dict=False): if not return_dict: return [i[1] for i in sorted(filenames.items(), key=lambda x: int(x[0]))] else: - return { - int(i[0]): i[1] for i in sorted(filenames.items(), key=lambda x: int(x[0])) - } + return {int(i[0]): i[1] for i in sorted(filenames.items(), key=lambda x: int(x[0]))} def load_URI(dset_path): @@ -289,8 +287,7 @@ def __init__( os.remove(file_path) else: raise IOError( - "Subfolder in traj folder; not deleting. Ensure folder is " - "correct and delete manually. " + "Subfolder in traj folder; not deleting. Ensure folder is " "correct and delete manually. " ) if check_exists: @@ -298,9 +295,7 @@ def __init__( for the_file in os.listdir(folder): for prefix in self.prefixes: if the_file.startswith(prefix): - raise RuntimeError( - f"folder {folder} is not empty: set check_exists=False to ignore" - ) + raise RuntimeError(f"folder {folder} is not empty: set check_exists=False to ignore") def continue_trajectory(self, continue_from=None, continue_max_delete=5): """ @@ -338,9 +333,7 @@ def continue_trajectory(self, continue_from=None, continue_max_delete=5): """ uris = list_URIs(self.folder, return_dict=True) - uri_inds = np.array( - list(uris.keys()) - ) # making a list of all URIs and filenames + uri_inds = np.array(list(uris.keys())) # making a list of all URIs and filenames uri_vals = np.array(list(uris.values())) uri_fnames = np.array([i.split("::")[0] for i in uris.values()]) if continue_from is None: @@ -349,29 +342,21 @@ def continue_trajectory(self, continue_from=None, continue_max_delete=5): if int(continue_from) not in uris: raise ValueError(f"block {continue_from} not in folder") - ind = np.nonzero(uri_inds == continue_from)[0][ - 0 - ] # position of a starting block in arrays + ind = np.nonzero(uri_inds == continue_from)[0][0] # position of a starting block in arrays newdata = load_URI(uri_vals[ind]) todelete = np.nonzero(uri_inds >= continue_from)[0] if len(todelete) > continue_max_delete: - raise ValueError( - "Refusing to delete {uris_delete} blocks - set continue_max_delete accordingly" - ) + raise ValueError("Refusing to delete {uris_delete} blocks - set continue_max_delete accordingly") fnames_delete = np.unique(uri_fnames[todelete]) inds_tosave = np.nonzero((uri_fnames == uri_fnames[ind]) * (uri_inds <= ind))[0] - for ( - saveind - ) in inds_tosave: # we are saving some data and deleting the whole last file + for saveind in inds_tosave: # we are saving some data and deleting the whole last file self.datas[uri_inds[saveind]] = load_URI(uri_vals[saveind]) self.counter["data"] = ind + 1 - files = os.listdir( - self.folder - ) # some heuristics to infer values of counters - not crucial but maybe useful + files = os.listdir(self.folder) # some heuristics to infer values of counters - not crucial but maybe useful for prefix in self.prefixes: if prefix != "data": myfiles = [i for i in files if i.startswith(prefix)] diff --git a/polychrom/legacy/contactmaps.py b/polychrom/legacy/contactmaps.py index 101a54a..d0e1bb4 100644 --- a/polychrom/legacy/contactmaps.py +++ b/polychrom/legacy/contactmaps.py @@ -184,9 +184,7 @@ def averageBinnedContactMap( Nbase = len(bins) - 1 if Nbase > 10000: - warnings.warn( - UserWarning("very large contact map" " may be difficult to visualize") - ) + warnings.warn(UserWarning("very large contact map" " may be difficult to visualize")) chromosomeStarts = np.cumsum(chainBinNums) chromosomeStarts = np.hstack((0, chromosomeStarts)) @@ -287,17 +285,10 @@ def myaction(values): # our worker receives some filenames if len(data) > 6000: warnings.warn( - UserWarning( - "very large contact map" - " may cause errors. these may be fixed with n=1 threads." - ) + UserWarning("very large contact map" " may cause errors. these may be fixed with n=1 threads.") ) if len(data) > 20000: - warnings.warn( - UserWarning( - "very large contact map" " may be difficult to visualize." - ) - ) + warnings.warn(UserWarning("very large contact map" " may be difficult to visualize.")) mysum = pureMap(data, cutoff) # create a map diff --git a/polychrom/legacy/forces.py b/polychrom/legacy/forces.py index 1312da2..0e52fe5 100644 --- a/polychrom/legacy/forces.py +++ b/polychrom/legacy/forces.py @@ -24,9 +24,7 @@ def minimizing_repulsive_Force(sim_object): nbCutOffDist = radius * 1.0 repul_energy = "1000* REPe * (1-r/REPr)^2 " - sim_object.force_dict["Nonbonded_minimizing_Force"] = openmm.CustomNonbondedForce( - repul_energy - ) + sim_object.force_dict["Nonbonded_minimizing_Force"] = openmm.CustomNonbondedForce(repul_energy) repulforceGr = sim_object.force_dict["Nonbonded_minimizing_Force"] repulforceGr.addGlobalParameter("REPe", sim_object.kT) repulforceGr.addGlobalParameter("REPr", sim_object.kT) @@ -35,9 +33,7 @@ def minimizing_repulsive_Force(sim_object): repulforceGr.setCutoffDistance(nbCutOffDist) -def fix_particles_Z_coordinate( - sim_object, particles, zCoordinates, k=0.3, useOtherAxis="z", mode="abs", gap=None -): +def fix_particles_Z_coordinate(sim_object, particles, zCoordinates, k=0.3, useOtherAxis="z", mode="abs", gap=None): """Limits position of a set of particles in z coordinate Parameters @@ -63,14 +59,10 @@ def fix_particles_Z_coordinate( start, stop = tuple(zCoordinates) zCoordinates = [] for par in particles: - zCoordinates.append( - start + float(stop - start) * (par / float(sim_object.N)) - ) + zCoordinates.append(start + float(stop - start) * (par / float(sim_object.N))) if (mode == "abs") and (gap is None): - zFixForce = openmm.CustomExternalForce( - "ZFIXk * (sqrt((%s - ZFIXr0)^2 + ZFIXa^2) - ZFIXa)" % (useOtherAxis,) - ) + zFixForce = openmm.CustomExternalForce("ZFIXk * (sqrt((%s - ZFIXr0)^2 + ZFIXa^2) - ZFIXa)" % (useOtherAxis,)) zFixForce.addGlobalParameter("ZFIXk", k * sim_object.kT / (sim_object.conlen)) elif (mode == "abs") and (gap is not None): zFixForce = openmm.CustomExternalForce( @@ -85,24 +77,17 @@ def fix_particles_Z_coordinate( zFixForce.addGlobalParameter("ZFIXgap", sim_object.conlen * gap) elif (mode == "quadratic") and (gap is None): - zFixForce = openmm.CustomExternalForce( - "ZFIXk * ((%s - ZFIXr0)^2)" % (useOtherAxis,) - ) - zFixForce.addGlobalParameter( - "ZFIXk", k * sim_object.kT / (sim_object.conlen * sim_object.conlen) - ) + zFixForce = openmm.CustomExternalForce("ZFIXk * ((%s - ZFIXr0)^2)" % (useOtherAxis,)) + zFixForce.addGlobalParameter("ZFIXk", k * sim_object.kT / (sim_object.conlen * sim_object.conlen)) elif (mode == "quadratic") and (gap is not None): zFixForce = openmm.CustomExternalForce( "ZFIXk * (step(%s - ZFIXr0 - ZFIXgap * 0.5) * " "(%s - ZFIXr0 - ZFIXgap * 0.5)^2 + " "step(-%s + ZFIXr0 - ZFIXgap * 0.5) * " - "(-%s + ZFIXr0 - ZFIXgap * 0.5)^2)" - % (useOtherAxis, useOtherAxis, useOtherAxis, useOtherAxis) + "(-%s + ZFIXr0 - ZFIXgap * 0.5)^2)" % (useOtherAxis, useOtherAxis, useOtherAxis, useOtherAxis) ) - zFixForce.addGlobalParameter( - "ZFIXk", k * sim_object.kT / (sim_object.conlen * sim_object.conlen) - ) + zFixForce.addGlobalParameter("ZFIXk", k * sim_object.kT / (sim_object.conlen * sim_object.conlen)) zFixForce.addGlobalParameter("ZFIXgap", sim_object.conlen * gap) else: raise RuntimeError("Not implemented") @@ -132,9 +117,7 @@ def lamina_attraction(sim_object, width=1, depth=1, r=None): from previously defined spherical potential. """ - sim_object.metadata["laminaAttraction"] = repr( - {"width": width, "depth": depth, "r": r} - ) + sim_object.metadata["laminaAttraction"] = repr({"width": width, "depth": depth, "r": r}) laminaForce = openmm.CustomExternalForce( "step(LAMr-LAMaa + LAMwidth) * step(LAMaa + LAMwidth - LAMr) " "* LAMdepth * (LAMr-LAMaa + LAMwidth) * (LAMaa + LAMwidth - LAMr) " @@ -151,10 +134,7 @@ def lamina_attraction(sim_object, width=1, depth=1, r=None): try: r = sim_object.sphericalConfinementRadius except: - raise ValueError( - "No spherical confinement radius defined" - " yet. Apply spherical confinement first!" - ) + raise ValueError("No spherical confinement radius defined" " yet. Apply spherical confinement first!") if sim_object.verbose == True: print("Lamina attraction added with r = %d" % r) @@ -190,9 +170,7 @@ def useDomains(sim_object, domains=None, filename=None): if len(sim_object.domains) != sim_object.N: sim_object._exitProgram("Wrong domain lengths") - pickle.dump( - sim_object.domains, open(os.path.join(sim_object.folder, "domains.dat"), "wb") - ) + pickle.dump(sim_object.domains, open(os.path.join(sim_object.folder, "domains.dat"), "wb")) def lennard_jones_force( @@ -249,10 +227,7 @@ def lennard_jones_force( ) if blindFraction > 0.99: - sim_object._exitProgram( - "why do you need this force without particles???" - " set blindFraction between 0 and 1" - ) + sim_object._exitProgram("why do you need this force without particles???" " set blindFraction between 0 and 1") if (sigmaRep is None) and (sigmaAttr is None): sigmaAttr = sigmaRep = sim_object.conlen else: @@ -331,9 +306,7 @@ def attractive_interaction(sim_object, i, j, epsilon, sigma=None, length=3): """ if type(sim_object.force_dict["Nonbonded"]) != openmm.NonbondedForce: - sim_object.exit( - "Cannot add interactions" " without Lennard-Jones nonbonded force" - ) + sim_object.exit("Cannot add interactions" " without Lennard-Jones nonbonded force") if sigma is None: sigma = 1.1 * sim_object.conlen @@ -342,24 +315,16 @@ def attractive_interaction(sim_object, i, j, epsilon, sigma=None, length=3): print("!!!!!!!!!bond with %d and %d is out of range!!!!!" % (i, j)) return repulforce = sim_object.force_dict["Nonbonded"] - for t1 in range( - int(np.ceil(i - length / 2)), int(np.ceil(i + (length - length / 2))) - ): - for t2 in range( - int(np.ceil(j - length / 2)), int(np.ceil(j + (length - length / 2))) - ): + for t1 in range(int(np.ceil(i - length / 2)), int(np.ceil(i + (length - length / 2)))): + for t2 in range(int(np.ceil(j - length / 2)), int(np.ceil(j + (length - length / 2)))): repulforce.addException(t1, t2, 0, sigma, epsilon, True) if sim_object.verbose == True: print("Exception added between" " particles %d and %d" % (t1, t2)) for tt in range(i - length, i + length): - repulforce.setParticleParameters( - tt, 0, sim_object.conlen, sim_object.epsilonRep - ) + repulforce.setParticleParameters(tt, 0, sim_object.conlen, sim_object.epsilonRep) for tt in range(j - length, j + length): - repulforce.setParticleParameters( - tt, 0, sim_object.conlen, sim_object.epsilonRep - ) + repulforce.setParticleParameters(tt, 0, sim_object.conlen, sim_object.epsilonRep) def gravity(sim_object, k=0.1, cutoff=None): @@ -411,9 +376,7 @@ def attraction_to_the_core(sim_object, k, r0, coreParticles=[]): attractForce = openmm.CustomExternalForce( " COREk * ((COREr - CORErn) ^ 2) ; " "COREr = sqrt(x^2 + y^2 + COREtt^2)" ) - attractForce.addGlobalParameter( - "COREk", k * sim_object.kT / (sim_object.conlen * sim_object.conlen) - ) + attractForce.addGlobalParameter("COREk", k * sim_object.kT / (sim_object.conlen * sim_object.conlen)) attractForce.addGlobalParameter("CORErn", r0 * sim_object.conlen) attractForce.addGlobalParameter("COREtt", 0.001 * sim_object.conlen) sim_object.force_dict["CoreAttraction"] = attractForce @@ -423,12 +386,9 @@ def attraction_to_the_core(sim_object, k, r0, coreParticles=[]): if r0 > 0.1: excludeForce = openmm.CustomExternalForce( - " CORE2k * ((CORE2r - CORE2rn) ^ 2) * step(CORE2rn - CORE2r) ;" - "CORE2r = sqrt(x^2 + y^2 + CORE2tt^2)" - ) - excludeForce.addGlobalParameter( - "CORE2k", k * sim_object.kT / (sim_object.conlen * sim_object.conlen) + " CORE2k * ((CORE2r - CORE2rn) ^ 2) * step(CORE2rn - CORE2r) ;" "CORE2r = sqrt(x^2 + y^2 + CORE2tt^2)" ) + excludeForce.addGlobalParameter("CORE2k", k * sim_object.kT / (sim_object.conlen * sim_object.conlen)) excludeForce.addGlobalParameter("CORE2rn", r0 * sim_object.conlen) excludeForce.addGlobalParameter("CORE2tt", 0.001 * sim_object.conlen) sim_object.force_dict["CoreExclusion"] = excludeForce @@ -482,10 +442,7 @@ def spherical_well(sim_object, r=10, depth=1): try: r = sim_object.sphericalConfinementRadius * 0.5 except: - exit( - "No spherical confinement radius defined yet." - " Apply spherical confinement first!" - ) + exit("No spherical confinement radius defined yet." " Apply spherical confinement first!") if sim_object.verbose == True: print("Well attraction added with r = %d" % r) @@ -549,10 +506,7 @@ def add_lamina_attraction(sim_object, width=1, depth=1, r=None, particles=None): try: r = sim_object.sphericalConfinementRadius except: - exit( - "No spherical confinement radius defined yet." - "Apply spherical confinement first!" - ) + exit("No spherical confinement radius defined yet." "Apply spherical confinement first!") if sim_object.verbose == True: print("Lamina attraction added with r = %d" % r) @@ -567,9 +521,7 @@ def add_lamina_attraction(sim_object, width=1, depth=1, r=None, particles=None): # old energy minimization -def old_energy_minimization( - sim_object, stepsPerIteration=100, maxIterations=1000, failNotConverged=True -): +def old_energy_minimization(sim_object, stepsPerIteration=100, maxIterations=1000, failNotConverged=True): """Runs system at smaller timestep and higher collision rate to resolve possible conflicts. @@ -580,13 +532,9 @@ def old_energy_minimization( print("Performing energy minimization") sim_object._apply_forces() if (maxIterations is True) or (maxIterations is False): - raise ValueError( - "Please stop using the old notation and read the new energy minimization code" - ) + raise ValueError("Please stop using the old notation and read the new energy minimization code") if (failNotConverged is not True) and (failNotConverged is not False): - raise ValueError( - "Please stop using the old notation and read the new energy minimization code" - ) + raise ValueError("Please stop using the old notation and read the new energy minimization code") def_step = sim_object.integrator.getStepSize() def_fric = sim_object.integrator.getFriction() @@ -639,9 +587,7 @@ def check_connectivity(sim_object, newcoords=None, maxBondSizeMultipler=10): """ if not hasattr(sim_object, "bondLengths"): - raise ValueError( - "must use either harmonic or abs bonds to use checkConnectivty" - ) + raise ValueError("must use either harmonic or abs bonds to use checkConnectivty") if newcoords == None: newcoords = sim_object.get_data() @@ -653,11 +599,7 @@ def check_connectivity(sim_object, newcoords=None, maxBondSizeMultipler=10): bondArray = np.array(sim_object.bondLengths) bondDists = np.sqrt( np.sum( - ( - newcoords[np.array(bondArray[:, 0], dtype=int)] - - newcoords[np.array(bondArray[:, 1], dtype=int)] - ) - ** 2, + (newcoords[np.array(bondArray[:, 0], dtype=int)] - newcoords[np.array(bondArray[:, 1], dtype=int)]) ** 2, axis=1, ) ) diff --git a/polychrom/lib/extrusion.py b/polychrom/lib/extrusion.py index 3ac319c..78df4d4 100644 --- a/polychrom/lib/extrusion.py +++ b/polychrom/lib/extrusion.py @@ -37,9 +37,7 @@ def setup(self, bondForce, blocks=100, smcStepsPerBlock=1): """ if len(self.allBonds) != 0: - raise ValueError( - "Not all bonds were used; {0} sets left".format(len(self.allBonds)) - ) + raise ValueError("Not all bonds were used; {0} sets left".format(len(self.allBonds))) self.bondForce = bondForce @@ -48,10 +46,7 @@ def setup(self, bondForce, blocks=100, smcStepsPerBlock=1): loaded_positions = self.LEFpositions[self.curtime : self.curtime + blocks] allBonds = [ - [ - (int(loaded_positions[i, j, 0]), int(loaded_positions[i, j, 1])) - for j in range(loaded_positions.shape[1]) - ] + [(int(loaded_positions[i, j, 0]), int(loaded_positions[i, j, 1])) for j in range(loaded_positions.shape[1])] for i in range(blocks) ] @@ -63,14 +58,8 @@ def setup(self, bondForce, blocks=100, smcStepsPerBlock=1): self.curBonds = allBonds.pop(0) for bond in self.uniqueBonds: - paramset = ( - self.activeParamDict - if (bond in self.curBonds) - else self.inactiveParamDict - ) - ind = bondForce.addBond( - bond[0], bond[1], **paramset - ) # changed from addBond + paramset = self.activeParamDict if (bond in self.curBonds) else self.inactiveParamDict + ind = bondForce.addBond(bond[0], bond[1], **paramset) # changed from addBond self.bondInds.append(ind) self.bondToInd = {i: j for i, j in zip(self.uniqueBonds, self.bondInds)} @@ -86,9 +75,7 @@ def step(self, context, verbose=False): :return: (current bonds, previous step bonds); just for reference """ if len(self.allBonds) == 0: - raise ValueError( - "No bonds left to run; you should restart simulation and run setup again" - ) + raise ValueError("No bonds left to run; you should restart simulation and run setup again") pastBonds = self.curBonds self.curBonds = self.allBonds.pop(0) # getting current bonds @@ -106,10 +93,6 @@ def step(self, context, verbose=False): for bond, isAdd in zip(bondsToChange, bondsIsAdd): ind = self.bondToInd[bond] paramset = self.activeParamDict if isAdd else self.inactiveParamDict - self.bondForce.setBondParameters( - ind, bond[0], bond[1], **paramset - ) # actually updating bonds - self.bondForce.updateParametersInContext( - context - ) # now run this to update things in the context + self.bondForce.setBondParameters(ind, bond[0], bond[1], **paramset) # actually updating bonds + self.bondForce.updateParametersInContext(context) # now run this to update things in the context return self.curBonds, pastBonds diff --git a/polychrom/polymer_analyses.py b/polychrom/polymer_analyses.py index 6c93625..954e00c 100644 --- a/polychrom/polymer_analyses.py +++ b/polychrom/polymer_analyses.py @@ -293,9 +293,7 @@ def R2_scaling(data, bins=None, ring=False): for i in range(len(bins)): length = bins[i] if ring: - rads[i] = np.mean( - (np.sum((data[:, :N] - data[:, length : length + N]) ** 2, 0)) - ) + rads[i] = np.mean((np.sum((data[:, :N] - data[:, length : length + N]) ** 2, 0))) else: rads[i] = np.mean((np.sum((data[:, :-length] - data[:, length:]) ** 2, 0))) return np.array(bins), rads @@ -372,10 +370,7 @@ def combine_values(in_df): """ average_arrs = pd.Series( index=ndarray_cols, - data=[ - ndarray_agg([np.asarray(j) for j in in_df[i].values]) - for i in ndarray_cols - ], + data=[ndarray_agg([np.asarray(j) for j in in_df[i].values]) for i in ndarray_cols], ) average_values = value_agg(in_df[value_cols]) sample_values = in_df[sample_cols].iloc[0] @@ -553,9 +548,7 @@ def getLinkingNumber(data1, data2, simplify=True, randomOffset=True, verbose=Fal return _polymer_math.getLinkingNumber(data1, data2, randomOffset=randomOffset) -def calculate_cistrans( - data, chains, chain_id=0, cutoff=5, pbc_box=False, box_size=None -): +def calculate_cistrans(data, chains, chain_id=0, cutoff=5, pbc_box=False, box_size=None): """ Analysis of the territoriality of polymer chains from simulations, using the cis/trans ratio. diff --git a/polychrom/polymerutils.py b/polychrom/polymerutils.py index ba95041..549a40a 100644 --- a/polychrom/polymerutils.py +++ b/polychrom/polymerutils.py @@ -167,9 +167,7 @@ def save(data, filename, mode="txt", pdbGroups=None): raise ValueError("Not sure what to do with filename {0}".format(filename)) elif mode == "pdb": - data = ( - data - np.minimum(np.min(data, axis=0), np.zeros(3, float) - 100)[None, :] - ) + data = data - np.minimum(np.min(data, axis=0), np.zeros(3, float) - 100)[None, :] retret = "" def add(st, n): @@ -219,13 +217,7 @@ def add(st, n): def rotation_matrix(rotate): """Calculates rotation matrix based on three rotation angles""" tx, ty, tz = rotate - Rx = np.array( - [[1, 0, 0], [0, np.cos(tx), -np.sin(tx)], [0, np.sin(tx), np.cos(tx)]] - ) - Ry = np.array( - [[np.cos(ty), 0, -np.sin(ty)], [0, 1, 0], [np.sin(ty), 0, np.cos(ty)]] - ) - Rz = np.array( - [[np.cos(tz), -np.sin(tz), 0], [np.sin(tz), np.cos(tz), 0], [0, 0, 1]] - ) + Rx = np.array([[1, 0, 0], [0, np.cos(tx), -np.sin(tx)], [0, np.sin(tx), np.cos(tx)]]) + Ry = np.array([[np.cos(ty), 0, -np.sin(ty)], [0, 1, 0], [np.sin(ty), 0, np.cos(ty)]]) + Rz = np.array([[np.cos(tz), -np.sin(tz), 0], [np.sin(tz), np.cos(tz), 0], [0, 0, 1]]) return np.dot(Rx, np.dot(Ry, Rz)) diff --git a/polychrom/pymol_show.py b/polychrom/pymol_show.py index 328fa17..ff07924 100644 --- a/polychrom/pymol_show.py +++ b/polychrom/pymol_show.py @@ -197,15 +197,11 @@ def getSelectionString(start, end): if seg1 == seg2: return "resi {atom1}-{atom2} and segi {seg1}".format(**locals()) elif np.abs(seg1 - seg2) == 1: - return "(resi {atom1}-{maxNum} and segi {seg1}) or (resi 0-{atom2} and segi {seg2})".format( - **locals() - ) + return "(resi {atom1}-{maxNum} and segi {seg1}) or (resi 0-{atom2} and segi {seg2})".format(**locals()) elif np.abs(seg1 - seg2) >= 2: - line = "(resi {atom1}-{maxNum} and segi {seg1}) or (resi 0-{atom2} and segi {seg2})".format( - **locals() - ) + line = "(resi {atom1}-{maxNum} and segi {seg1}) or (resi 0-{atom2} and segi {seg2})".format(**locals()) for i in range(seg1 + 1, seg2): line = line + " or (segi {0})".format(i) return line @@ -305,9 +301,7 @@ def getSelectionString(start, end): elif (showChain == "none") or (not showChain): pass else: - raise ValueError( - "please select showChain to be 'worm' or 'spheres' or 'none'" - ) + raise ValueError("please select showChain to be 'worm' or 'spheres' or 'none'") for i in spherePositions: out.write("select {0} and {1}\n".format(name, getSelectionString(i, i))) @@ -346,11 +340,7 @@ def getSelectionString(start, end): sleep(0.5) - print( - os.system( - "pymol {1} -u {0} {2}".format(out.name, tmpPdbFilename, miscArguments) - ) - ) + print(os.system("pymol {1} -u {0} {2}".format(out.name, tmpPdbFilename, miscArguments))) return script @@ -522,11 +512,7 @@ def new_coloring( sleep(0.5) - print( - os.system( - "pymol {1} -u {0} {2}".format(out.name, tmpPdbFilename, miscArguments) - ) - ) + print(os.system("pymol {1} -u {0} {2}".format(out.name, tmpPdbFilename, miscArguments))) return script @@ -538,9 +524,7 @@ def getTmpPath(folder=None, **kwargs): return tmpPath, tmpFilename -def show_chain( - data, showGui=True, saveTo=None, showChain="worm", chains=None, support="", **kwargs -): +def show_chain(data, showGui=True, saveTo=None, showChain="worm", chains=None, support="", **kwargs): """Shows a single rainbow-colored chain using PyMOL. Arguments: @@ -600,11 +584,7 @@ def show_chain( tmpScript.flush() - os.system( - "pymol {0} {1} -u {2}".format( - tmpPdbPath, "" if showGui else "-c", tmpScript.name - ) - ) + os.system("pymol {0} {1} -u {2}".format(tmpPdbPath, "" if showGui else "-c", tmpScript.name)) tmpScript.close() @@ -653,9 +633,7 @@ def makeMoviePymol( rotationCode = "" if rotationPeriod > 0: for i in range(numFrames // rotationPeriod + 1): - rotationCode += "util.mroll {0},{1},0\n".format( - i * rotationPeriod + 1, (i + 1) * rotationPeriod - ) + rotationCode += "util.mroll {0},{1},0\n".format(i * rotationPeriod + 1, (i + 1) * rotationPeriod) if pymolScript == None: script += textwrap.dedent( @@ -701,9 +679,7 @@ def makeMoviePymol( os.system("cd {0}; pymol -c -u {1}; cd -".format(imgFolder, tmpScriptPath)) _mencoder(imgFolder, fps, aviFilename) - shutil.move( - os.path.join(imgFolder, aviFilename), os.path.join(destFolder, aviFilename) - ) + shutil.move(os.path.join(imgFolder, aviFilename), os.path.join(destFolder, aviFilename)) def _mencoder(imgFolder, fps, aviFilename): diff --git a/polychrom/simulation.py b/polychrom/simulation.py index 517d865..977a3b6 100644 --- a/polychrom/simulation.py +++ b/polychrom/simulation.py @@ -208,7 +208,7 @@ def __init__(self, **kwargs): Shout out loud about every change. precision: str, optional (not recommended to change) - single is the default now, because mixed is much slower on 3080 and other new GPUs + single is the default now, because mixed is much slower on 3080 and other new GPUs If you are using double precision, it will be slower by a factor of 10 or so. save_decimals: int or False, optional @@ -244,16 +244,10 @@ def __init__(self, **kwargs): ] for i in kwargs.keys(): if i not in valid_names: - raise ValueError( - "incorrect argument provided: {0}. Allowed are {1}".format( - i, valid_names - ) - ) + raise ValueError("incorrect argument provided: {0}. Allowed are {1}".format(i, valid_names)) if None in kwargs.values(): - raise ValueError( - "None is not allowed in arguments due to HDF5 incompatiliblity. Use False instead." - ) + raise ValueError("None is not allowed in arguments due to HDF5 incompatiliblity. Use False instead.") default_args.update(kwargs) kwargs = default_args self.kwargs = kwargs @@ -306,9 +300,7 @@ def __init__(self, **kwargs): kwargs["timestep"] * simtk.unit.femtosecond, ) elif self.integrator_type.lower() == "verlet": - self.integrator = openmm.VariableVerletIntegrator( - kwargs["timestep"] * simtk.unit.femtosecond - ) + self.integrator = openmm.VariableVerletIntegrator(kwargs["timestep"] * simtk.unit.femtosecond) elif self.integrator_type.lower() == "variableverlet": self.integrator = openmm.VariableVerletIntegrator(kwargs["error_tol"]) @@ -347,8 +339,7 @@ def __init__(self, **kwargs): self.conlen = 1.0 * simtk.unit.nanometer * self.length_scale self.kbondScalingFactor = float( - (2 * self.kT / self.conlen**2) - / (simtk.unit.kilojoule_per_mole / simtk.unit.nanometer**2) + (2 * self.kT / self.conlen**2) / (simtk.unit.kilojoule_per_mole / simtk.unit.nanometer**2) ) self.system = openmm.System() @@ -413,11 +404,7 @@ def set_data(self, data, center=False, random_offset=1e-5, report=True): raise ValueError(f"length of data, {len(data)} does not match N, {self.N}") if data.shape[1] != 3: - raise ValueError( - "Data is not shaped correctly. Needs (N,3), provided: {0}".format( - data.shape - ) - ) + raise ValueError("Data is not shaped correctly. Needs (N,3), provided: {0}".format(data.shape)) if np.isnan(data).any(): raise ValueError("Data contains NANs") @@ -453,21 +440,13 @@ def set_velocities(self, v): v = np.asarray(v, dtype="float") if len(v) != self.N: - raise ValueError( - f"length of velocity array, {len(v)} does not match N, {self.N}" - ) + raise ValueError(f"length of velocity array, {len(v)} does not match N, {self.N}") if v.shape[1] != 3: - raise ValueError( - "Data is not shaped correctly. Needs (N,3), provided: {0}".format( - v.shape - ) - ) + raise ValueError("Data is not shaped correctly. Needs (N,3), provided: {0}".format(v.shape)) if np.isnan(v).any(): raise ValueError("Data contains NANs") - self.velocities = simtk.unit.Quantity( - v, simtk.unit.nanometer / simtk.unit.picosecond - ) + self.velocities = simtk.unit.Quantity(v, simtk.unit.nanometer / simtk.unit.picosecond) if hasattr(self, "context"): self.init_velocities() @@ -501,9 +480,7 @@ def add_force(self, force): self.add_force(f) else: if force.name in self.force_dict: - raise ValueError( - "A force named {} was added to the system twice!".format(force.name) - ) + raise ValueError("A force named {} was added to the system twice!".format(force.name)) forces._prepend_force_name_to_params(force) self.force_dict[force.name] = force @@ -545,9 +522,7 @@ def _apply_forces(self): else: force.setNonbondedMethod(force.CutoffNonPeriodic) - logging.info( - "adding force {} {}".format(i, self.system.addForce(self.force_dict[i])) - ) + logging.info("adding force {} {}".format(i, self.system.addForce(self.force_dict[i]))) for reporter in self.reporters: reporter.report( @@ -555,9 +530,7 @@ def _apply_forces(self): {i: j.__getstate__() for i, j in self.force_dict.items()}, ) - self.context = openmm.Context( - self.system, self.integrator, self.platform, self.properties - ) + self.context = openmm.Context(self.system, self.integrator, self.platform, self.properties) self.init_positions() self.init_velocities() self.forces_applied = True @@ -572,9 +545,7 @@ def initialize(self, **kwargs): self.masses = np.zeros(self.N, dtype=float) + self.kwargs["mass"] for mass in self.masses: self.system.addParticle(mass) - self.context = openmm.Context( - self.system, self.integrator, self.platform, self.properties - ) + self.context = openmm.Context(self.system, self.integrator, self.platform, self.properties) self.init_positions() self.init_velocities(**kwargs) @@ -591,9 +562,7 @@ def init_velocities(self, temperature="current"): try: self.context except: - raise ValueError( - "No context, cannot set velocs." "Initialize context before that" - ) + raise ValueError("No context, cannot set velocs." "Initialize context before that") if hasattr(self, "velocities"): self.context.setVelocities(self.velocities) @@ -611,15 +580,9 @@ def init_positions(self): try: self.context except: - raise ValueError( - "No context, cannot set positions." " Initialize context before that" - ) + raise ValueError("No context, cannot set positions." " Initialize context before that") self.context.setPositions(self.data) - eP = ( - self.context.getState(getEnergy=True).getPotentialEnergy() - / self.N - / self.kT - ) + eP = self.context.getState(getEnergy=True).getPotentialEnergy() / self.N / self.kT logging.info("Particles loaded. Potential energy is %lf" % eP) def reinitialize(self): @@ -632,9 +595,7 @@ def reinitialize(self): self.init_positions() self.init_velocities() - def local_energy_minimization( - self, tolerance=0.3, maxIterations=0, random_offset=0.02 - ): + def local_energy_minimization(self, tolerance=0.3, maxIterations=0, random_offset=0.02): """ A wrapper to the build-in OpenMM Local Energy Minimization @@ -694,9 +655,7 @@ def local_energy_minimization( eK = self.state.getKineticEnergy() / self.N / self.kT eP = self.state.getPotentialEnergy() / self.N / self.kT locTime = self.state.getTime() - logging.info( - "before minimization eK={0}, eP={1}, time={2}".format(eK, eP, locTime) - ) + logging.info("before minimization eK={0}, eP={1}, time={2}".format(eK, eP, locTime)) openmm.LocalEnergyMinimizer.minimize(self.context, tolerance, maxIterations) @@ -715,9 +674,7 @@ def local_energy_minimization( locTime = self.state.getTime() - logging.info( - "after minimization eK={0}, eP={1}, time={2}".format(eK, eP, locTime) - ) + logging.info("after minimization eK={0}, eP={1}, time={2}".format(eK, eP, locTime)) def do_block( self, @@ -749,9 +706,7 @@ def do_block( a = time.time() self.integrator.step(steps) # integrate! - self.state = self.context.getState( - getPositions=True, getVelocities=get_velocities, getEnergy=True - ) + self.state = self.context.getState(getPositions=True, getVelocities=get_velocities, getEnergy=True) b = time.time() coords = self.state.getPositions(asNumpy=True) newcoords = coords / simtk.unit.nanometer @@ -791,10 +746,7 @@ def do_block( msg += "Rg=%.3lf " % self.RG() msg += "SPS=%.0lf " % (steps / (float(b - a))) - if ( - self.integrator_type.lower() == "variablelangevin" - or self.integrator_type.lower() == "variableverlet" - ): + if self.integrator_type.lower() == "variablelangevin" or self.integrator_type.lower() == "variableverlet": dt = self.integrator.getStepSize() msg += "dt=%.1lffs " % (dt / simtk.unit.femtosecond) mass = self.system.getParticleMass(1) @@ -811,9 +763,7 @@ def do_block( "block": self.block, } if get_velocities: - result["vel"] = self.state.getVelocities() / ( - simtk.unit.nanometer / simtk.unit.picosecond - ) + result["vel"] = self.state.getVelocities() / (simtk.unit.nanometer / simtk.unit.picosecond) result.update(save_extras) if save: for reporter in self.reporters: @@ -828,9 +778,7 @@ def print_stats(self): """Prints detailed statistics of a system. Will be run every 50 steps """ - state = self.context.getState( - getPositions=True, getVelocities=True, getEnergy=True - ) + state = self.context.getState(getPositions=True, getVelocities=True, getEnergy=True) eP = state.getPotentialEnergy() pos = np.array(state.getPositions() / simtk.unit.nanometer) @@ -874,9 +822,7 @@ def print_stats(self): print(" z: %.2lf, %.2lf, %.2lf, %.2lf" % minmedmax(z)) print() print("Statistics for velocities:") - print( - " mean kinetic energy is: ", np.mean(EkPerParticle), "should be:", 1.5 - ) + print(" mean kinetic energy is: ", np.mean(EkPerParticle), "should be:", 1.5) print(" fastest particles are (in kT): ", np.sort(EkPerParticle)[-5:]) print() @@ -902,9 +848,7 @@ def show(self, shifts=[0.0, 0.2, 0.4, 0.6, 0.8], scale="auto"): return # determining the 95 percentile distance between particles, if scale == "auto": - meandist = np.percentile( - np.sqrt(np.sum(np.diff(data, axis=0) ** 2, axis=1)), 95 - ) + meandist = np.percentile(np.sqrt(np.sum(np.diff(data, axis=0) ** 2, axis=1)), 95) # rescaling the data, so that bonds are of the order of 1. # This is because rasmol spheres are of the fixed diameter. data /= meandist @@ -919,10 +863,7 @@ def show(self, shifts=[0.0, 0.2, 0.4, 0.6, 0.8], scale="auto"): if dist < 1.3: count += 1 if count > 100: - raise RuntimeError( - "Too many particles are close together. " - "This will cause rasmol to choke" - ) + raise RuntimeError("Too many particles are close together. " "This will cause rasmol to choke") rascript = tempfile.NamedTemporaryFile() # writing the rasmol script. Spacefill controls radius of the sphere. @@ -937,16 +878,12 @@ def show(self, shifts=[0.0, 0.2, 0.4, 0.6, 0.8], scale="auto"): # creating the array, linearly chanhing from -225 to 225 # to serve as an array of colors - colors = np.array( - [int((j * 450.0) / (len(data))) - 225 for j in range(len(data))] - ) + colors = np.array([int((j * 450.0) / (len(data))) - 225 for j in range(len(data))]) # creating spheres along the trajectory newData = np.zeros((len(data) * len(shifts) - (len(shifts) - 1), 4)) for i in range(len(shifts)): - newData[i : -1 : len(shifts), :3] = data[:-1] * shifts[i] + data[1:] * ( - 1 - shifts[i] - ) + newData[i : -1 : len(shifts), :3] = data[:-1] * shifts[i] + data[1:] * (1 - shifts[i]) newData[i : -1 : len(shifts), 3] = colors[:-1] newData[-1, :3] = data[-1] newData[-1, 3] = colors[-1] @@ -956,11 +893,7 @@ def show(self, shifts=[0.0, 0.2, 0.4, 0.6, 0.8], scale="auto"): # number of atoms and a blank line after is a requirement of rasmol for i in newData: - towrite.write( - ( - "CA\t{:f}\t{:f}\t{:f}\t{:d}\n".format(i[0], i[1], i[2], int(i[3])) - ).encode("utf-8") - ) + towrite.write(("CA\t{:f}\t{:f}\t{:f}\t{:d}\n".format(i[0], i[1], i[2], int(i[3]))).encode("utf-8")) towrite.flush() "TODO: rewrite using subprocess.popen" @@ -968,10 +901,7 @@ def show(self, shifts=[0.0, 0.2, 0.4, 0.6, 0.8], scale="auto"): if os.name == "posix": # if linux os.system("rasmol -xyz %s -script %s" % (towrite.name, rascript.name)) else: # if windows - os.system( - "C:/RasWin/raswin.exe -xyz %s -script %s" - % (towrite.name, rascript.name) - ) + os.system("C:/RasWin/raswin.exe -xyz %s -script %s" % (towrite.name, rascript.name)) rascript.close() towrite.close() diff --git a/polychrom/starting_conformations.py b/polychrom/starting_conformations.py index fb9cbea..dffe2c5 100644 --- a/polychrom/starting_conformations.py +++ b/polychrom/starting_conformations.py @@ -177,10 +177,7 @@ def create_constrained_random_walk( vec_to_rot = d[j] rot_axis = np.cross(past_displacement, np.array([0, 0, 1])) rot_axis = rot_axis / np.linalg.norm(rot_axis) - rot_angle = -np.arccos( - np.dot(past_displacement, np.array([0, 0, 1])) - / np.linalg.norm(past_displacement) - ) + rot_angle = -np.arccos(np.dot(past_displacement, np.array([0, 0, 1])) / np.linalg.norm(past_displacement)) # Rotating with the Rodriques' rotation formula # https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula next_displacement = ( diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..85c3b07 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,5 @@ +[tool.black] +line-length = 120 + +[tool.isort] +profile = "black" diff --git a/tests/test_configurations.py b/tests/test_configurations.py index 88c37cf..2f96bbb 100644 --- a/tests/test_configurations.py +++ b/tests/test_configurations.py @@ -3,9 +3,7 @@ from polychrom import starting_conformations -def Jun_30_2022_create_constrained_random_walk( - N, constraint_f, starting_point=(0, 0, 0), step_size=1.0 -): +def Jun_30_2022_create_constrained_random_walk(N, constraint_f, starting_point=(0, 0, 0), step_size=1.0): """ Creates a constrained freely joined chain of length N with step step_size. Each step of a random walk is tested with the constraint function and is @@ -58,9 +56,7 @@ def test_for_correct_run_all_allowed(self): def alwaystrue(new_p): return True - polymer = starting_conformations.create_constrained_random_walk( - N, alwaystrue, step_size=step - ) + polymer = starting_conformations.create_constrained_random_walk(N, alwaystrue, step_size=step) increments = [polymer[i + 1] - polymer[i] for i in range(len(polymer) - 1)] self.assertTrue( @@ -110,14 +106,9 @@ def test_for_correct_angle_fixing_always_true(self): def alwaystrue(new_p): return True - polymer = starting_conformations.create_constrained_random_walk( - N, alwaystrue, polar_fixed=polar_fixed - ) + polymer = starting_conformations.create_constrained_random_walk(N, alwaystrue, polar_fixed=polar_fixed) angles = np.arccos( - [ - np.dot(polymer[i + 2] - polymer[i + 1], polymer[i + 1] - polymer[i]) - for i in range(len(polymer) - 2) - ] + [np.dot(polymer[i + 2] - polymer[i + 1], polymer[i + 1] - polymer[i]) for i in range(len(polymer) - 2)] ) # Check that the angles are correct @@ -127,14 +118,9 @@ def test_for_correct_angle_fixing_with_confinement(self): def confined(new_p): return np.linalg.norm(new_p) < confinement - polymer = starting_conformations.create_constrained_random_walk( - N, confined, polar_fixed=polar_fixed - ) + polymer = starting_conformations.create_constrained_random_walk(N, confined, polar_fixed=polar_fixed) angles = np.arccos( - [ - np.dot(polymer[i + 2] - polymer[i + 1], polymer[i + 1] - polymer[i]) - for i in range(len(polymer) - 2) - ] + [np.dot(polymer[i + 2] - polymer[i + 1], polymer[i + 1] - polymer[i]) for i in range(len(polymer) - 2)] ) # Check that the angles are correct diff --git a/tests/test_contactmaps.py b/tests/test_contactmaps.py index 9ff8225..090809d 100644 --- a/tests/test_contactmaps.py +++ b/tests/test_contactmaps.py @@ -52,9 +52,7 @@ def test_contactmaps(): assert np.allclose(cmap1, cmap3) for n in [1, 4]: - cmap6 = monomerResolutionContactMap( - range(8), cutoff=1, loadFunction=lambda x: ars[x], n=n - ) + cmap6 = monomerResolutionContactMap(range(8), cutoff=1, loadFunction=lambda x: ars[x], n=n) cmap5 = cmapPureMap( range(8), cutoff=1, diff --git a/tests/test_forces.py b/tests/test_forces.py index a6e1b6b..e0cd0cb 100644 --- a/tests/test_forces.py +++ b/tests/test_forces.py @@ -30,9 +30,7 @@ def test_harmonic_bond_force(): data = np.concatenate([particles, particles2]) - sim.set_data( - data, center=False - ) # loads a polymer, puts a center of mass at zero + sim.set_data(data, center=False) # loads a polymer, puts a center of mass at zero sim.add_force(forces.harmonic_bonds(sim, bonds, bondLength=dist)) res = sim.do_block(0) pot = res["potentialEnergy"] * 2 # only one bond per pair diff --git a/tests/test_polymer_analyses.py b/tests/test_polymer_analyses.py index 1be6087..3c7cbda 100644 --- a/tests/test_polymer_analyses.py +++ b/tests/test_polymer_analyses.py @@ -27,13 +27,9 @@ def _testMutualSimplify(): b = np.dot(b, mat) a = a + np.random.random(a.shape) * 0.0001 b = b + np.random.random(b.shape) * 0.0001 - c1 = polychrom.polymer_analyses.getLinkingNumber( - a, b, simplify=False, randomOffset=False - ) + c1 = polychrom.polymer_analyses.getLinkingNumber(a, b, simplify=False, randomOffset=False) a, b = polychrom.polymer_analyses.mutualSimplify(a, b, verbose=False) - c2 = polychrom.polymer_analyses.getLinkingNumber( - a, b, simplify=False, randomOffset=False - ) + c2 = polychrom.polymer_analyses.getLinkingNumber(a, b, simplify=False, randomOffset=False) print("simplified from 200 to {0} and {1}".format(len(a), len(b))) print("Link before: {0}, link after: {1}".format(c1, c2)) assert c1 == c2 @@ -48,13 +44,9 @@ def _testMutualSimplify(): a = a + np.random.random(a.shape) * 0.0001 b = b + np.random.random(b.shape) * 0.0001 - c1 = polychrom.polymer_analyses.getLinkingNumber( - a, b, simplify=False, randomOffset=False - ) + c1 = polychrom.polymer_analyses.getLinkingNumber(a, b, simplify=False, randomOffset=False) a, b = polychrom.polymer_analyses.mutualSimplify(a, b, verbose=False) - c2 = polychrom.polymer_analyses.getLinkingNumber( - a, b, simplify=False, randomOffset=False - ) + c2 = polychrom.polymer_analyses.getLinkingNumber(a, b, simplify=False, randomOffset=False) print("simplified from 3000 and 1000 to {0} and {1}".format(len(a), len(b))) print("Link before: {0}, link after: {1}".format(c1, c2)) assert c1 == c2 @@ -64,9 +56,7 @@ def test_scalings(): import polychrom.starting_conformations, polychrom.polymer_analyses import numpy as np - datas = [ - polychrom.starting_conformations.create_random_walk(1, 80) for _ in range(100) - ] + datas = [polychrom.starting_conformations.create_random_walk(1, 80) for _ in range(100)] scals = [polychrom.polymer_analyses.R2_scaling(i) for i in datas] scals = np.mean(scals, axis=0) diff --git a/utilities/quickMovie.py b/utilities/quickMovie.py index 235a9ea..d1061f4 100644 --- a/utilities/quickMovie.py +++ b/utilities/quickMovie.py @@ -17,10 +17,7 @@ trajectoryPath = sys.argv[1] -files = [ - os.path.join(trajectoryPath, "block{0}.dat".format(i)) - for i in range(int(sys.argv[3]), int(sys.argv[4])) -] +files = [os.path.join(trajectoryPath, "block{0}.dat".format(i)) for i in range(int(sys.argv[3]), int(sys.argv[4]))] pymol_show.makeMoviePymol( files, diff --git a/utilities/showChainWithRasmol.py b/utilities/showChainWithRasmol.py index 90b10bc..7461e32 100644 --- a/utilities/showChainWithRasmol.py +++ b/utilities/showChainWithRasmol.py @@ -72,24 +72,19 @@ def convertData(data, colors): for i in range(len(shifts)): # filling in the array like 0,5,10,15; then 1,6,11,16; then 2,7,12,17, etc. # this is just very fast - newData[i : -1 : len(shifts), :3] = data[:-1] * shifts[i] + data[1:] * ( - 1 - shifts[i] - ) + newData[i : -1 : len(shifts), :3] = data[:-1] * shifts[i] + data[1:] * (1 - shifts[i]) newData[i : -1 : len(shifts), 3] = colors[:-1] newData[-1, :3] = data[-1] newData[-1, 3] = colors[-1] return newData newDatas = [ - convertData(data[breaks[i] : breaks[i + 1]], colors[breaks[i] : breaks[i + 1]]) - for i in range(len(breaks) - 1) + convertData(data[breaks[i] : breaks[i + 1]], colors[breaks[i] : breaks[i + 1]]) for i in range(len(breaks) - 1) ] newData = np.concatenate(newDatas) towrite = tempfile.NamedTemporaryFile(mode="w") - towrite.write( - "%d\n\n" % (len(newData)) - ) # number of atoms and a blank line after is a requirement of rasmol + towrite.write("%d\n\n" % (len(newData))) # number of atoms and a blank line after is a requirement of rasmol for i in newData: towrite.write("CA\t%lf\t%lf\t%lf\t%d\n" % tuple(i)) @@ -98,9 +93,7 @@ def convertData(data, colors): if os.name == "posix": # if linux os.system("rasmol -xyz %s -script %s" % (towrite.name, rascript.name)) else: # if windows - os.system( - "C:/RasWin/raswin.exe -xyz %s -script %s" % (towrite.name, rascript.name) - ) + os.system("C:/RasWin/raswin.exe -xyz %s -script %s" % (towrite.name, rascript.name)) exit()