source: trunk/G2shapes.py @ 4106

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

fixes so plots are correct for SHAPES fitting

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