#!/bin/python ############################################################################ # LICENSE INFO: # ############################################################################ # This file is part of CAMPARI. # # # # Version 5.0 # # # # Copyright (C) 2024, The CAMPARI development team (current and former # # contributors) # # Andreas Vitalis, Adam Steffen, Rohit Pappu, Hoang # # Tran, Albert Mao, Xiaoling Wang, Jose Pulido, # # Nicholas Lyle, Nicolas Bloechliger, Marco Bacci, # # Davide Garolini, Jiri Vymetal, Cassiano Langini, # # Alessia Del Panta, Francesco Cocina, Julian # # Widmer # # # # Website: http://sourceforge.net/projects/campari/ # # # # CAMPARI is free software: you can redistribute it and/or modify # # it under the terms of the GNU General Public License as published by # # the Free Software Foundation, either version 3 of the License, or # # (at your option) any later version. # # # # CAMPARI is distributed in the hope that it will be useful, # # but WITHOUT ANY WARRANTY; without even the implied warranty of # # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # # GNU General Public License for more details. # # # # You should have received a copy of the GNU General Public License # # along with CAMPARI. If not, see . # ############################################################################ # AUTHORSHIP INFO: # ############################################################################ # # # MAIN AUTHOR: Andreas Vitalis # # # ############################################################################ import numpy as np from scipy.spatial.transform import Rotation as spRot from scipy import version as spVersion ############################################################################################## # predefined classes ############################################################################################## # This is a general-purpose class used for static objects. # The dictionary (col_dic) can be used to see the available quantities, which are also listed in the documentation. class props_holder: def __init__(self,arg1,arg2): self.col_dic = {} colns = arg2.split(',') cnt = 0 for coln in colns: self.col_dic[coln.strip()] = cnt cnt = cnt + 1 self.theprops = arg1 self.thens = 0 if (len(colns) > 0): self.thens = len(self.theprops[:,1]) # minimal wrapper of the constructor def init_static_props(arg1,arg2): return props_holder(arg1,arg2) # This is a specialized class used only to encapsulate the simulation container. class box_holder: def __init__(self,arg1,arg2,arg3,arg4): self.shape = arg1 # as in CAMPARI (1-4) self.type = arg2 # as in CAMPARI (1-4) self.periodicity = arg3 # as in PSQL CAMPARI convention (0-7) self.volume = arg4[0] # A^3 if ((self.shape == 2) or (self.shape == 3)): self.origin = arg4[1:4] # A self.radius = arg4[4] # A if (self.shape == 3): self.height = arg4[6] # A self.softwall = arg4[7] # kcal/mol/A^2 self.bmatrix = np.zeros(9) self.ibmatrix = np.zeros(9) self.bmatrix[range(0,9,4)] = self.radius*2 self.ibmatrix[range(0,9,4)] = 0.5/self.radius if (self.shape == 3): self.bmatrix[8] = self.height self.ibmatrix[8] = 1.0/self.height self.minimage2 = arg4[29]*arg4[29] # A^2, faked else: self.origin = arg4[4:7] # A self.sidelengths = arg4[1:4] # A self.cosangles = arg4[8:11] self.bmatrix = arg4[11:20] # A self.ibmatrix = arg4[20:29] # 1/A self.minimage2 = arg4[29]*arg4[29] # A^2 # This is a class that is used as an aggregator for most of the properties that might not always be available, mostly analysis quantities. # Some of these are cumulative (and would need to be diffed to recover instantaneous quantities), others are instantaneous (and would need to be binned or averaged). # Not all quantities are transmitted in this way: some are easily regenerated because the raw data are already there (like Ramachandran maps), others are easily regenerated with the help of pre-built functions (like the data in POLYAVG.dat). # There is no forced synchronization between FMCSC_PYCALC and other step intervals (like FMCSC_DSSPCALC, etc.), so it is completely up to the user to set these meaningfully; similarly, there is no bookkeeping of the numbers one would require to average (relevant especially in case of lumped analyses (by group)). # Lastly, what these contain when passed to "campana_final" should not be relied upon. class analysis_holder: def __init__(self): self.holder = {} self.holder['contact_map'] = None # cumulative contact map [i,j] index i>j is minimum pairwise distance, index j>i is center-of-mass based, residues are all those from molecules with nonzero "SOLUTE_NR" in static_mol_ints self.holder['contact_clusters'] = None # instantaneous connectivity of molecules with nonzero "SOLUTE_NR" in static_mol_ints; column 0 is # links per clusters, column 1 is size per cluster, column 2 is map of solute index to cluster index self.holder['dssp_assignment'] = None # the list of instantaneous DSSP assignment per eligible residue (can be retrieved from the "PEP_RS" column in static_res_ints), meaning same as for "DSSP_RES.dat" self.holder['dssp_d_hb_list'] = None # the list of partners j for which i (position in the array) is the donor in an H-bond, -1 if none self.holder['dssp_a_hb_list'] = None # the list of partners j for which i (position in the array) is the acceptor in an H-bond, -1 if none self.holder['dssp_d_hb_energies'] = None # energies in kcal/mol matched to dssp_d_hb_list, 0.0 if not actually populated self.holder['dssp_a_hb_energies'] = None # energies in kcal/mol matched to dssp_a_hb_list, 0.0 if not actually populated self.holder['em_density'] = None # instantaneous spatial density 4D array (first 3 are xyz lattice, 4th is property, usually just of size 1), quantities are instantaneous and normalized to density self.holder['em_parameters'] = None # 9-vector of lattice parameters for spatial density (spacing, origin, size (as real)), each in xyz self.holder['pol_persistence'] = None # instantaneous angular correlation function across analysis groups and sequence spacing self.holder['pol_density_profile'] = None # instantaneous, unnormalized radial density profile across analysis groups and Rg-binning (FMCSC_RGBINSIZE and FMCSC_POLRGBINS), needs volume normalization self.holder['pol_rg_hist'] = None # cumulative radius of gyration histogram across analysis groups and Rg-binning (FMCSC_RGBINSIZE and FMCSC_POLRGBINS) self.holder['pol_ree_hist'] = None # cumulative end-to-end distance histogram across analysis groups and Rg-binning (FMCSC_RGBINSIZE and FMCSC_POLRGBINS) self.holder['pol_rg_delta_hist'] = None # cumulative 2D histogram of t (normalzied Rg) and asphericity delta* across analysis groups and fixed binning self.holder['pol_inverse_rh'] = None # instantaneous pairwise internal-distance based estimate of inverse radius per analysis group (1/A) self.holder['pol_internal_scaling'] = None # instantaneous internal distance scaling (A) across analysis groups and sequence spacing self.holder['pol_scattering'] = None # instantaneous scattering profile (Kratky plot) across analysis groups and wave vectors (FMCSC_SCATTERVECS and FMCSC_SCATTERRES) self.holder['seg_z_alpha'] = None # values of global secondary structure measure Z-alpha (alpha-content) per molecule (meaningless for non-peptides) self.holder['seg_z_beta'] = None # values of global secondary structure measure Z-beta (beta-content) per molecule (meaningless for non-peptides) self.holder['seg_sec_struct_perrs'] = None # cumulative statistics of residues across 8 segments types (col. 0), across eligible residues (col. 1), and across segment length (col. 2), see doc. for output file BB_SEGMENTS_RES.dat self.holder['pc_center_of_mass'] = None # cumulative histogram of centers of mass distances by pairs of analysis groups (col.s 0-1) and across distance bins (FMCSC_PCBINSIZE) in col. 2 self.holder['pc_user_selection'] = None # cumulative histogram across numbered distance groups (FMCSC_PCCODEFILE) in col. 0 and across distance bins (FMCSC_PCBINSIZE) in col. 1 self.holder['ia_energies'] = None # energies as requested; column 0: total, columns 1-X: individual terms as in ENERGY.dat from column 2, last column: frame weight (read-in, if not, then ELS, if not either, then just 1.0) self.holder['ia_cart_forces'] = None # Cartesian forces by atom and xyz, the same minor caveats apply as explained for FMCSC_CHECKGRAD self.holder['ia_imd_explain'] = None # Internal coordinate space data structure: column 0 is its index in a mol. list (1-3: translation, 4-6: rotation), column 1 is the number of atoms that move, columns 2-5 are the atom indices defining dihedral angles (0 for RB-coordinates) self.holder['ia_imd_forces'] = None # Vector of internal coordinate forces in kcal/mol/A or kcal/mol(/rad): explained by and matched to ia_imd_explain (lines up with first dimension) self.holder['ia_imd_masses'] = None # Vector of internal coordinate intertial masses in a.m.u. or a.m.u.*A^2: explained by and matched to ia_imd_explain (lines up with first dimension) self.holder['sol_accessible_volume'] = None # atom-wise solvent-accessible volume as defined by ABSINTH (responds to many atomic parameters and two globals: FMCSC_SAVPROBE, FMCSC_SAVMODE) self.holder['sol_solv_state_fos'] = None # atom-wise solvation state as mapped from solvent-accessible volumes for DMFI (see FMCSC_FOSFUNC to find parameters of relevance) self.holder['sol_solv_state_scr'] = None # atom-wise solvation state as mapped from solvent-accessible volumes for charge screening (see FMCSC_SCRFUNC to find parameters of relevance) self.holder['mpi_info'] = None # vector of MPI rank (starting at 0), total ranks, and number of outgoing transitions by rank (includes all ranks, even self), number of incoming transitions by rank (ditto); the latter nonzero only in actual replica exchange runs self.holder['clu_features'] = None # Current batch of clustering feature extraction: first dimension is feature set, second dimension is batch (usually 1 for serial runs or, per replica, for replica exchange-type runs; for MPI averaging / PIGS runs, batches do not come in at every step necessarily and are always ordered in the second dimension by replica (outer) and by batch (inner) # These are two utility functions for the analysis_holder class def init_analysis_prop(self,arg1,arg2): if (isinstance(arg2, np.ndarray)): ndims = len(arg2.shape) if (ndims == 1): print("campana.py: Analysis holder object intializes",arg1," to a 1D array of length ",len(arg2)," and data type ",type(arg2[0])) elif (ndims == 2): print("campana.py: Analysis holder object intializes",arg1," to a 2D array of shape ",arg2.shape," and data type ",type(arg2[0,0])) elif (ndims == 3): print("campana.py: Analysis holder object intializes",arg1," to a 3D array of shape ",arg2.shape," and data type ",type(arg2[0,0,0])) elif (ndims == 4): print("campana.py: Analysis holder object intializes",arg1," to a 4D array of shape ",arg2.shape," and data type ",type(arg2[0,0,0,0])) elif (ndims == 5): print("campana.py: Analysis holder object intializes",arg1," to a 5D array of shape ",arg2.shape," and data type ",type(arg2[0,0,0,0,0])) self.holder[arg1] = arg2 def update_analysis_prop(self,arg1,arg2): if (isinstance(self.holder[arg1], np.ndarray)): ndims = len(self.holder[arg1].shape) if (ndims == 1): self.holder[arg1][:] = arg2 elif (ndims == 2): self.holder[arg1][:,:] = arg2 elif (ndims == 3): self.holder[arg1][:,:,:] = arg2 elif (ndims == 4): self.holder[arg1][:,:,:,:] = arg2 elif (ndims == 5): self.holder[arg1][:,:,:,:,:] = arg2 else: self.holder[arg1][:] = arg2 else: self.holder[arg1] = arg2 # redefinition, different from writing to preexisting memory location ############################################################################################## # predefined utility functions ############################################################################################## # # calculation of distances (vector_calc_dis2) ############################################################################################## # Arguments: # bsh should be the container shape as in box.shape # bty should be the CAMPARI container type as in box.type # bmat and ibmat should be the regular and inverted matrices in vector form as in box.bmatrix and box.ibmatrix # (must be there but possibly unused) # bminim2 should be the square of the shortest possible image distance # s1x, s1y, s1z should be the scalars defining the 3-vector to use as reference # v1x, v1y, v1z should be vectors encoding positions to compute distances from # ignore_images can be used to not check for images (this is what CAMPARI does for atoms in the same molecule) # the return value is a NumPy nd-array of squared distances # if "svecx" "svecy" and "svecz" are passed as named arguments pointing to correctly sized nd-arrays arrays, # the shift vectors can also be retrieved def vector_calc_dis2(bsh, bty, bmat, ibmat, bminim2, s1x, s1y, s1z, v2x, v2y, v2z, ignore_images, **kwargs): dvx = v2x-s1x dvy = v2y-s1y dvz = v2z-s1z ndvx = None ndvy = None ndvz = None if ignore_images: d2 = dvx*dvx + dvy*dvy + dvz*dvz return d2 if ((bty == 1) and (bsh == 4)): # for the triclinic, fully periodic case, note that CAMPARI sanity-checks against geometries that would # require a second layer of images to consider the nearest one # the following is slow and inefficient but reasonably safe # first make sure to map to unit cell fractional coordinates fracx = dvx*ibmat[0] + dvy*ibmat[3] + dvz*ibmat[6] fracy = dvx*ibmat[1] + dvy*ibmat[4] + dvz*ibmat[7] fracz = dvx*ibmat[2] + dvy*ibmat[5] + dvz*ibmat[8] fracx = fracx - np.rint(fracx) fracy = fracy - np.rint(fracy) fracz = fracz - np.rint(fracz) ndvx = fracx*bmat[0] + fracy*bmat[3] + fracz*bmat[6] ndvy = fracx*bmat[1] + fracy*bmat[4] + fracz*bmat[7] ndvz = fracx*bmat[2] + fracy*bmat[5] + fracz*bmat[8] mdd2 = ndvx*ndvx + ndvy*ndvy + ndvz*ndvz # the value in mdd2 is the naive but sometimes wrong answer to finding nearest images in triclinic cells # it works safely if and only if mdd2 is smaller than half the shortest image distance (or if the box # has other simplifying properties like being rectangular) l1 = (mdd2 > bminim2) nrechecks = sum(l1) if (nrechecks > 0): svecx = np.zeros(nrechecks) svecy = np.zeros(nrechecks) svecz = np.zeros(nrechecks) for i in range(-1,2): for j in range(-1,2): for k in range(-1,2): shf = i*bmat[0:3] + j*bmat[3:6] + k*bmat[6:9] dvtmpx = ndvx[l1] + shf[0] dvtmpy = ndvy[l1] + shf[1] dvtmpz = ndvz[l1] + shf[2] dd2tmp = dvtmpx*dvtmpx + dvtmpy*dvtmpy + dvtmpz*dvtmpz l2 = (dd2tmp < mdd2[l1]) if (sum(l2) > 0): mdd2[l1][l2] = dd2tmp[l2] svecx[l2] = shf[0] svecy[l2] = shf[1] svecz[l2] = shf[2] ndvx[l1] = ndvx[l1] + svecx ndvy[l1] = ndvy[l1] + svecy ndvz[l1] = ndvz[l1] + svecz elif ((bty == 1) and (bsh == 1)): # in a rectangular box, minimizing fractional coordinates is enough bsides = bmat[range(0,9,4)] ndvx = dvx - bsides[0]*np.rint(dvx/bsides[0]) ndvy = dvy - bsides[1]*np.rint(dvy/bsides[1]) ndvz = dvz - bsides[2]*np.rint(dvz/bsides[2]) elif ((bty == 1) and (bsh == 3)): # periodic cylinder has only the z-axis to consider and behaves like a rectangular box ndvz = dvz - bmat[8]*np.rint(dvz/bmat[8]) elif ((bty == 2) or (bty == 3) or (bty == 4)): ok = True # all other cases require no image corrections else: raise Exception("Unsupported choice of container shape and/or boundary in vector_calc_dis2. Please check your Python code.") if (ndvx is None) and (ndvy is None) and (ndvz is None): d2 = dvx*dvx + dvy*dvy + dvz*dvz elif (ndvx is None) and (ndvy is None): d2 = dvx*dvx + dvy*dvy + ndvz*ndvz else: d2 = ndvx*ndvx + ndvy*ndvy + ndvz*ndvz if ("svecx" in kwargs) and (ndvx is not None): kwargs["svecx"][:] = ndvx-dvx if ("svecy" in kwargs) and (ndvy is not None): kwargs["svecy"][:] = ndvy-dvy if ("svecz" in kwargs) and (ndvz is not None): kwargs["svecz"][:] = ndvz-dvz return d2 ############################################################################################## # # calculation of dihedral angles without image correction (scalar_calc_dih_by_index) ############################################################################################## # Arguments: # i1-i4 should be the indeces of the 4 atoms (in sequence) # coords should be at least the first 3 columns (xyz) of dyn_atom as passed to campana_* def scalar_calc_dih_by_index(i1,i2,i3,i4,coords): d1 = coords[i2,0:3] - coords[i1,0:3] d2 = coords[i3,0:3] - coords[i2,0:3] d3 = coords[i4,0:3] - coords[i3,0:3] cp1 = np.cross(d1,d2) cp2 = np.cross(d2,d3) np1 = np.sqrt(np.sum(cp1*cp1)) np2 = np.sqrt(np.sum(cp2*cp2)) if (np1*np2 > 0): dp12 = np.sum(cp1*cp2) si = 1.0 if (dp12 < 0.0): si = -1.0 d4 = cp2/np2 - (si/np1)*cp1 np4 = 0.5*np.sqrt(np.sum(d4*d4)) rval = si*2.0*np.arcsin(np4) if (dp12 < 0.0): rval += np.pi si = np.sum(d1*cp2) if (si < 0.0): rval = -rval return np.rad2deg(rval) else: # ill-defined (aka arbitrary) return 0.0 ############################################################################################## # # calculation of dihedral angles without image correction (scalar_calc_dih_by_vector) ############################################################################################## # Arguments: # p1-p4 should be the 3-vectors of 4 positions (in sequence) # coords should be at least the first 3 columns (xyz) of dyn_atom as passed to campana_* def scalar_calc_dih_by_vector(p1,p2,p3,p4): d1 = p2 - p1 d2 = p3 - p2 d3 = p4 - p3 cp1 = np.cross(d1,d2) cp2 = np.cross(d2,d3) np1 = np.sqrt(np.sum(cp1*cp1)) np2 = np.sqrt(np.sum(cp2*cp2)) if (np1*np2 > 0): dp12 = np.sum(cp1*cp2) si = 1.0 if (dp12 < 0.0): si = -1.0 d4 = cp2/np2 - (si/np1)*cp1 np4 = 0.5*np.sqrt(np.sum(d4*d4)) rval = si*2.0*np.arcsin(np4) if (dp12 < 0.0): rval += np.pi si = np.sum(d1*cp2) if (si < 0.0): rval = -rval return np.rad2deg(rval) else: # ill-defined (aka arbitrary) return 0.0 ############################################################################################## # # calculation of bond angles without image correction (scalar_calc_bang_by_index) ############################################################################################## # Arguments: # p1-p3 should be the 3-vectors of 3 positions (in sequence) def scalar_calc_bang_by_index(i1,i2,i3,coords): d1 = coords[i1,0:3] - coords[i2,0:3] d2 = coords[i3,0:3] - coords[i2,0:3] nd1 = np.sqrt(np.sum(d1*d1)) nd2 = np.sqrt(np.sum(d2*d2)) if ((nd1 <= 0.0) or (nd2 <= 0.0)): # use a negative value to indicate superposed points return -1.0 dotp = np.sum(d1*d2) si = 1.0 if (dotp < 0): si = -1.0 d3 = d2/nd2 - si*d1/nd1 nd3 = np.sqrt(np.sum(d3*d3)) rval = si*2.0*np.arcsin(0.5*nd3) if (dotp < 0.0): rval += np.pi return np.rad2deg(rval) ############################################################################################## # # calculation of bond angles without image correction (scalar_calc_bang_by_vector) ############################################################################################## # Arguments: # p1-p3 should be the 3-vectors of 3 positions (in sequence) def scalar_calc_bang_by_vector(p1,p2,p3): d1 = p1 - p2 d2 = p3 - p2 nd1 = np.sqrt(np.sum(d1*d1)) nd2 = np.sqrt(np.sum(d2*d2)) if ((nd1 <= 0.0) or (nd2 <= 0.0)): # use a negative value to indicate superposed points return -1.0 dotp = np.sum(d1*d2) si = 1.0 if (dotp < 0): si = -1.0 d3 = d2/nd2 - si*d1/nd1 nd3 = np.sqrt(np.sum(d3*d3)) rval = si*2.0*np.arcsin(0.5*nd3) if (dotp < 0.0): rval += np.pi return np.rad2deg(rval) ############################################################################################## # # calculation of the diagonalized gyration tensor of a group of atoms (gyration_tensor_by_bool) ############################################################################################## # Arguments: # which should be the full-size True/False vector of which atoms to include # coords should be at least the first 3 columns (xyz) of dyn_atom as passed to campana_* # returns a 4-tuple of eigenvalues, the right eigenvectors (as 3x3 matrix), the radius of gyration, and the instantaneous asphericity delta* def gyration_tensor_by_bool(which,coords): if (np.sum(which) == 0): return None hlp1 = coords[which,0:3] nps = hlp1.shape[0] tmpcom = np.sum(hlp1,axis=0)/float(nps) tmpgyrten = np.zeros([3,3]) for i in range(0,nps): tmpgyrten += np.outer(hlp1[i,:]-tmpcom,hlp1[i,:]-tmpcom) tmpgyrten[:,:] /= (nps) evals, evecs = np.linalg.eig(tmpgyrten) rg = (evals[0] + evals[1] + evals[2]) deltastar = 1.0 - 3.0*(evals[0]*evals[1] + evals[0]*evals[2] + evals[1]*evals[2])/(rg**2) rg = np.sqrt(rg) return evals,evecs,rg,deltastar ############################################################################################## # # calculation of the principal axes of inertia a group of atoms (inertia_tensor_by_bool) ############################################################################################## # Arguments: # which should be the full-size True/False vector of which atoms to include # coords should be at least the first 3 columns (xyz) of dyn_atom as passed to campana_* # mass should be the mass column of static_atom_dbls # returns a 2-tuple of eigenvalues (principal moments of inertia) and the right eigenvectors (as 3x3 matrix) (body-coordinates) def inertia_tensor_by_bool(which,coords,mass): if (np.sum(which) == 0): return None hlp1 = coords[which,0:3] hlp2 = mass[which] nps = hlp1.shape[0] tmpcom = np.zeros(3) for i in range(0,nps): tmpcom += hlp2[i]*hlp1[i,:] tmpcom /= sum(hlp2) # note that this is extensive with mass and not a measure of size tmpinerten = np.zeros([3,3]) hlpx = hlp1[:,0]-tmpcom[0] hlpy = hlp1[:,1]-tmpcom[1] hlpz = hlp1[:,2]-tmpcom[2] tmpinerten[0,0] = np.sum(hlp2*(hlpy*hlpy+hlpz*hlpz)) tmpinerten[1,1] = np.sum(hlp2*(hlpx*hlpx+hlpz*hlpz)) tmpinerten[2,2] = np.sum(hlp2*(hlpx*hlpx+hlpy*hlpy)) tmpinerten[0,1] = -1.0*np.sum(hlp2*hlpx*hlpy) tmpinerten[0,2] = -1.0*np.sum(hlp2*hlpx*hlpz) tmpinerten[1,2] = -1.0*np.sum(hlp2*hlpy*hlpz) tmpinerten[1,0] = tmpinerten[0,1] tmpinerten[2,0] = tmpinerten[0,2] tmpinerten[2,1] = tmpinerten[1,2] evals, evecs = np.linalg.eig(tmpinerten) return evals,evecs ############################################################################################## # # calculation of the (fudged) dipole moment of a group of atoms (dipole_moment_by_bool) ############################################################################################## # Arguments: # which should be the full-size True/False vector of which atoms to include # coords should be at least the first 3 columns (xyz) of dyn_atom as passed to campana_* # pcharge should be the partial charge column of static_atom_dbls # returns a 3-tuple of the dipole moment, its norm, and the net charge def dipole_moment_by_bool(which,coords,pcharge): if (np.sum(which) == 0): return None todebye = 4.8033324 hlp1 = coords[which,0:3] hlp2 = pcharge[which] tc = sum(hlp2) # if there is a monopole present, our reported dipole won't be invariant -> fix nps = hlp1.shape[0] hlp2a = hlp2 - tc/float(nps) tmpcom = np.zeros(3) for i in range(0,nps): tmpcom += hlp2a[i]*hlp1[i,:] # if there is a monopole, remove it return todebye*tmpcom,todebye*np.sqrt(np.sum(tmpcom*tmpcom)),tc ############################################################################################## # # (weighted) alignment based on two sets of 3-space vectors (align_by_set) ############################################################################################## # Arguments: # oldcoords should be the Nx3 vector of xyz reference coordinates, remains unchanged # newcoords should be the exact matching set of new coordinates (to be aligned), remains unchanged # center should be a Boolean indicating whether the coordinates should be centered # optional: # "weights" can be an optional argument to weigh the alignment (must be vector of length N) # "coords" can be an optional argument to have the alignment operation be applied to (Mx3) # "subset" can be an optional argument to subset coords for this apply (Boolean, length M) # returns a 4 or 5-tuple of the original SciPy rotation object, the pre-shift (3-vector), the post-shift (3-vector), # and the quaternion (4-vector), and, possibly, the aligned version of "coords[subset]" (<=Mx3) or "coords" (Mx3) # # Warning 1: unfortunately, the results will depend slightly on the SciPy version since "match_vectors" # uses SVD while "align_vectors" uses the more appropriate Kabsch algorithm # Warning 2: this function assumes that periodic images have been taken care of beforehand if necessary; # no additional corrections are applied def align_by_set(oldcoords,newcoords,center=True,**kwargs): # match_vectors was replaced by align_vectors in 1.6 but align_vectors does not exist before ): # coordinates must be zero-meaned for rotation to make sense but "center" allows control versionAnnoyance = spVersion.version.split('.') noMatch = False if ((int(versionAnnoyance[0]) > 1) or (int(versionAnnoyance[1]) >= 6)): noMatch = True if ((len(oldcoords.shape) != 2) or (len(newcoords.shape) != 2)): return None nparts = oldcoords.shape[0] if (nparts != newcoords.shape[0]): return None if ((oldcoords.shape[1] != 3) or (newcoords.shape[1] != 3)): return None if ("weights" in kwargs): if (len(kwargs["weights"].shape) != 1): return None if (kwargs["weights"].shape[0] != nparts): return None theweights = kwargs["weights"] netweight = np.sum(theweights) minweight = np.min(theweights) if (netweight <= 0.0): return None if (minweight < 0.0): return None avgo = np.array([np.sum(theweights[:]*oldcoords[:,0]),np.sum(theweights[:]*oldcoords[:,1]),np.sum(theweights[:]*oldcoords[:,2])])/netweight avgn = -1.0*np.array([np.sum(theweights[:]*newcoords[:,0]),np.sum(theweights[:]*newcoords[:,1]),np.sum(theweights[:]*newcoords[:,2])])/netweight else: avgo = np.array([np.sum(oldcoords[:,0]),np.sum(oldcoords[:,1]),np.sum(oldcoords[:,2])])/float(nparts) avgn = -1.0*np.array([np.sum(newcoords[:,0]),np.sum(newcoords[:,1]),np.sum(newcoords[:,2])])/float(nparts) theweights = None # this is the default value for the SciPy fxn if (center == True): refcoords = np.array(oldcoords) testcoords = np.array(newcoords) for k in range(0,3): testcoords[:,k] += avgn[k] refcoords[:,k] -= avgo[k] if (noMatch == False): rot, smat = spRot.match_vectors(refcoords,testcoords,weights=theweights) else: rot, rmsd = spRot.align_vectors(refcoords,testcoords,weights=theweights,return_sensitivity=False) else: if (noMatch == False): rot, smat = spRot.match_vectors(oldcoords,newcoords,weights=theweights) else: rot, rmsd = spRot.align_vectors(oldcoords,newcoords,weights=theweights,return_sensitivity=False) thequat = rot.as_quat() if ("coords" in kwargs): if (len(kwargs["coords"].shape) != 2): return None if (kwargs["coords"].shape[1] < 3): return None nalls = kwargs["coords"].shape[0] if ("subset" in kwargs): if (len(kwargs["subset"].shape) != 1): return None if (kwargs["subset"].shape[0] != nalls): return None aligned = np.array(kwargs["coords"][kwargs["subset"],0:3]) else: aligned = np.array(kwargs["coords"][:,0:3]) for k in range(0,3): aligned[:,k] += avgn[k] aligned = rot.apply(aligned[:,0:3]) for k in range(0,3): aligned[:,k] += avgo[k] return rot,avgn,avgo,thequat,aligned else: return rot,avgn,avgo,thequat ############################################################################################## # # (weighted) RMSD based on two sets of 3-space vectors (rmsd_by_set) ############################################################################################## # Arguments: # refcoords should be the Nx3 vector of xyz reference coordinates, unchanged # newcoords should be the exact matching set of new coordinates (to be aligned), unchanged # optional: # "weights" can be an optional argument to weigh the RMSD (must be vector of length N) # "subset" can be an optional argument to subset coords for the RMSD sum (Boolean) # returns the RMSD # # Warning: this function assumes that periodic images have been taken care of beforehand if necessary; # no additional corrections are applied def rmsd_by_set(refcoords,newcoords,**kwargs): if ((len(refcoords.shape) != 2) or (len(newcoords.shape) != 2)): return None nits = refcoords.shape[0] if (nits != newcoords.shape[0]): return None if ((refcoords.shape[1] != 3) or (newcoords.shape[1] != 3)): return None if ("weights" in kwargs): if (len(kwargs["weights"].shape) != 1): return None if (kwargs["weights"].shape[0] != nits): return None theweights = kwargs["weights"] if ("subset" in kwargs): if (len(kwargs["subset"].shape) != 1): return None if (kwargs["subset"].shape[0] != nits): return None diffv = (refcoords[kwargs["subset"],:]-newcoords[kwargs["subset"],:]) if ("weights" in kwargs): netweight = np.sum(theweights[kwargs["subset"]]) minweight = np.min(theweights[kwargs["subset"]]) if (netweight <= 0.0): return None if (minweight < 0.0): return None normer = 1.0/netweight rmsd = np.sqrt(np.sum(theweights[kwargs["subset"]]*(diffv[:,0]*diffv[:,0]+diffv[:,1]*diffv[:,1]+diffv[:,2]*diffv[:,2]))*normer) else: normer = 1.0/float(diffv.shape[0]) rmsd = np.sqrt(np.sum(diffv[:,0]*diffv[:,0]+diffv[:,1]*diffv[:,1]+diffv[:,2]*diffv[:,2])*normer) else: diffv = (refcoords-newcoords) if ("weights" in kwargs): netweight = np.sum(theweights) minweight = np.min(theweights) if (netweight <= 0.0): return None if (minweight < 0.0): return None normer = 1.0/netweight rmsd = np.sqrt(np.sum(theweights[:]*(diffv[:,0]*diffv[:,0]+diffv[:,1]*diffv[:,1]+diffv[:,2]*diffv[:,2]))*normer) else: normer = 1.0/float(diffv.shape[0]) rmsd = np.sqrt(np.sum(diffv[:,0]*diffv[:,0]+diffv[:,1]*diffv[:,1]+diffv[:,2]*diffv[:,2])*normer) return rmsd ############################################################################################## # end of list of predefined utility functions ############################################################################################## ############################################################################################## # # ALL EDITS SHOULD BE BELOW THIS POINT IN GENERAL # ############################################################################################## # Finally, this is the class that should hold the permanent data structures a custom analysis uses, i.e., those written by you. # In addition, it contains a call counter and a static map of residue type to residue name. class counters: def __init__(self): # a copy of the data array in mod_aminos with manual offset self.restype_to_name = np.array(['???', 'GLY','ALA','VAL','LEU','ILE','SER','THR','CYS',\ 'PRO','PHE','TYR','TRP','HIP','HID','HIE','ASP',\ 'ASN','GLU','GLN','MET','LYS','ARG','ORN','AIB',\ 'PCA','UNK','ACE','FOR','NME','NH2','GAM','HYP',\ 'ABA','NVA','NLE','NA+','CL-','URE','SPC','T3P',\ 'NMF','NMA','ACA','PPA','T4P','CH4','FOA','DMA',\ 'MOH','PCR','DAB','K+ ','BR-','CS+','I- ','NH4',\ 'AC-','GDN','CYT','URA','THY','ADE','GUA','R5P',\ 'D5P','RPU','RPC','RPT','RPA','RPG','DPU','DPC',\ 'DPT','DPA','DPG','RIB','DIB','RIU','RIC','RIT',\ 'RIA','RIG','DIU','DIC','DIT','DIA','DIG','PRP',\ 'NBU','IBU','TOL','EOH','MSH','EMT','IMD','IME',\ 'MIN','1MN','2MN','LCP','NO3','T5P','T4E','BEN',\ 'NAP','PUR','O2 ','GLH','ASH','CYX','TYO','LYD',\ 'PTR','SEP','TPO','KAC','KM1','KM2','KM3']) # call counter self.ncalls = 0 # insert own code from here ###################################################################################################################################### # Both functions ("campana_step" and "campana_final") receive the same arguments. # Static properties: ################################################################################################################# # static_atom_ints, etc.: see documentation and Tutorial 18 of CAMPARI to see what these static structures contain # Dynamic properties that are wrapped: ############################################################################################### # box: see definition of box_holder class above # analyses: see definition of analysis_holder class above # user_counters: defined largely by you, see directly above # Dynamic properties that are raw nd-arrays: ######################################################################################### # dyn_atom: columns 0-2 are xyz(A), columns 3-5 are Z-matrix coordinates (length(A), angle(deg), dihedral or 2nd angle(deg)) # dyn_res: columns 0-2 are omega, phi, psi, columns 3-8 are nucleic acid backbone angles, columns 9-X are chi angles; # see documentation on FMCSC_SEQFILE for details on what these mean; # note that psi is the actual N-CA-C-N not, as often in CAMPARI, 180.0 + N-CA-C-O mod 2PI # dyn_mol: columns 0-2 are xyz(A) of center of mass, columns 3-5 are xyz(A) of center of geometry; # if the data come from an actual run in a (semi-)grand ensemble, there is also column 6 as either 1.0 or 0.0 to # distinguish actual from reservoir particles ###################################################################################################################################### # in all cases, if something is not defined (e.g., a residue has no phi-angle or an atom has no Z-matrix bond angle), the corresponding # values will (most likely) be 0.0 or 180.0 ###################################################################################################################################### def campana_step(static_atom_ints,static_atom_dbls,static_res_ints,static_res_dbls,static_mol_ints,static_mol_dbls,dyn_atom,dyn_res,dyn_mol,box,analyses,user_counters): user_counters.ncalls = user_counters.ncalls + 1 print("campana.py: This was invocation #",user_counters.ncalls, " of campana_step.") def campana_final(static_atom_ints,static_atom_dbls,static_res_ints,static_res_dbls,static_mol_ints,static_mol_dbls,dyn_atom,dyn_res,dyn_mol,box,analyses,user_counters): print("campana.py: This is the end of campana_final.")