TorsionDrive API


class torsiondrive.dihedral_scanner.DihedralScanner(engine, dihedrals, grid_spacing, init_coords_M=None, energy_decrease_thresh=None, dihedral_ranges=None, energy_upper_limit=None, extra_constraints=None, verbose=False)

DihedralScanner class is designed to create a dihedral grid, and fill in optimized geometries and energies into the grid, by running wavefront propagations of constrained optimizations

engine: QMEngine() instance

An QMEngine object, e.g. EnginePsi4, EngineQChem or EngineTerachem

dihedrals: List[(d1, d2, d3, d4), ..]

list of dihedral index tuples (d1, d2, d3, d4). The length of list determines the dimension of the grid i.e. dihedrals = [(0,1,2,3)] –> 1-D scan, dihedrals = [(0,1,2,3),(1,2,3,4)] –> 2-D Scan

grid_spacing: Int

Distance (in Degrees) between grid points, correspond to each dihedral, every value must be a divisor of 360

init_coords_M: geometric.molecule.Molecule() instance

A Molecule constains a series of initial geometries to start with

energy_decrease_thresh: Float

The threshold of the smallest energy decrease amount to trigger activating optimizations from grid point.

dihedral_ranges: List[(lower, upper), ..]

A list of dihedral range limits as a pair (lower, upper), each range corresponds to the dihedrals in input.

energy_upper_limit: Float or None

The threshold if the energy of a grid point that is higher than the current global minimum, to start new optimizations, in unit of a.u. i.e. if energy_upper_limit = 0.05, current global minimum energy is -9.9 , then a new task starting with energy -9.8 will be skipped.

extra_constraints: Dict

A nested dictionary specifing extra constraints in geomeTRIC format. Details in

verbose: bool

let methods print more information when running


Build a dihedral mask based on specified ranges

dihedral_ranges: List[(lower: Int, upper: Int), ..]

The range limits corresponding to each dihedral angle A full dihedral range is [-180, 180] The upper limit up to 360 is supported for the purpose of specifying range limits crossing the boundary, e.g. [80, 240], which effectively become [-180, 120] + [80, 180]

dihedral_mask: List[set(), ..]

The dihedral mask is a list of sets, each set contains all available values for one dihedral angle


This function should be called after self.setup_grid()


Create an empty tmp folder structure, save the paths for each grid point into self.tmp_folder_dict


self.tmp_folder_dict = {(30,-70): “opt_tmp/gid_+030_-070”, ..}


Return a string with ANSI colors showing current running status


Return a string of Ramachandran plot showing current running status


Write qdata.txt and file based on converged scan results

get_dihedral_id(molecule, check_grid_id=None)

Compute the closest grid ID for molecule (only first frame) If check_grid_id is given, will perform a check if the computed dihedral_values are close to the grid_id provided If the check is not passed, this function will return None


create a job scratch folder inside tmp_folder_dict[grid_id] name starting from ‘1’, and will use larger numbers if exist return the new folder name that’s been created


Take a center grid id, return all the neighboring grid ids, in all dimensions


Take a center grid id, return all the neighboring grid ids, in each dimension

launch_constrained_opt(molecule, grid_id)

Called by launch_opt_jobs() to launch one opt job in a new scr folder Return the new folder path


Launch constrained optimizations for molecules in opt_queue Tasks current opt_queue will be popped in order. If a task exist in self.task_cache, the cached result will be checked, then put into self.current_finished_job_results Else, the task will be launched by self.launch_constrained_opt, and information is saved as self.running_job_path_info[job_path] = m, from_grid_id, to_grid_id


The master function that calls all other functions. This function will run the following steps: 1. Launch a new set of jobs from self.opt_queue, add their job path to a dictionary 2. Check if any running job has finished 3. For each finished job, check if energy is lower than existing one, if so, add its neighbor grid points to opt_queue 4. Go back to the 1st step, loop until all jobs finished, indicated by opt_queue and running jobs both empty.


Push a set of initial tasks to self.opt_queue A task is defined as (m, from_grid_id, to_grid_id) tuple, where geometry is stored in m


Restore previous finished tasks from tmp folder. 1. Look into tmp folder and read scanner_settings.json, check if it matches current setting 2. Read the result pickle file from each leaf folder, into task_cache If successful, self.tmp_folder_dict will be initialized, same as self.create_tmp_folder(), and self.task_cache will be populated, with task caches, defined in this way:

self.task_cache = {(30,-60): {geo_key: (final_geo, final_energy, final_gradient, job_folder)}}

final_gradient will be None if it’s not available.

save_task_cache(job_path, m_init, m_final)

Save a file containing the finished job information to a pickle file on disk. The format should be consistent with self.restore_task_cache()


