Mercurial > repos > chemteam > mdanalysis_extract_rmsd
comparison extract_rmsd.py @ 2:589f8ef21e58 draft default tip
"planemo upload for repository https://github.com/galaxycomputationalchemistry/galaxy-tools-compchem/ commit f1c3c88c7395f2e84cbc533199406aadb79c5c07"
| author | chemteam |
|---|---|
| date | Fri, 13 Nov 2020 19:44:09 +0000 |
| parents | 8f6ad93973cb |
| children |
comparison
equal
deleted
inserted
replaced
| 1:8f6ad93973cb | 2:589f8ef21e58 |
|---|---|
| 34 if sum(1 for line in f1) != sum(1 for line in f2): | 34 if sum(1 for line in f1) != sum(1 for line in f2): |
| 35 raise IOError('Number of structure and trajectory files unequal.') | 35 raise IOError('Number of structure and trajectory files unequal.') |
| 36 | 36 |
| 37 no_t = len(traj_file_list) | 37 no_t = len(traj_file_list) |
| 38 | 38 |
| 39 data = np.zeros((no_t, no_t, | 39 # hard to find array size before loading files |
| 40 int((end - start)/step + ((end - start) % step > 0)))) | 40 universe_coordinate_data = [] |
| 41 | |
| 42 # load files | |
| 43 universes = {} | |
| 44 | 41 |
| 45 for traj in range(no_t): | 42 for traj in range(no_t): |
| 46 # We no longer align here, users should do this themselves. | 43 # We no longer align here, users should do this themselves. |
| 47 universes[traj] = m.Universe(str_file_list[traj], traj_file_list[traj], | 44 u = m.Universe(str_file_list[traj], traj_file_list[traj], |
| 48 format=traj_format, | 45 format=traj_format, topology_format=str_format) |
| 49 topology_format=str_format) | 46 u.transfer_to_memory() |
| 47 grp = u.select_atoms(group).universe | |
| 48 coordinates = grp.trajectory.coordinate_array[start:end:step] | |
| 49 universe_coordinate_data.append(coordinates) | |
| 50 | 50 |
| 51 universe_coordinate_data = np.array(universe_coordinate_data) | |
| 51 print("All trajs loaded by MDAnalysis") | 52 print("All trajs loaded by MDAnalysis") |
| 53 data = np.zeros((no_t, no_t, universe_coordinate_data.shape[1])) | |
| 52 | 54 |
| 53 # calculate differences | 55 # calculate differences |
| 54 for traj1 in range(no_t): | 56 for traj1 in range(no_t): |
| 55 print("Calculating differences for traj {}".format(traj1)) | 57 print("Calculating differences for traj {}".format(traj1)) |
| 56 for traj2 in range(traj1): | 58 for traj2 in range(traj1): |
| 57 for frame in range(data.shape[2]): | 59 for frame in range(data.shape[2]): |
| 58 universes[traj1].trajectory[frame] | 60 A = universe_coordinate_data[traj1, frame] |
| 59 universes[traj2].trajectory[frame] | 61 B = universe_coordinate_data[traj2, frame] |
| 60 A = universes[traj1].select_atoms(group).positions | |
| 61 B = universes[traj2].select_atoms(group).positions | |
| 62 r = rms.rmsd(A, B) | 62 r = rms.rmsd(A, B) |
| 63 data[traj1, traj2, frame] = r | 63 data[traj1, traj2, frame] = r |
| 64 data[traj2, traj1, frame] = r | 64 data[traj2, traj1, frame] = r |
| 65 | 65 |
| 66 with open(filepath_out, 'w') as f: | 66 with open(filepath_out, 'w') as f: |
| 88 parser.add_argument('--end', type=int, | 88 parser.add_argument('--end', type=int, |
| 89 help="Last trajectory frame to calculate RMSD") | 89 help="Last trajectory frame to calculate RMSD") |
| 90 parser.add_argument('--step', type=int, | 90 parser.add_argument('--step', type=int, |
| 91 help="Frame sampling frequency for RMSD calculation") | 91 help="Frame sampling frequency for RMSD calculation") |
| 92 args = parser.parse_args() | 92 args = parser.parse_args() |
| 93 | |
| 94 calc_rmsd(args.strs, args.trajs, args.str_format, | 93 calc_rmsd(args.strs, args.trajs, args.str_format, |
| 95 args.traj_format, args.outfile, | 94 args.traj_format, args.outfile, |
| 96 args.group, args.start, args.end, args.step) | 95 args.group, args.start, args.end, args.step) |
| 97 | 96 |
| 98 | 97 |
