Structural computations

centerofmass(ta::TrjArray; isweight=true, index::Union{UnitRange,Vector,BitArray}) -> com::TrjArray

Calculates the center of mass coordinates (COM) of the given trajectory, a TrjArray object ta. If isweight is true (default), coordinates are weighted by ta.mass (as long as ta.mass is not empty). Uses can also specify column indices for the COM calculation by giving index which should be a Vector of integers, a UnitRange, or a BiArray.

Returns the center of mass coordinates as a TrjArray object, which have virtual single "atom" whose coordinates are COMs.

Example

julia> ta = mdload("ak.pdb")
julia> com = centersofmass(ta)

decenter(ta::TrjArray; isweight=true, index::Union{UnitRange,Vector,BitArray}) -> ta_decentered::TrjArray, com::TrjArray

Calculates the centers of mass coordinates (COM) of the given trajectory, a TrjArray object ta. And translates the coordinates so that the COM of the output is identical to the origin (x=y=z=0). If isweight is true (default), coordinates are weighted by ta.mass (as long as ta.mass is not empty). Uses can also specify column indices for the COM calculation by giving index which should be a Vector of integers, a UnitRange, or a BiArray.

Returns a TrjArray object whose COM is the origin (x=y=z=0), and the COM of the given TrjArray ta.

Example

julia> ta = mdload("ak.pdb")
julia> ta_decentered = decenter(ta)

decenter!(ta::TrjArray; isweight=true, index::Union{UnitRange,Vector,BitArray}) -> com::TrjArray

Translates the coordinates of the given trajectory, a TrjArray ta, so that the COM of it become identical to the origin (x=y=z=0). If isweight is true (default), coordinates are weighted by ta.mass (as long as ta.mass is not empty). Uses can also specify column indices for the COM calculation by giving index which should be an Array of integers.

Returns the center of mass coordinates of the given trajectory.

Example

julia> ta = mdload("ak.pdb")
julia> decenter!(ta)

orient!(ta::TrjArray, index::Union{UnitRange,Vector,BitArray})

Orient the molecule using its principal axes of inertia.

Example

julia> ta = mdload("ak.pdb")
julia> orient!(ta)

rotate(ta::TrjArray, quaternion::Vector) -> ta_rotated::TrjArray

Performs rotation of coordinates of the given TrjArray object ta. Rotation matrix is constructed fromt the given quaternion vector.

Returns the trajecotry with rotated coordinates.

Example

julia> ta = mdload("ak.dcd")
julia> quaternion = [0.80951112,  0.10657412,  0.35146912,  0.45804312]
julia> ta_rotated = rotate(ta, quaternion)

rotate!(ta::TrjArray, quaternion::Vector)

Performs rotation of coordinates of the given TrjArray object ta. Rotation matrix is constructed fromt the given quaternion vector.

The coordinates of ta is replaced with the rotated ones.

Example

julia> ta = mdload("ak.dcd")
julia> quaternion = [0.80951112,  0.10657412,  0.35146912,  0.45804312]
julia> rotate!(ta, quaternion)

rotate(ta::TrjArray, quaternions::Matrix) -> ta_rotated::TrjArray

Performs multiple rotations of the coordinates of the 1st frame of the given TrjArray object ta. Rotation matrices are construcuted from the 2-dimensional Array quaternions. One set of quaternions for rotaton should be given in a row vector in quaternions.

Returns the trajecotry with rotated coordinates from the 1st frame of ta.

Example

julia> ta = mdload("ak.pdb")
julia> quaternions = [0.80951112  0.10657412  0.35146912  0.45804312; 
>                     0.10657412  0.80951112  -0.35146912  0.45804312]
julia> ta_rotated = rotate(ta, quaternions)

rotate_with_matrix(ta::TrjArray, R::Matrix) -> ta_rotated::TrjArray

Performs rotation of coordinates of the given TrjArray object ta with the given rotation matrix.

The coordinates of ta is replaced with the rotated ones.

Example

julia> ta = mdload("ak.dcd")
julia> R = [0.36 0.48 -0.8;
            -0.8 0.60 0.0;
            -.48 0.64 0.60]
julia> ra_rotated = rotate(ta, R)

superimpose(ref::TrjArray, ta::TrjArray; isweight=true, index::Union{UnitRange,Vector,BitArray}, isdecenter=false) -> ta_superimposed::TrjArray

Performs the least-squares fitting of the given trajectory (TrjArray object ta) to the reference structure (1st frame of TrjArray object ref). If isweight is true (default), coordinates are weighted by ta.mass (as long as ta.mass is not empty). Uses can specify column indices for the COM calculation by giving index which should be a Vector of integers, a UnitRange, or a BiArray. When isdecenter is true, the function assumes the COMs of both structures are located at the origin (default is isdecenter=false).

Returns the superimposed trajecotry.

The code of this function is licensed under the BSD license (Copyright (c) 2009-2016 Pu Liu and Douglas L. Theobald), see LICENSE.md

Example

julia> ref = mdload("ak.pdb")
julia> ta = mdload("ak.dcd", top=pdb)
julia> ta_superimposed = superimpose(ref, ta)

References

The algorithm of this function is based on 
P. Liu, D. K. Agrafiotis, and D. L. Theobald, 
J. Comput. Chem. 31, 1561-1563 (2010).

compute_rmsd(ref::TrjArray, ta::TrjArray; isweight=true, index=index::Union{UnitRange,Vector,BitArray}) -> rmsd

Calculates the root mean square deviations (RMSD) of the given trajectory (TrjArray object ta) from the reference structure (1st frame of TrjArray object ref). If isweight is true (default), coordinates are weighted by ta.mass (as long as ta.mass is not empty). Uses can specify column indices for the RMSD calculation by giving index which should be a Vector of integers, a UnitRange, or a BiArray.

Returns RMSD values.

Example

julia> ref = mdload("ak.pdb")
julia> ta = mdload("ak.nc", top=pdb)
julia> ta_superimposed = superimpose(ref, ta)
julia> rmsd = compute_rmsd(ref, ta_superimposed)

meanstructure(ta::TrjArray; isweight=true, index=index::Union{UnitRange,Vector,BitArray}) -> ta_mean::TrjArray, ta_superimposed::TrjArray

Calculates the mean structure of the given trajectory (TrjArray object ta) by iteratively fitting the trajecotry to tentative mean structures until the mean structure converges. If isweight is true (default), coordinates are weighted by ta.mass (as long as ta.mass is not empty). Uses can specify column indices for the fitting (superimpose function) by giving index which should be a Vector of integers, a UnitRange, or a BiArray.

Returns the mean structure as TrjArray object, and the supserimposed trajectory to the mean structure.

Example

julia> ta = mdload("ak.nc", top=pdb)
julia> ta_mean, ta_superimposed = meanstructure(ta)

compute_rmsf(ta::TrjArray{T, U}; isweight::Bool=true) -> rmsf

rmsf (root mean square fluctuation)

Calculates the root mean square fluctuations (RMSFs) of the atoms in the given trajectory (TrjArray object ta). from the average coordinates of the trajectory. If isweight is true (default), coordinates are weighted by ta.mass (as long as ta.mass is not empty).

Returns RMSF values.

Example

julia> ta = mdload("ak.nc", top=pdb)
julia> ta_mean, ta_superimposed = meanstructure(ta)
julia> rmsf = compute_rmsf(ta)

compute_distancemap(ta::TrjArray; kneighbor=3)

Calculates distance-map vectors from the given TrjArray ta. If atomid=j and atomid=i are closer to each other than the specified kneighbor, i.e., abs(i-j) < kneighbor, the calculation of the distance is ignored, considering that they are not nonbonded pairs. The default value of kneighbor is 3.

Returns distance-map vectors of all frames. The distance-map of each frame is contained in a row vector of the output.

Example

julia> ta = mdload("ak.pdb")
julia> ta = mdload("ak.dcd", top=ta)
julia> d = compute_distancemap(ta["atomname CA"])

compute_contactmap(ta::TrjArray; rcut=8.5, kneighbor=3)

Calculates contact-map vectors from the given TrjArray ta. Cut-off value for contact can be given in rcut (default is rcut=8.5). If atomid=j and atomid=i are closer to each other than the specified kneighbor, i.e., abs(i-j) < kneighbor, the calculation of the distance is ignored, considering that they are not nonbonded pairs. The default value of kneighbor is 3.

Returns contact-map vectors of all frames. The contact-map of each frame is contained in a row vector of the output.

Example

julia> ta = mdload("ak.pdb")
julia> ta = mdload("ak.dcd", top=ta)
julia> d = compute_distancemap(ta["atomname CA"])