Set up grid ids, each as a tuple with size corresponding to grid dimension. i.e. 1-D: grid_ids = ( (-165, ), (-150, ), … (180, ) ) 2-D: grid_ids = ( (-165,-165), (-165,-150), … (180,180) ) This function is called by the initializer.

self.grid_axes is also initialized, to be a full range of grid values for each dihedral, i.e., 1-D: grid_axes = [range(-165, 195, 15)] 2-D: grid_axes = [range(-165, 195, 15), range(-165, 195, 15)]


Validate a constrained optimization task before pushing to the queue. This is useful to limit the dihedrals into a range of interest.

task: (m, from_grid_id, to_grid_id)

A constrained optimization task

isValid: bool

True if the task is valid


Interface with engine to check if any job finished. Will wait infinitely here until at least one job finished. The finished job paths will be removed from self.running_job_path_info. The finished job results (m, grid_id) will be checked, if the result geometry is not close enough to target grid id, the result will be ignored. Results passed the check will be added to self.current_finished_job_results.

torsiondrive.dihedral_scanner.cross3(v1, v2)

Quick convenient function to compute cross product betwee two 3-element vectors cross3: 326 ns | np.cross: 35.8 us

torsiondrive.dihedral_scanner.dot3(v1, v2)

Quick convenient function to compute dot product betwee two 3-element vectors dot3: 231 ns | 745 ns


Convert an numpy array of xyz coordinate to a hashable object, keeping 0.001 precision This function has the limitation that 3.1999 and 3.2000 will produce different results due to the limitation of float point representation.

torsiondrive.dihedral_scanner.measure_dihedrals(molecule, dihedral_list, check_linear=True, check_bonded=True)

Measure dihedral values from molecule coordinates.

molecule: geometric.molecule.Molecule

The molecule object that contains atom coordinates. Only the first frame will be used.

dihedral_list: List[List[Int]]

A list of dihedrals to compute their value. Each diedral is represented by a list of tuple of four integers, each is a 0-based atom index.

check_linear: Bool

If True, will check if i-j-k or j-k-l angles in each dihedral is close to linear ( > 165 degree ), print a warning if found.

check_bonded: Bool

If True, will check if all i-j, j-k, k-l are bonded for each dihedral, print a warning if not.


Quick convenient function to get the norm of a 3-element vector norm3: 475 ns | np.linalg.norm: 4.31 us


Normalize any number to the range (-180, 180], including 180


class torsiondrive.qm_engine.EngineBlank(input_file=None, work_queue=None, native_opt=False, extra_constraints=None)

A blank engine only used in testing


torsiondrive.extra_constraints.build_geometric_constraint_string(constraints_dict, dihedral_idx_values=None)

Build the geomeTRIC constraint string with constraints_dict and a set of dihedral_idx_values

constraints_dict: Dict

constraints dict built by make_constraints_dict() function

dihedral_idx_values: List[List[d1, d2, d3, d4, v]]

A list containing the definition of dihedrals and their values Example: [(0,1,2,3,90.0), (1,2,3,4,100.0)]

constraints_string: string

A string with multiple lines, to be used as the geomeTRIC constraints.txt

torsiondrive.extra_constraints.build_terachem_constraint_string(constraints_dict, dihedral_idx_values=None)

Build the TeraChem constraint string with constraints_dict and a set of dihedral_idx_values

constraints_dict: Dict

constraints dict built by make_constraints_dict() function

dihedral_idx_values: List[List[d1, d2, d3, d4, v]]

A list containing the definition of dihedrals and their values Example: [(0,1,2,3,90), (1,2,3,4,100)]

constraints_string: string

A string with multiple lines, to be used as the TeraChem constraints format

torsiondrive.extra_constraints.check_conflict_constraints(constraints_dict, dihedral_idxs)

Utility function to check if any extra constraints in constraints_dict is conflict with the scanning dihedrals


Create an ordered dictionary with constraints specification, consistent with geomeTRIC

constraints_string: str

String-formatted constraint specification consistent with geomeTRIC constraints.txt

constraints_dict: dict

A dictionary contains the definition of the extra constraints. The format is consistant with the JSON interface of geomeTRIC.


  1. Only constraints of type “freeze” and “set” are supported, since extra “scan” is undefined with torsiondrive scan.
  2. Four attributes are allowed to be constrained: ‘distance’, ‘angle’, ‘dihedral’, ‘xyz’
  3. The input string is one-indexed, the output dictionary is zero-indexed.
  4. For “xyz”, dashed inputs like “1-3,7-9” (no space) is allowed, and will be converted to [0,1,2,6,7,8].


>>> make_constraints_dict(r"$freeze\nxyz 1-3\n$set\nangle 2 1 5 110.0")
    'freeze': [{
        'type': 'xyz',
        'indices': [0, 1, 2]
    'set': [{
        'type': 'angle',
        'indices': [1, 0, 4],
        'value': 110.0}]