source: trunk/G2shapes.py @ 4099

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

more speedup of SHAPES

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