ijddZddlZddlZddlZddlZddlZddlZdZdZ dZ dZ dZ d d Z dS) aStreamplots (3D) --- :mod:`MDAnalysis.visualization.streamlines_3D` =================================================================== :Authors: Tyler Reddy and Matthieu Chavent :Year: 2014 :Copyright: Lesser GNU Public License v2.1+ The :func:`generate_streamlines_3d` function can generate a 3D flow field from a MD trajectory, for instance, lipid molecules in a virus capsid. It can make use of multiple cores to perform the analyis in parallel (using :mod:`multiprocessing`). .. rubric: References .. footbibliography:: See Also -------- MDAnalysis.visualization.streamlines : streamplots in 2D .. autofunction:: generate_streamlines_3d Nctj||}|d}|j}|d|d|d|d|d|dg\}}}} } } ||z ||z||z | |z| |z | |zf} | S)aCalculate the extent of the atom coordinates + buffer. A function for the parent process which should take the input trajectory and calculate the limits of the container for the system and return these limits. Parameters ---------- topology_file_path : str topology file name trajectory_file_path : str trajectory file name buffer_value : float buffer value (padding) in +/- {x, y, z} all.r..) MDAnalysisUniverse select_atoms positionsminmax) topology_file_pathtrajectory_file_path buffer_valueuniverse_objectall_atom_selectionall_atom_coordinate_arrayx_minx_maxy_miny_maxz_minz_maxtuple_of_limitss q/srv/www/vhosts/g4struct/public_html/venv/lib/python3.11/site-packages/MDAnalysis/visualization/streamlines_3D.pydetermine_container_limitsr;s$!)0O)55 !3 <!&)--//!&)--//!&)--//!&)--//!&)--//!&)--// 0,E5%u        O cT|\}}}}}}tj|||||||||f}|S)uProduce a 3D grid for the simulation system. The partitioning is based on the tuple of Cartesian Coordinate limits calculated in an earlier step. Parameters ---------- tuple_of_limits : tuple ``x_min, x_max, y_min, y_max, z_min, z_max`` grid_spacing : float grid size in all directions in ångström Returns ------- grid : array ``numpy.mgrid[x_min:x_max:grid_spacing, y_min:y_max:grid_spacing, z_min:z_max:grid_spacing]`` )npmgrid) r grid_spacingrrrrrrgrids r produce_gridr%gsL&0?,E5%u 8 eL  eL  eL  " D Krc |\}}}|jd}|jd}td|jD}g}|D]+} | } || ,g} |D]+} | } | | ,g}|D]+}|}||,g}t || |D]7\}}}|t t |||8t j|}d}i}d}||dz kr||}||dz}d}||jd|z krQd}||dz kr0g}|||dzdf}|||dzdf}|||zd|z|zdf} |||zd|z|zdf}!||| |!fD]}"| |"t j|}#|#jdks Jdt j t j|d }$|$|d ||<|dz }|dz }|dz }||dz kr|dz }||dz k0||jd|z kQ|dz }||dz kt|}%t j d|%}&t j |&|}'g}(d})|'D]n}*i}+t|*},d}-|-|,kr9|\}.}/|+|.|/i|-dz }-|-|,k9|(|+|)dz })o|(|%||fS) aSplit the grid into blocks of vertices. Take the overall `grid` for the system and split it into lists of cube vertices that can be distributed to each core. Parameters ---------- grid : numpy.array 3D array num_cores : int number of partitions to generate Returns ------- list_dictionaries_for_cores : list of dict total_cubes : int num_sheets : int delta_array_shape : tuple rcg|]}|dz S)r).0ns r zsplit_grid..s   1Q   rrr .)z vertex_array has incorrect shapeaxis)centroid vertex_list)shapetupleflattenappendziplistr!arraytolistextendaveragelenarange array_splitpopitemupdate)0r$ num_coresxyz num_z_values num_sheetsdelta_array_shapeordered_list_per_sheet_x_valuesx_sheet array_all_x_values_current_sheetordered_list_per_sheet_y_values y_columns array_all_y_values_current_sheetordered_list_per_sheet_z_valuesz_slices array_all_z_values_current_sheet,ordered_list_cartesian_coordinates_per_sheetx_sheet_coordsy_sheet_coordsz_sheet_coords(array_ordered_cartesian_coords_per_sheetcurrent_base_sheet"dictionary_cubes_centroids_indices cube_countercurrent_base_sheet_arraycurrent_top_sheet_array current_indexcolumn_z_levelcurrent_list_cube_verticesfirst_two_vertices_base_sheetfirst_two_vertices_top_sheetnext_two_vertices_base_sheetnext_two_vertices_top_sheet vertex_set vertex_array cube_centroid total_cubespseudo_cube_indices sublist_of_cube_indices_per_corelist_dictionaries_for_coressubdictionary_counterarray_cube_indicescurrent_core_dictionary items_to_pop items_poppedkeyvalues0 r split_gridrqsR,GAq!72;LJ     ')#   ,3??+<+<('.. ,    ')#  +4+<+<+>+>('.. ,    ')#  +3+;+;+=+=('.. ,    460:=''';;  6 5;; ^^^DD E E    02x400, )+&L zA~ - -#K $  #K  "#  6 # #"$-..  \));CCEEJC # * *C< 8 8 8 A L\)) $**+BCCC"#  rcV g}g}d fd}d} |||| |||S)znRun the analysis on one core. The code to perform on a given core given the dictionary of cube data. c*tj|}tjj|ddz }tjj||tjddf}tj || dddtj tj ||}|d}tj||k}|d } tj| |k} |d } tj| |k} tj|d| d} tj| | d}|S) z7Determine if an array of coordinates are within a cube. euclidean@Nrga2U0*#?zJnot all cube vertex to centroid distances are the same, so not a true cube)rtolatolerr_msgrrr)r!r9scipyspatialdistancepdistrcdistnewaxistestingassert_allcloserabsolutesubtractwhere intersect1d)array_point_coordinateslist_cube_verticesrearray_cube_verticescube_half_side_length)array_cube_vertex_distances_from_centroidabsolute_delta_coordsabsolute_delta_x_coordsindices_delta_x_acceptableabsolute_delta_y_coordsindices_delta_y_acceptableabsolute_delta_z_coordsindices_delta_z_acceptable!intersection_xy_acceptable_arrays&overall_indices_points_in_current_cubes r point_in_cubez$per_core_work..point_in_cube8s !h'9:: M " ( (#[  cee   M " ( (#]2:qqq=%A   2 "" 5 9 9 ; ; 5 9 9 ; ;! #   !# K/ ? ?! ! #8"?%'X #'< <& & "#8"?%'X #'< <& & "#8"?%'X #'< <& & "-/N &q )+Ea+H- - )24 -/I!/L2 2 .65rcd}|D]b\}}||d|d}||d<t|dkr"tj||d}||d<nd|d<|dz }cdS) a=Basically update the cube dictionary objects assigned to this core to contain a new key/value pair corresponding to the indices of the relevant particles that fall within a given cube. Also, for a given cube, store a key/value pair for the centroid of the particles that fall within the cube. rr2r1start_frame_index_list_in_cuber/!centroid_of_particles_first_frameNr)itemsr=r!r<)%array_simulation_particle_coordinatesdictionary_cube_data_this_corerYrocubeindex_list_in_cubecentroid_particles_in_cubers r+update_dictionary_point_in_cube_start_framezBper_core_work..update_dictionary_point_in_cube_start_frameis 7==??  IC!.5]#Z ""  6HD1 2&''!++-/Z9:LM...* /899=A89 A LL'  rc"d}|D]w\}}|dd|d<d|d<d|d<nU||d}tj|d}||d <||dz }|d|d<|d |d<|d |d<|d z }xdS) zWUpdate the cube dictionary objects again as appropriate for the second and final frame.rrNdxdydzrr/ centroid_of_paticles_final_framerr )rr!r<)rrrYror8new_coordinate_array_for_particles_starting_in_this_cube0new_centroid_for_particles_starting_in_this_cubedelta_centroid_array_this_cubes rupdate_dictionary_end_framez2per_core_work..update_dictionary_end_frames 7==??  IC78@T T T :=>I DF:LDDD@ E78E>?@/T ;A>T ;A>T A LL5  rr)) start_frame_coord_arrayend_frame_coord_arrayr MDA_selection start_frame end_framelist_previous_frame_centroidslist_previous_frame_indicesrrrs @r per_core_workr(s%'!"$/6/6/6b:J0/!? = *)rctj||}||}|jD]5}|j|krn'|j|kr|j}"|j|kr|j} 56|| fS)aGenerate coordinate arrays. To reduce memory footprint produce only a single MDA selection and get desired coordinate arrays; can later send these coordinate arrays to all child processes rather than having each child process open a trajectory and waste memory. )r r r trajectoryframer ) rrrrrrrelevant_particlests2start_frame_relevant_particle_coordinate_array_xyz0end_frame_relevant_particle_coordinate_array_xyzs r(produce_coordinate_arrays_single_processrs!)0O)55mDD(   8i   E 8{ " "", ? >X " "", = < :8 rrumaximumc B"| dkrtj} n| } tjddi""fd}|||| | | f}t ||}t || \}}}}t |||||\}}tj| }|D]%}|t||||||f|&| | tt|t|d z z }tj|}tj|}tj|}d }d }d }d } td |D]r}!"|!d ||||f<"|!d ||||f<"|!d ||||f<|d z }| d z } ||dkrd }|d z }||d kr |d z }d }d }d } sd|t!|| k<d|t!|| k<d|t!|| k<|||fS)aProduce the x, y and z components of a 3D streamplot data set. Parameters ---------- topology_file_path : str Absolute path to the topology file trajectory_file_path : str Absolute path to the trajectory file. It will normally be desirable to filter the trajectory with a tool such as GROMACS :program:`g_filter` (see :footcite:p:`Chavent2014`) grid_spacing : float The spacing between grid lines (angstroms) MDA_selection : str MDAnalysis selection string start_frame : int First frame number to parse end_frame : int Last frame number to parse xmin : float Minimum coordinate boundary for x-axis (angstroms) xmax : float Maximum coordinate boundary for x-axis (angstroms) ymin : float Minimum coordinate boundary for y-axis (angstroms) ymax : float Maximum coordinate boundary for y-axis (angstroms) maximum_delta_magnitude : float Absolute value of the largest displacement tolerated for the centroid of a group of particles ( angstroms). Values above this displacement will not count in the streamplot (treated as excessively large displacements crossing the periodic boundary) num_cores : int or 'maximum' (optional) The number of cores to use. (Default 'maximum' uses all available cores) Returns ------- dx_array : array of floats An array object containing the displacements in the x direction dy_array : array of floats An array object containing the displacements in the y direction dz_array : array of floats An array object containing the displacements in the z direction Examples -------- Generate 3D streamlines and visualize in `mayavi`_:: import numpy as np import MDAnalysis import MDAnalysis.visualization.streamlines_3D import mayavi, mayavi.mlab # assign coordinate system limits and grid spacing: x_lower,x_upper = -8.73, 1225.96 y_lower,y_upper = -12.58, 1224.34 z_lower,z_upper = -300, 300 grid_spacing_value = 20 x1, y1, z1 = MDAnalysis.visualization.streamlines_3D.generate_streamlines_3d( 'testing.gro', 'testing_filtered.xtc', xmin=x_lower, xmax=x_upper, ymin=y_lower, ymax=y_upper, zmin=z_lower, zmax=z_upper, grid_spacing=grid_spacing_value, MDA_selection = 'name PO4', start_frame=2, end_frame=3, num_cores='maximum') x, y, z = np.mgrid[x_lower:x_upper:x1.shape[0]*1j, y_lower:y_upper:y1.shape[1]*1j, z_lower:z_upper:z1.shape[2]*1j] # plot with mayavi: fig = mayavi.mlab.figure(bgcolor=(1.0, 1.0, 1.0), size=(800, 800), fgcolor=(0, 0, 0)) for z_value in np.arange(z_lower, z_upper, grid_spacing_value): st = mayavi.mlab.flow(x, y, z, x1, y1, z1, line_width=1, seedtype='plane', integration_direction='both') st.streamline_type = 'tube' st.tube_filter.radius = 2 st.seed.widget.origin = np.array([ x_lower, y_upper, z_value]) st.seed.widget.point1 = np.array([ x_upper, y_upper, z_value]) st.seed.widget.point2 = np.array([ x_lower, y_lower, z_value]) st.seed.widget.resolution = int(x1.shape[0]) st.seed.widget.enabled = False mayavi.mlab.axes(extent = [0, 1200, 0, 1200, -300, 300]) fig.scene.z_plus_view() mayavi.mlab.savefig('test_streamplot_3D.png') # more compelling examples can be produced for vesicles and other spherical systems .. image:: test_streamplot_3D.png See Also -------- MDAnalysis.visualization.streamlines.generate_streamlines .. _mayavi: http://docs.enthought.com/mayavi/mayavi/ rwarnraise)roverc2|dS)N)rA) process_dictparent_cube_dictionarys rlog_result_to_parentz5generate_streamlines_3d..log_result_to_parentWs%%l33333r)rr#)r$rB)argscallbackrrrrrr g?)multiprocessing cpu_countr!seterrr%rqrPool apply_asyncrclosejoinintfloatzerosrangeabs)#rrr#rrrxminxmaxyminymaxzminzmaxmaximum_delta_magnituderBrrr$rirfrGrHrrpoolsub_dictionary_of_cube_datacubes_per_sheetdx_arraydy_arraydz_array current_sheety_index_current_sheetz_index_current_columntotal_cubes_current_sheet cube_numberrs# @rgenerate_streamlines_3drshI#-//  I&w''''44444 T4tT:O 'l   D  222Lj:K 1       32   * *D'B   #  '%+ *   JJLLLIIKKK % ,,uZ!^/D/DDEEOx)**Hx)**Hx)**HM !Q ,,..  #; / 5  02H H #; / 5  02H H #; / 5  02H H  !#!Q&! !%6q%9 9 9%& " !Q & !%):1)===" ()%)*&,-):=HS]]5 569rs02 )))X8bbbJI*I*I*X$$$h R*R*R*R*R*R*r