Mercurial > repos > chemteam > mdanalysis_rdf

comparison extract_rmsd.py @ 9:49076afc90b1 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/galaxycomputationalchemistry/galaxy-tools-compchem/ commit f1c3c88c7395f2e84cbc533199406aadb79c5c07"
author chemteam
date Fri, 13 Nov 2020 19:41:36 +0000
parents 4510baed86ee
children
comparison
equal deleted inserted replaced
8:4510baed86ee 9:49076afc90b1
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