compute_angle(ta1::TrjArray, ta2::TrjArray, ta3::TrjArray)

Calculates angles for the atom coordinates or the COMs of groups from the triplet of ta1, ta2 and ta3.

Returns angles.

Example

julia> ta = mdload("ak.pdb")
julia> ta = mdload("ak.dcd", top=ta)
julia> a = compute_distance(ta[:, "atomid 1"], ta[:, "atomid 9"], ta[:, "atomid 11"])

compute_dihedral(ta1::TrjArray, ta2::TrjArray, ta3::TrjArray, ta4::TrjArray)

Calculates dihedral angles for the atom coordinates or the COMs of groups from the quadruplet of ta1, ta2, ta3 and ta4.

Returns dihedral angles.

Example

julia> ta = mdload("ak.pdb")
julia> ta = mdload("ak.dcd", top=ta)
julia> d = compute_dihedral(ta[:, "atomid 1"], ta[:, "atomid 9"], ta[:, "atomid 11", ta[:, "atom 13"]])

compute_dihedral(ta1::TrjArray, ta2::TrjArray, ta3::TrjArray, ta4::TrjArray)

Calculates dihedral angles for the atom coordinates or the COMs of groups from the quadruplet of ta1, ta2, ta3 and ta4.

Returns dihedral angles.

Example

julia> ta = mdload("ak.pdb")
julia> ta = mdload("ak.dcd", top=ta)
julia> d = compute_dihedral(ta[:, "atomid 1"], ta[:, "atomid 9"], ta[:, "atomid 11", ta[:, "atom 13"]])

compute_qscore(native::TrjArray, ta::TrjArray) -> qscore::Array

Calculates all-atom based Q-score from given heavy atom coordinates.

Returns Q-scores.

Example

julia> native = mdload("ak.pdb")
julia> ta = mdload("ak.dcd", top=native)
julia> qscore = compute_qscore(native["not hydrogen"], ta["not hydrogen"])

References

Definition of Q-score from heavy atoms is given in 
R. B. Best, G. Hummer, and W. A. Eaton, PNAS 110, 17874 (2013).

compute_drms(native1::TrjArray, native2::TrjArray, ta1::TrjArray, ta2::TrjArray) -> drms

Calculates distance-based root mean square displacements (DRMS) from the given trajectories. First, native contacts between molecule1 and molecule2 are identified from the two native structures native1 and native2 where cutoff distance of 6.0 Angstrom is used. Then, the DRMS are calculated from the trajectories of molecule1 ta1 and molecule2 ta2.

Returns distance-based root mean square displacements.

Example

julia> native1 = mdload("protein.pdb")["not hydrogen"]
julia> native2 = mdload("ligand.pdb")["not hydrogen"]
julia> ta1 = mdload("protein.dcd")["not hydrogen"]
julia> ta2 = mdload("ligand.dcd")["not hydrogen"]
julia> drms = compute_drms(native1, native2, ta1, ta2)

References

Definition of DRMS is given in 
J. Domański, G. Hedger, R. B. Best, P. J. Stansfeld, and M. S. P. Sansom, 
Convergence and Sampling in Determining Free Energy Landscapes for Membrane Protein Association, 
J. Phys. Chem. B 121, 3364 (2017).

compute_pairlist(ta::TrjArray, rcut; iframe=1::Int) -> F

Make a pairlist for the given strcuture ta1 by searching pairs within a cutoff distance rcut. By default, the 1st frame of ta1 is used. Users can specify the frame by iframe.

Returns a NamedTuple object F which contains the pair lists in F.pair, and the distances of corresponding pairs in F.dist.

Example

julia> ta = mdload("ak.pdb")
julia> F = compute_pairlist(ta, 8.0)

References

The algorithm of this function is based on 
T.N. Heinz, and P.H. Hünenberger, 
A fast pairlist-construction algorithm for molecular simulations under periodic boundary conditions. 
J Comput Chem 25, 1474–1486 (2004). 

compute_pairlist_bruteforce(ta::TrjArray, rcut; iframe=1::Int) -> F

Bruteforce search version of compute_pairlist function. Mainly used for debugging.

Returns a NamedTuple object F which contains the pair lists in F.pair, and the distances of corresponding pairs in F.dist.

Example

julia> ta = mdload("ak.pdb")
julia> F = compute_pairlist(ta, 8.0)