source: trunk/G2shapes.py @ 4079

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

modify SHAPES for better output to GSAS-II
begin plotting capability for SHAPES output

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