source:trunk/G2shapes.py@4098

Last change on this file since 4098 was 4098, checked in by vondreele, 3 years ago

speedup of SHAPES

File size: 62.5 KB
Line
1#
3#
4# This program is free software: you can redistribute it and/or modify
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
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
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')
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
209
210        xmean = 0.0
211        ymean = 0.0
212        zmean = 0.0
214
215        file = open(pdbfile_in,'r')
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
240
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
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
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
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
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     # Evaluate local bead densities and point density on a notional grid - fast version
382
384                aList_box_x_all,aList_box_y_all,aList_box_z_all,\
385                aList_box_score,box_step,dmax,rsearch):
386
387        dmax_over2 = dmax/2.0
388        num_box = int(dmax/box_step)
389        search_sq = rsearch**2
390
391        XYZ = np.meshgrid(np.linspace(-dmax_over2,dmax_over2,num_box),
392            np.linspace(-dmax_over2,dmax_over2,num_box),
393            np.linspace(-dmax_over2,dmax_over2,num_box))
394        XYZ = np.array([XYZ[0].flatten(),XYZ[1].flatten(),XYZ[2].flatten()]).T
396        for XYZi in XYZ:
397            dsq = np.sum((xyz-XYZi)**2,axis=1)
398            count = int(np.sum(np.array([1 for dist in dsq if dist < search_sq])))
399            if count>1:
400                aList_box_x_all.append(XYZi[0])
401                aList_box_y_all.append(XYZi[1])
402                aList_box_z_all.append(XYZi[2])
403                aList_box_score.append(count)
404        return;
405
406   # Establish a volume
407
408    def set_vol(aList_box_x_all,aList_box_y_all,aList_box_z_all,aList_box_score,\
409                aList_box_x,aList_box_y,aList_box_z,vol_target,box_pt_vol):
410
411        num_pts = len(aList_box_score)
412        num_tries = int(max(aList_box_score))
413        density_thresh = max(aList_box_score) - 1.0
414        vol = vol_target + 1.0
415
416        i = 0
417        while i < num_tries:
418
419            density_thresh = density_thresh - 1.0
420            num_box_pts = 0.0
421
422            j = 0
423            while j < num_pts:
424                density = float(aList_box_score[j])
425                if density >= density_thresh:
426                    num_box_pts = num_box_pts + 1.0
427                j = j + 1
428
429            vol = num_box_pts*box_pt_vol
430            if vol < vol_target:
431                density_use = density_thresh
432
433            i = i + 1
434
435        #
436
437        num_box_pts1 = 0.0
438        num_box_pts2 = 0.0
439        density_thresh1 = density_use
440        density_thresh2 = density_use - 1.0
441        i = 0
442        while i < num_pts:
443            density_value = float(aList_box_score[i])
444            if density_value >= density_thresh1:
445                num_box_pts1 = num_box_pts1 + 1.0
446            if density_value >= density_thresh2:
447                num_box_pts2 = num_box_pts2 + 1.0
448            i = i + 1
449
450        vol1 = num_box_pts1*box_pt_vol
451        vol2 = num_box_pts2*box_pt_vol
452        delta1 = abs(vol1-vol_target)
453        delta2 = abs(vol2-vol_target)
454
455        if delta1 < delta2:
456            density_thresh = density_thresh1
457        else:
458            density_thresh = density_thresh2
459
460        i = 0
461        while i < num_pts:
462
463            density_value = float(aList_box_score[i])
464            if density_value >= density_thresh:
465                aList_box_x.append(aList_box_x_all[i])
466                aList_box_y.append(aList_box_y_all[i])
467                aList_box_z.append(aList_box_z_all[i])
468
469            i = i + 1
470
471        return;
472
473    # Find beads that are not in allowed volume
474
476                        aList_box_x,aList_box_y,aList_box_z,rsearch):
477
479        num_boxes = len(aList_box_x)
480        contact_limit_sq = rsearch**2
481
482        count = 0
484
488
489            inbox = 'no'
490            i = 0
491            while i < num_boxes:
492
493                x_box = float(aList_box_x[i])
494                y_box = float(aList_box_y[i])
495                z_box = float(aList_box_z[i])
497
498                if dsq < contact_limit_sq:
499                    inbox = 'yes'
500                    i = num_boxes
501
502                i = i + 1
503
504            if inbox == 'no':
505                aList_contacts.append(count)
506
507            count = count + 1
508
509        return;
510
511    # Compute a P(r)
512
514
515        num_hist = len(aList_pr_model)
516        count = 0
517        while count < num_hist:
518            aList_pr_model[count] = 0.0
519            count = count + 1
520
522        max_dr = (float(num_hist)-1.0)*hist_grid
523
524        i = 0
526
527            j = 0
528            while j < i:
529
532                dr = min(dr,max_dr)
533
534                # Find pointers and do un-interpolation
535
536                dr_grid = dr/hist_grid
537                int_dr_grid = int(dr_grid)
538
539                int_dr_grid = min(int_dr_grid,num_hist-2)
540                ip_low = int_dr_grid
541                ip_high = ip_low + 1
542
543                ip_high_frac = dr_grid - float(int_dr_grid)
544                ip_low_frac = 1.0 - ip_high_frac
545
546                aList_pr_model[ip_low] = float(aList_pr_model[ip_low]) + ip_low_frac
547                aList_pr_model[ip_high] = float(aList_pr_model[ip_high]) + ip_high_frac
548
549                j = j + 1
550            i = i + 1
551
552        return;
553
554    # Score for rms difference between observed and model histograms
555
556    def pr_dif(aList_pr,aList_pr_model,skip):
557
558        num_hist = len(aList_pr)
559        delta_hist_sum = 0.0
560        delta_hist_sum_sq = 0.0
561        hist_sum = 0.0
562
563        i = skip
564        while i < num_hist:
565
566            model = float(aList_pr_model[i])
567            data = float(aList_pr[i])
568            delta_hist = abs(data - model)
569            delta_hist_sum = delta_hist_sum + delta_hist
570            hist_sum = hist_sum + data
571
572            delta_hist_sum_sq = delta_hist_sum_sq + delta_hist*delta_hist
573
574            i = i + 1
575
576        mean_hist_sum = hist_sum/(num_hist - skip)
577        delta_hist_sum_sq = delta_hist_sum_sq/(num_hist - skip)
578        delta_hist_sum_sq = math.sqrt(delta_hist_sum_sq)/mean_hist_sum
579
580        return delta_hist_sum_sq;
581
582    # Statistics for fractional difference between observed and model histograms
583
584    def pr_rfactor(aList_pr,aList_pr_sd,aList_pr_model,skip):
585
586        num_hist = len(aList_pr)
587        delta_hist_sum = 0.0
588        hist_sum = 0.0
589
590        i = skip
591        while i < num_hist:
592
593            model = float(aList_pr_model[i])
594            exp = float(aList_pr[i])
595            delta_hist = exp - model
596            delta_hist_sum = delta_hist_sum + abs(delta_hist)
597            hist_sum = hist_sum + exp
598
599            i = i + 1
600
601        delta_hist_sum = delta_hist_sum/hist_sum
602
603        return delta_hist_sum;
604
605    # Compute the VDW energy for a interaction
606
608
610            vdw = 0.0
611        else:
612            dr_e6 = dr**6
613            dr_e12 = dr_e6**2
614            vdw = econ12/dr_e12 - 2.0*econ6/dr_e6
615            vdw = max(vdw,e_width)
616
617        return vdw;
618
619    # Set a random distribution of beads in a box with maximum extent dmax
620
623
624        half_side = 0.5*dmax
625
626        scale_xy = 1.0 - bias_z
627        scale_z = 1.0 + bias_z
628        x_range = scale_xy * half_side
629        y_range = scale_xy * half_side
630        z_range = scale_z * half_side
631
632        num_ops = len(aList_symm)
633
634        i = 0
636
637            triangle = random.triangular(-0.7,0.7,0.0)
638            x = triangle*x_range
639            triangle = random.triangular(-0.7,0.7,0.0)
640            y = triangle*y_range
641            triangle = random.triangular(-0.7,0.7,0.0)
642            z = triangle*z_range
643
647
648            j = 0
649            while j < num_ops:
650                aList_s = aList_symm[j]
651                m11 = float(aList_s[0])
652                m12 = float(aList_s[1])
653                m21 = float(aList_s[2])
654                m22 = float(aList_s[3])
655
656                xs = m11*x + m12*y
657                ys = m21*x + m22*y
658                zs = z
662
663                j = j + 1
664
665            i = i + num_symm
666
667        return;
668
669    # Read experimentalal P(r) from GNOM output file
670
672                aList_pr_model_test,aList_pr_model_test2,inFile):
673
674        angstrom_scale = 1.0
675        Bins,Dbins,BinMag = data['Pair']['Distribution']
676
677        aList_r += list(Bins)
678        aList_pr += list(BinMag)
679        aList_pr_sd += list(np.ones_like(Bins)/100.)
680        aList_pr_model += list(np.zeros_like(Bins))
681        aList_pr_model_test += list(np.zeros_like(Bins))
682        aList_pr_model_test2 += list(np.zeros_like(Bins))
683
684#        file = open(inFile,'r')
686#        file.close()
687#
688#        parse_data = 'no'
689#        for eachLine in allLines:
690#
691#            if parse_data == 'yes':
692#
693#                aList = eachLine.split()
694#                num_params = len(aList)
695#
696#                if num_params == 3:
697#                    r = float(aList[0])
698#                    pr = float(aList[1])
699#                    pr_sd = float(aList[2])
700#
701#                    aList_pr.append(pr)
702#                    aList_pr_sd.append(pr_sd)
703#                    aList_r.append(r)
704#                    aList_pr_model.append(0.0)
705#                    aList_pr_model_test.append(0.0)
706#                    aList_pr_model_test2.append(0.0)
707#
708#            if eachLine.find('R          P(R)      ERROR') > -1:
709#                parse_data = 'yes'
710#
711        num_hist = len(aList_r)
712        hist_grid = float(aList_r[1]) - float(aList_r[0])
713
714
715        # Heuristic for checking units
716#        test_r = max(aList_r)
717#        if test_r < 30.0:
718#
719#            aString = 'P(r)appears to be in nm. Converting to Angstrom units'
720#            print (aString)
721#            angstrom_scale = 10.0
722#
723#            i = 0
724#            while i < num_hist:
725#                aList_r[i] = angstrom_scale * aList_r[i]
726#                i = i + 1
727#
728#            hist_grid = angstrom_scale * hist_grid
729#
730#        i = 0
731#        while i < num_hist:
732#            r = float(aList_r[i])
733#            r_calc = float(i)*hist_grid
734#
735#            if abs(r - r_calc) > 0.15:
736#                aString = 'Input P(r) grid is irregular! Exiting'
737#                print (aString)
738#                time.sleep(4)
739#                sys.exit(1)
740#
741#            i = i + 1
742#
743        dmax = aList_r[num_hist-1]
744
745        # Pad histogram by 5 Angstrom
746
748
749        i = 0
751            r = dmax + float(i)*hist_grid
752            aList_pr.append(0.0)
753            aList_pr_sd.append(0.0)
754            aList_r.append(r)
755            aList_pr_model.append(0.0)
756            aList_pr_model_test.append(0.0)
757            aList_pr_model_test2.append(0.0)
758            i = i + 1
759
760        return (dmax,hist_grid,num_hist,angstrom_scale);
761
762    # Scale P(r) onto model P(r) assuming same grid
763
764    def scale_pr(aList_pr,aList_pr_sd,aList_pr_model):
765
766        num_hist = len(aList_pr)
767        total_dr = 0.0
768        total_pr = 0.0
769
770        i = 0
771        while i < num_hist:
772            total_dr = total_dr + float(aList_pr_model[i])
773            total_pr = total_pr + float(aList_pr[i])
774            i = i + 1
775
776        scale = total_dr/total_pr
777
778        i = 0
779        while i < num_hist:
780            aList_pr[i] = scale*float(aList_pr[i])
781            aList_pr_sd[i] = scale * float(aList_pr_sd[i])
782            i = i + 1
783
784        return;
785
786    # Return a non-zero distance between two coordinates
787
788    def get_dr(x1,y1,z1,x2,y2,z2):
789
790        x1 = float(x1)
791        y1 = float(y1)
792        z1 = float(z1)
793        x2 = float(x2)
794        y2 = float(y2)
795        z2 = float(z2)
796        dr = (x1 - x2)**2 + (y1-y2)**2 + (z1-z2)**2
797        dr = max(dr,0.1)
798        dr = math.sqrt(dr)
799
800        return dr;
801
803
805
807
808        xsum = 0.0
809        ysum = 0.0
810        zsum = 0.0
812
813        i = 0
815
817
820                xsum = xsum + float(aList_beads_x[i])
821                ysum = ysum + float(aList_beads_y[i])
822                zsum = zsum + float(aList_beads_z[i])
823
824            i = i + 1
825
830            delta = (xsum - x)**2 + (ysum - y)**2 + (zsum - z)**2
831            delta = math.sqrt(delta)
832        else:
833            delta = 0.0
834
835        return delta;
836
837    # Obtain mean of total VDW energy
838
840
842        vdw_all = 0.0
843
844        i = 0
846            j = 0
847            while j < i:
851                vdw_all = vdw_all + vdw
852                j = j + 1
853            i = i + 1
854
856
857        return vdw_all;
858
859    # Energy minimize
860
862
864        eps12 = eps**12
865        eps6 = eps**6
867        scale = 0.0
868        icount = -1
869
871        num_ops = len(aList_symm)
872        num_cycles = 51
873
874        i = 0
875        while i < num_cycles:
876
877            icount = icount + 1
878
882
883            sum_forces_scale = 0.0
884
885            k = 0
886            while k < nbeads - num_ops:
887
891
892                fx = 0.0
893                fy = 0.0
894                fz = 0.0
895                j = 0
897
901
902                    dr = get_dr(xold,yold,zold,xj,yj,zj)
903
904                    # Truncate very steep
905                    dr = min(eps,dr)
906
908                        dr_sq = dr*dr
909                        dr12 = dr_sq**6
910                        dr6 = dr_sq**3
911
912                        dx = xold - xj
913                        dy = yold - yj
914                        dz = zold - zj
915
916                        force = (1.0/dr_sq)*(eps12/dr12 - 0.5*eps6/dr6)
917                        fx = fx + force*dx
918                        fy = fy + force*dy
919                        fz = fz + force*dz
920
921                        sum_forces_scale = sum_forces_scale + abs(fx) + abs(fy) + abs(fz)
922
923                    j = j + 1
924
925                #
926                xstep = scale*fx
927                ystep = scale*fy
928                zstep = scale*fz
929
930                if xstep > 0.0:
931                    xstep = min(xstep,step_max)
932                else:
933                    xstep = max(xstep,-step_max)
934
935                if ystep > 0.0:
936                    ystep = min(ystep,step_max)
937                else:
938                    ystep = max(ystep,-step_max)
939
940                if zstep > 0.0:
941                    zstep = min(zstep,step_max)
942                else:
943                    zstep = max(zstep,-step_max)
944
945                xtest = xold + xstep
946                ytest = yold + ystep
947                ztest = zold + zstep
951
952                # Apply shifs to symm positions
953                l = 0
954                while l < num_ops:
955                    aList_s = aList_symm[l]
956                    m11 = float(aList_s[0])
957                    m12 = float(aList_s[1])
958                    m21 = float(aList_s[2])
959                    m22 = float(aList_s[3])
960
961                    xs = m11*xtest + m12*ytest
962                    ys = m21*xtest + m22*ytest
963                    zs = ztest
964
968
969                    l = l + 1
970
971                #
972
973                k = k + num_ops + 1
974
975            # Apply shifted positions after first cycle
976            if i > 0:
977
978                m = 0
983                    m = m + 1
984
985            #
986
989
991
992            if icount == 0:
993                aString = 'Emin cycle: ' + str(i) + ' Energy: ' + str('%4.2f'%(vdw_all))
994                print (aString)
995                icount = -10
996
997            if vdw_all < 0.0:
998                i = num_cycles
999
1000            i = i + 1
1001
1002        return;
1003
1004    # Set up symmetry operators for rotational symmetry
1005
1006    def make_symm(aList_symm,num_symm):
1007
1008        angle_step = 360.0/float(num_symm)
1009
1010        i = 1
1011        while i < num_symm:
1012            theta = float(i) * angle_step
1014            cos_theta = math.cos(theta)
1015            sin_theta = math.sin(theta)
1016            aList_s = [cos_theta,sin_theta,-sin_theta,cos_theta]
1017            aList_symm.append(aList_s)
1018            i = i + 1
1019
1020        return aList_symm;
1021
1022    # Set up a shift vector in P(r) for a change in bead position
1023
1024    def pr_shift_atom(aList_pr_model_test2,x1,y1,z1,\
1026
1027        num_hist = len(aList_r)
1028        max_dr = (float(num_hist)-1.0)*hist_grid
1030
1031        i = 0
1032        while i < num_hist:
1033            aList_pr_model_test2[i] = 0.0
1034            i = i + 1
1035
1036        i = 0
1038
1039            if i != ii:
1043                dr = get_dr(x1,y1,z1,x2,y2,z2)
1044                dr = min(dr,max_dr)
1045                dr_grid = dr/hist_grid
1046                int_dr_grid = int(dr_grid)
1047                int_dr_grid = max(int_dr_grid,0)
1048                int_dr_grid = min(int_dr_grid,num_hist-2)
1049                ip_low = int_dr_grid
1050                ip_high = ip_low + 1
1051                ip_high_frac = dr_grid - float(int_dr_grid)
1052                ip_low_frac = 1.0 - ip_high_frac
1053                aList_pr_model_test2[ip_low] = float(aList_pr_model_test2[ip_low]) + ip_low_frac
1054                aList_pr_model_test2[ip_high] = float(aList_pr_model_test2[ip_high]) + ip_high_frac
1055
1056            i = i + 1
1057
1058        return;
1059
1060    # Recenter set of beads to origin
1061
1063
1065        xsum = 0.0
1066        ysum = 0.0
1067        zsum = 0.0
1068
1069        i = 0
1071            xsum = xsum + float(aList_beads_x[i])
1072            ysum = ysum + float(aList_beads_y[i])
1073            zsum = zsum + float(aList_beads_z[i])
1074            i = i + 1
1075
1079
1080        i = 0
1085            i = i + 1
1086
1087        return;
1088
1089    #############
1090    # EXECUTION #
1091    #############
1092
1093    #profiling start
1094    pr = cProfile.Profile()
1095    pr.enable()
1096    time0 = time.time()
1097
1098    version_aString = 'Program: SHAPES version 1.3'
1099
1100    print (version_aString)
1101    aString = 'Author: John Badger'
1102    print (aString)
1104    print (aString)
1105    aString = 'License: GNU GPLv3'
1106    print (aString)
1107
1108    localtime = time.asctime( time.localtime(time.time()) )
1109    aString = 'Starting time: ' + str(localtime) + '\n'
1110    print (aString)
1111
1112#    aList_summary = []
1113#    aList_summary.append(version_aString)
1114#    aList_summary.append(str(localtime))
1115
1116    ######################
1117    # Start up parmeters #
1118    ######################
1119#        data['Shapes'] = {'outName':'','NumAA':100,'Niter':1,'AAscale':1.0,'Symm':1,'bias-z':0.0,
1120#            'inflateV':1.0,'AAglue':0.0}
1121
1123    num_sols = 1
1124    num_aa = 1.0
1125    num_symm = 1
1126    bias_z = 0.0
1127    inflate = 1.0
1128    prefix = ''
1129    surface_scale = 0.0
1130    starting_pdb = 'no'
1131    inFile = 'none'
1132    pdbfile_in  = 'none'
1133    shapeDict = data['Shapes']
1134    prefix = shapeDict['outName']
1136    num_sols = shapeDict['Niter']
1137    num_aa = shapeDict['AAscale']
1138    num_symm = shapeDict['Symm']
1139    bias_z = shapeDict['bias-z']
1140    inflate = shapeDict['inflateV']
1141    surface_scale = shapeDict['AAglue']
1142    pdbOut = shapeDict['pdbOut']
1143    box_step = shapeDict.get('boxStep',4.0)
1144    Phases = []
1145    Patterns = []
1146    PRcalc = []
1147
1148#    # Parse
1149#
1150#    if os.path.exists('shapes_ip.txt'):
1151#        file = open('shapes_ip.txt','r')
1153#        file.close()
1154#    else:
1155#        aString = 'The local parameter file shapes_ip.txt was not found ! Exiting'
1156#        print (aString)
1157#        time.sleep(4)
1158#        sys.exit(1)
1159#
1160#    for eachLine in allLines:
1161#
1162#        aList = eachLine.split()
1163#
1164#        num_params = len(aList)
1165#        if num_params > 0:
1166#
1167#            if aList[0] == 'num_amino_acids':
1169##            elif aList[0] == 'input_pr':
1170##                inFile = str(aList[1])
1171#            elif aList[0] == 'num_solns':
1172#                num_sols = int(aList[1])
1173#            elif aList[0] == 'num_aa_scale':
1174#                num_aa = float(aList[1])
1175#            elif aList[0] == 'symm':
1176#                num_symm = int(aList[1])
1177#            elif aList[0] == 'bias_z':
1178#                bias_z = float(aList[1])
1179#            elif aList[0] == 'inflate_vol':
1180#                inflate = float(aList[1])
1181#            elif aList[0] == 'pdb_start':
1182#                pdbfile_in = str(aList[1])
1183#            elif aList[0] == 'id':
1184#                prefix = str(aList[1]) + '_'
1185#            elif aList[0] == 'glue':
1186#                surface_scale = float(aList[1])
1187
1188
1189    # Check inputs
1190
1191    if num_sols > 0:
1192        aString = 'Number of runs: ' + str(num_sols)
1193        print (aString)
1194    else:
1195        aString = 'Zero reconstruction runs specified! Exiting'
1196        print (aString)
1197        time.sleep(4)
1198        sys.exit(1)
1199
1200    #
1202        if os.path.exists(pdbfile_in):
1203            aString = 'Will use CA atoms from ' + str(pdbfile_in) + ' as the initial bead distribution.'
1204            print (aString)
1205            starting_pdb = 'yes'
1206        else:
1207            aString = 'Zero amino acid count specified and no starting file found. Exiting'
1208            print (aString)
1209            time.sleep(4)
1210            sys.exit(1)
1211    else:
1212        aString = 'Number of amino acids: ' + str(nbeads)
1213        print (aString)
1214
1215    #
1216#    if os.path.exists(inFile):
1217#        aString = 'Input P(r) file name: ' + str(inFile)
1218#        print (aString)
1219#    else:
1221#        print (aString)
1222#        time.sleep(4)
1223#        sys.exit(1)
1224
1225    #
1226    if num_aa == 0.0:
1227        aString = 'Scale for amino acid count to particle number cannot be zero! Exiting'
1228        print (aString)
1229        time.sleep(4)
1230        sys.exit(1)
1231    else:
1232        aString = 'Scale aa to bead count: ' + str(num_aa)
1233        print (aString)
1234
1235    #
1236    if num_symm == 0:
1237        aString = 'Rotational symmetry cannot be zero! Set to 1 for no symmetry. Exiting'
1238        print (aString)
1239        time.sleep(4)
1240        sys.exit(1)
1241    else:
1242        aString = 'Point symmetry: ' + str(num_symm)
1243        print (aString)
1244
1245    #
1246    if bias_z > 0.2:
1247        aString = 'Max bias on Z axis for initial particle distribution is 0.2 (rods). Reset to 0.2.'
1248        print (aString)
1249        bias_z = 0.2
1250    elif bias_z < -0.2:
1251        aString = 'Min bias on Z axis for initial particle distribution is -0.2 (disks). Reset to -0.2.'
1252        print (aString)
1253        bias_z = -0.2
1254    else:
1255        aString = 'Z-axis bias: ' + str(bias_z)
1256        print (aString)
1257
1258    #
1259    if inflate < 0.0:
1260        aString = 'Inflation of PSV cannot be less than zero! Exiting'
1261        print (aString)
1262        time.sleep(4)
1263        sys.exit(1)
1264    elif inflate > 2.0:
1265        aString = 'Inflation of PSV cannt be greater than 2.0! Exiting'
1266        print (aString)
1267        time.sleep(4)
1268        sys.exit(1)
1269    else:
1270        aString = 'PSV inflation factor: ' + str(inflate)
1271        print (aString)
1272
1273    #
1274    if surface_scale > 0.0:
1275        aString = 'Cavity weight: ' + str(surface_scale)
1276        print (aString)
1277
1278    ########## UNIVERSAL CONSTANTS ######################
1279
1280    # No of macrocycles (gives extra cycles at constant volume after num_contract)
1281    niter = 160
1282
1283    # No of contraction cycles
1284    num_contract = 140
1285
1286    # Number of cycles at each fixed volume
1287    num_micro_cyc = 10
1288
1289    # Final quench
1290    num_sa_max = niter - num_micro_cyc
1291
1292    # Initial scale for P(r) shifts versus E shifts
1293    hscale = 3000.0
1294
1295    # Standard deviation for annealing acceptance (cf well-depth of -1 unit for two beads)
1296    sd_mc = float(num_symm) * 2.0
1297
1298    # Fiddle factor for keeping the accessible, molecular volume larger than PSV
1299    scale_vol = 1.15
1300
1301    # Standard amino acid volume MW = 110.0 x 1.21 i.e. mean mw x mw-to-vol-scale
1303
1304    # Bead separation for best packing ~5.6 (I think)
1305    #- 75% better than rectangular grid 5.1 for this amino acid vol
1307
1308    # Usually num_aa is unity. Adjust parameters otherwise
1309    if num_aa != 1 and nbeads != 0:
1313
1314    # Increase bead separation for inflated volumes
1316
1317    # Partial specific volumes at start and end
1318
1319    if starting_pdb == 'yes':
1320        nmols_vol_start = 1.1 * inflate
1321    else:
1322        nmols_vol_start = 2.0 * inflate
1323
1324    nmols_vol_end = 1.0 * inflate
1325    nmols_vol_subtract = nmols_vol_start - nmols_vol_end
1326
1327    # Box parametere
1328#    box_step = 4.0      #5.0
1329    box_pt_vol = box_step*box_step*box_step
1330
1331    # Energy parameters - flat bottomed VDW (2.0A for a 5.6A bead separation)
1332
1336    r_width = bead_sep + well_width
1337    r_width6 = r_width**6
1338    r_width12 = r_width6**2
1339    e_width = econ12/r_width12 - 2.0*econ6/r_width6
1341    abs_e_width = abs(e_width)
1342
1343    # Range for box identification (might need to increase for poor data)
1344    rsearch = (bead_sep + 0.5*well_width)*1.5
1345
1346    # Search range for optional cavity inhibition energy term
1349
1350    # Setup symmetry operators
1351
1352    aList_symm = []
1353    aList_symm = make_symm(aList_symm,num_symm)
1354    num_ops = len(aList_symm)
1355
1357
1358    aList_r = []
1359    aList_pr = []
1360    aList_pr_sd = []
1361    aList_pr_model = []
1362    aList_pr_model_test = []
1363    aList_pr_model_test2 = []
1364
1366                                            aList_pr_model,aList_pr_model_test,\
1367                                            aList_pr_model_test2,inFile)
1368
1369#    dmax_over2 = dmax/2.0
1370    num_hist = len(aList_r)
1371
1372    aString = 'Number of points read from P(r): ' + str(num_hist_in)
1373    print (aString)
1374    aString = 'Grid sampling: ' + str(hist_grid) + ' Dmax: ' + str(dmax)
1375    print (aString)
1376
1377    # Skip over initial points in scoring
1378
1379    skip = r_width/float(num_hist)
1380    skip = int(skip) + 2
1381
1382    # Read intensity data that was used for P(r)
1383
1384    aList_q = []
1385    aList_i = []
1386    aList_i_sd = []
1387
1389
1390    num_intensities = len(aList_q)
1391    aString = 'Number of intensity data points read: ' + str(num_intensities)
1392    print (aString)
1393
1394    #########################
1395    # CYCLE OVER SOLUTIONS  #
1396    #########################
1397
1398    i_soln = 0
1399    while i_soln < num_sols:
1400
1401        file_no = str(i_soln + 1)
1402
1403        aString = '\nReconstruction trial: ' + str(file_no)
1404        print (aString)
1405
1406        aString  = 'Trial:' + file_no
1407#        aList_summary.append(aString)
1408
1410#        file_pr = prefix + 'pr_calc_' + file_no + '.dat'
1411        file_psv = prefix + 'psv_shape_' + file_no + '.pdb'
1412#        file_intensity = prefix + 'intensity_' + file_no + '.dat'
1413
1414        # Setup initial bead distribution
1415
1419
1420        # Re-initialize standard deviation for annealing acceptance
1421        sd_mc = float(num_symm) * 2.0
1422
1423        # Set random bead distribution
1424
1425        if starting_pdb == 'yes':
1428            num_symm = 1
1430            print (aString)
1431            aString = 'Symmetry set to 1 (required)'
1432            print (aString)
1433            aString = 'Input center was shifted to the origin'
1434            print (aString)
1435        else:
1439            print (aString)
1440
1441        # Protein partial specific volume
1443
1444        # Histogram of inter-bead distance
1445
1447
1448        # Scale experimental P(r) and model histogram
1449
1450        scale_pr(aList_pr,aList_pr_sd,aList_pr_model)
1451
1452        # Minimize initial energy using expanded VDW
1453
1454        if starting_pdb != 'yes':
1455            aString = 'Minimize energy of initial positions'
1456            print (aString)
1460        else:
1461            aString = 'Skipping energy minimization of initial positions'
1462            print (aString)
1463
1464        # Get the initial score between observed and calculated P(r)
1465
1466        hist_score_best = pr_dif(aList_pr,aList_pr_model,skip)
1467        aString = 'Initial rms P(r): ' + str('%4.3f'%(hist_score_best))
1468        print (aString)
1469
1470        ###########################
1471        # Iterate                 #
1472        ###########################
1473
1474        num_boxes = 0
1475        count_boxing = 0
1476        fraction_psv = 0
1477        success_rate_all = 0.0
1478        box_iter = num_micro_cyc - 1
1479
1480        sum_delta_pack = 0.0
1481
1482        count_it = 0
1483        while count_it < niter:
1484
1485            success = 0
1486            count_hist_yes = 0
1487            sum_e = 0.0
1488            sum_h = 0.0
1489
1490            # Find populated volumes and fix solution every 10 macrocycles
1491
1492            box_iter = box_iter + 1
1493
1494            if box_iter == num_micro_cyc:
1495
1496                box_iter = 0
1497                count_boxing = count_boxing + 1
1498
1499                if count_it < num_contract - 1:
1500                    scale = float(count_it)/float(num_contract)
1501                else:
1502                    scale = 1.0
1503
1504                # Establish confinement volume using a local average
1505
1506                aList_box_x_all = []
1507                aList_box_y_all = []
1508                aList_box_z_all = []
1509                aList_box_score = []
1510
1512
1514                # rsearch_use = (2.0 - scale)*rsearch
1515
1517                        aList_box_x_all,aList_box_y_all,aList_box_z_all,\
1518                        aList_box_score,box_step,dmax,rsearch)
1519
1520                aList_box_x = []
1521                aList_box_y = []
1522                aList_box_z = []
1523
1524                psv_ratio = nmols_vol_start - scale*nmols_vol_subtract
1525                vol_target = scale_vol*(psv_ratio*psv_vol)
1526
1527                set_vol(aList_box_x_all,aList_box_y_all,aList_box_z_all,\
1528                        aList_box_score,aList_box_x,aList_box_y,aList_box_z,\
1529                        vol_target,box_pt_vol)
1530
1531                num_boxes = len(aList_box_x)
1532                fraction_psv = float(num_boxes)*box_pt_vol/psv_vol
1533
1534                # Find beads that are ouside the allowed volume
1535
1536                aList_contacts = []
1538                        aList_box_x,aList_box_y,aList_box_z,rsearch)
1539                num_outof_box = len(aList_contacts)
1540
1541                aString = 'Target volume: ' + str('%4.2f'%(scale_vol*psv_ratio)) + ' Actual volume: ' + \
1542                          str('%4.2f'%(fraction_psv)) + ' Beads outside volume: ' + str(num_outof_box)
1543                print (aString)
1544
1545                # Recalculate P(r) and rescore for reliability
1546
1548                hist_score_best = pr_dif(aList_pr,aList_pr_model,skip)
1549
1550                # Reset SA deviation if mean success rate over last trials is under 0.1
1551
1552                mean_success_rate = float(success_rate_all)/float(num_micro_cyc)
1553
1554                if count_it < num_contract and count_boxing != 1:
1555
1556                    if mean_success_rate < 0.1:
1557                        sd_mc = 1.3*sd_mc
1558                        aString = 'Raising allowed energy deviation to ' + str('%4.2f'%(sd_mc))
1559                        print (aString)
1560
1561                    if mean_success_rate > 0.2:
1562                        sd_mc = 0.7*sd_mc
1563                        aString = 'Reducing allowed energy deviation to ' + str('%4.2f'%(sd_mc))
1564                        print (aString)
1565
1566                success_rate_all = 0.0
1567
1568            # Make one macrocycle that is a run over the nbeads
1569
1570            ii = 0
1572
1573                # Initialize
1574
1575                energy_old = 0.0
1576                energy_new = 0.0
1577
1578                i = 0
1579                while i < num_hist:
1580                    aList_pr_model_test[i]  = 0.0
1581                    i = i + 1
1582
1583                # Select a target bead and make trial shift
1584
1588
1589                ibox = random.randint(0,num_boxes-1)
1590                xtest = float(aList_box_x[ibox]) + random.uniform(-rsearch,rsearch)
1591                ytest = float(aList_box_y[ibox]) + random.uniform(-rsearch,rsearch)
1592                ztest = float(aList_box_z[ibox]) + random.uniform(-rsearch,rsearch)
1593
1594                # Calculate and capture symmetry mates
1595
1596                aList_temp_save_x = []
1597                aList_temp_save_y = []
1598                aList_temp_save_z = []
1599                aList_symm_x = []
1600                aList_symm_y = []
1601                aList_symm_z = []
1602
1603                l = 0
1604                while l < num_ops:
1605                    aList_s = aList_symm[l]
1606                    m11 = float(aList_s[0])
1607                    m12 = float(aList_s[1])
1608                    m21 = float(aList_s[2])
1609                    m22 = float(aList_s[3])
1610
1611                    xs = m11*xtest + m12*ytest
1612                    ys = m21*xtest + m22*ytest
1613                    zs = ztest
1614
1615                    aList_symm_x.append(xs)
1616                    aList_symm_y.append(ys)
1617                    aList_symm_z.append(zs)
1618
1619                    ipt = ii + l + 1
1623
1624                    l = l + 1
1625
1626                # Get initial VDW energy for interactions of this bead with all others
1627
1628                i = 0
1630                    if i != ii:
1634                        dr_old = get_dr(xold,yold,zold,x,y,z)
1636                        energy_old = energy_old + num_symm*vdw_old
1637                    i = i + 1
1638
1639                # Get initial contributions to P(r)
1640
1643
1644                i = 0
1645                while i < num_hist:
1646                    aList_pr_model_test[i] = aList_pr_model_test2[i]
1647                    i = i + 1
1648
1649                # Get VDW energy for interactions of the shifted bead with all others
1650
1651                l = 0
1652                while l < num_ops:
1653                    ipt = ii + l + 1
1657                    l = l + 1
1658
1659                i = 0
1661                    if i != ii:
1665                        dr_new = get_dr(xtest,ytest,ztest,x,y,z)
1667                        energy_new = energy_new + num_symm*vdw_new
1668                    i = i + 1
1669
1670                # Get cavity energy difference
1671
1676
1678                sum_delta_pack = sum_delta_pack + abs(delta_pack)
1679
1680                # Get shifted contributions to P(r)
1681
1684
1685                # Get net shift in contribution to P(r)
1686
1687                i = 0
1688                while i < num_hist:
1689                    aList_pr_model_test[i] = aList_pr_model_test2[i] - aList_pr_model_test[i]
1690                    i = i + 1
1691
1692                # Get statistic for agreement with P(r) after accumulating shift vectors
1693
1694                i = 0
1695                while i < num_hist:
1696                    aList_pr_model_test[i] = float(aList_pr_model[i]) + num_symm*float(aList_pr_model_test[i])
1697                    i = i + 1
1698
1699                hist_score = pr_dif(aList_pr,aList_pr_model_test,skip)
1700
1701                # Scoring shifts
1702
1703                delta_h = hist_score - hist_score_best
1704                delta_e = energy_new - energy_old + delta_pack
1705
1706                # Recalibrate scale so impact of energy and P(r) is equal on plausible shifts
1707
1708                if delta_e < abs_e_width:
1709                    sum_e = sum_e + abs(delta_e)
1710                    sum_h = sum_h + abs(delta_h)
1711
1712                # Monitor potential moves based of P(r)
1713
1714                if delta_h < 0.0:
1715                    count_hist_yes = count_hist_yes + 1.0
1716
1717                # Acceptance and update
1718
1719                score = delta_e + delta_h*hscale
1720
1721                if count_it < num_sa_max:
1722                    barrier = abs(random.gauss(0.0,sd_mc))
1723                else:
1724                    barrier = 0.0
1725
1726                if score < barrier:
1727
1728                    success = success + 1.0
1729                    hist_score_best = hist_score
1730
1731                    # Update model but symmetry positions that were already put in
1732
1736
1737                    # Update P(r)
1738
1739                    i = 0
1740                    while i < num_hist:
1741                        aList_pr_model[i] = aList_pr_model_test[i]
1742                        i = i + 1
1743
1744                else:
1745
1746                    # Revert to original model
1747
1751
1752                    l = 0
1753                    while l < num_ops:
1754                        ipt = ii + l + 1
1758                        l = l + 1
1759                #
1760
1761                ii = ii + num_symm
1762
1763            # Get energy statistics at end of macrocycle
1764
1766
1767            # Rescale and convergence statistics
1768
1769            if sum_h > 0.0:
1770                hscale = sum_e/sum_h
1771
1774            success_rate_all = success_rate_all + success_rate
1775
1776            aString = 'Cycle ' + str(count_it+1) + ' Moves ' + str('%.2f'%(success_rate)) + \
1777                      ' Possibles ' + str('%.2f'%(count_hist_yes)) + ' rms P(r) '+ str('%4.3f'%(hist_score)) + \
1778                      ' Energy ' + str('%4.2f'%(vdw_all))
1779            print (aString)
1780
1781            # Debug statitics. Weight of 10 gives about 1.0
1783            #print (sum_delta_pack)
1784            #
1785
1786            count_it = count_it + 1
1787
1788        #####################
1789        # ANALYZE AND WRITE #
1790        #####################
1791
1792        aString = '\nFinal model statistics'
1793        print (aString)
1794
1796        hist_score_best = pr_dif(aList_pr,aList_pr_model,skip)
1797
1798        # P(r) fitting statistics
1799        delta_hist_sum = pr_rfactor(aList_pr,aList_pr_sd,aList_pr_model_test,skip)
1800
1801        aString = 'Delta P(r): ' + str('%4.3f'%(delta_hist_sum))
1802        print (aString)
1803
1804        # Get final energy
1806
1807        aString = 'VDW energy: ' + str('%4.2f'%(vdw_all))
1808        print (aString)
1809
1811        # Write out beads as pseudo a PDB file
1812        if pdbOut:
1814
1815        # Calculate and write final PSV shape
1816
1817        aList_box_x_all = []
1818        aList_box_y_all = []
1819        aList_box_z_all = []
1820        aList_box_score = []
1821
1823                aList_box_x_all,aList_box_y_all,aList_box_z_all,\
1824                aList_box_score,box_step,dmax,rsearch)
1825
1826        aList_box_x = []
1827        aList_box_y = []
1828        aList_box_z = []
1829        psv_vol_use = psv_vol*inflate
1830
1831        set_vol(aList_box_x_all,aList_box_y_all,aList_box_z_all,\
1832                aList_box_score,aList_box_x,aList_box_y,aList_box_z,\
1833                psv_vol_use,box_pt_vol)
1834
1835        num_boxes = len(aList_box_x)
1836        fraction_psv = num_boxes*box_pt_vol/psv_vol
1837
1838        # Correct final volume if initial estimate is too small
1839
1840        fraction_psv_use = num_boxes*box_pt_vol/psv_vol_use
1841
1842        if fraction_psv_use < 1.0:
1843            aList_box_x = []
1844            aList_box_y = []
1845            aList_box_z = []
1846            vol_use = 1.05*psv_vol_use
1847            set_vol(aList_box_x_all,aList_box_y_all,aList_box_z_all,\
1848                    aList_box_score,aList_box_x,aList_box_y,aList_box_z,\
1849                    vol_use,box_pt_vol)
1850
1851            num_boxes = len(aList_box_x)
1852            fraction_psv = float(num_boxes)*box_pt_vol/psv_vol
1853
1854        aString = 'Final PSV of protein envelope: ' + str('%4.2f'%(fraction_psv))
1855        print (aString)
1856
1857        # Write input and model P(r)
1858#        pr_writer(aList_pr,aList_r,aList_pr_model,file_pr)
1859        PRcalc.append([aList_r,aList_pr,copy.copy(aList_pr_model),delta_hist_sum])
1860
1861        # Calculate comparison versus intensities
1862
1863        if num_intensities > 0:
1864
1865            # calculate intensity
1866            aList_i_calc = []
1869
1870            # Scale and obtain statistics
1871            (chi_sq,rvalue) = score_Ic(aList_q,aList_i,aList_i_sd,aList_i_calc)
1872
1873            aString = 'Rvalue: ' + str('%4.3f'%(rvalue)) + ' CHI-squared: ' + str('%4.3f'%(chi_sq))
1874            print (aString)
1875
1876            # Write output intensity file
1877            Patterns.append([aList_q,aList_i,aList_i_calc,rvalue])
1878#            write_all_data(file_intensity,aList_q,aList_i,aList_i_calc,aString)
1879
1880#            aString = 'Output intensity file: ' + str(file_intensity)
1881#            print (aString)
1882
1883        else:
1884
1885            chi_sq = 'NA'
1886
1887        # Write final volume
1888
1889        delta_hist_sum = '%4.3f'%(delta_hist_sum)
1890        vdw_all = '%4.2f'%(vdw_all)
1891        fraction_psv = '%4.2f'%(fraction_psv)
1892        chi_sq = '%4.3f'%(chi_sq)
1893
1894        aString = 'REMARK     P(r) dif:'+ str(delta_hist_sum) + ' E:'\
1895                  + str(vdw_all) + ' CHIsq:' + str(chi_sq) + \
1896                  ' PSV:' + str(fraction_psv)
1897
1898        Phases.append([file_psv.split('.')[0],aList_box_x,aList_box_y,aList_box_z])
1899        if pdbOut:
1900            pdb_writer(aList_box_x,aList_box_y,aList_box_z,file_psv,aString)
1901
1902        # Write Summary
1903
1904        aString = 'P(r) dif:' + str(delta_hist_sum) + ' E:' \
1905                         + str(vdw_all) + ' CHISQ:' + str(chi_sq) + ' PSV:' + str(fraction_psv)
1906#        aList_summary.append(aString)
1907
1908        i_soln = i_soln + 1
1909
1910    #########################################
1911    #### END OF LOOP OVER MULTI-SOLUTIONS ###
1912    #########################################
1913
1914    #end profiling
1915    pr.disable()
1916    s = StringIO.StringIO()
1917    sortby = 'tottime'
1918    ps = pstats.Stats(pr, stream=s).strip_dirs().sort_stats(sortby)
1919    print('Profiler of function calculation; top 50% of routines:')
1920    ps.print_stats(.5)
1921    print(s.getvalue())
1922    print('%s%.3f'%('Run time = ',time.time()-time0))
1923
1924    localtime = time.asctime( time.localtime(time.time()) )
1925
1926    aString = '\nCompletion time: ' + str(localtime)
1927    print (aString)
1928
1929#    aList_summary.append(str(localtime))
1930#
1931#    # Create summary
1932#
1933#    aFile_log = prefix + 'shapes_summary.log'
1934#    num_lines = len(aList_summary)
1935#
1936#    file = open(aFile_log,'w')
1937#    i = 0
1938#    while i < num_lines:
1939#        aString = str(aList_summary[i])
1940#        file.write(aString)
1941#        file.write('\n')
1942#        i = i + 1
1943#    file.close()
1944
1945    return Phases,Patterns,PRcalc
1946
1947
1948
Note: See TracBrowser for help on using the repository browser.