source: trunk/G2shapes.py @ 4084

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

fixes to SHAPES results table
fix remaining plot amnesia problems

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