Context Navigation

source: trunk/G2shapes.py @ 4650

Last change on this file since 4650 was 4339, checked in by toby, 5 years ago
set svn flags; improve compare F-test & add info & covariance plot
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 64.6 KB

Line
1	#
2	# Copyright v1.0, 1.2, 1.3: 2019, John Badger.
3	#
4	# This program is free software: you can redistribute it and/or modify
5	# it under the terms of the GNU General Public License as published by
6	# the Free Software Foundation, either version 3 of the License, or
7	# (at your option) any later version.
8	#
9	# This program is distributed in the hope that it will be useful,
10	# but WITHOUT ANY WARRANTY; without even the implied warranty of
11	# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12	# GNU General Public License for more details.
13	#
14	# You should have received a copy of the GNU General Public License
15	# along with this program. If not, see <https://www.gnu.org/licenses/>.
16	#
17
18	# Version 1.2 is intended to be runnable under Python2 and Python3
19	# Version 1.3 includes an optional 'glue' term for extended structures
20	#
21	# this version modified by R B. Von Dreele for inclusion in GSAS-II
22
23	from __future__ import division, print_function
24	import math
25	import sys
26	import os
27	import copy
28	import random
29	import time
30	import cProfile,pstats
31	import io as StringIO
32	import numpy as np
33	nxs = np.newaxis
34
35	def G2shapes(Profile,ProfDict,Limits,data,dlg=None):
36
37	########## FUNCTIONS ########################
38
39	# NEW 1.1 Calculate intensity from P(r) dropping all scaling factors
40
41	def ft_to_intensity(aList_q,aList_i_calc,aList_r,aList_pr_model,nbeads):
42
43	num_q = len(aList_q)
44	num_r = len(aList_r)
45
46	count_q = 0
47	while count_q < num_q:
48
49	q = float(aList_q[count_q])
50
51	# Sets P(r=0) =0.0. Later scaling includes a baseline term.
52	integral = 0.0
53
54	count_r = 1
55	while count_r < num_r:
56
57	r = float(aList_r[count_r])
58	pr = float(aList_pr_model[count_r])
59	qr = q*r
60	integral = integral + pr*math.sin(qr)/qr
61
62	count_r = count_r + 1
63
64	aList_i_calc.append(integral)
65
66	count_q = count_q + 1
67
68	return;
69
70	# NEW 1.1 Scale and Compare I and Ic. Includes a baseline correction
71
72	def score_Ic(aList_q,aList_i,aList_i_sd,aList_i_calc):
73
74	num_q = len(aList_q)
75
76	idif = 0.0
77	isum = 0.0
78	sd_sq = 0.0
79	chi_sq = 0.0
80
81	# Least squares scale for calculated I onto observed I
82
83	S = 0.0
84	Sx = 0.0
85	Sy = 0.0
86	Sxx = 0.0
87	Sxy = 0.0
88
89	count = 0
90	while count < num_q:
91	x = float(aList_i_calc[count])
92	y = float(aList_i[count])
93
94	S = S + 1.0
95	Sx = Sx + x
96	Sy = Sy + y
97	Sxx = Sxx + x*x
98	Sxy = Sxy + x*y
99
100	count = count + 1
101
102	delta = SSxx - SxSx
103	a = (SxxSy - SxSxy)/delta
104	b = (SSxy - SxSy)/delta
105
106	# Apply scale and find statistics
107
108	i = 0
109	while i < num_q:
110	iobs = float(aList_i[i])
111	sd = float(aList_i_sd[i])
112
113	aList_i_calc[i] = b*float(aList_i_calc[i]) + a
114	icalc = aList_i_calc[i]
115
116	idif = idif + abs(iobs - icalc)
117	isum = isum + iobs + icalc
118
119	dif = iobs - icalc
120	dif_sq = dif*dif
121	sd_sq = sd*sd
122
123	chi_sq = chi_sq + dif_sq/sd_sq
124
125	i = i + 1
126
127	rvalue = 2.0*idif/isum
128	chi_sq = chi_sq/num_q
129
130	return (chi_sq,rvalue);
131
132	# NEW 1.1 Write original and calculated data.
133
134	def write_all_data(file_intensity,aList_q,aList_i,aList_i_calc,aString):
135
136	num_q = len(aList_q)
137
138	file = open(file_intensity,'w',)
139	aString = '# ' + aString + '\n'
140	file.write(aString)
141
142	i = 0
143	while i < num_q:
144
145	q = aList_q[i]
146	intensity = aList_i[i]
147	intensity_calc = aList_i_calc[i]
148
149	aString = str(q) + ' ' + str(intensity) + ' ' + str(intensity_calc) + '\n'
150
151	file.write(aString)
152
153	i = i + 1
154
155	file.close()
156
157	return;
158
159	# NEW 1.1 Read intensity data from GNOM output file
160
161	def read_i(aList_q,aList_i,aList_i_sd,inFile,angstrom_scale):
162
163	scale_units = 1.0/angstrom_scale
164
165	Q,Io,wt,Ic,Ib,Ifb = Profile[:6]
166	Qmin = Limits[1][0]
167	Qmax = Limits[1][1]
168	wtFactor = ProfDict['wtFactor']
169	Back,ifBack = data['Back']
170	Ibeg = np.searchsorted(Q,Qmin)
171	Ifin = np.searchsorted(Q,Qmax)+1 #include last point
172	aList_q += list(Q[Ibeg:Ifin]*scale_units)
173	aList_i += list(Io[Ibeg:Ifin])
174	aList_i_sd += list(1./np.sqrt(wtFactor*wt[Ibeg:Ifin]))
175
176	# file = open(inFile,'r')
177	# allLines = file.readlines()
178	# file.close()
179	#
180	# parse_data = 'no'
181	# for eachLine in allLines:
182	#
183	# if parse_data == 'yes':
184	#
185	# aList = eachLine.split()
186	# num_params = len(aList)
187	#
188	# if num_params == 5:
189	#
190	# q = float(aList[0]) * scale_units
191	# if q > 0.0:
192	# i = float(aList[1])
193	# i_sd = float(aList[2])
194	# aList_q.append(q)
195	# aList_i.append(i)
196	# aList_i_sd.append(i_sd)
197	#
198	# else:
199	# if num_params == 0 and len(aList_q) > 0:
200	# parse_data = 'no'
201	#
202	# if eachLine.find('S J EXP ERROR J REG I REG') > -1:
203	# parse_data = 'yes'
204	#
205	return;
206
207	# Read PDB for starting point
208
209	def read_pdb(aList_beads_x,aList_beads_y,aList_beads_z,pdbfile_in):
210
211	xmean = 0.0
212	ymean = 0.0
213	zmean = 0.0
214	nbeads = 0
215
216	file = open(pdbfile_in,'r')
217	allLines = file.readlines()
218	file.close()
219
220	for eachLine in allLines:
221
222	tag = eachLine[0:6]
223	tag = tag.strip()
224
225	if tag == 'ATOM':
226
227	atom_name = eachLine[13:16]
228	atom_name = atom_name.strip()
229
230	if atom_name == 'CA':
231
232	x = float(eachLine[30:38])
233	y = float(eachLine[38:46])
234	z = float(eachLine[46:54])
235
236	xmean = xmean + x
237	ymean = ymean + y
238	zmean = zmean + z
239
240	nbeads = nbeads + 1
241
242	xmean = xmean/float(nbeads)
243	ymean = ymean/float(nbeads)
244	zmean = zmean/float(nbeads)
245
246	for eachLine in allLines:
247
248	tag = eachLine[0:6]
249	tag = tag.strip()
250
251	if tag == 'ATOM':
252
253	atom_name = eachLine[13:16]
254	atom_name = atom_name.strip()
255
256	if atom_name == 'CA':
257
258	x = float(eachLine[30:38]) - xmean
259	y = float(eachLine[38:46]) - ymean
260	z = float(eachLine[46:54]) - zmean
261	aList_beads_x.append(x)
262	aList_beads_y.append(y)
263	aList_beads_z.append(z)
264
265	return;
266
267	# # Write P(r) with SD and calculated value.
268	#
269	# def pr_writer(aList_pr,aList_r,aList_pr_model,outfile_pr):
270	#
271	# num_pr = len(aList_pr)
272	#
273	# file = open(outfile_pr,'w')
274	# file.write('#\n')
275	#
276	# i = 0
277	# while i < num_pr:
278	#
279	# r = float(aList_r[i])
280	# pr = float(aList_pr[i])
281	# pr_calc = float(aList_pr_model[i])
282	# aString = str(r) + ' ' + str(pr) + ' ' + str(pr_calc) + '\n'
283	# file.write(aString)
284	#
285	# i = i + 1
286	#
287	# file.close()
288	#
289	# return;
290
291	# Write a set of points as a pseudo-PDB file
292
293	def pdb_writer(aList_x_write,aList_y_write,aList_z_write,out_file,aString):
294
295	atom_number = 0
296	res_number = 0
297	dummy_b = 0
298	num_atoms = len(aList_x_write)
299
300	file = open(out_file,'w')
301	file.write(aString)
302	file.write('\n')
303
304	i = 0
305	while i < num_atoms:
306
307	x = float(aList_x_write[i])
308	y = float(aList_y_write[i])
309	z = float(aList_z_write[i])
310	x = '%.3f'%(x)
311	y = '%.3f'%(y)
312	z = '%.3f'%(z)
313	x = x.rjust(8)
314	y = y.rjust(8)
315	z = z.rjust(8)
316	atom_number = atom_number + 1
317	res_number = res_number + 1
318	atom_number_str = str(atom_number)
319	atom_number_str = atom_number_str.rjust(5)
320	res_number_str = str(res_number)
321	res_number_str = res_number_str.rjust(4)
322	bfactor = str(dummy_b) + '.00'
323	bfactor = bfactor.rjust(6)
324
325	if res_number < 10000:
326	aLine_out = 'ATOM ' + atom_number_str + ' CA ALA A' + res_number_str + ' ' + x + y + z + ' 1.00' + bfactor + '\n'
327	else:
328	res_number_str = str(res_number - 9999)
329	aLine_out = 'ATOM ' + atom_number_str + ' CA ALA B' + res_number_str + ' ' + x + y + z + ' 1.00' + bfactor + '\n'
330	file.write(aLine_out)
331
332	i = i + 1
333
334	file.close()
335
336	return;
337
338	# Evaluate local bead densities and point density on a notional grid
339
340	def set_box(aList_beads_x,aList_beads_y,aList_beads_z,\
341	aList_box_x_all,aList_box_y_all,aList_box_z_all,\
342	aList_box_score,box_step,dmax,rsearch):
343
344	dmax_over2 = dmax/2.0
345	search_sq = rsearch**2
346	num_beads = len(aList_beads_x)
347
348	count_x = -dmax_over2
349	while count_x < dmax_over2:
350	count_y = -dmax_over2
351	while count_y < dmax_over2:
352	count_z = -dmax_over2
353	while count_z < dmax_over2:
354
355	count_local = 0
356	i = 0
357	while i < num_beads:
358
359	dx = float(aList_beads_x[i]) - count_x
360	dy = float(aList_beads_y[i]) - count_y
361	dz = float(aList_beads_z[i]) - count_z
362
363	dsq = dxdx + dydy + dz*dz
364
365	if dsq < search_sq:
366	count_local = count_local + 1
367
368	i = i + 1
369
370	if count_local > 1:
371	aList_box_x_all.append(count_x)
372	aList_box_y_all.append(count_y)
373	aList_box_z_all.append(count_z)
374	aList_box_score.append(count_local)
375
376	count_z = count_z + box_step
377	count_y = count_y + box_step
378	count_x = count_x + box_step
379
380	return;
381
382	# Evaluate local bead densities and point density on a notional grid - fast version
383
384	def set_box_fast(aList_beads_x,aList_beads_y,aList_beads_z,\
385	aList_box_x_all,aList_box_y_all,aList_box_z_all,\
386	aList_box_score,box_step,dmax,rsearch):
387
388	dmax_over2 = dmax/2.0
389	num_box = int(dmax/box_step)
390	search_sq = rsearch**2
391
392	XYZ = np.meshgrid(np.linspace(-dmax_over2,dmax_over2,num_box),
393	np.linspace(-dmax_over2,dmax_over2,num_box),
394	np.linspace(-dmax_over2,dmax_over2,num_box))
395	XYZ = np.array([XYZ[0].flatten(),XYZ[1].flatten(),XYZ[2].flatten()]).T
396	xyz = np.array((aList_beads_y,aList_beads_x,aList_beads_z)).T
397	for XYZi in XYZ:
398	dsq = np.sum((xyz-XYZi)**2,axis=1)
399	count = int(np.sum(np.where(dsq<search_sq,1,0)))
400	if count>1:
401	aList_box_x_all.append(XYZi[0])
402	aList_box_y_all.append(XYZi[1])
403	aList_box_z_all.append(XYZi[2])
404	aList_box_score.append(count)
405	return;
406
407	# Establish a volume
408
409	def set_vol(aList_box_x_all,aList_box_y_all,aList_box_z_all,aList_box_score,\
410	aList_box_x,aList_box_y,aList_box_z,vol_target,box_pt_vol):
411
412	num_pts = len(aList_box_score)
413	num_tries = int(max(aList_box_score))
414	density_thresh = max(aList_box_score) - 1.0
415	vol = vol_target + 1.0
416
417	i = 0
418	while i < num_tries:
419
420	density_thresh = density_thresh - 1.0
421	num_box_pts = 0.0
422
423	j = 0
424	while j < num_pts:
425	density = float(aList_box_score[j])
426	if density >= density_thresh:
427	num_box_pts = num_box_pts + 1.0
428	j = j + 1
429
430	vol = num_box_pts*box_pt_vol
431	if vol < vol_target:
432	density_use = density_thresh
433
434	i = i + 1
435
436	#
437
438	num_box_pts1 = 0.0
439	num_box_pts2 = 0.0
440	density_thresh1 = density_use
441	density_thresh2 = density_use - 1.0
442	i = 0
443	while i < num_pts:
444	density_value = float(aList_box_score[i])
445	if density_value >= density_thresh1:
446	num_box_pts1 = num_box_pts1 + 1.0
447	if density_value >= density_thresh2:
448	num_box_pts2 = num_box_pts2 + 1.0
449	i = i + 1
450
451	vol1 = num_box_pts1*box_pt_vol
452	vol2 = num_box_pts2*box_pt_vol
453	delta1 = abs(vol1-vol_target)
454	delta2 = abs(vol2-vol_target)
455
456	if delta1 < delta2:
457	density_thresh = density_thresh1
458	else:
459	density_thresh = density_thresh2
460
461	i = 0
462	while i < num_pts:
463
464	density_value = float(aList_box_score[i])
465	if density_value >= density_thresh:
466	aList_box_x.append(aList_box_x_all[i])
467	aList_box_y.append(aList_box_y_all[i])
468	aList_box_z.append(aList_box_z_all[i])
469
470	i = i + 1
471
472	return;
473
474	# Find beads that are not in allowed volume
475
476	def disallowed_beads(aList_beads_x,aList_beads_y,aList_beads_z,aList_contacts,\
477	aList_box_x,aList_box_y,aList_box_z,rsearch):
478
479	num_beads = len(aList_beads_x)
480	num_boxes = len(aList_box_x)
481	contact_limit_sq = rsearch**2
482
483	count = 0
484	while count < num_beads:
485
486	x_bead = float(aList_beads_x[count])
487	y_bead = float(aList_beads_y[count])
488	z_bead = float(aList_beads_z[count])
489
490	inbox = 'no'
491	i = 0
492	while i < num_boxes:
493
494	x_box = float(aList_box_x[i])
495	y_box = float(aList_box_y[i])
496	z_box = float(aList_box_z[i])
497	dsq = (x_bead - x_box)2 + (y_bead - y_box)2 + (z_bead - z_box)**2
498
499	if dsq < contact_limit_sq:
500	inbox = 'yes'
501	i = num_boxes
502
503	i = i + 1
504
505	if inbox == 'no':
506	aList_contacts.append(count)
507
508	count = count + 1
509
510	return;
511
512	# Compute a P(r)
513
514	def calc_pr(aList_beads_x,aList_beads_y,aList_beads_z,aList_pr_model,hist_grid):
515
516	num_hist = len(aList_pr_model)
517	count = 0
518	while count < num_hist:
519	aList_pr_model[count] = 0.0
520	count = count + 1
521
522	nbeads = len(aList_beads_x)
523	max_dr = (float(num_hist)-1.0)*hist_grid
524
525	i = 0
526	while i < nbeads:
527
528	j = 0
529	while j < i:
530
531	dr = get_dr(aList_beads_x[i],aList_beads_y[i],aList_beads_z[i],\
532	aList_beads_x[j],aList_beads_y[j],aList_beads_z[j])
533	dr = min(dr,max_dr)
534
535	# Find pointers and do un-interpolation
536
537	dr_grid = dr/hist_grid
538	int_dr_grid = int(dr_grid)
539
540	int_dr_grid = min(int_dr_grid,num_hist-2)
541	ip_low = int_dr_grid
542	ip_high = ip_low + 1
543
544	ip_high_frac = dr_grid - float(int_dr_grid)
545	ip_low_frac = 1.0 - ip_high_frac
546
547	aList_pr_model[ip_low] = float(aList_pr_model[ip_low]) + ip_low_frac
548	aList_pr_model[ip_high] = float(aList_pr_model[ip_high]) + ip_high_frac
549
550	j = j + 1
551	i = i + 1
552
553	return;
554
555	# Score for rms difference between observed and model histograms
556
557	def pr_dif(aList_pr,aList_pr_model,skip):
558
559	num_hist = len(aList_pr)
560	delta_hist_sum = 0.0
561	delta_hist_sum_sq = 0.0
562	hist_sum = 0.0
563
564	i = skip
565	while i < num_hist:
566
567	model = float(aList_pr_model[i])
568	data = float(aList_pr[i])
569	delta_hist = abs(data - model)
570	delta_hist_sum = delta_hist_sum + delta_hist
571	hist_sum = hist_sum + data
572
573	delta_hist_sum_sq = delta_hist_sum_sq + delta_hist*delta_hist
574
575	i = i + 1
576
577	mean_hist_sum = hist_sum/(num_hist - skip)
578	delta_hist_sum_sq = delta_hist_sum_sq/(num_hist - skip)
579	delta_hist_sum_sq = math.sqrt(delta_hist_sum_sq)/mean_hist_sum
580
581	return delta_hist_sum_sq;
582
583	# Statistics for fractional difference between observed and model histograms
584
585	def pr_rfactor(aList_pr,aList_pr_sd,aList_pr_model,skip):
586
587	num_hist = len(aList_pr)
588	delta_hist_sum = 0.0
589	hist_sum = 0.0
590
591	i = skip
592	while i < num_hist:
593
594	model = float(aList_pr_model[i])
595	exp = float(aList_pr[i])
596	delta_hist = exp - model
597	delta_hist_sum = delta_hist_sum + abs(delta_hist)
598	hist_sum = hist_sum + exp
599
600	i = i + 1
601
602	delta_hist_sum = delta_hist_sum/hist_sum
603
604	return delta_hist_sum;
605
606	# Compute the VDW energy for a interaction
607
608	def vdw_energy(econ12,econ6,e_width,dr,bead_sep3):
609
610	if dr > bead_sep3:
611	vdw = 0.0
612	else:
613	dr_e6 = dr**6
614	dr_e12 = dr_e6**2
615	vdw = econ12/dr_e12 - 2.0*econ6/dr_e6
616	vdw = max(vdw,e_width)
617
618	return vdw;
619
620	def vdw_energies(econ12,econ6,e_width,drs,bead_sep3):
621
622	drs_e6 = drs**6
623	drs_e12 = drs_e6**2
624	vdws = econ12/drs_e12 - 2.0*econ6/drs_e6
625	vdws = np.where(drs>bead_sep3,0.,vdws)
626	vdws = np.where(vdws>e_width,vdws,e_width)
627	return vdws
628
629	# Set a random distribution of beads in a box with maximum extent dmax
630
631	def random_beads(aList_beads_x,aList_beads_y,aList_beads_z,\
632	nbeads,dmax,aList_symm,bias_z):
633
634	half_side = 0.5*dmax
635
636	scale_xy = 1.0 - bias_z
637	scale_z = 1.0 + bias_z
638	x_range = scale_xy * half_side
639	y_range = scale_xy * half_side
640	z_range = scale_z * half_side
641
642	num_ops = len(aList_symm)
643
644	i = 0
645	while i < nbeads:
646
647	triangle = random.triangular(-0.7,0.7,0.0)
648	x = triangle*x_range
649	triangle = random.triangular(-0.7,0.7,0.0)
650	y = triangle*y_range
651	triangle = random.triangular(-0.7,0.7,0.0)
652	z = triangle*z_range
653
654	aList_beads_x.append(x)
655	aList_beads_y.append(y)
656	aList_beads_z.append(z)
657
658	j = 0
659	while j < num_ops:
660	aList_s = aList_symm[j]
661	m11 = float(aList_s[0])
662	m12 = float(aList_s[1])
663	m21 = float(aList_s[2])
664	m22 = float(aList_s[3])
665
666	xs = m11x + m12y
667	ys = m21x + m22y
668	zs = z
669	aList_beads_x.append(xs)
670	aList_beads_y.append(ys)
671	aList_beads_z.append(zs)
672
673	j = j + 1
674
675	i = i + num_symm
676
677	return;
678
679	# Read experimentalal P(r) from GNOM output file
680
681	def read_pr(aList_r,aList_pr,aList_pr_sd,aList_pr_model,\
682	aList_pr_model_test,aList_pr_model_test2,inFile):
683
684	angstrom_scale = 1.0
685	Bins,Dbins,BinMag = data['Pair']['Distribution']
686
687	aList_r += list(Bins)
688	aList_pr += list(BinMag)
689	aList_pr_sd += list(np.ones_like(Bins)/100.)
690	aList_pr_model += list(np.zeros_like(Bins))
691	aList_pr_model_test += list(np.zeros_like(Bins))
692	aList_pr_model_test2 += list(np.zeros_like(Bins))
693
694	# file = open(inFile,'r')
695	# allLines = file.readlines()
696	# file.close()
697	#
698	# parse_data = 'no'
699	# for eachLine in allLines:
700	#
701	# if parse_data == 'yes':
702	#
703	# aList = eachLine.split()
704	# num_params = len(aList)
705	#
706	# if num_params == 3:
707	# r = float(aList[0])
708	# pr = float(aList[1])
709	# pr_sd = float(aList[2])
710	#
711	# aList_pr.append(pr)
712	# aList_pr_sd.append(pr_sd)
713	# aList_r.append(r)
714	# aList_pr_model.append(0.0)
715	# aList_pr_model_test.append(0.0)
716	# aList_pr_model_test2.append(0.0)
717	#
718	# if eachLine.find('R P(R) ERROR') > -1:
719	# parse_data = 'yes'
720	#
721	num_hist = len(aList_r)
722	hist_grid = float(aList_r[1]) - float(aList_r[0])
723
724
725	# Heuristic for checking units
726	# test_r = max(aList_r)
727	# if test_r < 30.0:
728	#
729	# aString = 'P(r)appears to be in nm. Converting to Angstrom units'
730	# print (aString)
731	# angstrom_scale = 10.0
732	#
733	# i = 0
734	# while i < num_hist:
735	# aList_r[i] = angstrom_scale * aList_r[i]
736	# i = i + 1
737	#
738	# hist_grid = angstrom_scale * hist_grid
739	#
740	# i = 0
741	# while i < num_hist:
742	# r = float(aList_r[i])
743	# r_calc = float(i)*hist_grid
744	#
745	# if abs(r - r_calc) > 0.15:
746	# aString = 'Input P(r) grid is irregular! Exiting'
747	# print (aString)
748	# time.sleep(4)
749	# sys.exit(1)
750	#
751	# i = i + 1
752	#
753	dmax = aList_r[num_hist-1]
754
755	# Pad histogram by 5 Angstrom
756
757	ipad = int(5.0/hist_grid)
758
759	i = 0
760	while i < ipad:
761	r = dmax + float(i)*hist_grid
762	aList_pr.append(0.0)
763	aList_pr_sd.append(0.0)
764	aList_r.append(r)
765	aList_pr_model.append(0.0)
766	aList_pr_model_test.append(0.0)
767	aList_pr_model_test2.append(0.0)
768	i = i + 1
769
770	return (dmax,hist_grid,num_hist,angstrom_scale);
771
772	# Scale P(r) onto model P(r) assuming same grid
773
774	def scale_pr(aList_pr,aList_pr_sd,aList_pr_model):
775
776	num_hist = len(aList_pr)
777	total_dr = 0.0
778	total_pr = 0.0
779
780	i = 0
781	while i < num_hist:
782	total_dr = total_dr + float(aList_pr_model[i])
783	total_pr = total_pr + float(aList_pr[i])
784	i = i + 1
785
786	scale = total_dr/total_pr
787
788	i = 0
789	while i < num_hist:
790	aList_pr[i] = scale*float(aList_pr[i])
791	aList_pr_sd[i] = scale * float(aList_pr_sd[i])
792	i = i + 1
793
794	return;
795
796	# Return a non-zero distance between two coordinates
797
798	def get_dr(x1,y1,z1,x2,y2,z2):
799
800	x1 = float(x1)
801	y1 = float(y1)
802	z1 = float(z1)
803	x2 = float(x2)
804	y2 = float(y2)
805	z2 = float(z2)
806	dr = (x1 - x2)2 + (y1-y2)2 + (z1-z2)**2
807	dr = max(dr,0.1)
808	dr = math.sqrt(dr)
809
810	return dr;
811
812	# Return non-zero distances one coordinate & the rest
813
814	def get_drs(xyz,XYZ):
815
816	return np.sqrt(np.sum((XYZ-xyz)**2,axis=1))
817
818	# Find center of beads within a radii
819
820	def center_beads(x,y,z,aList_beads_x,aList_beads_y,aList_beads_z,radii_1,radii_2):
821
822	num_beads = len(aList_beads_x)
823
824	# xsum = 0.0
825	# ysum = 0.0
826	# zsum = 0.0
827	# count_beads = 0.0
828	#
829	# i = 0
830	# while i < num_beads:
831	#
832	# dr = get_dr(x,y,z,aList_beads_x[i],aList_beads_y[i],aList_beads_z[i])
833	#
834	# if dr > radii_1 and dr < radii_2:
835	# count_beads = count_beads + 1.0
836	# xsum = xsum + float(aList_beads_x[i])
837	# ysum = ysum + float(aList_beads_y[i])
838	# zsum = zsum + float(aList_beads_z[i])
839	#
840	# i = i + 1
841	#
842	# if count_beads > 0.0:
843	# xsum = xsum/count_beads
844	# ysum = ysum/count_beads
845	# zsum = zsum/count_beads
846	# delta = (xsum - x)2 + (ysum - y)2 + (zsum - z)**2
847	# delta = math.sqrt(delta)
848	# else:
849	# delta = 0.0
850
851	XYZ = np.array([aList_beads_x,aList_beads_y,aList_beads_z]).T
852	xyz = np.array([x,y,z])
853	drs = get_drs(xyz,XYZ)
854	sumXYZ = np.array([XYZ[i] for i in range(num_beads) if radii_1<drs[i]<radii_2])
855	count_beads = sumXYZ.shape[0]
856
857	delta = 0.0
858	if count_beads:
859	delta = np.sqrt(np.sum(((np.sum(sumXYZ,axis=0)/count_beads)-xyz)**2))
860
861	return delta;
862
863	# Obtain mean of total VDW energy
864
865	def get_total_energy(aList_beads_x,aList_beads_y,aList_beads_z,econ12,econ6,bead_sep3):
866
867	nbeads = len(aList_beads_x)
868	vdw_all = 0.0
869
870
871	i = 0
872	while i < nbeads:
873	xyz = np.array([aList_beads_x[i],aList_beads_y[i],aList_beads_z[i]])
874	XYZ = np.array([aList_beads_x[:i],aList_beads_y[:i],aList_beads_z[:i]]).T
875	drs = get_drs(xyz,XYZ)
876	vdws = vdw_energies(econ12,econ6,e_width,drs,bead_sep3)
877	vdw_all += np.sum(vdws)
878	i += 1
879
880	# i = 0
881	# while i < nbeads:
882	# j = 0
883	# while j < i:
884	# dr = get_dr(aList_beads_x[i],aList_beads_y[i],aList_beads_z[i],\
885	# aList_beads_x[j],aList_beads_y[j],aList_beads_z[j])
886	# vdw = vdw_energy(econ12,econ6,e_width,dr,bead_sep3)
887	# vdw_all = vdw_all + vdw
888	# j = j + 1
889	# i = i + 1
890
891	vdw_all = vdw_all/float(nbeads)
892
893	return vdw_all;
894
895	# Energy minimize
896
897	def e_min(aList_beads_x,aList_beads_y,aList_beads_z,bead_sep,bead_sep3,aList_symm):
898
899	eps = bead_sep/(2.0**(1.0/6.0))
900	eps12 = eps**12
901	eps6 = eps**6
902	step_max = bead_sep
903	scale = 0.0
904	icount = -1
905
906	nbeads = len(aList_beads_x)
907	num_ops = len(aList_symm)
908	num_cycles = 51
909
910	i = 0
911	while i < num_cycles:
912
913	icount = icount + 1
914
915	aList_beads_x_new = []
916	aList_beads_y_new = []
917	aList_beads_z_new = []
918
919	sum_forces_scale = 0.0
920
921	k = 0
922	while k < nbeads - num_ops:
923
924	xold = float(aList_beads_x[k])
925	yold = float(aList_beads_y[k])
926	zold = float(aList_beads_z[k])
927
928
929	# fxyz = np.zeros(3)
930	# XYZ = np.array([aList_beads_x,aList_beads_y,aList_beads_z]).T
931	# xyz = np.array([xold,yold,zold])
932	# drs = get_drs(xyz,XYZ)
933	# drs = np.where(drs>eps,drs,eps)
934	# drs_sq = drs*drs
935	# drs12 = drs_sq**6
936	# drs6 = drs_sq**3
937	# dxs = xyz-XYZ
938	# forces = np.where(drs<bead_sep3,(1.0/drs_sq)(eps12/drs12 - 0.5eps6/drs6),0.0)
939	# fxyz = np.sum(forces[:,nxs]*dxs,axis=0)
940	# sum_forces_scale = np.sum(np.abs(fxyz))
941	#
942	# xstep = scale*fxyz[0]
943	# ystep = scale*fxyz[1]
944	# zstep = scale*fxyz[2]
945
946
947
948	fx = 0.0
949	fy = 0.0
950	fz = 0.0
951
952	j = 0
953	while j < nbeads:
954
955	xj = aList_beads_x[j]
956	yj = aList_beads_y[j]
957	zj = aList_beads_z[j]
958
959	dr = get_dr(xold,yold,zold,xj,yj,zj)
960
961	# Truncate very steep
962	dr = min(eps,dr)
963
964	if dr < bead_sep3:
965	dr_sq = dr*dr
966	dr12 = dr_sq**6
967	dr6 = dr_sq**3
968
969	dx = xold - xj
970	dy = yold - yj
971	dz = zold - zj
972
973	force = (1.0/dr_sq)(eps12/dr12 - 0.5eps6/dr6)
974	fx = fx + force*dx
975	fy = fy + force*dy
976	fz = fz + force*dz
977
978	sum_forces_scale = sum_forces_scale + abs(fx) + abs(fy) + abs(fz)
979
980	j = j + 1
981
982	#
983	xstep = scale*fx
984	ystep = scale*fy
985	zstep = scale*fz
986
987	if xstep > 0.0:
988	xstep = min(xstep,step_max)
989	else:
990	xstep = max(xstep,-step_max)
991
992	if ystep > 0.0:
993	ystep = min(ystep,step_max)
994	else:
995	ystep = max(ystep,-step_max)
996
997	if zstep > 0.0:
998	zstep = min(zstep,step_max)
999	else:
1000	zstep = max(zstep,-step_max)
1001
1002	xtest = xold + xstep
1003	ytest = yold + ystep
1004	ztest = zold + zstep
1005	aList_beads_x_new.append(xtest)
1006	aList_beads_y_new.append(ytest)
1007	aList_beads_z_new.append(ztest)
1008
1009	# Apply shifs to symm positions
1010	l = 0
1011	while l < num_ops:
1012	aList_s = aList_symm[l]
1013	m11 = float(aList_s[0])
1014	m12 = float(aList_s[1])
1015	m21 = float(aList_s[2])
1016	m22 = float(aList_s[3])
1017
1018	xs = m11xtest + m12ytest
1019	ys = m21xtest + m22ytest
1020	zs = ztest
1021
1022	aList_beads_x_new.append(xs)
1023	aList_beads_y_new.append(ys)
1024	aList_beads_z_new.append(zs)
1025
1026	l = l + 1
1027
1028	#
1029
1030	k = k + num_ops + 1
1031
1032	# Apply shifted positions after first cycle
1033	if i > 0:
1034
1035	m = 0
1036	while m < nbeads:
1037	aList_beads_x[m] = aList_beads_x_new[m]
1038	aList_beads_y[m] = aList_beads_y_new[m]
1039	aList_beads_z[m] = aList_beads_z_new[m]
1040	m = m + 1
1041
1042	#
1043
1044	mean_force = (num_ops+1)sum_forces_scale/(nbeads3.0)
1045	scale = bead_sep/mean_force
1046
1047	vdw_all = get_total_energy(aList_beads_x,aList_beads_y,aList_beads_z,econ12,econ6,bead_sep3)
1048
1049	if icount == 0:
1050	aString = 'Emin cycle: ' + str(i) + ' Energy: ' + str('%4.2f'%(vdw_all))
1051	print (aString)
1052	icount = -10
1053
1054	if vdw_all < 0.0:
1055	i = num_cycles
1056
1057	i = i + 1
1058
1059	return;
1060
1061	# Set up symmetry operators for rotational symmetry
1062
1063	def make_symm(aList_symm,num_symm):
1064
1065	angle_step = 360.0/float(num_symm)
1066
1067	i = 1
1068	while i < num_symm:
1069	theta = float(i) * angle_step
1070	theta = math.radians(theta)
1071	cos_theta = math.cos(theta)
1072	sin_theta = math.sin(theta)
1073	aList_s = [cos_theta,sin_theta,-sin_theta,cos_theta]
1074	aList_symm.append(aList_s)
1075	i = i + 1
1076
1077	return aList_symm;
1078
1079	# Set up a shift vector in P(r) for a change in bead position
1080
1081	def pr_shift_atom(aList_pr_model_test2,x1,y1,z1,\
1082	aList_beads_x,aList_beads_y,aList_beads_z,hist_grid,ii):
1083
1084	num_hist = len(aList_r)
1085	max_dr = (float(num_hist)-1.0)*hist_grid
1086	# num_beads = len(aList_beads_x)
1087
1088	# aList_pr_model_test2 = num_hist*[0.0,]
1089	#
1090	# i = 0
1091	# while i < num_hist:
1092	# aList_pr_model_test2[i] = 0.0
1093	# i = i + 1
1094
1095	XYZ = np.array([aList_beads_x,aList_beads_y,aList_beads_z]).T
1096	xyz = np.array([x1,y1,z1])
1097	drs = get_drs(xyz,XYZ)
1098	drs_grid = np.where(drs>max_dr,max_dr,drs)/hist_grid
1099	int_drs_grid = np.array(drs_grid,dtype=np.int)
1100	int_drs_grid = np.where(int_drs_grid>num_hist-2,num_hist-2,int_drs_grid)
1101	ip_lows = int_drs_grid
1102	ip_highs = ip_lows + 1
1103	ip_high_fracs = drs_grid - int_drs_grid
1104	ip_low_fracs = 1.0 - ip_high_fracs
1105	for ip_low in ip_lows:
1106	aList_pr_model_test2[ip_low] += ip_low_fracs[ip_low]
1107	for ip_high in ip_highs:
1108	aList_pr_model_test2[ip_high] += ip_high_fracs[ip_high]
1109
1110
1111	# i = 0
1112	# while i < num_beads:
1113	#
1114	# if i != ii:
1115	# x2 = float(aList_beads_x[i])
1116	# y2 = float(aList_beads_y[i])
1117	# z2 = float(aList_beads_z[i])
1118	# dr = get_dr(x1,y1,z1,x2,y2,z2)
1119	# dr = min(dr,max_dr)
1120	# dr_grid = dr/hist_grid
1121	# int_dr_grid = int(dr_grid)
1122	# int_dr_grid = max(int_dr_grid,0)
1123	# int_dr_grid = min(int_dr_grid,num_hist-2)
1124	# ip_low = int_dr_grid
1125	# ip_high = ip_low + 1
1126	# ip_high_frac = dr_grid - float(int_dr_grid)
1127	# ip_low_frac = 1.0 - ip_high_frac
1128	# aList_pr_model_test2[ip_low] = float(aList_pr_model_test2[ip_low]) + ip_low_frac
1129	# aList_pr_model_test2[ip_high] = float(aList_pr_model_test2[ip_high]) + ip_high_frac
1130	#
1131	# i = i + 1
1132	#
1133	return;
1134
1135	# Recenter set of beads to origin
1136
1137	def recenter_pdb(aList_beads_x,aList_beads_y,aList_beads_z):
1138
1139	nbeads = len(aList_beads_x)
1140	xsum = 0.0
1141	ysum = 0.0
1142	zsum = 0.0
1143
1144	i = 0
1145	while i < nbeads:
1146	xsum = xsum + float(aList_beads_x[i])
1147	ysum = ysum + float(aList_beads_y[i])
1148	zsum = zsum + float(aList_beads_z[i])
1149	i = i + 1
1150
1151	xmean = xsum/float(nbeads)
1152	ymean = ysum/float(nbeads)
1153	zmean = zsum/float(nbeads)
1154
1155	i = 0
1156	while i < nbeads:
1157	aList_beads_x[i] = float(aList_beads_x[i]) - xmean
1158	aList_beads_y[i] = float(aList_beads_y[i]) - ymean
1159	aList_beads_z[i] = float(aList_beads_z[i]) - zmean
1160	i = i + 1
1161
1162	return;
1163
1164	#############
1165	# EXECUTION #
1166	#############
1167
1168	#profiling start
1169	pr = cProfile.Profile()
1170	pr.enable()
1171	time0 = time.time()
1172
1173	version_aString = 'Program: SHAPES version 1.3'
1174
1175	print (version_aString)
1176	aString = 'Author: John Badger'
1177	print (aString)
1178	aString = 'Copyright: 2019, John Badger'
1179	print (aString)
1180	aString = 'License: GNU GPLv3'
1181	print (aString)
1182
1183	localtime = time.asctime( time.localtime(time.time()) )
1184	aString = 'Starting time: ' + str(localtime) + '\n'
1185	print (aString)
1186
1187	# aList_summary = []
1188	# aList_summary.append(version_aString)
1189	# aList_summary.append(str(localtime))
1190
1191	######################
1192	# Start up parmeters #
1193	######################
1194	# data['Shapes'] = {'outName':'','NumAA':100,'Niter':1,'AAscale':1.0,'Symm':1,'bias-z':0.0,
1195	# 'inflateV':1.0,'AAglue':0.0}
1196
1197	nbeads = 0
1198	num_sols = 1
1199	num_aa = 1.0
1200	num_symm = 1
1201	bias_z = 0.0
1202	inflate = 1.0
1203	prefix = ''
1204	surface_scale = 0.0
1205	starting_pdb = 'no'
1206	inFile = 'none'
1207	pdbfile_in = 'none'
1208	shapeDict = data['Shapes']
1209	prefix = shapeDict['outName']
1210	nbeads = shapeDict['NumAA']
1211	num_sols = shapeDict['Niter']
1212	num_aa = shapeDict['AAscale']
1213	num_symm = shapeDict['Symm']
1214	bias_z = shapeDict['bias-z']
1215	inflate = shapeDict['inflateV']
1216	surface_scale = shapeDict['AAglue']
1217	pdbOut = shapeDict['pdbOut']
1218	box_step = shapeDict.get('boxStep',4.0)
1219	Phases = []
1220	Patterns = []
1221	PRcalc = []
1222
1223	# # Parse
1224	#
1225	# if os.path.exists('shapes_ip.txt'):
1226	# file = open('shapes_ip.txt','r')
1227	# allLines = file.readlines()
1228	# file.close()
1229	# else:
1230	# aString = 'The local parameter file shapes_ip.txt was not found ! Exiting'
1231	# print (aString)
1232	# time.sleep(4)
1233	# sys.exit(1)
1234	#
1235	# for eachLine in allLines:
1236	#
1237	# aList = eachLine.split()
1238	#
1239	# num_params = len(aList)
1240	# if num_params > 0:
1241	#
1242	# if aList[0] == 'num_amino_acids':
1243	# nbeads = int(aList[1])
1244	## elif aList[0] == 'input_pr':
1245	## inFile = str(aList[1])
1246	# elif aList[0] == 'num_solns':
1247	# num_sols = int(aList[1])
1248	# elif aList[0] == 'num_aa_scale':
1249	# num_aa = float(aList[1])
1250	# elif aList[0] == 'symm':
1251	# num_symm = int(aList[1])
1252	# elif aList[0] == 'bias_z':
1253	# bias_z = float(aList[1])
1254	# elif aList[0] == 'inflate_vol':
1255	# inflate = float(aList[1])
1256	# elif aList[0] == 'pdb_start':
1257	# pdbfile_in = str(aList[1])
1258	# elif aList[0] == 'id':
1259	# prefix = str(aList[1]) + '_'
1260	# elif aList[0] == 'glue':
1261	# surface_scale = float(aList[1])
1262
1263
1264	# Check inputs
1265
1266	if num_sols > 0:
1267	aString = 'Number of runs: ' + str(num_sols)
1268	print (aString)
1269	else:
1270	aString = 'Zero reconstruction runs specified! Exiting'
1271	print (aString)
1272	time.sleep(4)
1273	sys.exit(1)
1274
1275	#
1276	if nbeads == 0:
1277	if os.path.exists(pdbfile_in):
1278	aString = 'Will use CA atoms from ' + str(pdbfile_in) + ' as the initial bead distribution.'
1279	print (aString)
1280	starting_pdb = 'yes'
1281	else:
1282	aString = 'Zero amino acid count specified and no starting file found. Exiting'
1283	print (aString)
1284	time.sleep(4)
1285	sys.exit(1)
1286	else:
1287	aString = 'Number of amino acids: ' + str(nbeads)
1288	print (aString)
1289
1290	#
1291	# if os.path.exists(inFile):
1292	# aString = 'Input P(r) file name: ' + str(inFile)
1293	# print (aString)
1294	# else:
1295	# aString = 'P(r) input file not found. Exiting'
1296	# print (aString)
1297	# time.sleep(4)
1298	# sys.exit(1)
1299
1300	#
1301	if num_aa == 0.0:
1302	aString = 'Scale for amino acid count to particle number cannot be zero! Exiting'
1303	print (aString)
1304	time.sleep(4)
1305	sys.exit(1)
1306	else:
1307	aString = 'Scale aa to bead count: ' + str(num_aa)
1308	print (aString)
1309
1310	#
1311	if num_symm == 0:
1312	aString = 'Rotational symmetry cannot be zero! Set to 1 for no symmetry. Exiting'
1313	print (aString)
1314	time.sleep(4)
1315	sys.exit(1)
1316	else:
1317	aString = 'Point symmetry: ' + str(num_symm)
1318	print (aString)
1319
1320	#
1321	if bias_z > 0.2:
1322	aString = 'Max bias on Z axis for initial particle distribution is 0.2 (rods). Reset to 0.2.'
1323	print (aString)
1324	bias_z = 0.2
1325	elif bias_z < -0.2:
1326	aString = 'Min bias on Z axis for initial particle distribution is -0.2 (disks). Reset to -0.2.'
1327	print (aString)
1328	bias_z = -0.2
1329	else:
1330	aString = 'Z-axis bias: ' + str(bias_z)
1331	print (aString)
1332
1333	#
1334	if inflate < 0.0:
1335	aString = 'Inflation of PSV cannot be less than zero! Exiting'
1336	print (aString)
1337	time.sleep(4)
1338	sys.exit(1)
1339	elif inflate > 2.0:
1340	aString = 'Inflation of PSV cannt be greater than 2.0! Exiting'
1341	print (aString)
1342	time.sleep(4)
1343	sys.exit(1)
1344	else:
1345	aString = 'PSV inflation factor: ' + str(inflate)
1346	print (aString)
1347
1348	#
1349	if surface_scale > 0.0:
1350	aString = 'Cavity weight: ' + str(surface_scale)
1351	print (aString)
1352
1353	########## UNIVERSAL CONSTANTS ######################
1354
1355	# No of macrocycles (gives extra cycles at constant volume after num_contract)
1356	niter = 160
1357
1358	# No of contraction cycles
1359	num_contract = 140
1360
1361	# Number of cycles at each fixed volume
1362	num_micro_cyc = 10
1363
1364	# Final quench
1365	num_sa_max = niter - num_micro_cyc
1366
1367	# Initial scale for P(r) shifts versus E shifts
1368	hscale = 3000.0
1369
1370	# Standard deviation for annealing acceptance (cf well-depth of -1 unit for two beads)
1371	sd_mc = float(num_symm) * 2.0
1372
1373	# Fiddle factor for keeping the accessible, molecular volume larger than PSV
1374	scale_vol = 1.15
1375
1376	# Standard amino acid volume MW = 110.0 x 1.21 i.e. mean mw x mw-to-vol-scale
1377	vol_bead = 133.1
1378
1379	# Bead separation for best packing ~5.6 (I think)
1380	#- 75% better than rectangular grid 5.1 for this amino acid vol
1381	bead_sep = 5.6
1382
1383	# Usually num_aa is unity. Adjust parameters otherwise
1384	if num_aa != 1 and nbeads != 0:
1385	nbeads = int(nbeads*num_aa)
1386	vol_bead = vol_bead / num_aa
1387	bead_sep = (vol_bead * 4/3)**(1.0/3.0)
1388
1389	# Increase bead separation for inflated volumes
1390	bead_sep = bead_sep * inflate**(1.0/3.0)
1391
1392	# Partial specific volumes at start and end
1393
1394	if starting_pdb == 'yes':
1395	nmols_vol_start = 1.1 * inflate
1396	else:
1397	nmols_vol_start = 2.0 * inflate
1398
1399	nmols_vol_end = 1.0 * inflate
1400	nmols_vol_subtract = nmols_vol_start - nmols_vol_end
1401
1402	# Box parametere
1403	# box_step = 4.0 #5.0
1404	box_pt_vol = box_stepbox_stepbox_step
1405
1406	# Energy parameters - flat bottomed VDW (2.0A for a 5.6A bead separation)
1407
1408	well_width = 0.36*bead_sep
1409	econ12 = bead_sep**12
1410	econ6 = bead_sep**6
1411	r_width = bead_sep + well_width
1412	r_width6 = r_width**6
1413	r_width12 = r_width6**2
1414	e_width = econ12/r_width12 - 2.0*econ6/r_width6
1415	bead_sep3 = 3.0*bead_sep
1416	abs_e_width = abs(e_width)
1417
1418	# Range for box identification (might need to increase for poor data)
1419	rsearch = (bead_sep + 0.5well_width)1.5
1420
1421	# Search range for optional cavity inhibition energy term
1422	radii_1 = 1.5*bead_sep
1423	radii_2 = 4.0*bead_sep
1424
1425	# Setup symmetry operators
1426
1427	aList_symm = []
1428	aList_symm = make_symm(aList_symm,num_symm)
1429	num_ops = len(aList_symm)
1430
1431	# Read experimental histogram
1432
1433	aList_r = []
1434	aList_pr = []
1435	aList_pr_sd = []
1436	aList_pr_model = []
1437	aList_pr_model_test = []
1438	aList_pr_model_test2 = []
1439
1440	(dmax,hist_grid,num_hist_in,angstrom_scale) = read_pr(aList_r,aList_pr,aList_pr_sd,\
1441	aList_pr_model,aList_pr_model_test,\
1442	aList_pr_model_test2,inFile)
1443
1444	# dmax_over2 = dmax/2.0
1445	num_hist = len(aList_r)
1446
1447	aString = 'Number of points read from P(r): ' + str(num_hist_in)
1448	print (aString)
1449	aString = 'Grid sampling: ' + str(hist_grid) + ' Dmax: ' + str(dmax)
1450	print (aString)
1451
1452	# Skip over initial points in scoring
1453
1454	skip = r_width/float(num_hist)
1455	skip = int(skip) + 2
1456
1457	# Read intensity data that was used for P(r)
1458
1459	aList_q = []
1460	aList_i = []
1461	aList_i_sd = []
1462
1463	read_i(aList_q,aList_i,aList_i_sd,inFile,angstrom_scale)
1464
1465	num_intensities = len(aList_q)
1466	aString = 'Number of intensity data points read: ' + str(num_intensities)
1467	print (aString)
1468
1469	#########################
1470	# CYCLE OVER SOLUTIONS #
1471	#########################
1472
1473	i_soln = 0
1474	while i_soln < num_sols:
1475
1476	file_no = str(i_soln + 1)
1477
1478	aString = '\nReconstruction trial: ' + str(file_no)
1479	print (aString)
1480
1481	aString = 'Trial:' + file_no
1482	# aList_summary.append(aString)
1483
1484	file_beads = prefix + 'beads_' + file_no + '.pdb'
1485	# file_pr = prefix + 'pr_calc_' + file_no + '.dat'
1486	file_psv = prefix + 'psv_shape_' + file_no + '.pdb'
1487	# file_intensity = prefix + 'intensity_' + file_no + '.dat'
1488
1489	# Setup initial bead distribution
1490
1491	aList_beads_x = []
1492	aList_beads_y = []
1493	aList_beads_z = []
1494
1495	# Re-initialize standard deviation for annealing acceptance
1496	sd_mc = float(num_symm) * 2.0
1497
1498	# Set random bead distribution
1499
1500	if starting_pdb == 'yes':
1501	read_pdb(aList_beads_x,aList_beads_y,aList_beads_z,pdbfile_in)
1502	nbeads = len(aList_beads_x)
1503	num_symm = 1
1504	aString = 'Number of CA sites read: ' + str(nbeads)
1505	print (aString)
1506	aString = 'Symmetry set to 1 (required)'
1507	print (aString)
1508	aString = 'Input center was shifted to the origin'
1509	print (aString)
1510	else:
1511	random_beads(aList_beads_x,aList_beads_y,aList_beads_z,nbeads,dmax,aList_symm,bias_z)
1512	nbeads = len(aList_beads_x)
1513	aString = 'Number of beads randomly placed: ' + str(nbeads)
1514	print (aString)
1515
1516	# Protein partial specific volume
1517	psv_vol = float(nbeads)*vol_bead
1518
1519	# Histogram of inter-bead distance
1520
1521	calc_pr(aList_beads_x,aList_beads_y,aList_beads_z,aList_pr_model,hist_grid)
1522
1523	# Scale experimental P(r) and model histogram
1524
1525	scale_pr(aList_pr,aList_pr_sd,aList_pr_model)
1526
1527	# Minimize initial energy using expanded VDW
1528
1529	if starting_pdb != 'yes':
1530	aString = 'Minimize energy of initial positions'
1531	print (aString)
1532	bead_sep_e = 1.35*bead_sep
1533	bead_sep3_e = 3.0*bead_sep_e
1534	e_min(aList_beads_x,aList_beads_y,aList_beads_z,bead_sep_e,bead_sep3_e,aList_symm)
1535	else:
1536	aString = 'Skipping energy minimization of initial positions'
1537	print (aString)
1538
1539	# Get the initial score between observed and calculated P(r)
1540
1541	hist_score_best = pr_dif(aList_pr,aList_pr_model,skip)
1542	aString = 'Initial rms P(r): ' + str('%4.3f'%(hist_score_best))
1543	print (aString)
1544
1545	aList_i_calc = []
1546	ft_to_intensity(aList_q,aList_i_calc,aList_r,aList_pr_model,nbeads)
1547	(chi_sq,rvalue) = score_Ic(aList_q,aList_i,aList_i_sd,aList_i_calc)
1548	aString = 'Initial Rvalue: ' + str('%4.3f'%(rvalue)) + ' CHI-squared: ' + str('%4.3f'%(chi_sq))
1549	print (aString)
1550
1551	###########################
1552	# Iterate #
1553	###########################
1554
1555	num_boxes = 0
1556	count_boxing = 0
1557	fraction_psv = 0
1558	success_rate_all = 0.0
1559	box_iter = num_micro_cyc - 1
1560
1561	sum_delta_pack = 0.0
1562
1563	count_it = 0
1564	while count_it < niter:
1565
1566	success = 0
1567	count_hist_yes = 0
1568	sum_e = 0.0
1569	sum_h = 0.0
1570
1571	# Find populated volumes and fix solution every 10 macrocycles
1572
1573	box_iter = box_iter + 1
1574
1575	if box_iter == num_micro_cyc:
1576
1577	box_iter = 0
1578	count_boxing = count_boxing + 1
1579
1580	if count_it < num_contract - 1:
1581	scale = float(count_it)/float(num_contract)
1582	else:
1583	scale = 1.0
1584
1585	# Establish confinement volume using a local average
1586
1587	aList_box_x_all = []
1588	aList_box_y_all = []
1589	aList_box_z_all = []
1590	aList_box_score = []
1591
1592	recenter_pdb(aList_beads_x,aList_beads_y,aList_beads_z)
1593
1594	# Adaptive masking was not helpful
1595	# rsearch_use = (2.0 - scale)*rsearch
1596
1597	set_box_fast(aList_beads_x,aList_beads_y,aList_beads_z,\
1598	aList_box_x_all,aList_box_y_all,aList_box_z_all,\
1599	aList_box_score,box_step,dmax,rsearch)
1600
1601	aList_box_x = []
1602	aList_box_y = []
1603	aList_box_z = []
1604
1605	psv_ratio = nmols_vol_start - scale*nmols_vol_subtract
1606	vol_target = scale_vol(psv_ratiopsv_vol)
1607
1608	set_vol(aList_box_x_all,aList_box_y_all,aList_box_z_all,\
1609	aList_box_score,aList_box_x,aList_box_y,aList_box_z,\
1610	vol_target,box_pt_vol)
1611
1612	num_boxes = len(aList_box_x)
1613	fraction_psv = float(num_boxes)*box_pt_vol/psv_vol
1614
1615	# Find beads that are ouside the allowed volume
1616
1617	aList_contacts = []
1618	disallowed_beads(aList_beads_x,aList_beads_y,aList_beads_z,aList_contacts,\
1619	aList_box_x,aList_box_y,aList_box_z,rsearch)
1620	num_outof_box = len(aList_contacts)
1621
1622	aString = 'Target volume: ' + str('%4.2f'%(scale_vol*psv_ratio)) + ' Actual volume: ' + \
1623	str('%4.2f'%(fraction_psv)) + ' Beads outside volume: ' + str(num_outof_box)
1624	print (aString)
1625
1626	# Recalculate P(r) and rescore for reliability
1627
1628	calc_pr(aList_beads_x,aList_beads_y,aList_beads_z,aList_pr_model,hist_grid)
1629	hist_score_best = pr_dif(aList_pr,aList_pr_model,skip)
1630
1631	# aList_i_calc = []
1632	# ft_to_intensity(aList_q,aList_i_calc,aList_r,aList_pr_model,nbeads)
1633	# (chi_sq,rvalue) = score_Ic(aList_q,aList_i,aList_i_sd,aList_i_calc)
1634
1635	# Reset SA deviation if mean success rate over last trials is under 0.1
1636
1637	mean_success_rate = float(success_rate_all)/float(num_micro_cyc)
1638
1639	if count_it < num_contract and count_boxing != 1:
1640
1641	if mean_success_rate < 0.1:
1642	sd_mc = 1.3*sd_mc
1643	aString = 'Raising allowed energy deviation to ' + str('%4.2f'%(sd_mc))
1644	print (aString)
1645
1646	if mean_success_rate > 0.2:
1647	sd_mc = 0.7*sd_mc
1648	aString = 'Reducing allowed energy deviation to ' + str('%4.2f'%(sd_mc))
1649	print (aString)
1650
1651	success_rate_all = 0.0
1652
1653	# Make one macrocycle that is a run over the nbeads
1654
1655	ii = 0
1656	while ii < nbeads:
1657
1658	# Initialize
1659
1660	energy_old = 0.0
1661	energy_new = 0.0
1662
1663	i = 0
1664	while i < num_hist:
1665	aList_pr_model_test[i] = 0.0
1666	i = i + 1
1667
1668	# Select a target bead and make trial shift
1669
1670	xold = float(aList_beads_x[ii])
1671	yold = float(aList_beads_y[ii])
1672	zold = float(aList_beads_z[ii])
1673
1674	ibox = random.randint(0,num_boxes-1)
1675	xtest = float(aList_box_x[ibox]) + random.uniform(-rsearch,rsearch)
1676	ytest = float(aList_box_y[ibox]) + random.uniform(-rsearch,rsearch)
1677	ztest = float(aList_box_z[ibox]) + random.uniform(-rsearch,rsearch)
1678
1679	# Calculate and capture symmetry mates
1680
1681	aList_temp_save_x = []
1682	aList_temp_save_y = []
1683	aList_temp_save_z = []
1684	aList_symm_x = []
1685	aList_symm_y = []
1686	aList_symm_z = []
1687
1688	l = 0
1689	while l < num_ops:
1690	aList_s = aList_symm[l]
1691	m11 = float(aList_s[0])
1692	m12 = float(aList_s[1])
1693	m21 = float(aList_s[2])
1694	m22 = float(aList_s[3])
1695
1696	xs = m11xtest + m12ytest
1697	ys = m21xtest + m22ytest
1698	zs = ztest
1699
1700	aList_symm_x.append(xs)
1701	aList_symm_y.append(ys)
1702	aList_symm_z.append(zs)
1703
1704	ipt = ii + l + 1
1705	aList_temp_save_x.append(aList_beads_x[ipt])
1706	aList_temp_save_y.append(aList_beads_y[ipt])
1707	aList_temp_save_z.append(aList_beads_z[ipt])
1708
1709	l = l + 1
1710
1711	# Get initial VDW energy for interactions of this bead with all others
1712
1713	i = 0
1714	while i < nbeads:
1715	if i != ii:
1716	x = float(aList_beads_x[i])
1717	y = float(aList_beads_y[i])
1718	z = float(aList_beads_z[i])
1719	dr_old = get_dr(xold,yold,zold,x,y,z)
1720	vdw_old = vdw_energy(econ12,econ6,e_width,dr_old,bead_sep3)
1721	energy_old = energy_old + num_symm*vdw_old
1722	i = i + 1
1723
1724	# Get initial contributions to P(r)
1725
1726	aList_pr_model_test2 = num_hist*[0.0,]
1727	pr_shift_atom(aList_pr_model_test2,xold,yold,zold,aList_beads_x,\
1728	aList_beads_y,aList_beads_z,hist_grid,ii)
1729
1730	i = 0
1731	while i < num_hist:
1732	aList_pr_model_test[i] = aList_pr_model_test2[i]
1733	i = i + 1
1734
1735	# Get VDW energy for interactions of the shifted bead with all others
1736
1737	l = 0
1738	while l < num_ops:
1739	ipt = ii + l + 1
1740	aList_beads_x[ipt] = aList_symm_x[l]
1741	aList_beads_y[ipt] = aList_symm_y[l]
1742	aList_beads_z[ipt] = aList_symm_z[l]
1743	l = l + 1
1744
1745	i = 0
1746	while i < nbeads:
1747	if i != ii:
1748	x = float(aList_beads_x[i])
1749	y = float(aList_beads_y[i])
1750	z = float(aList_beads_z[i])
1751	dr_new = get_dr(xtest,ytest,ztest,x,y,z)
1752	vdw_new = vdw_energy(econ12,econ6,e_width,dr_new,bead_sep3)
1753	energy_new = energy_new + num_symm*vdw_new
1754	i = i + 1
1755
1756	# Get cavity energy difference
1757
1758	delta_old = center_beads(xold,yold,zold,aList_beads_x,\
1759	aList_beads_y,aList_beads_z,radii_1,radii_2)
1760	delta_new = center_beads(xtest,ytest,ztest,aList_beads_x,\
1761	aList_beads_y,aList_beads_z,radii_1,radii_2)
1762
1763	delta_pack = num_symmsurface_scale(delta_new - delta_old)/(radii_1 + radii_2)
1764	sum_delta_pack = sum_delta_pack + abs(delta_pack)
1765
1766	# Get shifted contributions to P(r)
1767
1768	aList_pr_model_test2 = num_hist*[0.0,]
1769	pr_shift_atom(aList_pr_model_test2,xtest,ytest,ztest,aList_beads_x,\
1770	aList_beads_y,aList_beads_z,hist_grid,ii)
1771
1772	# Get net shift in contribution to P(r)
1773
1774	i = 0
1775	while i < num_hist:
1776	aList_pr_model_test[i] = aList_pr_model_test2[i] - aList_pr_model_test[i]
1777	i = i + 1
1778
1779	# Get statistic for agreement with P(r) after accumulating shift vectors
1780
1781	i = 0
1782	while i < num_hist:
1783	aList_pr_model_test[i] = float(aList_pr_model[i]) + num_symm*float(aList_pr_model_test[i])
1784	i = i + 1
1785
1786	hist_score = pr_dif(aList_pr,aList_pr_model_test,skip)
1787
1788	# aList_i_calc = []
1789	# ft_to_intensity(aList_q,aList_i_calc,aList_r,aList_pr_model,nbeads)
1790	# (chi_sq,rvalue) = score_Ic(aList_q,aList_i,aList_i_sd,aList_i_calc)
1791
1792	# Scoring shifts
1793
1794	delta_h = hist_score - hist_score_best
1795	delta_e = energy_new - energy_old + delta_pack
1796
1797	# Recalibrate scale so impact of energy and P(r) is equal on plausible shifts
1798
1799	if delta_e < abs_e_width:
1800	sum_e = sum_e + abs(delta_e)
1801	sum_h = sum_h + abs(delta_h)
1802
1803	# Monitor potential moves based of P(r)
1804
1805	if delta_h < 0.0:
1806	count_hist_yes = count_hist_yes + 1.0
1807
1808	# Acceptance and update
1809
1810	score = delta_e + delta_h*hscale
1811
1812	if count_it < num_sa_max:
1813	barrier = abs(random.gauss(0.0,sd_mc))
1814	else:
1815	barrier = 0.0
1816
1817	if score < barrier:
1818
1819	success = success + 1.0
1820	hist_score_best = hist_score
1821
1822	# Update model but symmetry positions that were already put in
1823
1824	aList_beads_x[ii] = xtest
1825	aList_beads_y[ii] = ytest
1826	aList_beads_z[ii] = ztest
1827
1828	# Update P(r)
1829
1830	i = 0
1831	while i < num_hist:
1832	aList_pr_model[i] = aList_pr_model_test[i]
1833	i = i + 1
1834
1835	else:
1836
1837	# Revert to original model
1838
1839	aList_beads_x[ii] = xold
1840	aList_beads_y[ii] = yold
1841	aList_beads_z[ii] = zold
1842
1843	l = 0
1844	while l < num_ops:
1845	ipt = ii + l + 1
1846	aList_beads_x[ipt] = aList_temp_save_x[l]
1847	aList_beads_y[ipt] = aList_temp_save_y[l]
1848	aList_beads_z[ipt] = aList_temp_save_z[l]
1849	l = l + 1
1850	#
1851
1852	ii = ii + num_symm
1853
1854	# Get energy statistics at end of macrocycle
1855
1856	vdw_all = get_total_energy(aList_beads_x,aList_beads_y,aList_beads_z,econ12,econ6,bead_sep3)
1857
1858	# Rescale and convergence statistics
1859
1860	if sum_h > 0.0:
1861	hscale = sum_e/sum_h
1862
1863	count_hist_yes = count_hist_yes*float(num_symm)/float(nbeads)
1864	success_rate = success*float(num_symm)/float(nbeads)
1865	success_rate_all = success_rate_all + success_rate
1866
1867	if not (count_it+1)%10:
1868	aString = 'Cycle ' + str(count_it+1) + ' Moves ' + str('%.2f'%(success_rate)) + \
1869	' Possibles ' + str('%.2f'%(count_hist_yes)) + ' rms P(r) '+ str('%4.3f'%(hist_score)) + \
1870	' Energy ' + str('%4.2f'%(vdw_all))
1871	print (aString)
1872	# print('Rvalue: %4.3f CHI-squared: %4.3f'%(rvalue,chi_sq))
1873
1874	if dlg:
1875	dlg.Update(count_it+1,newmsg='Cycle no.: '+str(count_it)+' of 160')
1876
1877	# Debug statitics. Weight of 10 gives about 1.0
1878	#sum_delta_pack = sum_delta_pack/float(nbeads)
1879	#print (sum_delta_pack)
1880	#
1881
1882	count_it = count_it + 1
1883
1884	#####################
1885	# ANALYZE AND WRITE #
1886	#####################
1887
1888	aString = '\nFinal model statistics'
1889	print (aString)
1890
1891	calc_pr(aList_beads_x,aList_beads_y,aList_beads_z,aList_pr_model,hist_grid)
1892	hist_score_best = pr_dif(aList_pr,aList_pr_model,skip)
1893
1894	# P(r) fitting statistics
1895	delta_hist_sum = pr_rfactor(aList_pr,aList_pr_sd,aList_pr_model_test,skip)
1896
1897	aString = 'Delta P(r): ' + str('%4.3f'%(delta_hist_sum))
1898	print (aString)
1899
1900	# Get final energy
1901	vdw_all = get_total_energy(aList_beads_x,aList_beads_y,aList_beads_z,econ12,econ6,bead_sep3)
1902
1903	aString = 'VDW energy: ' + str('%4.2f'%(vdw_all))
1904	print (aString)
1905
1906	Phases.append([file_beads.split('.')[0],aList_beads_x,aList_beads_y,aList_beads_z])
1907	# Write out beads as pseudo a PDB file
1908	if pdbOut:
1909	pdb_writer(aList_beads_x,aList_beads_y,aList_beads_z,file_beads,aString)
1910
1911	# Calculate and write final PSV shape
1912
1913	aList_box_x_all = []
1914	aList_box_y_all = []
1915	aList_box_z_all = []
1916	aList_box_score = []
1917
1918	set_box_fast(aList_beads_x,aList_beads_y,aList_beads_z,\
1919	aList_box_x_all,aList_box_y_all,aList_box_z_all,\
1920	aList_box_score,box_step,dmax,rsearch)
1921
1922	aList_box_x = []
1923	aList_box_y = []
1924	aList_box_z = []
1925	psv_vol_use = psv_vol*inflate
1926
1927	set_vol(aList_box_x_all,aList_box_y_all,aList_box_z_all,\
1928	aList_box_score,aList_box_x,aList_box_y,aList_box_z,\
1929	psv_vol_use,box_pt_vol)
1930
1931	num_boxes = len(aList_box_x)
1932	fraction_psv = num_boxes*box_pt_vol/psv_vol
1933
1934	# Correct final volume if initial estimate is too small
1935
1936	fraction_psv_use = num_boxes*box_pt_vol/psv_vol_use
1937
1938	if fraction_psv_use < 1.0:
1939	aList_box_x = []
1940	aList_box_y = []
1941	aList_box_z = []
1942	vol_use = 1.05*psv_vol_use
1943	set_vol(aList_box_x_all,aList_box_y_all,aList_box_z_all,\
1944	aList_box_score,aList_box_x,aList_box_y,aList_box_z,\
1945	vol_use,box_pt_vol)
1946
1947	num_boxes = len(aList_box_x)
1948	fraction_psv = float(num_boxes)*box_pt_vol/psv_vol
1949
1950	aString = 'Final PSV of protein envelope: ' + str('%4.2f'%(fraction_psv))
1951	print (aString)
1952
1953	# Write input and model P(r)
1954	# pr_writer(aList_pr,aList_r,aList_pr_model,file_pr)
1955	PRcalc.append([aList_r,aList_pr,copy.copy(aList_pr_model),delta_hist_sum])
1956
1957	# Calculate comparison versus intensities
1958
1959	if num_intensities > 0:
1960
1961	# calculate intensity
1962	aList_i_calc = []
1963	# num_beads = len(aList_box_x)
1964	ft_to_intensity(aList_q,aList_i_calc,aList_r,aList_pr_model,nbeads)
1965
1966	# Scale and obtain statistics
1967	(chi_sq,rvalue) = score_Ic(aList_q,aList_i,aList_i_sd,aList_i_calc)
1968
1969	aString = 'Rvalue: ' + str('%4.3f'%(rvalue)) + ' CHI-squared: ' + str('%4.3f'%(chi_sq))
1970	print (aString)
1971
1972	# Write output intensity file
1973	Patterns.append([aList_q,aList_i,aList_i_calc,rvalue])
1974	# write_all_data(file_intensity,aList_q,aList_i,aList_i_calc,aString)
1975
1976	# aString = 'Output intensity file: ' + str(file_intensity)
1977	# print (aString)
1978
1979	else:
1980
1981	chi_sq = 'NA'
1982
1983	# Write final volume
1984
1985	delta_hist_sum = '%4.3f'%(delta_hist_sum)
1986	vdw_all = '%4.2f'%(vdw_all)
1987	fraction_psv = '%4.2f'%(fraction_psv)
1988	chi_sq = '%4.3f'%(chi_sq)
1989
1990	aString = 'REMARK P(r) dif:'+ str(delta_hist_sum) + ' E:'\
1991	+ str(vdw_all) + ' CHIsq:' + str(chi_sq) + \
1992	' PSV:' + str(fraction_psv)
1993
1994	Phases.append([file_psv.split('.')[0],aList_box_x,aList_box_y,aList_box_z])
1995	if pdbOut:
1996	pdb_writer(aList_box_x,aList_box_y,aList_box_z,file_psv,aString)
1997
1998	# Write Summary
1999
2000	aString = 'P(r) dif:' + str(delta_hist_sum) + ' E:' \
2001	+ str(vdw_all) + ' CHISQ:' + str(chi_sq) + ' PSV:' + str(fraction_psv)
2002	# aList_summary.append(aString)
2003
2004	i_soln = i_soln + 1
2005
2006	#########################################
2007	#### END OF LOOP OVER MULTI-SOLUTIONS ###
2008	#########################################
2009
2010	#end profiling
2011	pr.disable()
2012	s = StringIO.StringIO()
2013	sortby = 'tottime'
2014	ps = pstats.Stats(pr, stream=s).strip_dirs().sort_stats(sortby)
2015	print('Profiler of function calculation; top 50% of routines:')
2016	ps.print_stats(.5)
2017	print(s.getvalue())
2018	print('%s%.3f'%('Run time = ',time.time()-time0))
2019
2020	localtime = time.asctime( time.localtime(time.time()) )
2021
2022	aString = '\nCompletion time: ' + str(localtime)
2023	print (aString)
2024
2025	# aList_summary.append(str(localtime))
2026	#
2027	# # Create summary
2028	#
2029	# aFile_log = prefix + 'shapes_summary.log'
2030	# num_lines = len(aList_summary)
2031	#
2032	# file = open(aFile_log,'w')
2033	# i = 0
2034	# while i < num_lines:
2035	# aString = str(aList_summary[i])
2036	# file.write(aString)
2037	# file.write('\n')
2038	# i = i + 1
2039	# file.close()
2040
2041	return Phases,Patterns,PRcalc
2042
2043
2044

Note: See TracBrowser for help on using the repository browser.

Download in other formats: