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 7322173f..28606c08 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) +from blues import utils try: import openeye.oechem as oechem if not oechem.OEChemIsLicensed(): @@ -42,6 +49,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 +196,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 @@ -198,6 +207,8 @@ 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._calculateProperties() def getAtomIndices(self, structure, resname): @@ -433,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. @@ -462,15 +477,17 @@ class SideChainMove(Move): """ - def __init__(self, structure, residue_list, 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 + 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 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.""" @@ -535,8 +552,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 +564,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 +639,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 +681,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(numpy.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 +763,69 @@ 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 - - - Returns - ------- - theta_ran : - targetatoms : + # Read in yaml file + 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(): + 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] + 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] + rotfile.close() + return rot_atoms, rot_bonds, qry_atoms - 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']] + 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:"%(len(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,7 +848,18 @@ 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 move(self, context, verbose=False): + def beforeMove(self, context): + self.selected_bond = 1 + self.start_pos = context.getState(getPositions=True).getPositions(asNumpy=True) + 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): """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 +878,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,42 +901,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 + #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 - #calculate the rotation matrix - rot_matrix = self.rotation_matrix(rot_axis, theta) + #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) - # apply the rotation matrix to the target atoms - for idx, atom in enumerate(target_atoms): + target_atoms = self.selected_bond[0]['atms2mv'] - my_position = positions[atom] + # apply the rotation matrix to the target atoms + for idx, atom in enumerate(target_atoms): - if self.verbose: - print('The current position for %i is: %s' % (atom, my_position)) + my_position = positions[atom] - # 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 current position for %i is: %s' % (atom, my_position)) - if self.verbose: print("The new position should be:", new_position) + # 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] - 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 new position should be:", new_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 + 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 updated position for this atom is:', model.positions[atom]) + #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): + + my_position = positions[atom] + + if self.verbose: + print('The current position for %i is: %s' % (atom, my_position)) + + # 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]) # update the actual ncmc context object with the new positions context.setPositions(nc_positions) @@ -838,9 +998,34 @@ def move(self, context, verbose=False): 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. + + """ + 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 + 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..5dbfffc8 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: @@ -1287,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']) 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/blues/tests/test_sidechain.py b/blues/tests/test_sidechain.py index 75ea7146..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. """ @@ -60,13 +141,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) @@ -79,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() diff --git a/blues/tests/test_simulation.py b/blues/tests/test_simulation.py index caf742b3..7b55565c 100644 --- a/blues/tests/test_simulation.py +++ b/blues/tests/test_simulation.py @@ -78,6 +78,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 @@ -493,7 +508,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': { 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 340b0c87..b3518a86 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 0f30c431..fd86452e 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