Gromos Trajectory Examples¶
Example file for the evaluation of GROMOS trajectory files in pygromos¶
[3]:
#specific imports from pygromos for trc and tre file support
import pygromos.files.trajectory.trc as traj_trc
import pygromos.files.trajectory.tre as traj_tre
1) TRC¶
1.1) TRC import¶
[4]:
# import the trajectory file into a Trc class
trc = traj_trc.Trc(traj_path='../example_files/Traj_files/b_emin_vacuum_1.trc',
in_cnf='../example_files/Traj_files/b_emin_vacuum.cnf')
The TRC class bridges the GROMOS TRC file based strucutres and mdtraj.Trajectories.
One can read in TRCs and use them just like mdtraj.Trajectories. One can than write them out both as .h5 as well as GROMOS .trc/trc.gz files.
If you have a function that’s generally useful, please contact the developers to possibly add it to the pygromos code to help other people :)
[5]:
[x for x in dir(trc) if not x.startswith("_")]
[5]:
['TITLE',
'atom_slice',
'center_coordinates',
'distances',
'generate_TITLE_entry',
'generate_entry_for_frame',
'get_dummy_cnf',
'image_molecules',
'join',
'load',
'make_molecules_whole',
'n_atoms',
'n_chains',
'n_frames',
'n_residues',
'openmm_boxes',
'openmm_positions',
'parse_trc_efficiently',
'path',
'recreate_view',
'remove_solvent',
'restrict_atoms',
'rmsd',
'save',
'save_amberrst7',
'save_binpos',
'save_dcd',
'save_dtr',
'save_gro',
'save_gsd',
'save_hdf5',
'save_lammpstrj',
'save_lh5',
'save_mdcrd',
'save_netcdf',
'save_netcdfrst',
'save_pdb',
'save_tng',
'save_trr',
'save_xtc',
'save_xyz',
'slice',
'smooth',
'stack',
'step',
'superpose',
'time',
'timestep',
'to_cnf',
'top',
'topology',
'unitcell_angles',
'unitcell_lengths',
'unitcell_vectors',
'unitcell_volumes',
'view',
'write',
'write_trc',
'xyz']
1.2) TRC file handling¶
Save as .h5¶
[6]:
trc.save('./traj.h5')
[6]:
'./traj.h5'
Save as trc¶
[7]:
trc.save('./traj.trc.gz')
export last frame as cnf¶
[8]:
trc[-1].to_cnf(base_cnf='../example_files/Traj_files/b_emin_vacuum.cnf')
[8]:
TITLE
THIS IS THE FRAME AT TIMESTEP: 4.0 OF THE TRAJECTORY WITH TITLE: MD++steepest descent energy minimization of the peptide in vacuum>>> Generated with PyGromosTools (riniker group) <<<position trajectoryEND
POSITION
#
1 VAL H1 0 0.858298182 1.499708533 1.397253275
1 VAL H2 1 0.953678370 1.440333724 1.515184283
1 VAL N 2 0.913451612 1.421260476 1.425633550
1 VAL H3 3 0.855719924 1.339777827 1.431015491
1 VAL CA 4 1.022031069 1.399157882 1.329032302
1 VAL CB 5 0.952420831 1.374197721 1.195075631
1 VAL CG1 6 0.913190961 1.499370813 1.116295457
1 VAL CG2 7 1.026286483 1.267775059 1.113663197
1 VAL C 8 1.119535685 1.516167402 1.343593240
1 VAL O 9 1.121832848 1.571724653 1.453307033
2 TYR N 10 1.202524424 1.547034740 1.244336247
2 TYR H 11 1.207148552 1.490973473 1.161649466
2 TYR CA 12 1.303517222 1.653589725 1.236642957
2 TYR CB 13 1.247369647 1.791285038 1.200619698
2 TYR CG 14 1.197796106 1.877276540 1.317051053
2 TYR CD1 15 1.071256876 1.857888699 1.371185541
2 TYR HD1 16 0.993085802 1.811082721 1.311375499
2 TYR CD2 17 1.288784027 1.961568475 1.379797101
2 TYR HD2 18 1.383159995 1.985753894 1.330927134
2 TYR CE1 19 1.040328622 1.914669991 1.494227052
2 TYR HE1 20 0.939020157 1.905765176 1.533425450
2 TYR CE2 21 1.257037044 2.019920826 1.501881838
2 TYR HE2 22 1.325480461 2.093687534 1.543774128
2 TYR CZ 23 1.134425282 1.992075920 1.561123490
2 TYR OH 24 1.123193622 2.005527496 1.695985913
2 TYR HH 25 1.175849319 2.085860014 1.723783255
2 TYR C 26 1.418586850 1.648926616 1.337371469
2 TYR O 27 1.523679972 1.595995188 1.301521301
3 ARG N 28 1.396972656 1.675651789 1.465864182
3 ARG H 29 1.309179425 1.712227583 1.496792197
3 ARG CA 30 1.504508853 1.665650368 1.565602183
3 ARG CB 31 1.486977339 1.767935872 1.678033829
3 ARG CG 32 1.373603702 1.738026977 1.776329041
3 ARG CD 33 1.415657640 1.794736385 1.912081718
3 ARG NE 34 1.386569738 1.698562503 2.019374609
3 ARG HE 35 1.292315245 1.701340795 2.052666664
3 ARG CZ 36 1.472650051 1.610763788 2.072641850
3 ARG NH1 37 1.596361995 1.592568755 2.024506330
3 ARG HH11 38 1.661001921 1.530819535 2.069304466
3 ARG HH12 39 1.621459961 1.630685449 1.935536146
3 ARG NH2 40 1.436868072 1.541736007 2.181771994
3 ARG HH21 41 1.500163198 1.475605249 2.222011805
3 ARG HH22 42 1.348556042 1.557563543 2.225923061
3 ARG C 43 1.542207122 1.526291847 1.616297126
3 ARG O 44 1.574849486 1.506344676 1.733203053
4 LYSH N 45 1.538578033 1.431676269 1.522896171
4 LYSH H 46 1.539155364 1.463887453 1.428223729
4 LYSH CA 47 1.564581990 1.287459135 1.534663439
4 LYSH CB 48 1.482068777 1.220921040 1.644995928
4 LYSH CG 49 1.553891420 1.094621658 1.692990303
4 LYSH CD 50 1.455659389 0.988943756 1.743898392
4 LYSH CE 51 1.532733679 0.865419090 1.790902019
4 LYSH NZ 52 1.442433715 0.751384020 1.812048912
4 LYSH HZ1 53 1.408901930 0.719389319 1.723440170
4 LYSH HZ2 54 1.364803195 0.780083120 1.868170381
4 LYSH HZ3 55 1.491788745 0.677637756 1.858152866
4 LYSH C 56 1.540888786 1.215819478 1.401547551
4 LYSH O 57 1.445617437 1.139712334 1.385369897
5 GLN N 58 1.620790362 1.257511973 1.303739190
5 GLN H 59 1.691180825 1.327280164 1.317106128
5 GLN CA 60 1.639090419 1.212603927 1.164958239
5 GLN CB 61 1.516473532 1.236391783 1.076570988
5 GLN CG 62 1.487887502 1.382966638 1.043351531
5 GLN CD 63 1.366745234 1.400402308 0.951538980
5 GLN OE1 64 1.254561663 1.426969886 0.994376242
5 GLN NE2 65 1.389542460 1.397161841 0.820551813
5 GLN HE21 66 1.312573314 1.408466220 0.757729173
5 GLN HE22 67 1.483796239 1.384176731 0.789786577
5 GLN C 68 1.759001017 1.293522954 1.115122199
5 GLN O1 69 1.792673707 1.280218840 0.995480478
5 GLN O2 70 1.810874224 1.373489738 1.196003079
END
GENBOX
0
0.000000000 0.000000000 0.000000000
0.000000000 0.000000000 0.000000000
0.000000000 0.000000000 0.000000000
0.000000000 0.000000000 0.000000000
END
Get selection of trajectories¶
[9]:
trc[::10].n_frames
[9]:
9
Save out selection¶
[10]:
trc[::10].save('9_frames.trc')
1.3) TRC Calculate¶
[11]:
# Calculate the rmsd to the initial frame (0th frame).
# Alternatively a different trajectory can be provide as argument to the rmsd function.
# The accepted arguments are integer or single trajectory frame.
rmsd = trc.rmsd(0)
[12]:
# Which returns the rmsd for every time frame to the initial frame.
# It can be seen how the rmsd slowly gets larger as the simulations get farther away from the initial setup.
rmsd
[12]:
| rmsd | |
|---|---|
| time | |
| 0.00 | 0.000000 |
| 0.05 | 0.013725 |
| 0.10 | 0.022836 |
| 0.15 | 0.030384 |
| 0.20 | 0.036549 |
| ... | ... |
| 3.80 | 0.177028 |
| 3.85 | 0.177924 |
| 3.90 | 0.178814 |
| 3.95 | 0.179699 |
| 4.00 | 0.180578 |
81 rows × 1 columns
[13]:
# The mean over all frames can be easily taken with the pandas function mean()
rmsd.mean()
[13]:
rmsd 0.121639
dtype: float32
1.4) TRC Visualization¶
[14]:
view = trc.view
view
Do some changes and see how the visualization changes¶
To jump to a certain frame, use the slider in the widget or select frames with the following cell:
[15]:
trc.view.frame = 70
let’s find hydrogenbonds¶
[16]:
trc.view.add_contact(hydrogenBond=True)
add surface¶
[17]:
trc.view.add_surface(color="lightgrey",opacity=0.3)
remove surface¶
[18]:
trc.view.remove_surface()
add transperent cartoon¶
[19]:
trc.view.add_representation("cartoon",selection="protein",color='lightblue',opacity=0.3)
show distance between atom pair¶
[20]:
trc.view.add_distance(atom_pair=[[10,16]], label_color="black",color="darkgrey")
once you are happy with the result render and download image¶
[21]:
trc.view.download_image('first_protein.png')
Go crazy and make a movie to show off¶
[ ]:
### uncomment this for fast installation of NOT provided package
# import sys
# !{sys.executable} -m pip install moviepy
import moviepy.editor as mpy
from IPython import display
from time import sleep
import tqdm
import os
# Setup frames
frames = len(trc.step)
template = '0image{}.png'
download_dir = os.getcwd()
imagefiles = [download_dir + template.format(str(i)) for i in range(0, frames, 1)]
# render frames
for frame in tqdm.tqdm(range(0, frames)):
trc.view.frame = frame
sleep(0.2) # depending on the speed of the computer this might need to be increased
trc.view.download_image(filename=template.format(frame))
sleep(1) # depending on the speed of the computer this might need to be increased
# set frames per second and generate object
frame_per_second = 8
im = mpy.ImageSequenceClip(imagefiles, fps=frame_per_second)
im.write_gif('protein.gif', fps=frame_per_second)
# show movie
display.HTML("<img src='protein.gif'></img>")
2) TRE¶
2.1) Tre import and structure¶
[22]:
# import the trajectory file into a Tre class
from pygromos.files.trajectory.tre_field_libs.ene_fields import gromos_2015_tre_block_names_table
tre = traj_tre.Tre(input_value="../example_files/Traj_files/test_CHE_H2O_bilayer.tre", _ene_ana_names=gromos_2015_tre_block_names_table)
/home/bschroed/Documents/projects/PyGromosTools/pygromos/files/trajectory/_general_trajectory.py:296: PerformanceWarning:
your performance may suffer as PyTables will pickle object types that it cannot
map directly to c-types [inferred_type->mixed,key->block2_values] [items->Index(['totals', 'baths', 'bonded', 'nonbonded', 'special', 'eds', 'mass',
'temperature', 'volume', 'pressure'],
dtype='object')]
self.database.to_hdf(
[23]:
tre.database
[23]:
| step | time | totals | baths | bonded | nonbonded | special | eds | mass | temperature | volume | pressure | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 0.0 | [-408806.2436, 393521.852, -802328.0956, 13840... | [[99739.83003, 29519.49412, 70220.33591], [293... | [[0.0, 49931.35526, 0.0, 88472.52195, 0.0]] | [[984602.6158, -1925334.589, 0.0, 0.0]] | [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,... | [0.0] | [1277287.736] | [[296.2008302, 350.6569429, 278.0450574, 1.000... | [2425.07286, 8.463596169, 0.0, 0.0, 0.0, 8.463... | [29.98293036, 94773.2037, 131128.599, 30.37611... |
| 1 | 1000 | 2.0 | [-411576.4224, 392980.7347, -804557.157, 13955... | [[99960.45845, 29481.77196, 70478.6865], [2930... | [[0.0, 50409.38587, 0.0, 89141.04865, 0.0]] | [[988489.6284, -1932597.22, 0.0, 0.0]] | [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,... | [0.0] | [1277287.736] | [[296.8560381, 350.2088478, 279.0680247, 1.000... | [2425.07286, 8.463596169, 0.0, 0.0, 0.0, 8.463... | [32.10481199, 92112.53213, 131040.7862, 26.183... |
| 2 | 2000 | 4.0 | [-409167.9172, 392370.1191, -801538.0363, 1389... | [[100096.9564, 29163.29693, 70933.65948], [292... | [[0.0, 50155.81989, 0.0, 88764.59053, 0.0]] | [[983688.8106, -1924147.257, 0.0, 0.0]] | [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,... | [0.0] | [1277287.736] | [[297.2614008, 346.4257383, 280.8695397, 0.999... | [2425.07286, 8.463596169, 0.0, 0.0, 0.0, 8.463... | [28.23515616, 96520.43393, 130756.5894, 30.172... |
| 3 | 3000 | 6.0 | [-408468.813, 392911.3266, -801380.1396, 13890... | [[100572.3051, 29097.77284, 71474.53225], [292... | [[0.0, 50205.91055, 0.0, 88694.38663, 0.0]] | [[984891.3076, -1925171.744, 0.0, 0.0]] | [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,... | [0.0] | [1277287.736] | [[298.6730603, 345.647389, 283.0111843, 1.0000... | [2425.07286, 8.463596169, 0.0, 0.0, 0.0, 8.463... | [28.79320593, 96009.21744, 130922.0286, 24.551... |
| 4 | 4000 | 8.0 | [-409322.7261, 391652.3813, -800975.1074, 1398... | [[100279.2719, 28656.33854, 71622.93335], [291... | [[0.0, 50608.00218, 0.0, 89216.65927, 0.0]] | [[980597.3587, -1921397.128, 0.0, 0.0]] | [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,... | [0.0] | [1277287.736] | [[297.8028294, 340.4036676, 283.5987946, 0.999... | [2425.07286, 8.463596169, 0.0, 0.0, 0.0, 8.463... | [24.7954137, 100472.9563, 130538.2987, 27.3050... |
| 5 | 5000 | 10.0 | [-409060.5321, 391503.6168, -800564.1489, 1401... | [[99621.54627, 28465.00254, 71156.54373], [291... | [[0.0, 50668.20147, 0.0, 89507.19113, 0.0]] | [[980518.6316, -1921258.173, 0.0, 0.0]] | [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,... | [0.0] | [1277287.736] | [[295.8495588, 338.1308206, 281.7520741, 1.000... | [2425.07286, 8.463596169, 0.0, 0.0, 0.0, 8.463... | [24.44895665, 100856.1354, 130501.386, 20.8121... |
| 6 | 6000 | 12.0 | [-407724.5301, 393376.3395, -801100.8696, 1399... | [[100618.1056, 29001.93804, 71616.16754], [292... | [[0.0, 50612.96859, 0.0, 89338.70979, 0.0]] | [[980975.7104, -1922028.258, 0.0, 0.0]] | [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,... | [0.0] | [1277287.736] | [[298.8090756, 344.5089841, 283.5720046, 0.999... | [2425.07286, 8.463596169, 0.0, 0.0, 0.0, 8.463... | [24.88400774, 100946.7969, 131119.5628, 21.646... |
| 7 | 7000 | 14.0 | [-406045.9433, 393472.141, -799518.0842, 14009... | [[101200.4734, 29288.55636, 71911.91706], [292... | [[0.0, 50652.38951, 0.0, 89438.21007, 0.0]] | [[979207.7071, -1918816.391, 0.0, 0.0]] | [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,... | [0.0] | [1277287.736] | [[300.5385535, 347.9136733, 284.743057, 1.0000... | [2425.07286, 8.463596169, 0.0, 0.0, 0.0, 8.463... | [22.9661595, 103318.3573, 131165.6623, 27.4988... |
| 8 | 8000 | 16.0 | [-408499.0398, 391678.3093, -800177.3491, 1397... | [[100402.1573, 28648.71805, 71753.43928], [291... | [[0.0, 50465.54724, 0.0, 89294.44028, 0.0]] | [[977305.7837, -1917243.12, 0.0, 0.0]] | [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,... | [0.0] | [1277287.736] | [[298.1677666, 340.3131452, 284.115547, 0.9999... | [2425.07286, 8.463596169, 0.0, 0.0, 0.0, 8.463... | [19.87564957, 106444.6562, 130544.6054, 14.226... |
| 9 | 9000 | 18.0 | [-406991.4101, 392482.0724, -799473.4824, 1399... | [[100192.0289, 28084.21204, 72107.81681], [292... | [[0.0, 50693.74147, 0.0, 89209.25341, 0.0]] | [[976722.0947, -1916098.572, 0.0, 0.0]] | [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,... | [0.0] | [1277287.736] | [[297.5437407, 333.6074763, 285.5187434, 1.000... | [2425.07286, 8.463596169, 0.0, 0.0, 0.0, 8.463... | [18.24226614, 108700.5608, 130819.9731, 21.832... |
[24]:
[x for x in dir(tre) if not x.startswith("_")]
[24]:
['ENEVERSION',
'TITLE',
'add_traj',
'database',
'get_Hvap',
'get_baths',
'get_bondedContributions',
'get_density',
'get_eds',
'get_mass',
'get_nonbondedContributions',
'get_precalclam',
'get_specialContributions',
'get_temperature',
'get_temperature_Info',
'get_time_step',
'get_totals',
'get_totangle',
'get_totbonded',
'get_totcov',
'get_totcrf',
'get_totdihedral',
'get_totene',
'get_totkin',
'get_totlj',
'get_totnonbonded',
'get_totpot',
'path',
'step',
'time',
'tre_block_name_table',
'write']
Tre files contain all energy related data (like split up energy terms, temperature, pressure, …..). In PyGromos they generally share the same block structure as other files, but all the data inside the specific timesteps is stored efficiently inside a pandas DataFrame, here called tre.database . This database offers manipulation with all pandas functions. Alternatively many common functions are provided inside the Tre class.
This class should in principle replace further usage of the gromos++ ene_ana function, since all these operation can be done efficiently on the pandas DataFrame.
We are currently working on adding more common functions to the Tre class. If you find a useful function please contact the developers so the function can be added for general usage :)
2.2) Common Tre functions¶
[25]:
# calculate the average density over all timesteps
tre.get_density().mean()
[25]:
874.6072007296449
[26]:
# calculate the mean temperature over all frames for all baths in the system. In this example two baths with slightly different temperatures.
tre.get_temperature().mean()
[26]:
bath1 297.770285
bath2 297.713431
dtype: float64
Tables and lists inside the database are stored in numpy arrays. For example the two temperatures from the previous example are stored in a numpy array of size 2 since it has two temperature baths
Specific values inside a tre file can also be directly accessed with numpy and pandas syntax
[27]:
tre.database.iloc[2]
[27]:
step 2000
time 4.0
totals [-409167.9172, 392370.1191, -801538.0363, 1389...
baths [[100096.9564, 29163.29693, 70933.65948], [292...
bonded [[0.0, 50155.81989, 0.0, 88764.59053, 0.0]]
nonbonded [[983688.8106, -1924147.257, 0.0, 0.0]]
special [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,...
eds [0.0]
mass [1277287.736]
temperature [[297.2614008, 346.4257383, 280.8695397, 0.999...
volume [2425.07286, 8.463596169, 0.0, 0.0, 0.0, 8.463...
pressure [28.23515616, 96520.43393, 130756.5894, 30.172...
Name: 2, dtype: object
[28]:
# select the first nonbonded energy value for the first force group over all time frames
tre.database["nonbonded"].apply(lambda x: x[0][0])
[28]:
0 984602.6158
1 988489.6284
2 983688.8106
3 984891.3076
4 980597.3587
5 980518.6316
6 980975.7104
7 979207.7071
8 977305.7837
9 976722.0947
Name: nonbonded, dtype: float64
[29]:
tre.get_totals()
[29]:
| totene | totkin | totpot | totcov | totbond | totangle | totimproper | totdihedral | totcrossdihedral | totnonbonded | ... | totjval | totxray | totle | totorder | totsymm | eds_vr,entropy | totqm | totbsleus | totrdc | wip1 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| time | |||||||||||||||||||||
| 0.0 | -408806.2436 | 393521.8520 | -802328.0956 | 138403.8772 | 0.0 | 49931.35526 | 0.0 | 88472.52195 | 0.0 | -940731.9728 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 2.0 | -411576.4224 | 392980.7347 | -804557.1570 | 139550.4345 | 0.0 | 50409.38587 | 0.0 | 89141.04865 | 0.0 | -944107.5915 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 4.0 | -409167.9172 | 392370.1191 | -801538.0363 | 138920.4104 | 0.0 | 50155.81989 | 0.0 | 88764.59053 | 0.0 | -940458.4467 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 6.0 | -408468.8130 | 392911.3266 | -801380.1396 | 138900.2972 | 0.0 | 50205.91055 | 0.0 | 88694.38663 | 0.0 | -940280.4368 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 8.0 | -409322.7261 | 391652.3813 | -800975.1074 | 139824.6615 | 0.0 | 50608.00218 | 0.0 | 89216.65927 | 0.0 | -940799.7689 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 10.0 | -409060.5321 | 391503.6168 | -800564.1489 | 140175.3926 | 0.0 | 50668.20147 | 0.0 | 89507.19113 | 0.0 | -940739.5415 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 12.0 | -407724.5301 | 393376.3395 | -801100.8696 | 139951.6784 | 0.0 | 50612.96859 | 0.0 | 89338.70979 | 0.0 | -941052.5480 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 14.0 | -406045.9433 | 393472.1410 | -799518.0842 | 140090.5996 | 0.0 | 50652.38951 | 0.0 | 89438.21007 | 0.0 | -939608.6838 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 16.0 | -408499.0398 | 391678.3093 | -800177.3491 | 139759.9875 | 0.0 | 50465.54724 | 0.0 | 89294.44028 | 0.0 | -939937.3367 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 18.0 | -406991.4101 | 392482.0724 | -799473.4824 | 139902.9949 | 0.0 | 50693.74147 | 0.0 | 89209.25341 | 0.0 | -939376.4773 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
10 rows × 38 columns
\(\lambda\)-Sampling & TREs¶
[30]:
# import the trajectory file into a Tre class
tre = traj_tre.Tre(input_value="../example_files/Traj_files/RAFE_TI_l0_5.tre")
tre.get_precalclam()
/home/bschroed/Documents/projects/PyGromosTools/pygromos/files/trajectory/_general_trajectory.py:296: PerformanceWarning:
your performance may suffer as PyTables will pickle object types that it cannot
map directly to c-types [inferred_type->mixed,key->block2_values] [items->Index(['totals', 'baths', 'bonded', 'nonbonded', 'special', 'eds',
'precalclam', 'mass', 'temperature', 'volume', 'pressure'],
dtype='object')]
self.database.to_hdf(
[30]:
| nr_lambdas | A_e_lj_1 | B_e_lj_1 | A_e_crf_1 | B_e_crf_1 | AB_kinetic_1 | AB_bond_1 | AB_angle_1 | AB_improper_1 | AB_disres_1 | ... | B_e_crf_2 | AB_kinetic_2 | AB_bond_2 | AB_angle_2 | AB_improper_2 | AB_disres_2 | AB_dihres_2 | AB_disfld_2 | A_dihedral | B_dihedral | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| time | |||||||||||||||||||||
| 0.0 | 2.0 | -4.469426 | -53.669555 | -79.373179 | -52.887373 | 36.091797 | 0.0 | 0.0 | 0.0 | 0.0 | ... | -122.394320 | 36.091797 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 1.0 | 2.0 | -55.493090 | -58.228689 | -86.711496 | -52.410062 | 44.886645 | 0.0 | 0.0 | 0.0 | 0.0 | ... | -122.457590 | 44.886645 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 2.0 | 2.0 | -3.935587 | -54.492379 | -79.515189 | -54.142546 | 46.199981 | 0.0 | 0.0 | 0.0 | 0.0 | ... | -139.035309 | 46.199981 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 3.0 | 2.0 | -44.047130 | -54.945185 | -83.150923 | -53.673066 | 41.268221 | 0.0 | 0.0 | 0.0 | 0.0 | ... | -130.358158 | 41.268221 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 4.0 | 2.0 | -23.340313 | -52.961388 | -75.743516 | -53.920133 | 36.680211 | 0.0 | 0.0 | 0.0 | 0.0 | ... | -134.272063 | 36.680211 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 495.0 | 2.0 | 39.941549 | -64.577818 | -71.993096 | -53.485889 | 35.829309 | 0.0 | 0.0 | 0.0 | 0.0 | ... | -130.534976 | 35.829309 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 496.0 | 2.0 | -8.487283 | -54.840254 | -73.625420 | -54.112768 | 33.260161 | 0.0 | 0.0 | 0.0 | 0.0 | ... | -117.262832 | 33.260161 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 497.0 | 2.0 | 29.934914 | -54.557127 | -76.326046 | -51.811587 | 41.735086 | 0.0 | 0.0 | 0.0 | 0.0 | ... | -104.326673 | 41.735086 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 498.0 | 2.0 | -35.691879 | -60.330236 | -77.859153 | -54.295433 | 44.944847 | 0.0 | 0.0 | 0.0 | 0.0 | ... | -123.231158 | 44.944847 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 499.0 | 2.0 | 21.308545 | -58.470225 | -72.376205 | -53.554977 | 41.740483 | 0.0 | 0.0 | 0.0 | 0.0 | ... | -125.539239 | 41.740483 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
500 rows × 25 columns
EDS in TREs¶
[31]:
# import the trajectory file into a Tre class
tre = traj_tre.Tre(input_value="../example_files/Traj_files/RAFE_eds.tre")
tre.get_eds()
/home/bschroed/Documents/projects/PyGromosTools/pygromos/files/trajectory/_general_trajectory.py:296: PerformanceWarning:
your performance may suffer as PyTables will pickle object types that it cannot
map directly to c-types [inferred_type->mixed,key->block2_values] [items->Index(['totals', 'baths', 'bonded', 'nonbonded', 'special', 'eds',
'precalclam', 'mass', 'temperature', 'volume', 'pressure'],
dtype='object')]
self.database.to_hdf(
[31]:
| numstates | total_1 | nonbonded_1 | special_1 | offset_1 | total_2 | nonbonded_2 | special_2 | offset_2 | total_3 | ... | special_7 | offset_7 | total_8 | nonbonded_8 | special_8 | offset_8 | total_9 | nonbonded_9 | special_9 | offset_9 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| time | |||||||||||||||||||||
| 0.00 | 9.0 | -304.241276 | -304.241276 | 0.0 | 0.0 | -309.738568 | -309.738568 | 0.0 | 0.0 | -286.397004 | ... | 0.0 | 0.0 | -306.757290 | -306.757290 | 0.0 | 0.0 | 1.661115e+05 | 1.661115e+05 | 0.0 | 0.0 |
| 0.04 | 9.0 | -327.432243 | -327.432243 | 0.0 | 0.0 | -343.990987 | -343.990987 | 0.0 | 0.0 | -330.156886 | ... | 0.0 | 0.0 | -337.505471 | -337.505471 | 0.0 | 0.0 | 2.114028e+05 | 2.114028e+05 | 0.0 | 0.0 |
| 0.08 | 9.0 | -342.969913 | -342.969913 | 0.0 | 0.0 | -348.338017 | -348.338017 | 0.0 | 0.0 | -342.678812 | ... | 0.0 | 0.0 | -343.085575 | -343.085575 | 0.0 | 0.0 | 2.335587e+05 | 2.335587e+05 | 0.0 | 0.0 |
| 0.12 | 9.0 | -307.787722 | -307.787722 | 0.0 | 0.0 | -320.642022 | -320.642022 | 0.0 | 0.0 | -275.892803 | ... | 0.0 | 0.0 | -287.546356 | -287.546356 | 0.0 | 0.0 | 1.263949e+05 | 1.263949e+05 | 0.0 | 0.0 |
| 0.16 | 9.0 | -325.101329 | -325.101329 | 0.0 | 0.0 | -333.832261 | -333.832261 | 0.0 | 0.0 | -309.882401 | ... | 0.0 | 0.0 | -328.133697 | -328.133697 | 0.0 | 0.0 | 3.637117e+05 | 3.637117e+05 | 0.0 | 0.0 |
| 0.20 | 9.0 | -341.014069 | -341.014069 | 0.0 | 0.0 | -338.126642 | -338.126642 | 0.0 | 0.0 | -319.108916 | ... | 0.0 | 0.0 | -331.591241 | -331.591241 | 0.0 | 0.0 | 3.426647e+06 | 3.426647e+06 | 0.0 | 0.0 |
| 0.24 | 9.0 | -332.591891 | -332.591891 | 0.0 | 0.0 | -313.860613 | -313.860613 | 0.0 | 0.0 | -270.145527 | ... | 0.0 | 0.0 | -295.307547 | -295.307547 | 0.0 | 0.0 | 6.720166e+05 | 6.720166e+05 | 0.0 | 0.0 |
| 0.28 | 9.0 | -384.741762 | -384.741762 | 0.0 | 0.0 | -379.723111 | -379.723111 | 0.0 | 0.0 | -340.573094 | ... | 0.0 | 0.0 | -363.097220 | -363.097220 | 0.0 | 0.0 | 6.147874e+04 | 6.147874e+04 | 0.0 | 0.0 |
| 0.32 | 9.0 | -341.408526 | -341.408526 | 0.0 | 0.0 | -352.793311 | -352.793311 | 0.0 | 0.0 | -268.495614 | ... | 0.0 | 0.0 | -327.075935 | -327.075935 | 0.0 | 0.0 | 2.863084e+04 | 2.863084e+04 | 0.0 | 0.0 |
| 0.36 | 9.0 | -336.573404 | -336.573404 | 0.0 | 0.0 | -340.350887 | -340.350887 | 0.0 | 0.0 | -238.003947 | ... | 0.0 | 0.0 | -325.957742 | -325.957742 | 0.0 | 0.0 | 3.048352e+04 | 3.048352e+04 | 0.0 | 0.0 |
10 rows × 37 columns
Concatenate and Copy multiple Trajectories¶
Trajectories offer a wide range of additional file manipulations. Trajectory classes can be copied (deep) and added to each other to concatenate multiple small simulation pieces into one large trajectory.
[32]:
tre_copy = traj_tre.Tre(input_value=tre)
[33]:
tre_copy.database.shape
[33]:
(10, 13)
[34]:
tre_combined = tre + tre_copy
/home/bschroed/Documents/projects/PyGromosTools/pygromos/files/trajectory/_general_trajectory.py:180: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
new_traj.database = self.database.append(new_data, ignore_index=True)
[35]:
tre_combined.database.shape
[35]:
(20, 13)
In the new combined trajectory we have one long trajectory made from the two smaller ones. The length is one element shorter, since normally the last element of the first trajectory and the first element of the second trajectory is the same element. This can be controlled via the option “skip_new_0=True” in the add_traj() function which is the core of the “+” operator for trajectories. In the following line the default behavior can be seen as a smooth numbering in the TIMESTEPs.
[36]:
tre_combined.database.time
[36]:
0 0.00
1 0.04
2 0.08
3 0.12
4 0.16
5 0.20
6 0.24
7 0.28
8 0.32
9 0.36
10 0.00
11 0.04
12 0.08
13 0.12
14 0.16
15 0.20
16 0.24
17 0.28
18 0.32
19 0.36
Name: time, dtype: float64
[37]:
print(len(tre_combined.database), len(tre.database))
20 10
[ ]: