From 8066aee98df99e10ee635292ae47d0ea63c07d7f Mon Sep 17 00:00:00 2001 From: Kalistyn Burley Date: Tue, 12 Feb 2019 12:11:47 -0800 Subject: [PATCH 01/11] working on up-to-date bias branch --- blues/moves.py | 291 +++++++++++++++++++++++++++++---------- blues/simulation.py | 14 +- blues/tests/rotpref.yaml | 58 ++++++++ 3 files changed, 283 insertions(+), 80 deletions(-) create mode 100644 blues/tests/rotpref.yaml diff --git a/blues/moves.py b/blues/moves.py index 7322173f..4c13b5dd 100644 --- a/blues/moves.py +++ b/blues/moves.py @@ -21,6 +21,13 @@ import parmed from simtk import unit +import yaml +from mdtraj.geometry.dihedral import (indices_chi2, + indices_chi3, + indices_chi1, + indices_chi4, + indices_chi5) + try: import openeye.oechem as oechem if not oechem.OEChemIsLicensed(): @@ -198,6 +205,10 @@ def __init__(self, structure, resname='LIG', random_state=None): self.center_of_mass = None self.positions = structure[self.atom_indices].positions + ## adding placeholder values for random ligand moves (always true for these moves) vs biased sidechain KB + self.make_NCMC_move = True + self.bin_boolean = True + self._calculateProperties() def getAtomIndices(self, structure, resname): @@ -462,10 +473,11 @@ class SideChainMove(Move): """ - def __init__(self, structure, residue_list, verbose=False, write_move=False): + def __init__(self, structure, residue_list, bias_range=15, verbose=False, write_move=False): self.structure = structure self.molecule = self._pmdStructureToOEMol() self.residue_list = residue_list + self.bias_range = bias_range self.all_atoms = [atom.index for atom in self.structure.topology.atoms()] self.rot_atoms, self.rot_bonds, self.qry_atoms = self.getRotBondAtoms() self.atom_indices = self.rot_atoms @@ -535,8 +547,6 @@ def getTargetAtoms(self, molecule, backbone_atoms, residue_list): qry_atoms = {} qry_atoms.clear() - reslib = [] - #print('Searching residue list for atoms...') # loop through all the atoms in the PDB OEGraphMol structure for atom in molecule.GetAtoms(): @@ -549,8 +559,6 @@ def getTargetAtoms(self, molecule, backbone_atoms, residue_list): # store the atom location in a query atom dict keyed by its atom index qry_atoms.update({atom: atom.GetIdx()}) #print('Found atom %s in residue number %i %s'%(atom,myres.GetResidueNumber(),myres.GetName())) - if myres not in reslib: - reslib.append(myres) return qry_atoms, backbone_atoms @@ -626,33 +634,41 @@ def getRotAtoms(self, rotbonds, molecule, backbone_atoms): idx_list.clear() query_list.clear() resnum = (rotbonds[bond]) - thisbond = bond ax1 = bond.GetBgn() ax2 = bond.GetEnd() if resnum in rot_atom_dict.keys(): - rot_atom_dict[resnum].update({thisbond: []}) + chi = len(rot_atom_dict[resnum]['chis'].keys())+1 + rot_atom_dict[resnum]['chis'].update({chi:{'bond_ptr':bond,'atms2mv':[],\ + 'dihed_atms':[],'bin_pref':[]}}) else: - rot_atom_dict.update({resnum: {thisbond: []}}) - - idx_list.append(ax1.GetIdx()) - idx_list.append(ax2.GetIdx()) + chi = 1 + res = oechem.OEAtomGetResidue(ax2) + residue_name = res.GetName() + rot_atom_dict.update({resnum : {'res_name': residue_name,'chis':{chi:\ + {'bond_ptr':bond,'atms2mv':[],'dihed_atms':[],'bin_pref':[]}}}}) if ax1 not in query_list and ax1.GetIdx() not in backbone_atoms: query_list.append(ax1) if ax2 not in query_list and ax2.GetIdx() not in backbone_atoms: query_list.append(ax2) + if ax2.GetIdx() > ax1.GetIdx(): + idx_list.append(ax1.GetIdx()) + idx_list.append(ax2.GetIdx()) + else: + idx_list.append(ax2.GetIdx()) + idx_list.append(ax1.GetIdx()) + for atom in query_list: checklist = atom.GetAtoms() for candidate in checklist: - if candidate not in query_list and candidate.GetIdx() not in backbone and candidate != ax2: + if candidate not in query_list and candidate.GetIdx() not in backbone: query_list.append(candidate) if candidate.GetAtomicNum() > 1: can_nbors = candidate.GetAtoms() for can_nbor in can_nbors: - if can_nbor not in query_list and candidate.GetIdx( - ) not in backbone and candidate != ax2: + if can_nbor not in query_list and candidate.GetIdx() not in backbone: query_list.append(can_nbor) for atm in query_list: @@ -660,11 +676,58 @@ def getRotAtoms(self, rotbonds, molecule, backbone_atoms): if y not in idx_list: idx_list.append(y) - rot_atom_dict[resnum].update({thisbond: list(idx_list)}) - #print("Moving these atoms:", idx_list) + rot_atom_dict[resnum]['chis'][chi]['atms2mv'] = list(idx_list) return rot_atom_dict + def getDihedral(self, positions, atomlist): + """This function computes the dihedral angle for the given atoms at the given positions""" + top = self.structure.topology + traj = mdtraj.Trajectory(np.asarray(positions),top) + angle = mdtraj.compute_dihedrals(traj, atomlist) + return angle + + def normalize_angles(self,ang_rad): + '''This function takes in an angle in radians, wraps it so it falls between 0 and 2pi and returns the + angle.''' + normalized_angle = (ang_rad+2*math.pi)%(2*math.pi) + return(normalized_angle) + + def is_in_bin(self,test_angle,bin_center): + '''This function takes in a test_angle, a bin_center and a bias_range and checks if the test_angle + falls within the bounds of the bin_center +/- the bias_window. All values are input in degrees.''' + #define upper and lower bin boundaries + lower = bin_center/180*math.pi - self.bias_range/180*math.pi + upper = bin_center/180*math.pi + self.bias_range/180*math.pi + + #normalize all input values to 0 to 2pi + n = self.normalize_angles(test_angle) + upper = self.normalize_angles(upper) + lower = self.normalize_angles(lower) + #if verbose: print(n,upper,lower) + + # check if value in bin + if n <= upper and n > lower: + outcome = True + else: outcome = False + #if verbose: print(outcome) + return outcome + + def getDihedralIndices(self,axis1,axis2): + """This function takes in a PDB filename (as a string) and list of residue numbers. It returns + a nested dictionary of rotatable bonds (containing only heavy atoms), that are keyed by residue number, + then keyed by bond pointer, containing values of atom indicies [axis1, axis2, atoms to be rotated] + **Note: The atom indicies start at 0, and are offset by -1 from the PDB file indicies""" + top = mdtraj.Topology.from_openmm(self.structure.topology) + indices = indices_chi1(top).tolist()+indices_chi2(top).tolist()+indices_chi3(top).tolist()\ + +indices_chi4(top).tolist()+indices_chi5(top).tolist() + dihed = [] + for list in indices: + if axis1==list[1] and axis2==list[2]: + dihed = list + break + return dihed + def getRotBondAtoms(self): """This function is called on class initialization. @@ -695,35 +758,71 @@ def getRotBondAtoms(self): # Generate dictionary of residues, bonds and atoms to be rotated rot_atoms = self.getRotAtoms(rot_bonds, self.molecule, backbone_atoms) - return rot_atoms, rot_bonds, qry_atoms - - def chooseBondandTheta(self): - """This function is called on class initialization. - Takes a dictionary containing nested dictionary, keyed by res#, - then keyed by bond_ptrs, containing a list of atoms to move, randomly selects a bond, - and generates a random angle (radians). It returns the atoms associated with the - the selected bond, the pointer for the selected bond and the randomly generated angle + # Read in yaml file + rot_pref = yaml.load(open('rotpref.yaml')) + # Restructure dictionary and populate atms2mv and bin_pref + for residx in rot_atoms.keys(): + resname = rot_atoms[residx]['res_name'] + for chi in rot_atoms[residx]['chis']: + atom1=rot_atoms[residx]['chis'][chi]['atms2mv'][0] + atom2=rot_atoms[residx]['chis'][chi]['atms2mv'][1] + print(atom1,atom2) + atoms = self.getDihedralIndices(atom1,atom2) + if chi == 5 and resname == 'ARG': + print('Arginine chi5 detected. Atom indices assigned based on numbers') + atoms = [atom1-3, atom1, atom2, atom2+1] + rot_atoms[residx]['chis'][chi]['dihed_atms'] = atoms + rot_atoms[residx]['chis'][chi]['bin_pref'] = rot_pref[resname][chi] + return rot_atoms, rot_bonds, qry_atoms - Returns - ------- - theta_ran : - - targetatoms : - - res_choice : + def chooseBond(self,positions): + bonds_in_bins = [] - bond_choice : + for residx in self.rot_atoms.keys(): + resname = self.rot_atoms[residx]['res_name'] + for chi in self.rot_atoms[residx]['chis']: + dihed_atoms = [self.rot_atoms[residx]['chis'][chi]['dihed_atms']] + print(dihed_atoms) + curr_angle = self.getDihedral(positions,dihed_atoms) + bin_ct=0 + for bin in self.rot_atoms[residx]['chis'][chi]['bin_pref']: + bin_ct+=1 + if self.is_in_bin(curr_angle,bin): + bin_idx = bin_ct-1 + bonds_in_bins.append([(self.rot_atoms[residx]['chis'][chi]),curr_angle[0][0],bin_idx]) + break - """ + print("\nThere are %i bonds that fall within the biasing bins: %s\n"%(len(bonds_in_bins),bonds_in_bins)) - res_choice = random.choice(list(self.rot_atoms.keys())) - bond_choice = random.choice(list(self.rot_atoms[res_choice].keys())) - targetatoms = self.rot_atoms[res_choice][bond_choice] - theta_ran = random.random() * 2 * math.pi - return theta_ran, targetatoms, res_choice, bond_choice + try: + bond_choice = random.choice(bonds_in_bins) + print("The bond randomly chosen to rotate is:\n",bond_choice) + except: + print("No rotatable bonds are within range for rotation for biased NCMC sidechain\ + \move proposal. Proceed with additional round of MD.") + bond_choice = False + + return bond_choice + + def chooseTheta(self,bond_choice): + moveok = False + while moveok == False: + theta_ran = random.random()*2*math.pi-math.pi + new_ang = theta_ran + bond_choice[1] + for bin in bond_choice[0]['bin_pref']: + if bin == bond_choice[0]['bin_pref'][bond_choice[2]]: + continue + elif self.is_in_bin(new_ang,bin): + conv = new_ang*180/math.pi + print("Your theta is %.2f and your new angle is %.2f or %.2f."%(theta_ran,new_ang,conv)) + print("You started in bin [%i] and are proposing move to bin [%i]"%(bond_choice[0]['bin_pref'][bond_choice[2]],bin)) + new_bin = bin + moveok = True + break + return theta_ran, new_bin def rotation_matrix(self, axis, theta): """Function returns the rotation matrix associated with counterclockwise rotation @@ -746,6 +845,12 @@ def rotation_matrix(self, axis, theta): 2 * (bd - ac)], [2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)], [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]]) + def beforeMove(self, context): + + self.start_pos = context.getState(getPositions=True).getPositions(asNumpy=True) + + return context + def move(self, context, verbose=False): """Rotates the target atoms around a selected bond by angle theta and updates the atom coordinates in the parmed structure as well as the ncmc context object @@ -765,10 +870,6 @@ def move(self, context, verbose=False): """ - # determine the axis, theta, residue, and bond + atoms to be rotated - theta, target_atoms, res, bond = self.chooseBondandTheta() - print('Rotating bond: %s in resnum: %s by %.2f radians' % (bond, res, theta)) - #retrieve the current positions initial_positions = context.getState(getPositions=True).getPositions(asNumpy=True) nc_positions = copy.deepcopy(initial_positions) @@ -792,55 +893,93 @@ def move(self, context, verbose=False): positions = model.positions - # find the rotation axis using the updated positions - axis1 = target_atoms[0] - axis2 = target_atoms[1] - rot_axis = (positions[axis1] - positions[axis2]) / positions.unit + #Pick a bond to rotate (by identifying which bond dihedrals fall within biasing regions) + selected_bond = self.chooseBond(self.start_pos) - #calculate the rotation matrix - rot_matrix = self.rotation_matrix(rot_axis, theta) + if not selected_bond: + self.make_NCMC_move = False + else: + self.make_NCMC_move = True + new_theta, self.target_bin = self.chooseTheta(selected_bond) + self.dihed_atoms = [selected_bond[0]['dihed_atms']] + axis1 = self.dihed_atoms[0][1] + axis2 = self.dihed_atoms[0][2] - # apply the rotation matrix to the target atoms - for idx, atom in enumerate(target_atoms): + rot_axis = (positions[axis1] - positions[axis2])/positions.unit - my_position = positions[atom] + #Adjust theta to account for repositioning during NCMC relaxation + post_relax_dihed = self.getDihedral(initial_positions,self.dihed_atoms) + theta_adj = selected_bond[1] - post_relax_dihed + new_theta - if self.verbose: - print('The current position for %i is: %s' % (atom, my_position)) + #calculate the rotation matrix + rot_matrix = self.rotation_matrix(rot_axis, theta_adj) - # find the reduced position (substract out axis) - red_position = (my_position - model.positions[axis2])._value - # find the new positions by multiplying by rot matrix - new_position = numpy.dot(rot_matrix, red_position) * positions.unit + positions[axis2] + # apply the rotation matrix to the target atoms + for idx, atom in enumerate(target_atoms): - if self.verbose: print("The new position should be:", new_position) + my_position = positions[atom] - positions[atom] = new_position - # Update the parmed model with the new positions - model.atoms[atom].xx = new_position[0] / positions.unit - model.atoms[atom].xy = new_position[1] / positions.unit - model.atoms[atom].xz = new_position[2] / positions.unit + if self.verbose: + print('The current position for %i is: %s' % (atom, my_position)) - #update the copied ncmc context array with the new positions - nc_positions[atom][0] = model.atoms[atom].xx * nc_positions.unit / 10 - nc_positions[atom][1] = model.atoms[atom].xy * nc_positions.unit / 10 - nc_positions[atom][2] = model.atoms[atom].xz * nc_positions.unit / 10 + # find the reduced position (substract out axis) + red_position = (my_position - model.positions[axis2])._value + # find the new positions by multiplying by rot matrix + new_position = numpy.dot(rot_matrix, red_position) * positions.unit + positions[axis2] - if self.verbose: - print('The updated position for this atom is:', model.positions[atom]) + if self.verbose: print("The new position should be:", new_position) + + positions[atom] = new_position + # Update the parmed model with the new positions + model.atoms[atom].xx = new_position[0] / positions.unit + model.atoms[atom].xy = new_position[1] / positions.unit + model.atoms[atom].xz = new_position[2] / positions.unit + + #update the copied ncmc context array with the new positions + nc_positions[atom][0] = model.atoms[atom].xx * nc_positions.unit / 10 + nc_positions[atom][1] = model.atoms[atom].xy * nc_positions.unit / 10 + nc_positions[atom][2] = model.atoms[atom].xz * nc_positions.unit / 10 + + if self.verbose: + print('The updated position for this atom is:', model.positions[atom]) + + # update the actual ncmc context object with the new positions + context.setPositions(nc_positions) - # update the actual ncmc context object with the new positions - context.setPositions(nc_positions) + # update the class structure positions + self.structure.positions = model.positions - # update the class structure positions - self.structure.positions = model.positions + if self.write_move: + filename = 'sc_move_%s_%s_%s.pdb' % (res, axis1, axis2) + mod_prot = model.save(filename, overwrite=True) - if self.write_move: - filename = 'sc_move_%s_%s_%s.pdb' % (res, axis1, axis2) - mod_prot = model.save(filename, overwrite=True) return context + def afterMove(self, context): + """This method is called at the end of the NCMC portion if the + context needs to be checked or modified before performing the move + at the halfway point. + + Parameters + ---------- + context: simtk.openmm.Context object: + Context containing the positions to be moved. + + Returns + ------- + context: simtk.openmm.Context object + The same input context, but whose context were changed by this function. + + """ + post_pos = context.getState(getPositions=True).getPositions(asNumpy=True) + final_angle = self.getDihedral(post_pos,self.dihed_atoms) + if self.is_in_bin(final_angle,self.target_bin): + self.bin_boolean = True + else: self.bin_boolean = False + + return context + class SmartDartMove(RandomLigandRotationMove): """**WARNING:** This class has not been completely tested. Use at your own risk. diff --git a/blues/simulation.py b/blues/simulation.py index ddab123d..a9e41741 100644 --- a/blues/simulation.py +++ b/blues/simulation.py @@ -1068,7 +1068,8 @@ def _stepNCMC(self, nstepsNC, moveStep, move_engine=None): #Attempt anything related to the move before protocol is performed if not step: self._ncmc_sim.context = move_engine.selected_move.beforeMove(self._ncmc_sim.context) - + if self.move_engine.moves[self.move_engine.selected_move].make_NCMC_move == False: + break # Attempt selected MoveEngine Move at the halfway point #to ensure protocol is symmetric if step == moveStep: @@ -1134,7 +1135,11 @@ def _acceptRejectMove(self, write_move=False): 'NCMCLogAcceptanceProbability = %.6f + Alchemical Correction = %.6f' % (work_ncmc, correction_factor)) work_ncmc = work_ncmc + correction_factor - if work_ncmc > randnum: + #added for sidechain biasing + if self.move_engine.moves[self.move_engine.selected_move].bin_boolean == False: + print("Bin Boolean false: dihedral angle departed from target bin during NCMC move") + + if work_ncmc > randnum and self.move_engine.moves[self.move_engine.selected_move].bin_boolean: self.accept += 1 logger.info('NCMC MOVE ACCEPTED: work_ncmc {} > randnum {}'.format(work_ncmc, randnum)) @@ -1244,8 +1249,9 @@ def run(self, nIter=0, nstepsNC=0, moveStep=0, nstepsMD=0, temperature=300, writ logger.info('BLUES Iteration: %s' % N) self._syncStatesMDtoNCMC() self._stepNCMC(nstepsNC, moveStep) - self._acceptRejectMove(write_move) - self._resetSimulations(temperature) + if self.move_engine.moves[self.move_engine.selected_move].make_NCMC_move == True: + self._acceptRejectMove(write_move) + self._resetSimulations(temperature) self._stepMD(nstepsMD) # END OF NITER diff --git a/blues/tests/rotpref.yaml b/blues/tests/rotpref.yaml new file mode 100644 index 00000000..47b3353e --- /dev/null +++ b/blues/tests/rotpref.yaml @@ -0,0 +1,58 @@ +ARG: + 1: [-68, 56, -171] + 2: [-68, 67, 180] + 3: [-63, 68, -177] + 4: [-89, 86, -174] + 5: [-60, 60, 180] +ASN: + 1: [-68, 59, 180] + 2: [-107, 53] +ASP: + 1: [-75,74,-178] + 2: [107,62,-171] +CYS: + 1: [-68, 56, -171] +GLN: + 1: [-68, 56, -171] + 2: [-68, 56, -171] + 3: [-68, 56, -171] +GLU: + 1: [-68, 56, -171] + 2: [-68, 56, -171] + 3: [-68, 56, -171] +HIS: + 1: [-68, 56, -171] + 2: [-68, 56, -171] +ILE: + 1: [-68, 56, -171] + 2: [-68, 56, -171] +LEU: + 1: [-68, 56, -171] + 2: [-68, 56, -171] +LYS: + 1: [-68, 56, -171] + 2: [-68, 56, -171] + 3: [-68, 56, -171] + 4: [-68, 56, -171] + 5: [-68, 56, -171] +MET: + 1: [-60, 60, 180] + 2: [-60, 60, 180] + 3: [-70, 70, 175] +PHE: + 1: [-79, 52, -177] + 2: [-90,90] +PRO: + 1: [-68, 56, -171] +SER: + 1: [-68, 56, -171] +THR: + 1: [-60, 60, 180] +TRP: + 1: [-68, 56, -171] + 2: [-68, 56, -171] +TYR: + 1: [-68, 56, -171] + 2: [-68, 56, -171] +VAL: + 1: [-72,57,181] From e4e9ef190a6a3a1e080b3a1b26a69c070f0ad051 Mon Sep 17 00:00:00 2001 From: Kalistyn Burley Date: Tue, 12 Feb 2019 12:39:21 -0800 Subject: [PATCH 02/11] working version of bias on updated BLUES --- blues/moves.py | 6 +++--- blues/tests/test_sidechain.py | 6 ++++-- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/blues/moves.py b/blues/moves.py index 4c13b5dd..b7c07bee 100644 --- a/blues/moves.py +++ b/blues/moves.py @@ -683,7 +683,7 @@ def getRotAtoms(self, rotbonds, molecule, backbone_atoms): def getDihedral(self, positions, atomlist): """This function computes the dihedral angle for the given atoms at the given positions""" top = self.structure.topology - traj = mdtraj.Trajectory(np.asarray(positions),top) + traj = mdtraj.Trajectory(numpy.asarray(positions),top) angle = mdtraj.compute_dihedrals(traj, atomlist) return angle @@ -767,7 +767,6 @@ def getRotBondAtoms(self): for chi in rot_atoms[residx]['chis']: atom1=rot_atoms[residx]['chis'][chi]['atms2mv'][0] atom2=rot_atoms[residx]['chis'][chi]['atms2mv'][1] - print(atom1,atom2) atoms = self.getDihedralIndices(atom1,atom2) if chi == 5 and resname == 'ARG': print('Arginine chi5 detected. Atom indices assigned based on numbers') @@ -784,7 +783,6 @@ def chooseBond(self,positions): resname = self.rot_atoms[residx]['res_name'] for chi in self.rot_atoms[residx]['chis']: dihed_atoms = [self.rot_atoms[residx]['chis'][chi]['dihed_atms']] - print(dihed_atoms) curr_angle = self.getDihedral(positions,dihed_atoms) bin_ct=0 for bin in self.rot_atoms[residx]['chis'][chi]['bin_pref']: @@ -914,6 +912,8 @@ def move(self, context, verbose=False): #calculate the rotation matrix rot_matrix = self.rotation_matrix(rot_axis, theta_adj) + target_atoms = selected_bond[0]['atms2mv'] + # apply the rotation matrix to the target atoms for idx, atom in enumerate(target_atoms): diff --git a/blues/tests/test_sidechain.py b/blues/tests/test_sidechain.py index 75ea7146..83543083 100644 --- a/blues/tests/test_sidechain.py +++ b/blues/tests/test_sidechain.py @@ -60,13 +60,15 @@ def setUp(self): self.simulations = SimulationFactory(self.systems, self.engine, self.cfg) def test_getRotBondAtoms(self): - vals = [v for v in self.sidechain.rot_atoms[1].values()][0] + vals = [v for v in self.sidechain.rot_atoms[1]['chis'][1]['atms2mv']] assert len(vals) == 11 #Ensure it selects 1 rotatable bond in Valine assert len(self.sidechain.rot_bonds) == 1 def test_sidechain_move(self): - atom_indices = [v for v in self.sidechain.rot_atoms[1].values()][0] + atom_indices = [v for v in self.sidechain.rot_atoms[1]['chis'][1]['atms2mv']] + before_context = self.simulations.ncmc.context + dobefore = self.sidechain.beforeMove(before_context) before_move = self.simulations.ncmc.context.getState(getPositions=True).getPositions( asNumpy=True)[atom_indices, :] self.simulations.ncmc.context = self.engine.runEngine(self.simulations.ncmc.context) From 5ca0f1786f2de30effa68b938a82a73f1e7b7b3d Mon Sep 17 00:00:00 2001 From: Kalistyn Burley Date: Wed, 13 Feb 2019 12:32:37 -0800 Subject: [PATCH 03/11] Fixed issues with bin_boolean and make_NCMC_move variables --- blues/moves.py | 41 ++++++++++++++++++++++++----------------- blues/simulation.py | 9 +++++---- 2 files changed, 29 insertions(+), 21 deletions(-) diff --git a/blues/moves.py b/blues/moves.py index b7c07bee..644e1500 100644 --- a/blues/moves.py +++ b/blues/moves.py @@ -799,8 +799,7 @@ def chooseBond(self,positions): bond_choice = random.choice(bonds_in_bins) print("The bond randomly chosen to rotate is:\n",bond_choice) except: - print("No rotatable bonds are within range for rotation for biased NCMC sidechain\ - \move proposal. Proceed with additional round of MD.") + print("No rotatable bonds are within range for rotation for biased NCMC sidechain move proposal. Proceed with additional round of MD.") bond_choice = False return bond_choice @@ -846,10 +845,15 @@ def rotation_matrix(self, axis, theta): def beforeMove(self, context): self.start_pos = context.getState(getPositions=True).getPositions(asNumpy=True) + self.selected_bond = self.chooseBond(self.start_pos) + if self.selected_bond: + self.make_NCMC_move = True + else: + self.make_NCMC_move = False return context - def move(self, context, verbose=False): + def move(self, context, verbose=True): """Rotates the target atoms around a selected bond by angle theta and updates the atom coordinates in the parmed structure as well as the ncmc context object @@ -891,28 +895,27 @@ def move(self, context, verbose=False): positions = model.positions - #Pick a bond to rotate (by identifying which bond dihedrals fall within biasing regions) - selected_bond = self.chooseBond(self.start_pos) - - if not selected_bond: - self.make_NCMC_move = False - else: + if self.selected_bond: self.make_NCMC_move = True - new_theta, self.target_bin = self.chooseTheta(selected_bond) - self.dihed_atoms = [selected_bond[0]['dihed_atms']] + print("choosing new bond") + new_theta, self.target_bin = self.chooseTheta(self.selected_bond) + self.dihed_atoms = [self.selected_bond[0]['dihed_atms']] + print("Dihed atoms:",self.dihed_atoms) + print("Dihed before NCMC relax:",self.selected_bond[1]) axis1 = self.dihed_atoms[0][1] axis2 = self.dihed_atoms[0][2] - + print("Rot Axis:",axis1,axis2) rot_axis = (positions[axis1] - positions[axis2])/positions.unit #Adjust theta to account for repositioning during NCMC relaxation post_relax_dihed = self.getDihedral(initial_positions,self.dihed_atoms) - theta_adj = selected_bond[1] - post_relax_dihed + new_theta - + theta_adj = self.selected_bond[1] - post_relax_dihed + new_theta + print("Dihed at halfway before move:",post_relax_dihed) + print("Adjusted theta:",theta_adj) #calculate the rotation matrix - rot_matrix = self.rotation_matrix(rot_axis, theta_adj) + rot_matrix = self.rotation_matrix(rot_axis, -theta_adj) - target_atoms = selected_bond[0]['atms2mv'] + target_atoms = self.selected_bond[0]['atms2mv'] # apply the rotation matrix to the target atoms for idx, atom in enumerate(target_atoms): @@ -943,6 +946,7 @@ def move(self, context, verbose=False): if self.verbose: print('The updated position for this atom is:', model.positions[atom]) + print("Dihed after rot:",self.getDihedral(nc_positions,self.dihed_atoms)) # update the actual ncmc context object with the new positions context.setPositions(nc_positions) @@ -976,7 +980,10 @@ def afterMove(self, context): final_angle = self.getDihedral(post_pos,self.dihed_atoms) if self.is_in_bin(final_angle,self.target_bin): self.bin_boolean = True - else: self.bin_boolean = False + else: + self.bin_boolean = False + print(final_angle) + print(self.target_bin) return context diff --git a/blues/simulation.py b/blues/simulation.py index a9e41741..1da6349c 100644 --- a/blues/simulation.py +++ b/blues/simulation.py @@ -1068,7 +1068,8 @@ def _stepNCMC(self, nstepsNC, moveStep, move_engine=None): #Attempt anything related to the move before protocol is performed if not step: self._ncmc_sim.context = move_engine.selected_move.beforeMove(self._ncmc_sim.context) - if self.move_engine.moves[self.move_engine.selected_move].make_NCMC_move == False: + if self._move_engine.moves[0].make_NCMC_move == False: + print("NCMC move false") break # Attempt selected MoveEngine Move at the halfway point #to ensure protocol is symmetric @@ -1136,10 +1137,10 @@ def _acceptRejectMove(self, write_move=False): work_ncmc = work_ncmc + correction_factor #added for sidechain biasing - if self.move_engine.moves[self.move_engine.selected_move].bin_boolean == False: + if self._move_engine.moves[0].bin_boolean == False: print("Bin Boolean false: dihedral angle departed from target bin during NCMC move") - if work_ncmc > randnum and self.move_engine.moves[self.move_engine.selected_move].bin_boolean: + if work_ncmc > randnum and self._move_engine.moves[0].bin_boolean == True: self.accept += 1 logger.info('NCMC MOVE ACCEPTED: work_ncmc {} > randnum {}'.format(work_ncmc, randnum)) @@ -1249,7 +1250,7 @@ def run(self, nIter=0, nstepsNC=0, moveStep=0, nstepsMD=0, temperature=300, writ logger.info('BLUES Iteration: %s' % N) self._syncStatesMDtoNCMC() self._stepNCMC(nstepsNC, moveStep) - if self.move_engine.moves[self.move_engine.selected_move].make_NCMC_move == True: + if self._move_engine.moves[0].make_NCMC_move == True: self._acceptRejectMove(write_move) self._resetSimulations(temperature) self._stepMD(nstepsMD) From 9d4a93c76e3455a7f5acee56959a49e9aa367576 Mon Sep 17 00:00:00 2001 From: Kalistyn Burley Date: Wed, 13 Feb 2019 12:51:05 -0800 Subject: [PATCH 04/11] cleaned up print statements. closed yaml --- blues/moves.py | 18 +++++------------- 1 file changed, 5 insertions(+), 13 deletions(-) diff --git a/blues/moves.py b/blues/moves.py index 644e1500..43052b5c 100644 --- a/blues/moves.py +++ b/blues/moves.py @@ -760,7 +760,9 @@ def getRotBondAtoms(self): rot_atoms = self.getRotAtoms(rot_bonds, self.molecule, backbone_atoms) # Read in yaml file - rot_pref = yaml.load(open('rotpref.yaml')) + rotfile = open('rotpref.yaml',"r") + rot_pref = yaml.load(rotfile) + #rot_pref = yaml.load(open('rotpref.yaml')) # Restructure dictionary and populate atms2mv and bin_pref for residx in rot_atoms.keys(): resname = rot_atoms[residx]['res_name'] @@ -773,7 +775,7 @@ def getRotBondAtoms(self): atoms = [atom1-3, atom1, atom2, atom2+1] rot_atoms[residx]['chis'][chi]['dihed_atms'] = atoms rot_atoms[residx]['chis'][chi]['bin_pref'] = rot_pref[resname][chi] - + rotfile.close() return rot_atoms, rot_bonds, qry_atoms def chooseBond(self,positions): @@ -792,7 +794,7 @@ def chooseBond(self,positions): bonds_in_bins.append([(self.rot_atoms[residx]['chis'][chi]),curr_angle[0][0],bin_idx]) break - print("\nThere are %i bonds that fall within the biasing bins: %s\n"%(len(bonds_in_bins),bonds_in_bins)) + print("\nThere are %i bonds that fall within the biasing bins:"%(len(bonds_in_bins))) try: @@ -896,22 +898,15 @@ def move(self, context, verbose=True): positions = model.positions if self.selected_bond: - self.make_NCMC_move = True - print("choosing new bond") new_theta, self.target_bin = self.chooseTheta(self.selected_bond) self.dihed_atoms = [self.selected_bond[0]['dihed_atms']] - print("Dihed atoms:",self.dihed_atoms) - print("Dihed before NCMC relax:",self.selected_bond[1]) axis1 = self.dihed_atoms[0][1] axis2 = self.dihed_atoms[0][2] - print("Rot Axis:",axis1,axis2) rot_axis = (positions[axis1] - positions[axis2])/positions.unit #Adjust theta to account for repositioning during NCMC relaxation post_relax_dihed = self.getDihedral(initial_positions,self.dihed_atoms) theta_adj = self.selected_bond[1] - post_relax_dihed + new_theta - print("Dihed at halfway before move:",post_relax_dihed) - print("Adjusted theta:",theta_adj) #calculate the rotation matrix rot_matrix = self.rotation_matrix(rot_axis, -theta_adj) @@ -946,7 +941,6 @@ def move(self, context, verbose=True): if self.verbose: print('The updated position for this atom is:', model.positions[atom]) - print("Dihed after rot:",self.getDihedral(nc_positions,self.dihed_atoms)) # update the actual ncmc context object with the new positions context.setPositions(nc_positions) @@ -982,8 +976,6 @@ def afterMove(self, context): self.bin_boolean = True else: self.bin_boolean = False - print(final_angle) - print(self.target_bin) return context From 05f532db59046fd3cac8738f7fe336e547e31189 Mon Sep 17 00:00:00 2001 From: sgill2 Date: Fri, 15 Feb 2019 14:03:45 -0800 Subject: [PATCH 05/11] added support for acceptance_ratio option to adjust/reject move acceptance --- blues/moves.py | 2 + blues/simulation.py | 88 ++++++++++++++++++++-------------- blues/tests/test_simulation.py | 27 ++++++++++- 3 files changed, 81 insertions(+), 36 deletions(-) diff --git a/blues/moves.py b/blues/moves.py index 7322173f..da6db3bc 100644 --- a/blues/moves.py +++ b/blues/moves.py @@ -42,6 +42,7 @@ def __init__(self): """Initialize the Move object Currently empy. """ + self.acceptance_ratio = 1 def initializeSystem(self, system, integrator): """If the system or integrator needs to be modified to perform the move @@ -188,6 +189,7 @@ class RandomLigandRotationMove(Move): """ def __init__(self, structure, resname='LIG', random_state=None): + super(RandomLigandRotationMove, self).__init__() self.structure = structure self.resname = resname self.random_state = random_state diff --git a/blues/simulation.py b/blues/simulation.py index ddab123d..72e23eed 100644 --- a/blues/simulation.py +++ b/blues/simulation.py @@ -1062,36 +1062,41 @@ def _stepNCMC(self, nstepsNC, moveStep, move_engine=None): self._ncmc_sim.currentIter = self.currentIter move_engine.selectMove() - lastStep = nstepsNC - 1 - for step in range(int(nstepsNC)): - try: - #Attempt anything related to the move before protocol is performed - if not step: - self._ncmc_sim.context = move_engine.selected_move.beforeMove(self._ncmc_sim.context) - - # Attempt selected MoveEngine Move at the halfway point - #to ensure protocol is symmetric - if step == moveStep: - if hasattr(logger, 'report'): - logger.info = logger.report - #Do move - logger.info('Performing %s...' % move_engine.move_name) - self._ncmc_sim.context = move_engine.runEngine(self._ncmc_sim.context) - - # Do 1 NCMC step with the integrator - self._ncmc_sim.step(1) - - #Attempt anything related to the move after protocol is performed - if step == lastStep: - self._ncmc_sim.context = move_engine.selected_move.afterMove(self._ncmc_sim.context) + #Attempt anything related to the move before protocol is performed + self._ncmc_sim.context = move_engine.selected_move.beforeMove(self._ncmc_sim.context) + acceptance_ratio = move_engine.selected_move.acceptance_ratio + #check if acceptance_ratio is non-zero, if so perform integration + if not np.isclose(0, acceptance_ratio): + + for step in range(int(nstepsNC)): + try: + # Attempt selected MoveEngine Move at the halfway point + #to ensure protocol is symmetric + if step == moveStep: + if hasattr(logger, 'report'): + logger.info = logger.report + #Do move + logger.info('Performing %s...' % move_engine.move_name) + self._ncmc_sim.context = move_engine.runEngine(self._ncmc_sim.context) + if np.isclose(0, acceptance_ratio): + break + + # Do 1 NCMC step with the integrator + self._ncmc_sim.step(1) + + except Exception as e: + logger.error(e) + move_engine.selected_move._error(self._ncmc_sim.context) + break + #Attempt anything related to the move after protocol is performed + self._ncmc_sim.context = move_engine.selected_move.afterMove(self._ncmc_sim.context) - except Exception as e: - logger.error(e) - move_engine.selected_move._error(self._ncmc_sim.context) - break + else: + logger.info('move acceptance_ratio is 0; skipping integration') # ncmc_state1 stores the state AFTER a proposed move. ncmc_state1 = self.getStateFromContext(self._ncmc_sim.context, self._state_keys) + self._setStateTable('ncmc', 'state1', ncmc_state1) def _computeAlchemicalCorrection(self): @@ -1114,6 +1119,18 @@ def _computeAlchemicalCorrection(self): -1.0 / self._ncmc_sim.context._integrator.kT) return correction_factor + def _checkPEConsistencyOnRejection(self): + """Checks wheter the potential energy is the same as the last MD iteration if a move is rejected + (it should be) + """ + md_state0 = self.stateTable['md']['state0'] + md_PE = self._md_sim.context.getState(getEnergy=True).getPotentialEnergy() + if not math.isclose(md_state0['potential_energy']._value, md_PE._value, rel_tol=float('1e-%s' % rtol)): + logger.error( + 'Last MD potential energy %s != Current MD potential energy %s. Potential energy should match the prior state.' + % (md_state0['potential_energy'], md_PE)) + sys.exit(1) + def _acceptRejectMove(self, write_move=False): """Choose to accept or reject the proposed move based @@ -1126,6 +1143,7 @@ def _acceptRejectMove(self, write_move=False): """ work_ncmc = self._ncmc_sim.context._integrator.getLogAcceptanceProbability(self._ncmc_sim.context) randnum = math.log(np.random.random()) + acceptance_ratio = self._move_engine.selected_move.acceptance_ratio # Compute correction if work_ncmc is not NaN if not np.isnan(work_ncmc): @@ -1134,7 +1152,14 @@ def _acceptRejectMove(self, write_move=False): 'NCMCLogAcceptanceProbability = %.6f + Alchemical Correction = %.6f' % (work_ncmc, correction_factor)) work_ncmc = work_ncmc + correction_factor - if work_ncmc > randnum: + if np.isclose(0, acceptance_ratio): + print('accept probablility', acceptance_ratio) + logger.info('NCMC MOVE REJECTED: acceptance_ratio for the selected move is 0; work_ncmc {}'.format(work_ncmc)) + self._checkPEConsistencyOnRejection() + + + + elif work_ncmc + math.log(acceptance_ratio) > randnum: self.accept += 1 logger.info('NCMC MOVE ACCEPTED: work_ncmc {} > randnum {}'.format(work_ncmc, randnum)) @@ -1152,15 +1177,8 @@ def _acceptRejectMove(self, write_move=False): # If reject move, do nothing, # NCMC simulation be updated from MD Simulation next iteration. - # Potential energy should be from last MD step in the previous iteration - md_state0 = self.stateTable['md']['state0'] - md_PE = self._md_sim.context.getState(getEnergy=True).getPotentialEnergy() - if not math.isclose(md_state0['potential_energy']._value, md_PE._value, rel_tol=float('1e-%s' % rtol)): - logger.error( - 'Last MD potential energy %s != Current MD potential energy %s. Potential energy should match the prior state.' - % (md_state0['potential_energy'], md_PE)) - sys.exit(1) + self._checkPEConsistencyOnRejection() def _resetSimulations(self, temperature=None): """At the end of each iteration: diff --git a/blues/tests/test_simulation.py b/blues/tests/test_simulation.py index 8bf4a13f..48b2495e 100644 --- a/blues/tests/test_simulation.py +++ b/blues/tests/test_simulation.py @@ -75,6 +75,21 @@ def system(structure, system_cfg): class NoRandomLigandRotation(RandomLigandRotationMove): + def __init__(self, structure, resname='LIG'): + print('type structure', type(structure), structure) + super(NoRandomLigandRotation, self).__init__(structure, resname) + self.before = False + self.after= False + self.acceptance_ratio = 1 + def beforeMove(self, context): + if self.before: + self.acceptance_ratio = 0 + return context + def afterMove(self, context): + if self.after: + self.acceptance_ratio = 0 + return context + def move(self, context): return context @@ -490,7 +505,13 @@ def test_blues_simulationRunYAML(self, tmpdir, structure, tol_atom_indices, syst pos_compare = np.not_equal(before_iter, after_iter).all() assert pos_compare - def test_blues_simulationRunPython(self, systems, simulations, engine, tmpdir, sim_cfg): + @pytest.mark.parametrize('before,after', [ + (True, False), + (False, True), + (True, True), + (False, False) + ]) + def test_blues_simulationRunPython(self, systems, simulations, engine, tmpdir, sim_cfg, before, after): print('Testing BLUESSimulation.run() from pure python') md_rep_cfg = { 'stream': { @@ -520,6 +541,10 @@ def test_blues_simulationRunPython(self, systems, simulations, engine, tmpdir, s md_reporters = ReporterConfig(tmpdir.join('tol-test'), md_rep_cfg).makeReporters() ncmc_reporters = ReporterConfig(tmpdir.join('tol-test-ncmc'), ncmc_rep_cfg).makeReporters() + engine.moves[0].acceptance_ratio = 1 + engine.moves[0].before = before + engine.moves[0].after = after + simulations = SimulationFactory( systems, engine, sim_cfg, md_reporters=md_reporters, ncmc_reporters=ncmc_reporters) From 9ec3cd636136559ab2065d5ab51d0efdaa4714bb Mon Sep 17 00:00:00 2001 From: sgill2 Date: Fri, 15 Feb 2019 14:27:35 -0800 Subject: [PATCH 06/11] added acceptance_ratio handling for MC --- blues/simulation.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/blues/simulation.py b/blues/simulation.py index 72e23eed..5dbfffc8 100644 --- a/blues/simulation.py +++ b/blues/simulation.py @@ -1305,8 +1305,14 @@ def _acceptRejectMove(self, temperature=None): work_mc = (md_state1['potential_energy'] - md_state0['potential_energy']) * ( -1.0 / self._ncmc_sim.context._integrator.kT) randnum = math.log(np.random.random()) + acceptance_ratio = self._move_engine.selected_move.acceptance_ratio + + if np.isclose(0, acceptance_ratio): + self.reject += 1 + logger.info('MC MOVE REJECTED: work_mc {} < {}'.format(work_mc, randnum)) + self._md_sim.context.setPositions(md_state0['positions']) - if work_mc > randnum: + elif work_mc > randnum: self.accept += 1 logger.info('MC MOVE ACCEPTED: work_mc {} > randnum {}'.format(work_mc, randnum)) self._md_sim.context.setPositions(md_state1['positions']) From d1d246afa1c5a138f86c87301045a08115e2fdb0 Mon Sep 17 00:00:00 2001 From: Kalistyn Burley Date: Tue, 19 Feb 2019 13:07:27 -0800 Subject: [PATCH 07/11] updated yaml --- blues/example.py | 6 ++--- blues/moves.py | 3 +-- blues/tests/rotpref.yaml | 58 ---------------------------------------- 3 files changed, 4 insertions(+), 63 deletions(-) delete mode 100644 blues/tests/rotpref.yaml diff --git a/blues/example.py b/blues/example.py index 190a916a..72f8da53 100644 --- a/blues/example.py +++ b/blues/example.py @@ -35,7 +35,7 @@ def sidechain_example(yaml_file): structure = cfg['Structure'] #Select move type - sidechain = SideChainMove(structure, [1]) + sidechain = SideChainMove(structure, [1], verbose=False) #Iniitialize object that selects movestep sidechain_mover = MoveEngine(sidechain) @@ -62,5 +62,5 @@ def sidechain_example(yaml_file): output.write("%s\n" % str(value)[1:-1]) -ligrot_example(get_data_filename('blues', '../examples/rotmove_cuda.yml')) -#sidechain_example(get_data_filename('blues', '../examples/sidechain_cuda.yml')) +#ligrot_example(get_data_filename('blues', '../examples/rotmove_cuda.yml')) +sidechain_example(get_data_filename('blues', '../examples/sidechain_cuda.yml')) diff --git a/blues/moves.py b/blues/moves.py index 43052b5c..17c54810 100644 --- a/blues/moves.py +++ b/blues/moves.py @@ -760,9 +760,8 @@ def getRotBondAtoms(self): rot_atoms = self.getRotAtoms(rot_bonds, self.molecule, backbone_atoms) # Read in yaml file - rotfile = open('rotpref.yaml',"r") + rotfile = open('rotpref.yml',"r") rot_pref = yaml.load(rotfile) - #rot_pref = yaml.load(open('rotpref.yaml')) # Restructure dictionary and populate atms2mv and bin_pref for residx in rot_atoms.keys(): resname = rot_atoms[residx]['res_name'] diff --git a/blues/tests/rotpref.yaml b/blues/tests/rotpref.yaml deleted file mode 100644 index 47b3353e..00000000 --- a/blues/tests/rotpref.yaml +++ /dev/null @@ -1,58 +0,0 @@ -ARG: - 1: [-68, 56, -171] - 2: [-68, 67, 180] - 3: [-63, 68, -177] - 4: [-89, 86, -174] - 5: [-60, 60, 180] -ASN: - 1: [-68, 59, 180] - 2: [-107, 53] -ASP: - 1: [-75,74,-178] - 2: [107,62,-171] -CYS: - 1: [-68, 56, -171] -GLN: - 1: [-68, 56, -171] - 2: [-68, 56, -171] - 3: [-68, 56, -171] -GLU: - 1: [-68, 56, -171] - 2: [-68, 56, -171] - 3: [-68, 56, -171] -HIS: - 1: [-68, 56, -171] - 2: [-68, 56, -171] -ILE: - 1: [-68, 56, -171] - 2: [-68, 56, -171] -LEU: - 1: [-68, 56, -171] - 2: [-68, 56, -171] -LYS: - 1: [-68, 56, -171] - 2: [-68, 56, -171] - 3: [-68, 56, -171] - 4: [-68, 56, -171] - 5: [-68, 56, -171] -MET: - 1: [-60, 60, 180] - 2: [-60, 60, 180] - 3: [-70, 70, 175] -PHE: - 1: [-79, 52, -177] - 2: [-90,90] -PRO: - 1: [-68, 56, -171] -SER: - 1: [-68, 56, -171] -THR: - 1: [-60, 60, 180] -TRP: - 1: [-68, 56, -171] - 2: [-68, 56, -171] -TYR: - 1: [-68, 56, -171] - 2: [-68, 56, -171] -VAL: - 1: [-72,57,181] From 859d219a9828f16cfacb21c7571d6aec92c16fb0 Mon Sep 17 00:00:00 2001 From: Kalistyn Burley Date: Tue, 19 Feb 2019 14:10:39 -0800 Subject: [PATCH 08/11] reattempting commit with accept-rat --- blues/moves.py | 14 ++++---------- blues/simulation.py | 8 -------- 2 files changed, 4 insertions(+), 18 deletions(-) diff --git a/blues/moves.py b/blues/moves.py index bdc0aad0..8c2821a3 100644 --- a/blues/moves.py +++ b/blues/moves.py @@ -208,8 +208,6 @@ def __init__(self, structure, resname='LIG', random_state=None): self.center_of_mass = None self.positions = structure[self.atom_indices].positions ## adding placeholder values for random ligand moves (always true for these moves) vs biased sidechain KB - self.make_NCMC_move = True - self.bin_boolean = True self._calculateProperties() @@ -849,10 +847,8 @@ def beforeMove(self, context): self.start_pos = context.getState(getPositions=True).getPositions(asNumpy=True) self.selected_bond = self.chooseBond(self.start_pos) - if self.selected_bond: - self.make_NCMC_move = True - else: - self.make_NCMC_move = False + if not self.selected_bond: + self.acceptance_ratio = 0 return context @@ -973,10 +969,8 @@ def afterMove(self, context): """ post_pos = context.getState(getPositions=True).getPositions(asNumpy=True) final_angle = self.getDihedral(post_pos,self.dihed_atoms) - if self.is_in_bin(final_angle,self.target_bin): - self.bin_boolean = True - else: - self.bin_boolean = False + if not self.is_in_bin(final_angle,self.target_bin): + self.acceptance_ratio = False return context diff --git a/blues/simulation.py b/blues/simulation.py index 8e647b0f..b8f3523a 100644 --- a/blues/simulation.py +++ b/blues/simulation.py @@ -1152,13 +1152,6 @@ def _acceptRejectMove(self, write_move=False): 'NCMCLogAcceptanceProbability = %.6f + Alchemical Correction = %.6f' % (work_ncmc, correction_factor)) work_ncmc = work_ncmc + correction_factor -<<<<<<< HEAD - #added for sidechain biasing - if self._move_engine.moves[0].bin_boolean == False: - print("Bin Boolean false: dihedral angle departed from target bin during NCMC move") - - if work_ncmc > randnum and self._move_engine.moves[0].bin_boolean == True: -======= if np.isclose(0, acceptance_ratio): print('accept probablility', acceptance_ratio) logger.info('NCMC MOVE REJECTED: acceptance_ratio for the selected move is 0; work_ncmc {}'.format(work_ncmc)) @@ -1167,7 +1160,6 @@ def _acceptRejectMove(self, write_move=False): elif work_ncmc + math.log(acceptance_ratio) > randnum: ->>>>>>> 9ec3cd636136559ab2065d5ab51d0efdaa4714bb self.accept += 1 logger.info('NCMC MOVE ACCEPTED: work_ncmc {} > randnum {}'.format(work_ncmc, randnum)) From 263594a7f84efd4803fb2b46ae9644e40d0d354a Mon Sep 17 00:00:00 2001 From: Kalistyn Burley Date: Wed, 8 May 2019 10:52:06 -0700 Subject: [PATCH 09/11] added/updated rotpref --- blues/moves.py | 3 +- blues/simulation.py | 5 ++- blues/tests/rotpref.yml | 58 +++++++++++++++++++++++++++++++++++ examples/example_rotmove.py | 4 +-- examples/example_sidechain.py | 6 ++-- examples/rotmove_cuda.yml | 4 +-- examples/rotpref.yml | 58 +++++++++++++++++++++++++++++++++++ examples/sidechain_cuda.yml | 4 +-- 8 files changed, 129 insertions(+), 13 deletions(-) create mode 100644 blues/tests/rotpref.yml create mode 100644 examples/rotpref.yml diff --git a/blues/moves.py b/blues/moves.py index 8c2821a3..20e60572 100644 --- a/blues/moves.py +++ b/blues/moves.py @@ -849,7 +849,8 @@ def beforeMove(self, context): self.selected_bond = self.chooseBond(self.start_pos) if not self.selected_bond: self.acceptance_ratio = 0 - + else: + self.acceptance_ratio = 1 return context def move(self, context, verbose=True): diff --git a/blues/simulation.py b/blues/simulation.py index b8f3523a..5dbfffc8 100644 --- a/blues/simulation.py +++ b/blues/simulation.py @@ -1262,9 +1262,8 @@ def run(self, nIter=0, nstepsNC=0, moveStep=0, nstepsMD=0, temperature=300, writ logger.info('BLUES Iteration: %s' % N) self._syncStatesMDtoNCMC() self._stepNCMC(nstepsNC, moveStep) - if self._move_engine.moves[0].make_NCMC_move == True: - self._acceptRejectMove(write_move) - self._resetSimulations(temperature) + self._acceptRejectMove(write_move) + self._resetSimulations(temperature) self._stepMD(nstepsMD) # END OF NITER diff --git a/blues/tests/rotpref.yml b/blues/tests/rotpref.yml new file mode 100644 index 00000000..9d7a6697 --- /dev/null +++ b/blues/tests/rotpref.yml @@ -0,0 +1,58 @@ +ARG: + 1: [-68, 60, -171] + 2: [-67, 70, 182] + 3: [-66, 68, -178] + 4: [-88, 86, -179] + 5: [-60, 60, 180] +ASN: + 1: [-68, 59, 180] + 2: [-107, 53] +ASP: + 1: [-75, 74, -182] + 2: [107, 62, -171] +CYS: + 1: [-62, 55, -176] +GLN: + 1: [-67, 56, -171] + 2: [-66, 82, -175] + 3: [-85, 90] +GLU: + 1: [-179,-68,60] + 2: [-179,-70,85] + 3: [-68,88] +HIS: + 1: [-60, 58, -179] + 2: [-80, 100] +ILE: + 1: [-60, 50, 170] + 2: [-60, 70, -178] +LEU: + 1: [-60, 180] + 2: [68, 180] +LYS: + 1: [-68, 60, -178] + 2: [-68, 56, -178] + 3: [-68, 56, -178] + 4: [-68, 60, -178] + 5: [-68, 56, -178] +MET: + 1: [-70, 60, 180] + 2: [-60, 60, 180] + 3: [-70, 70, 175] +PHE: + 1: [-79, 52, -177] + 2: [-90, 90] +PRO: + 1: [-30, 30] +SER: + 1: [-60, 60, -171] +THR: + 1: [-60, 55, 170] +TRP: + 1: [-68, 68, 178] + 2: [-90, 0, 90] +TYR: + 1: [-68, 68, 178] + 2: [-90, 90] +VAL: + 1: [-65,65,178] diff --git a/examples/example_rotmove.py b/examples/example_rotmove.py index cd1d0fce..b9808fbd 100644 --- a/examples/example_rotmove.py +++ b/examples/example_rotmove.py @@ -6,7 +6,7 @@ def rotmove_cuda(yaml_file): # Parse a YAML configuration, return as Dict - cfg = Settings('rotmove_cuda.yaml').asDict() + cfg = Settings('rotmove_cuda.yml').asDict() structure = cfg['Structure'] #Select move type @@ -30,4 +30,4 @@ def rotmove_cuda(yaml_file): if __name__ == "__main__": - rotmove_cuda('rotmove_cuda.yaml') + rotmove_cuda('rotmove_cuda.yml') diff --git a/examples/example_sidechain.py b/examples/example_sidechain.py index a7f980b2..7f3aa343 100644 --- a/examples/example_sidechain.py +++ b/examples/example_sidechain.py @@ -1,11 +1,11 @@ from blues.moves import SideChainMove -from blues.engine import MoveEngine +from blues.moves import MoveEngine from blues.simulation import * import json from blues.settings import * # Parse a YAML configuration, return as Dict -cfg = Settings('sidechain_cuda.yaml').asDict() +cfg = Settings('sidechain_cuda.yml').asDict() structure = cfg['Structure'] #Select move type @@ -28,7 +28,7 @@ import mdtraj as md import numpy as np -traj = md.load_netcdf('vacDivaline-test/vacDivaline.nc', top='tests/data/vacDivaline.prmtop') +traj = md.load_netcdf('vacDivaline-test/vacDivaline.nc', top='../blues/tests/data/vacDivaline.prmtop') indicies = np.array([[0, 4, 6, 8]]) dihedraldata = md.compute_dihedrals(traj, indicies) with open("vacDivaline-test/dihedrals.txt", 'w') as output: diff --git a/examples/rotmove_cuda.yml b/examples/rotmove_cuda.yml index d15f62ec..9c394e50 100644 --- a/examples/rotmove_cuda.yml +++ b/examples/rotmove_cuda.yml @@ -13,8 +13,8 @@ logger: stream: True structure: - filename: tests/data/eqToluene.prmtop - xyz: tests/data/eqToluene.inpcrd + filename: ../blues/tests/data/eqToluene.prmtop + xyz: ../blues/tests/data/eqToluene.inpcrd system: nonbondedMethod: PME diff --git a/examples/rotpref.yml b/examples/rotpref.yml new file mode 100644 index 00000000..9d7a6697 --- /dev/null +++ b/examples/rotpref.yml @@ -0,0 +1,58 @@ +ARG: + 1: [-68, 60, -171] + 2: [-67, 70, 182] + 3: [-66, 68, -178] + 4: [-88, 86, -179] + 5: [-60, 60, 180] +ASN: + 1: [-68, 59, 180] + 2: [-107, 53] +ASP: + 1: [-75, 74, -182] + 2: [107, 62, -171] +CYS: + 1: [-62, 55, -176] +GLN: + 1: [-67, 56, -171] + 2: [-66, 82, -175] + 3: [-85, 90] +GLU: + 1: [-179,-68,60] + 2: [-179,-70,85] + 3: [-68,88] +HIS: + 1: [-60, 58, -179] + 2: [-80, 100] +ILE: + 1: [-60, 50, 170] + 2: [-60, 70, -178] +LEU: + 1: [-60, 180] + 2: [68, 180] +LYS: + 1: [-68, 60, -178] + 2: [-68, 56, -178] + 3: [-68, 56, -178] + 4: [-68, 60, -178] + 5: [-68, 56, -178] +MET: + 1: [-70, 60, 180] + 2: [-60, 60, 180] + 3: [-70, 70, 175] +PHE: + 1: [-79, 52, -177] + 2: [-90, 90] +PRO: + 1: [-30, 30] +SER: + 1: [-60, 60, -171] +THR: + 1: [-60, 55, 170] +TRP: + 1: [-68, 68, 178] + 2: [-90, 0, 90] +TYR: + 1: [-68, 68, 178] + 2: [-90, 90] +VAL: + 1: [-65,65,178] diff --git a/examples/sidechain_cuda.yml b/examples/sidechain_cuda.yml index eac7cf58..c4ea6870 100644 --- a/examples/sidechain_cuda.yml +++ b/examples/sidechain_cuda.yml @@ -12,8 +12,8 @@ logger: stream: True structure: - filename: tests/data/vacDivaline.prmtop - xyz: tests/data/vacDivaline.inpcrd + filename: ../blues/tests/data/vacDivaline.prmtop + xyz: ../blues/tests/data/vacDivaline.inpcrd system: nonbondedMethod: PME From ef57e4d390036f6e43455b3e7330f9be5d75cb96 Mon Sep 17 00:00:00 2001 From: sgill2 Date: Wed, 2 Sep 2020 17:27:32 -0700 Subject: [PATCH 10/11] Updated to handle bias and unbiased sidechain rotation attempts --- blues/moves.py | 115 +++++++++++++++++++++++++++++++++++-------------- 1 file changed, 83 insertions(+), 32 deletions(-) diff --git a/blues/moves.py b/blues/moves.py index 20e60572..28606c08 100644 --- a/blues/moves.py +++ b/blues/moves.py @@ -27,7 +27,7 @@ indices_chi1, indices_chi4, indices_chi5) - +from blues import utils try: import openeye.oechem as oechem if not oechem.OEChemIsLicensed(): @@ -444,6 +444,10 @@ class SideChainMove(Move): List of the residue numbers of the sidechains to be rotated. verbose : bool, default=False Enable verbosity to print out detailed information of the rotation. + bias_range: float, default=15.0 + The range around each Dunbrak minima to be considered for rotation. + If 0, then the Dunbrak minimas are not considered and a random rotation + of the sidechain takes place instead. write_move : bool, default=False If True, writes a PDB of the system after rotation. @@ -473,7 +477,7 @@ class SideChainMove(Move): """ - def __init__(self, structure, residue_list, bias_range=15, verbose=False, write_move=False): + def __init__(self, structure, residue_list, bias_range=15.0, verbose=False, write_move=False): self.structure = structure self.molecule = self._pmdStructureToOEMol() self.residue_list = residue_list @@ -483,6 +487,7 @@ def __init__(self, structure, residue_list, bias_range=15, verbose=False, write_ self.atom_indices = self.rot_atoms self.verbose = verbose self.write_move = write_move + self.bias = bool(bias_range) def _pmdStructureToOEMol(self): """Helper function for converting the parmed structure into an OEMolecule.""" @@ -760,7 +765,7 @@ def getRotBondAtoms(self): rot_atoms = self.getRotAtoms(rot_bonds, self.molecule, backbone_atoms) # Read in yaml file - rotfile = open('rotpref.yml',"r") + rotfile = open(utils.get_data_filename('blues', 'tests/rotpref.yml'),"r") rot_pref = yaml.load(rotfile) # Restructure dictionary and populate atms2mv and bin_pref for residx in rot_atoms.keys(): @@ -844,13 +849,14 @@ def rotation_matrix(self, axis, theta): [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]]) def beforeMove(self, context): - + self.selected_bond = 1 self.start_pos = context.getState(getPositions=True).getPositions(asNumpy=True) - self.selected_bond = self.chooseBond(self.start_pos) - if not self.selected_bond: - self.acceptance_ratio = 0 - else: - self.acceptance_ratio = 1 + if self.bias: + self.selected_bond = self.chooseBond(self.start_pos) + if not self.selected_bond: + self.acceptance_ratio = 0 + else: + self.acceptance_ratio = 1 return context def move(self, context, verbose=True): @@ -895,20 +901,64 @@ def move(self, context, verbose=True): positions = model.positions - if self.selected_bond: - new_theta, self.target_bin = self.chooseTheta(self.selected_bond) - self.dihed_atoms = [self.selected_bond[0]['dihed_atms']] - axis1 = self.dihed_atoms[0][1] - axis2 = self.dihed_atoms[0][2] - rot_axis = (positions[axis1] - positions[axis2])/positions.unit + #if bias is being used + if self.bias : + #check if it's within + if self.selected_bond: + new_theta, self.target_bin = self.chooseTheta(self.selected_bond) + self.dihed_atoms = [self.selected_bond[0]['dihed_atms']] + axis1 = self.dihed_atoms[0][1] + axis2 = self.dihed_atoms[0][2] + rot_axis = (positions[axis1] - positions[axis2])/positions.unit - #Adjust theta to account for repositioning during NCMC relaxation - post_relax_dihed = self.getDihedral(initial_positions,self.dihed_atoms) - theta_adj = self.selected_bond[1] - post_relax_dihed + new_theta - #calculate the rotation matrix - rot_matrix = self.rotation_matrix(rot_axis, -theta_adj) + #Adjust theta to account for repositioning during NCMC relaxation + post_relax_dihed = self.getDihedral(initial_positions,self.dihed_atoms) + theta_adj = self.selected_bond[1] - post_relax_dihed + new_theta + #calculate the rotation matrix + rot_matrix = self.rotation_matrix(rot_axis, -theta_adj) + + target_atoms = self.selected_bond[0]['atms2mv'] + + # apply the rotation matrix to the target atoms + for idx, atom in enumerate(target_atoms): + + my_position = positions[atom] + + if self.verbose: + print('The current position for %i is: %s' % (atom, my_position)) - target_atoms = self.selected_bond[0]['atms2mv'] + # find the reduced position (substract out axis) + red_position = (my_position - model.positions[axis2])._value + # find the new positions by multiplying by rot matrix + new_position = numpy.dot(rot_matrix, red_position) * positions.unit + positions[axis2] + + if self.verbose: print("The new position should be:", new_position) + + positions[atom] = new_position + # Update the parmed model with the new positions + model.atoms[atom].xx = new_position[0] / positions.unit + model.atoms[atom].xy = new_position[1] / positions.unit + model.atoms[atom].xz = new_position[2] / positions.unit + + #update the copied ncmc context array with the new positions + nc_positions[atom][0] = model.atoms[atom].xx * nc_positions.unit / 10 + nc_positions[atom][1] = model.atoms[atom].xy * nc_positions.unit / 10 + nc_positions[atom][2] = model.atoms[atom].xz * nc_positions.unit / 10 + + if self.verbose: + print('The updated position for this atom is:', model.positions[atom]) + + #if bias is not being used we can rotate randomly + else: + positions = model.positions + + # find the rotation axis using the updated positions + axis1 = target_atoms[0] + axis2 = target_atoms[1] + rot_axis = (positions[axis1] - positions[axis2]) / positions.unit + + #calculate the rotation matrix + rot_matrix = self.rotation_matrix(rot_axis, theta) # apply the rotation matrix to the target atoms for idx, atom in enumerate(target_atoms): @@ -939,15 +989,15 @@ def move(self, context, verbose=True): if self.verbose: print('The updated position for this atom is:', model.positions[atom]) - # update the actual ncmc context object with the new positions - context.setPositions(nc_positions) + # update the actual ncmc context object with the new positions + context.setPositions(nc_positions) - # update the class structure positions - self.structure.positions = model.positions + # update the class structure positions + self.structure.positions = model.positions - if self.write_move: - filename = 'sc_move_%s_%s_%s.pdb' % (res, axis1, axis2) - mod_prot = model.save(filename, overwrite=True) + if self.write_move: + filename = 'sc_move_%s_%s_%s.pdb' % (res, axis1, axis2) + mod_prot = model.save(filename, overwrite=True) return context @@ -968,10 +1018,11 @@ def afterMove(self, context): The same input context, but whose context were changed by this function. """ - post_pos = context.getState(getPositions=True).getPositions(asNumpy=True) - final_angle = self.getDihedral(post_pos,self.dihed_atoms) - if not self.is_in_bin(final_angle,self.target_bin): - self.acceptance_ratio = False + if self.bias: + post_pos = context.getState(getPositions=True).getPositions(asNumpy=True) + final_angle = self.getDihedral(post_pos,self.dihed_atoms) + if not self.is_in_bin(final_angle,self.target_bin): + self.acceptance_ratio = 0 return context From 492ff74f19822eeca4b633776be79497db41e00d Mon Sep 17 00:00:00 2001 From: sgill2 Date: Wed, 2 Sep 2020 17:28:22 -0700 Subject: [PATCH 11/11] update test to handle bias and unbiased rotation case --- blues/tests/test_sidechain.py | 84 ++++++++++++++++++++++++++++++++++- 1 file changed, 82 insertions(+), 2 deletions(-) diff --git a/blues/tests/test_sidechain.py b/blues/tests/test_sidechain.py index 83543083..b0a17fa2 100644 --- a/blues/tests/test_sidechain.py +++ b/blues/tests/test_sidechain.py @@ -24,7 +24,88 @@ @skipUnless(has_openeye, 'Cannot test SideChainMove without openeye-toolkits and oeommtools.') -class SideChainTester(unittest.TestCase): +class SideChainTesterDunbrak(unittest.TestCase): + """ + Test the SmartDartMove.move() function. + """ + + def setUp(self): + # Obtain topologies/positions + prmtop = utils.get_data_filename('blues', 'tests/data/vacDivaline.prmtop') + inpcrd = utils.get_data_filename('blues', 'tests/data/vacDivaline.inpcrd') + self.struct = parmed.load_file(prmtop, xyz=inpcrd) + + self.sidechain = SideChainMove(self.struct, [1], bias_range=15.0) + self.engine = MoveEngine(self.sidechain) + self.engine.selectMove() + + self.system_cfg = {'nonbondedMethod': app.NoCutoff, 'constraints': app.HBonds} + self.systems = SystemFactory(self.struct, self.sidechain.atom_indices, self.system_cfg) + + self.cfg = { + 'dt': 0.002 * unit.picoseconds, + 'friction': 1 * 1 / unit.picoseconds, + 'temperature': 300 * unit.kelvin, + 'nIter': 1, + 'nstepsMD': 1, + 'nstepsNC': 4, + 'alchemical_functions': { + 'lambda_sterics': + 'step(0.199999-lambda) + step(lambda-0.2)*step(0.8-lambda)*abs(lambda-0.5)*1/0.3 + step(lambda-0.800001)', + 'lambda_electrostatics': + 'step(0.2-lambda)- 1/0.2*lambda*step(0.2-lambda) + 1/0.2*(lambda-0.8)*step(lambda-0.8)' + } + } + + self.simulations = SimulationFactory(self.systems, self.engine, self.cfg) + + def test_getRotBondAtoms(self): + vals = [v for v in self.sidechain.rot_atoms[1]['chis'][1]['atms2mv']] + assert len(vals) == 11 + #Ensure it selects 1 rotatable bond in Valine + assert len(self.sidechain.rot_bonds) == 1 + + def test_sidechain_move(self): + #check to make sure the dart stays within the specified ranges during the simulation + atom_indices = [v for v in self.sidechain.rot_atoms[1]['chis'][1]['atms2mv']] + + for i in range(20): + before_context = self.simulations.ncmc.context + dobefore = self.sidechain.beforeMove(before_context) + before_move = self.simulations.ncmc.context.getState(getPositions=True).getPositions( + asNumpy=True)[atom_indices, :] + self.simulations.ncmc.context = self.engine.runEngine(self.simulations.ncmc.context) + after_move = self.simulations.ncmc.context.getState(getPositions=True).getPositions( + asNumpy=True)[atom_indices, :] + + #Check that our system has performed the move correctly + # Integrator must step for context to update positions + # Remove the first two atoms in check as these are the anchor atoms and are not rotated. + pos_compare = np.not_equal(before_move, after_move)[2:, :].all() + positions = self.simulations.ncmc.context.getState(getPositions=True).getPositions(asNumpy=True) + bond_check = False + bonds_in_bins = [] + #check the current angle after the move and make sure it's within range of the dihedral bins specified + for residx in self.sidechain.rot_atoms.keys(): + resname = self.sidechain.rot_atoms[residx]['res_name'] + for chi in self.sidechain.rot_atoms[residx]['chis']: + dihed_atoms = [self.sidechain.rot_atoms[residx]['chis'][chi]['dihed_atms']] + curr_angle = self.sidechain.getDihedral(positions,dihed_atoms) + bin_ct=0 + for bin in self.sidechain.rot_atoms[residx]['chis'][chi]['bin_pref']: + bin_ct+=1 + if self.sidechain.is_in_bin(curr_angle,bin): + + bin_idx = bin_ct-1 + bonds_in_bins.append([(self.sidechain.rot_atoms[residx]['chis'][chi]),curr_angle[0][0],bin_idx]) + bond_check = True + break + + assert bond_check == True + + assert pos_compare + +class SideChainTesterRandom(unittest.TestCase): """ Test the SmartDartMove.move() function. """ @@ -81,6 +162,5 @@ def test_sidechain_move(self): pos_compare = np.not_equal(before_move, after_move)[2:, :].all() assert pos_compare - if __name__ == "__main__": unittest.main()