source: trunk/G2shapes.py @ 4095

Last change on this file since 4095 was 4095, checked in by vondreele, 2 years ago

improve dmax estimation - use proper selected weights
ditto P(R) calculation
remove Seq fit from all but particle fits - not appropriate for rest.
Add option for bead separation in SHAPES results - improves appearance in small proteins

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