pyPRISM.trajectory.Debyer module

class pyPRISM.trajectory.Debyer.Debyer(domain, nthreads=None)

Bases: object

Parallelized Debye method to calculate \(\hat{\omega}_{\alpha,\beta}(k)\)

Warning

This pyPRISM functionality is still under testing. Use with caution.

Mathematical Definition

\[\hat{\omega}_{\alpha,\beta}(k) = \delta_{\alpha,\beta} + \frac{C_{\alpha,\beta}}{N_{frame}} \sum_f \sum_{i,j} \frac{\sin(kr^{f}_{ij})}{kr^{f}_{ij}}\]

Variable Definitions

  • \(\hat{\omega}_{\alpha,\beta}(k)\)
    Intra-molecular correlation function in Fourier-space
  • \(\delta_{\alpha,\beta}\)
    Kronecker delta for when considering a self (\(\alpha == \beta\)) versus not-self (\(\alpha /= \beta\)) site type pair.
  • \(N_{frame}\)
    Number of frames in simulation trajectory
  • \(C_{\alpha,\beta}\)
    Scaling coefficient. If \(\alpha == \beta\), then \(C_{\alpha,\beta} = 1/N_{\alpha}\) else \(C_{\alpha,\beta} = 1/(N_{\alpha} + N_{\beta})\)
  • \(r^{f}_{i,j}\)
    At frame \(f\) of simulation, the scalar distance between sites with index i and j.
  • \(k\)
    Wavenumber of calculation.

Description

One of the most powerful uses of PRISM is to combine it with modern simulation techniques to calculate the intra-molecular correlation functions \(\hat{\omega}(k)\). This allows PRISM to be used for systems which do not have analytical descriptions for their \(\hat{\omega}(k)\) and, furthermore, allows PRISM to predic the structure of non-ideal systems.

Unfortunately, this calculation is extremely computationally intensive and requires care to calculate correctly. Here we provide a parallelized implementation of the Debye Method which can be used to calculate \(\hat{\omega}(k)\) for small to medium sized simulations.

This method works by allowing the user to selectively provide two sets of coordinate trajectories. To calculate \(\hat{\omega}(k)\) between site-types \(\alpha\) and \(\beta\), the user should pass one trajectory of all sites of type \(\alpha\) and the other where all sites are of type \(\beta\). To calculate a self \(\hat{\omega}(k)\) where \(\alpha == \beta\) (selfOmega=*True*), both trajectories should be the same. Note that selfOmega should be correctly set in either case.

The calculate method below takes six arguments: positions1, positions2, molecules1, molecules2, box, selfOmega.

The positions1 and positions2 arguments are numpy array containing multiple frames of coordinates. Each of the two arrays should only contain the coordinates of the site-type pair being considered. In other words, positions1 should have a trajectory of positions for site type \(\alpha\) and positions2 for \(\beta\) when calculating \(\hat{\omega}_{\alpha,\beta}(k)\). If \(\alpha == \beta\), the positions1 and positions2 should be the same array.

The molecules1 and molecules2 arguments specify to which molecule each site belongs. These arrays are necessary for calculations using trajectories that contain multiple molecules. For each site in the simulation of a site-type, these lists have an integer index which specify which molecule the site belongs to. All sites which share the same molecular index in this array belong to the same molecue. For single molecule simulations, an array of zeros or ones can be specified.

See the example below for an example of this classes use.

Note

Note that the num_chunks argument does not set the number of threads used in the calculation. This is governed by the OMP_NUM_THREADS (at least on *nix and OSX). If you have many cores on your machine, you may actually see improved performance if you reduce the value of OMP_NUM_THREADS from its maximum. Setting this variable is also useful if you’re sharing a machine or node. If running a script, you can set this variables globally:

$ export OMP_NUM_THREADS=4

or locally for a single execution

$ OMP_NUM_THREADS=4 python pyPRISM_script.py

Note

As currently written, this class assumes periodicity in all dimensions.

Example

Below we consider an arbitrary system that has at least two types of sites in it.

import pyPRISM
import numpy as np

# load position, type, and molecule information using a users method of
# choice. For the array sizes, F = number of frames from simulation, N =
# total number of atoms/beads/sites
positions = ...  # 3D Array of size (F,N,3)
types = ...      # 1D Array of size (N)
molecules = ...  # 1D Array of size (N)
box = ...        # 2D Array of size (F,3)

# create pyPRISM domain
domain = pyPRISM.Domain(dr = 0.25, length = 1024)

# create Debyer object
debyer = pyPRISM.trajectory.Debyer(domain=domain,nthreads=4)

# calculate omega_1_1
mask1 = (types == 1)
positions1 =  positions[:,mask1,:]
molecules1 =  molecules[mask1]
mask2 = (types == 1)
positions2 =  positions[:,mask2,:]
molecules2 =  molecules[mask2]
selfOmega = True
omega_1_1  = debyer.calculate(positions1, positions2, molecules1, molecules2, box, selfOmega)

# calculate omega_1_2
mask1 = (types == 1)
positions1 =  positions[:,mask1,:]
molecules1 =  molecules[mask1]
mask2 = (types == 2)
positions2 =  positions[:,mask2,:]
molecules2 =  molecules[mask2]
selfOmega = False
omega_1_2  = debyer.calculate(positions1, positions2, molecules1, molecules2, box, selfOmega)
calculate(self, positions1, positions2, molecules1, molecules2, box, selfOmega)

Calculate omega from two site type trajectories

Warning

The first five arguments of this method (positions1, positions2, molecules1, molecules2, box) must be numpy arrays (not Python lists) for this method to work.

Parameters:
  • positions1 (np.ndarray, float, size (\(F,N_{\alpha}\),3)) – 3-D numpy array containing coordinate trajectory from a simulation for site type \(\alpha\). See above for description.
  • positions2 (np.ndarray, float, size (\(F,N_{\beta}\),3)) – 3-D numpy array containing coordinate trajectory from a simulation for site type \(\beta\). See above for description.
  • molecules1 (np.ndarray, float, size (\(N_{\alpha}\))) – 1-D numpy arrays which specify the molecular identity of a site for sites of type \(\alpha\). See above for a description.
  • molecules2 (np.ndarray, float, size (\(N_{\beta}\))) – 1-D numpy arrays which specify the molecular identity of a site for sites of type \(\beta\). See above for a description.
  • box (np.ndarray, float, size(3)) – 1-D array of box dimensions, \(l_{x}, l_{y}, l_{z}\)
  • selfOmega (bool) – Set to True if \(\alpha == \beta\) and False otherwise.
Returns:

omega – Intra-molecular correlation function
Return type:

np.ndarray