-
Notifications
You must be signed in to change notification settings - Fork 0
/
parse_nmr_g16
executable file
·157 lines (125 loc) · 4.67 KB
/
parse_nmr_g16
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
import sys
import re
import csv
from math import exp
try:
import typer
except ImportError:
sys.exit("typer not installed, install it with 'pip install typer[all]")
# constants
BOLTZMANN = 0.0019872041
def open_file(file_name):
with open(file_name) as f:
opened_file = f.readlines()
return opened_file
def get_nmr_ranges(file):
line_ranges = []
for line in range(0, len(file)):
if "Total nuclear spin-spin coupling J" in file[line]:
for i in range(line, len(file)):
if "End of Minotr F.D" in file[i]:
break
line_ranges.append((line + 2, i))
# check if coupling data was found
if len(line_ranges) == 0:
sys.exit("No coupling data found in the NMR file provided, exiting...")
else:
print("Found {} conformers in the NMR file.".format(len(line_ranges)))
return line_ranges
def parse_tables(file):
nmr_file = open_file(file)
line_ranges = get_nmr_ranges(nmr_file)
tables = []
for line_range in line_ranges: # line_range is a tuple (start,end)
table = {}
formatted_lines = []
for line in nmr_file[line_range[0]:line_range[1]]:
if line.startswith(" "):
continue
else:
nline = re.sub("D", "E", re.sub(" ", " ", line.strip()))
formatted_lines.append([float(a) for a in nline.split(" ")])
for formatted_line in formatted_lines:
if formatted_line[0] in table.keys():
table[formatted_line[0]] = table[formatted_line[0]] + \
formatted_line[1::]
else:
table[formatted_line[0]] = formatted_line[1::]
tables.append(table)
return tables
def write_conformers(tables, filename):
"""
From a dict to a .csv file with all the conformers
"""
# first save all the conformers
with open("conformers.csv", "w") as f:
write = csv.writer(f)
for table in tables:
atom_nrs = list(table.keys())
write.writerow(["atom nr"] + [int(a) for a in atom_nrs])
for row in atom_nrs:
write.writerow([int(row)] + table[row])
print("Conformer data saved to conformers.csv")
def get_energies(energy_file):
energies = []
for line in energy_file:
if "Sum of electronic and thermal Free Energies" in line:
energy = line.split("=")[1].lstrip().rstrip()
energies.append(float(energy))
# check if freq data was found
if len(energies) == 0:
sys.exit("No freq data found in the energy file provided, exiting...")
else:
print("Found {} conformers in the energy file.".format(len(energies)))
return energies
def get_populations(energy_file, temperature=298.15):
energy_lines = open_file(energy_file)
populations = []
energies = get_energies(energy_lines)
minima = min(energies)
for i in range(0, len(energies)):
energies[i] = (energies[i] - minima) * 627.503
for energy in energies:
pop = exp(-energy/(BOLTZMANN*temperature))
populations.append(pop)
total_sum = sum(populations)
for i in range(0, len(populations)):
populations[i] = populations[i]/total_sum
return populations
def average_coupling(tables, populations, save=True):
final_table = {}
print(len(populations), len(tables))
if len(populations) != len(tables):
sys.exit(
"Number of conformers in energy file is different from NMR file, exiting...")
for i in range(0, len(tables)):
keys = tables[i].keys()
if i == 0:
for line in keys:
final_table[line] = [coupling*populations[i]
for coupling in tables[i][line]]
else:
for line in keys:
final_table[line] = [final_table[line][j] +
tables[i][line][j]*populations[i]
for j in range(len(final_table[line]))]
if save:
with open("averaged.csv", "w") as f:
write = csv.writer(f)
atom_nrs = list(final_table.keys())
write.writerow(["atom nr"] + [int(a) for a in atom_nrs])
for row in atom_nrs:
write.writerow([int(row)] + final_table[row])
def main(
nmr_file: str,
energy_file: str,
temperature: float = 298.15,
save_conformers: bool = True
):
populations = get_populations(energy_file)
tables = parse_tables(nmr_file)
average_coupling(tables, populations)
if save_conformers:
write_conformers(tables, nmr_file)
if __name__ == "__main__":
typer.run(main)