i2dZddlZddlZddlZddlZddlmZddl m Z m Z ddl m Z ddlmZmZejeZe d d ZGd d eZd ZdS)a Polymer analysis --- :mod:`MDAnalysis.analysis.polymer` ======================================================= :Author: Richard J. Gowers :Year: 2015, 2018 :Copyright: Lesser GNU Public License v2.1+ This module contains various commonly used tools in analysing polymers. N) NoDataError)requires AtomGroup) calc_bonds) AnalysisBase ResultsGroupbondscxfdD}dt|D}dt|D}t|dkr?tddd|Dt|dkr?td dd |Dt |dg}t td z D]}|d }|jz|z }||z }|S) aRearrange a linear AtomGroup into backbone order Requires that the backbone has bond information, and that only backbone atoms are provided (ie no side chains or hydrogens). Parameters ---------- backbone : AtomGroup the backbone atoms, not necessarily in order Returns ------- sorted_backbone : AtomGroup backbone in order, so `sorted_backbone[i]` is bonded to `sorted_backbone[i - 1]` and `sorted_backbone[i + 1]` .. versionadded:: 0.20.0 c>g|]}t|jzS)len bonded_atoms).0atombackbones e/srv/www/vhosts/g4struct/public_html/venv/lib/python3.11/site-packages/MDAnalysis/analysis/polymer.py z!sort_backbone..Fs)FFFTs4$x/00FFFc$g|] \}}|dk |S)rrrrds rrz!sort_backbone..Gs!GGG74Q$rc g|] \}}|dv | S))rrrrs rrz!sort_backbone..Hs(qavooooorrzBackbone contains atoms with connectivity degree not equal to 1 or 2. This suggests branches or isolated atoms. Problematic atoms: {}.z, c34K|]}t|VdSNstrras r z sort_backbone..Ps(<<1A<<<<<.Ys(;;1A;;;;;;rr)zipr ValueErrorformatjoinrranger)rdegrees deg1_atoms wrong_atomssorted_backbone_end_atom next_atoms` r sort_backboner00su,GFFFXFFFGGGc(G&<&<GGGJ'22K ;! vdii<< <<<<<==    :! vdii;; ;;;;;<<     A00O 3x==1$ % % % %"2&*X5H  9$ rceZdZdZdZedZfdZdZdZ dZ e dZ e d Z e d Zd Zd Zd ZddZxZS)PersistenceLengthaCalculate the persistence length for polymer chains The persistence length is the length at which two points on the polymer chain become decorrelated. This is determined by first measuring the autocorrelation (:math:`C(n)`) of two bond vectors (:math:`\mathbf{a}_i, \mathbf{a}_{i + n}`) separated by :math:`n` bonds .. math:: C(n) = \langle \cos\theta_{i, i+n} \rangle = \langle \mathbf{a_i} \cdot \mathbf{a_{i+n}} \rangle where :math:`a_i` and :math:`a_{i+n}` are unit vectors along the bonds. An exponential decay is then fitted to this, which yields the persistence length .. math:: C(n) \approx \exp\left( - \frac{n \bar{l_B}}{l_P} \right) where :math:`\bar{l_B}` is the average bond length, and :math:`l_P` is the persistence length which is fitted Parameters ---------- atomgroups : iterable List of AtomGroups. Each should represent a single polymer chain, ordered in the correct order. verbose : bool, optional Show detailed progress of the calculation if set to ``True``. Attributes ---------- results.bond_autocorrelation : numpy.ndarray the measured bond autocorrelation results.lb : float the average bond length .. versionadded:: 2.0.0 lb : float Alias to the :attr:`results.lb`. .. deprecated:: 2.0.0 Will be removed in MDAnalysis 3.0.0. Please use :attr:`results.lb` instead. results.x : numpy.ndarray length of the decorrelation predicted by *lp* .. versionadded:: 2.0.0 results.lp : float calculated persistence length .. versionadded:: 2.0.0 lp : float Alias to the :attr:`results.lp`. .. deprecated:: 2.0.0 Will be removed in MDAnalysis 3.0.0. Please use :attr:`results.lp` instead. results.fit : numpy.ndarray the modelled backbone decorrelation predicted by *lp* .. versionadded:: 2.0.0 fit : float Alias to the :attr:`results.fit`. .. deprecated:: 2.0.0 Will be removed in MDAnalysis 3.0.0. Please use :attr:`results.fit` instead. See Also -------- :func:`sort_backbone` for producing the sorted AtomGroup required for input. Example ------- .. code-block:: python from MDAnalysis.tests.datafiles import TRZ_psf, TRZ import MDAnalysis as mda from MDAnalysis.analysis import polymer u = mda.Universe(TRZ_psf, TRZ) # this system is a pure polymer melt of polyamide, # so we can select the chains by using the .fragments attribute chains = u.atoms.fragments # select only the backbone atoms for each chain backbones = [chain.select_atoms('not name O* H*') for chain in chains] # sort the chains, removing any non-backbone atoms sorted_backbones = [polymer.sort_backbone(bb) for bb in backbones] persistence_length = polymer.PersistenceLength(sorted_backbones) # Run the analysis, this will average over all polymer chains # and all timesteps in trajectory persistence_length = persistence_length.run() print(f'The persistence length is: {persistence_length.results.lp}') # always check the visualisation of this: persistence_length.plot() .. versionadded:: 0.13.0 .. versionchanged:: 0.20.0 The run method now automatically performs the exponential fit .. versionchanged:: 1.0.0 Deprecated :meth:`PersistenceLength.perform_fit` has now been removed. .. versionchanged:: 2.0.0 Former ``results`` are now stored as ``results.bond_autocorrelation``. :attr:`lb`, :attr:`lp`, :attr:`fit` are now stored in a :class:`MDAnalysis.analysis.base.Results` instance. .. versionchanged:: 2.10.0 introduced :meth:`get_supported_backends` allowing for parallel execution on ``multiprocessing`` and ``dask`` backends. TcdS)N)serialmultiprocessingdaskr)clss rget_supported_backendsz(PersistenceLength.get_supported_backendss44rc tt|j|djjfi|||_d|D}t |dtfd|Dstd|_ dS)Nrc,g|]}t|Sr)r)rags rrz.PersistenceLength.__init__..s---BB---rc3$K|] }|kV dSrr)rl chainlengths rr!z-PersistenceLength.__init__..s(221 #222222rz%Not all AtomGroups were the same size) superr2__init__universe trajectory _atomgroupsrallr%r>)self atomgroupskwargslensr> __class__s @rr@zPersistenceLength.__init__s/&&/ qM " -  17   &.-*---*Q-(( 2222T22222 FDEE E&rcjtj|jdz tj|j_dS)Nr)dtype)npzerosr>float32resultsraw_bond_autocorrrEs r_preparezPersistenceLength._prepares2)+  q  * * *  &&&rc |jD]}|jdd|jddz }|tj||zddddfz}tj||}t |jdz D]3}|jj d|jdz |z xx|||dfz cc<4dS)Nrr#)axis) rC positionsrLsqrtsuminnerr(r>rOrP)rEchainvecsinner_pris r _single_framezPersistenceLength._single_frames % % %E?122&")==D BGTD[--1-5566qqq$w? ?Dxd++H4+a/00 % % .0t'!+q00ae_% % % %rc:tdtjiS)NrP)lookup)r ndarray_sumrQs r_get_aggregatorz!PersistenceLength._get_aggregators&#\%=    rcRd}tj|t|jjS)NzThe `lb` attribute was deprecated in MDAnalysis 2.0.0 and will be removed in MDAnalysis 3.0.0. Please use `results.variance` instead.)warningswarnDeprecationWarningrOlbrEwmsgs rrfzPersistenceLength.lb + 5  d.///|rcRd}tj|t|jjS)NzThe `lp` attribute was deprecated in MDAnalysis 2.0.0 and will be removed in MDAnalysis 3.0.0. Please use `results.variance` instead.)rcrdrerOlprgs rrkzPersistenceLength.lp*rircRd}tj|t|jjS)NzThe `fit` attribute was deprecated in MDAnalysis 2.0.0 and will be removed in MDAnalysis 3.0.0. Please use `results.variance` instead.)rcrdrerOfitrgs rrmzPersistenceLength.fit4s, 5  d.///|rctj|jdz d|jdz }|t|j|jjzz}|jj|z |j_ | | dS)Nr) rLlinspacer>rrC _trajectoryn_framesrOrPbond_autocorrelation_calc_bond_length _perform_fit)rEnorms r _concludezPersistenceLength._conclude>s{4+a/D4Dq4HII D$%%(8(AAA L *T 1 )     rcg}|jD]P}|j}t|dd|dd}||Qt j||j_dS)zcalculate average bond lengthNr#r)rCrUrmeanappendrLrOrf)rEbsr;posbs rrsz#PersistenceLength._calc_bond_lengthGsq "  B,C3ss8SW--2244A IIaLLLL'"++ rc |jjn#t$rtddwxYw|jjt jt|jjz|j_t|jj|jj|j_ t j |jj |jj z |j_ dS)z'Fit the results to an exponential decayzUse the run method firstN) rOrrAttributeErrorrrfrLarangerxfit_exponential_decayrkexprmrQs rrtzPersistenceLength._perform_fitPs D L - - - D D D899t C D29  1 2 2, ,   0 LNDL=   64<>/DLO"CDD s *Ncddlm}||\}}||jj|jjdd||jj|jjd|d| d| d d |jj z| d |S) aVisualize the results and fit Parameters ---------- ax : matplotlib.Axes, optional if provided, the graph is plotted on this axis Returns ------- ax : the axis that the graph was plotted on rNroResult)labelFitrz$C(x)$g(best)loc) matplotlib.pyplotpyplotsubplotsplotrOrrrrm set_xlabel set_ylabelset_xlimrflegend)rEaxpltr-s rrzPersistenceLength.plot`s (''''' :LLNNEAr  LN L -       0>>> d i    Cdlo-... f  rr)__name__ __module__ __qualname____doc__%_analysis_algorithm_is_parallelizable classmethodr8r@rRr]rapropertyrfrkrmrvrsrtr __classcell__)rIs@rr2r2ms(B-1)55[5 ' ' ' ' '   %%%$   XX  X &&&EEE rr2cfd}tj|||dd}|S)aZFit a function to an exponential decay .. math:: y = \exp\left(- \frac{x}{a}\right) Parameters ---------- x, y : array_like The two arrays of data Returns ------- a : float The coefficient *a* for this decay Notes ----- This function assumes that data starts at 1.0 and decays to 0.0 c2tj| |z Sr)rLr)rr s rexpfuncz&fit_exponential_decay..expfuncsvqb1f~~rr)scipyoptimize curve_fit)ryrr s rrrs;*   !Q//215A Hr)rnumpyrLscipy.optimizerrcloggingr core.groupsrr lib.distancesrbaser r getLoggerrloggerr0r2rrrrrs0  --------&&&&&&,,,,,,,,  8 $ $ '999xPPPPP PPPf     r