#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Thu Nov 14 14:38:45 2019 @author: fujisawa """ # # Copyright v1.0, 1.2, 1.3: 2019, John Badger. # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program. If not, see . # # Version 1.2 is intended to be runnable under Python2 and Python3 # Version 1.3 includes an optional 'glue' term for extended structures import math import sys import os import random import time ########## FUNCTIONS ######################## # NEW 1.1 Calculate intensity from P(r) dropping all scaling factors def ft_to_intensity(aList_q,aList_i_calc,aList_r,aList_pr_model,nbeads): num_q = len(aList_q) num_r = len(aList_r) count_q = 0 while count_q < num_q: q = float(aList_q[count_q]) # Sets P(r=0) =0.0. Later scaling includes a baseline term. integral = 0.0 count_r = 1 while count_r < num_r: r = float(aList_r[count_r]) pr = float(aList_pr_model[count_r]) qr = q*r integral = integral + pr*math.sin(qr)/qr count_r = count_r + 1 aList_i_calc.append(integral) count_q = count_q + 1 return; # NEW 1.1 Scale and Compare I and Ic. Includes a baseline correction def score_Ic(aList_q,aList_i,aList_i_sd,aList_i_calc): num_q = len(aList_q) idif = 0.0 isum = 0.0 sd_sq = 0.0 idif_sq = 0.0 chi_sq = 0.0 # Least squares scale for calculated I onto observed I S = 0.0 Sx = 0.0 Sy = 0.0 Sxx = 0.0 Sxy = 0.0 count = 0 while count < num_q: x = float(aList_i_calc[count]) y = float(aList_i[count]) S = S + 1.0 Sx = Sx + x Sy = Sy + y Sxx = Sxx + x*x Sxy = Sxy + x*y count = count + 1 delta = S*Sxx - Sx*Sx a = (Sxx*Sy - Sx*Sxy)/delta b = (S*Sxy - Sx*Sy)/delta # Apply scale and find statistics i = 0 while i < num_q: q = float(aList_q[i]) iobs = float(aList_i[i]) sd = float(aList_i_sd[i]) aList_i_calc[i] = b*float(aList_i_calc[i]) + a icalc = aList_i_calc[i] idif = idif + abs(iobs - icalc) isum = isum + iobs + icalc dif = iobs - icalc dif_sq = dif*dif sd_sq = sd*sd chi_sq = chi_sq + dif_sq/sd_sq i = i + 1 rvalue = 2.0*idif/isum chi_sq = chi_sq/num_q return (chi_sq,rvalue); # NEW 1.1 Write original and calculated data. def write_all_data(file_intensity,aList_q,aList_i,aList_i_calc,aString): num_q = len(aList_q) file = open(file_intensity,'w',) aString = '# ' + aString + '\n' file.write(aString) i = 0 while i < num_q: q = aList_q[i] intensity = aList_i[i] intensity_calc = aList_i_calc[i] aString = str(q) + ' ' + str(intensity) + ' ' + str(intensity_calc) + '\n' file.write(aString) i = i + 1 file.close() return; # NEW 1.1 Read intensity data from GNOM output file def read_i(aList_q,aList_i,aList_i_sd,aList_i_reg,inFile,angstrom_scale): scale_units = 1.0/angstrom_scale file = open(inFile,'r') allLines = file.readlines() file.close() parse_data = 'no' for eachLine in allLines: if parse_data == 'yes': aList = eachLine.split() num_params = len(aList) if num_params == 5: q = float(aList[0]) * scale_units if q > 0.0: i = float(aList[1]) i_sd = float(aList[2]) i_reg = float(aList[3]) aList_q.append(q) aList_i.append(i) aList_i_sd.append(i_sd) aList_i_reg.append(i_reg) else: if num_params == 0 and len(aList_q) > 0: parse_data = 'no' if eachLine.find('S J EXP ERROR J REG I REG') > -1: parse_data = 'yes' return; # Read PDB for starting point def read_pdb(aList_beads_x,aList_beads_y,aList_beads_z,pdbfile_in): xmean = 0.0 ymean = 0.0 zmean = 0.0 nbeads = 0 file = open(pdbfile_in,'r') allLines = file.readlines() file.close() for eachLine in allLines: tag = eachLine[0:6] tag = tag.strip() if tag == 'ATOM': atom_name = eachLine[13:16] atom_name = atom_name.strip() if atom_name == 'CA': x = float(eachLine[30:38]) y = float(eachLine[38:46]) z = float(eachLine[46:54]) xmean = xmean + x ymean = ymean + y zmean = zmean + z nbeads = nbeads + 1 xmean = xmean/float(nbeads) ymean = ymean/float(nbeads) zmean = zmean/float(nbeads) for eachLine in allLines: tag = eachLine[0:6] tag = tag.strip() if tag == 'ATOM': atom_name = eachLine[13:16] atom_name = atom_name.strip() if atom_name == 'CA': x = float(eachLine[30:38]) - xmean y = float(eachLine[38:46]) - ymean z = float(eachLine[46:54]) - zmean aList_beads_x.append(x) aList_beads_y.append(y) aList_beads_z.append(z) return; # Write P(r) with SD and calculated value. def pr_writer(aList_pr,aList_r,aList_pr_model,outfile_pr): num_pr = len(aList_pr) file = open(outfile_pr,'w') file.write('#\n') i = 0 while i < num_pr: r = float(aList_r[i]) pr = float(aList_pr[i]) pr_calc = float(aList_pr_model[i]) aString = str(r) + ' ' + str(pr) + ' ' + str(pr_calc) + '\n' file.write(aString) i = i + 1 file.close() return; # Write a set of points as a pseudo-PDB file def pdb_writer(aList_x_write,aList_y_write,aList_z_write,out_file,aString): atom_number = 0 res_number = 0 dummy_b = 0 num_atoms = len(aList_x_write) file = open(out_file,'w') file.write(aString) file.write('\n') i = 0 while i < num_atoms: x = float(aList_x_write[i]) y = float(aList_y_write[i]) z = float(aList_z_write[i]) x = '%.3f'%(x) y = '%.3f'%(y) z = '%.3f'%(z) x = x.rjust(8) y = y.rjust(8) z = z.rjust(8) atom_number = atom_number + 1 res_number = res_number + 1 atom_number_str = str(atom_number) atom_number_str = atom_number_str.rjust(5) res_number_str = str(res_number) res_number_str = res_number_str.rjust(4) bfactor = str(dummy_b) + '.00' bfactor = bfactor.rjust(6) if res_number < 10000: aLine_out = 'ATOM ' + atom_number_str + ' CA ALA A' + res_number_str + ' ' + x + y + z + ' 1.00' + bfactor + '\n' else: res_number_str = str(res_number - 9999) aLine_out = 'ATOM ' + atom_number_str + ' CA ALA B' + res_number_str + ' ' + x + y + z + ' 1.00' + bfactor + '\n' file.write(aLine_out) i = i + 1 file.close() return; # Evaluate local bead densities and point density on a notional grid def set_box(aList_beads_x,aList_beads_y,aList_beads_z,\ aList_box_x_all,aList_box_y_all,aList_box_z_all,\ aList_box_score,box_step,dmax,rsearch): dmax_over2 = dmax/2.0 search_sq = rsearch**2 num_beads = len(aList_beads_x) count_x = -dmax_over2 while count_x < dmax_over2: count_y = -dmax_over2 while count_y < dmax_over2: count_z = -dmax_over2 while count_z < dmax_over2: count_local = 0 i = 0 while i < num_beads: dx = float(aList_beads_x[i]) - count_x dy = float(aList_beads_y[i]) - count_y dz = float(aList_beads_z[i]) - count_z dsq = dx*dx + dy*dy + dz*dz if dsq < search_sq: count_local = count_local + 1 i = i + 1 if count_local > 1: aList_box_x_all.append(count_x) aList_box_y_all.append(count_y) aList_box_z_all.append(count_z) aList_box_score.append(count_local) count_z = count_z + box_step count_y = count_y + box_step count_x = count_x + box_step return; # Establish a volume def set_vol(aList_box_x_all,aList_box_y_all,aList_box_z_all,aList_box_score,\ aList_box_x,aList_box_y,aList_box_z,vol_target,box_pt_vol): num_pts = len(aList_box_score) num_tries = int(max(aList_box_score)) density_thresh = max(aList_box_score) - 1.0 vol = vol_target + 1.0 i = 0 while i < num_tries: density_thresh = density_thresh - 1.0 num_box_pts = 0.0 j = 0 while j < num_pts: density = float(aList_box_score[j]) if density >= density_thresh: num_box_pts = num_box_pts + 1.0 j = j + 1 vol = num_box_pts*box_pt_vol if vol < vol_target: density_use = density_thresh i = i + 1 # num_box_pts1 = 0.0 num_box_pts2 = 0.0 density_thresh1 = density_use density_thresh2 = density_use - 1.0 i = 0 while i < num_pts: density_value = float(aList_box_score[i]) if density_value >= density_thresh1: num_box_pts1 = num_box_pts1 + 1.0 if density_value >= density_thresh2: num_box_pts2 = num_box_pts2 + 1.0 i = i + 1 vol1 = num_box_pts1*box_pt_vol vol2 = num_box_pts2*box_pt_vol delta1 = abs(vol1-vol_target) delta2 = abs(vol2-vol_target) if delta1 < delta2: density_thresh = density_thresh1 else: density_thresh = density_thresh2 i = 0 while i < num_pts: density_value = float(aList_box_score[i]) if density_value >= density_thresh: aList_box_x.append(aList_box_x_all[i]) aList_box_y.append(aList_box_y_all[i]) aList_box_z.append(aList_box_z_all[i]) i = i + 1 return; # Find beads that are not in allowed volume def disallowed_beads(aList_beads_x,aList_beads_y,aList_beads_z,aList_contacts,\ aList_box_x,aList_box_y,aList_box_z,rsearch): num_beads = len(aList_beads_x) num_boxes = len(aList_box_x) contact_limit_sq = rsearch**2 count = 0 while count < num_beads: x_bead = float(aList_beads_x[count]) y_bead = float(aList_beads_y[count]) z_bead = float(aList_beads_z[count]) inbox = 'no' i = 0 while i < num_boxes: x_box = float(aList_box_x[i]) y_box = float(aList_box_y[i]) z_box = float(aList_box_z[i]) dsq = (x_bead - x_box)**2 + (y_bead - y_box)**2 + (z_bead - z_box)**2 if dsq < contact_limit_sq: inbox = 'yes' i = num_boxes i = i + 1 if inbox == 'no': aList_contacts.append(count) count = count + 1 return; # Compute a P(r) def calc_pr(aList_beads_x,aList_beads_y,aList_beads_z,aList_pr_model,hist_grid): num_hist = len(aList_pr_model) count = 0 while count < num_hist: aList_pr_model[count] = 0.0 count = count + 1 nbeads = len(aList_beads_x) max_dr = (float(num_hist)-1.0)*hist_grid i = 0 while i < nbeads: j = 0 while j < i: dr = get_dr(aList_beads_x[i],aList_beads_y[i],aList_beads_z[i],\ aList_beads_x[j],aList_beads_y[j],aList_beads_z[j]) dr = min(dr,max_dr) # Find pointers and do un-interpolation dr_grid = dr/hist_grid int_dr_grid = int(dr_grid) int_dr_grid = min(int_dr_grid,num_hist-2) ip_low = int_dr_grid ip_high = ip_low + 1 ip_high_frac = dr_grid - float(int_dr_grid) ip_low_frac = 1.0 - ip_high_frac aList_pr_model[ip_low] = float(aList_pr_model[ip_low]) + ip_low_frac aList_pr_model[ip_high] = float(aList_pr_model[ip_high]) + ip_high_frac j = j + 1 i = i + 1 return; # Score for rms difference between observed and model histograms def pr_dif(aList_pr,aList_pr_model,skip): num_hist = len(aList_pr) delta_hist_sum = 0.0 delta_hist_sum_sq = 0.0 hist_sum = 0.0 i = skip while i < num_hist: model = float(aList_pr_model[i]) data = float(aList_pr[i]) delta_hist = abs(data - model) delta_hist_sum = delta_hist_sum + delta_hist hist_sum = hist_sum + data delta_hist_sum_sq = delta_hist_sum_sq + delta_hist*delta_hist i = i + 1 mean_hist_sum = hist_sum/(num_hist - skip) delta_hist_sum_sq = delta_hist_sum_sq/(num_hist - skip) delta_hist_sum_sq = math.sqrt(delta_hist_sum_sq)/mean_hist_sum return delta_hist_sum_sq; # Statistics for fractional difference between observed and model histograms def pr_rfactor(aList_pr,aList_pr_sd,aList_pr_model,skip): num_hist = len(aList_pr) delta_hist_sum = 0.0 hist_sum = 0.0 i = skip while i < num_hist: model = float(aList_pr_model[i]) exp = float(aList_pr[i]) delta_hist = exp - model delta_hist_sum = delta_hist_sum + abs(delta_hist) hist_sum = hist_sum + exp i = i + 1 delta_hist_sum = delta_hist_sum/hist_sum return delta_hist_sum; # Compute the VDW energy for a interaction def vdw_energy(econ12,econ6,e_width,dr,bead_sep3): if dr > bead_sep3: vdw = 0.0 else: dr_e6 = dr**6 dr_e12 = dr_e6**2 vdw = econ12/dr_e12 - 2.0*econ6/dr_e6 vdw = max(vdw,e_width) return vdw; # Set a random distribution of beads in a box with maximum extent dmax def random_beads(aList_beads_x,aList_beads_y,aList_beads_z,\ nbeads,dmax,aList_symm,bias_z,num_symm): half_side = 0.5*dmax half_side_sq = half_side*half_side scale_xy = 1.0 - bias_z scale_z = 1.0 + bias_z x_range = scale_xy * half_side y_range = scale_xy * half_side z_range = scale_z * half_side num_ops = len(aList_symm) i = 0 while i < nbeads: triangle = random.triangular(-0.7,0.7,0.0) x = triangle*x_range triangle = random.triangular(-0.7,0.7,0.0) y = triangle*y_range triangle = random.triangular(-0.7,0.7,0.0) z = triangle*z_range aList_beads_x.append(x) aList_beads_y.append(y) aList_beads_z.append(z) j = 0 while j < num_ops: aList_s = aList_symm[j] m11 = float(aList_s[0]) m12 = float(aList_s[1]) m21 = float(aList_s[2]) m22 = float(aList_s[3]) xs = m11*x + m12*y ys = m21*x + m22*y zs = z aList_beads_x.append(xs) aList_beads_y.append(ys) aList_beads_z.append(zs) j = j + 1 i = i + num_symm return; # Read experimentalal P(r) from GNOM output file def read_pr(aList_r,aList_pr,aList_pr_sd,aList_pr_model,\ aList_pr_model_test,aList_pr_model_test2,inFile): angstrom_scale = 1.0 file = open(inFile,'r') allLines = file.readlines() file.close() parse_data = 'no' for eachLine in allLines: if parse_data == 'yes': aList = eachLine.split() num_params = len(aList) if num_params == 3: r = float(aList[0]) pr = float(aList[1]) pr_sd = float(aList[2]) aList_pr.append(pr) aList_pr_sd.append(pr_sd) aList_r.append(r) aList_pr_model.append(0.0) aList_pr_model_test.append(0.0) aList_pr_model_test2.append(0.0) if eachLine.find('R P(R) ERROR') > -1: parse_data = 'yes' num_hist = len(aList_r) hist_grid = float(aList_r[1]) - float(aList_r[0]) # Heuristic for checking units test_r = max(aList_r) if test_r < 30.0: aString = 'P(r)appears to be in nm. Converting to Angstrom units' print (aString) angstrom_scale = 10.0 i = 0 while i < num_hist: aList_r[i] = angstrom_scale * aList_r[i] i = i + 1 hist_grid = angstrom_scale * hist_grid i = 0 while i < num_hist: r = float(aList_r[i]) r_calc = float(i)*hist_grid if abs(r - r_calc) > 0.15: aString = 'Input P(r) grid is irregular! Exiting' print (aString) time.sleep(4) sys.exit(1) i = i + 1 dmax = aList_r[num_hist-1] # Pad histogram by 5 Angstrom ipad = int(5.0/hist_grid) i = 0 while i < ipad: r = dmax + float(i)*hist_grid aList_pr.append(0.0) aList_pr_sd.append(0.0) aList_r.append(r) aList_pr_model.append(0.0) aList_pr_model_test.append(0.0) aList_pr_model_test2.append(0.0) i = i + 1 return (dmax,hist_grid,num_hist,angstrom_scale); # Scale P(r) onto model P(r) assuming same grid def scale_pr(aList_pr,aList_pr_sd,aList_pr_model): num_hist = len(aList_pr) total_dr = 0.0 total_pr = 0.0 i = 0 while i < num_hist: total_dr = total_dr + float(aList_pr_model[i]) total_pr = total_pr + float(aList_pr[i]) i = i + 1 scale = total_dr/total_pr i = 0 while i < num_hist: aList_pr[i] = scale*float(aList_pr[i]) aList_pr_sd[i] = scale * float(aList_pr_sd[i]) i = i + 1 return; # Return a non-zero distance between two coordinates def get_dr(x1,y1,z1,x2,y2,z2): x1 = float(x1) y1 = float(y1) z1 = float(z1) x2 = float(x2) y2 = float(y2) z2 = float(z2) dr = (x1 - x2)**2 + (y1-y2)**2 + (z1-z2)**2 dr = max(dr,0.1) dr = math.sqrt(dr) return dr; # Find center of beads within a radii def center_beads(x,y,z,aList_beads_x,aList_beads_y,aList_beads_z,radii_1,radii_2): num_beads = len(aList_beads_x) xsum = 0.0 ysum = 0.0 zsum = 0.0 count_beads = 0.0 i = 0 while i < num_beads: dr = get_dr(x,y,z,aList_beads_x[i],aList_beads_y[i],aList_beads_z[i]) if dr > radii_1 and dr < radii_2: count_beads = count_beads + 1.0 xsum = xsum + float(aList_beads_x[i]) ysum = ysum + float(aList_beads_y[i]) zsum = zsum + float(aList_beads_z[i]) i = i + 1 if count_beads > 0.0: xsum = xsum/count_beads ysum = ysum/count_beads zsum = zsum/count_beads delta = (xsum - x)**2 + (ysum - y)**2 + (zsum - z)**2 delta = math.sqrt(delta) else: delta = 0.0 return delta; # Obtain mean of total VDW energy def get_total_energy(aList_beads_x,aList_beads_y,aList_beads_z,econ12,econ6,bead_sep3,e_width): nbeads = len(aList_beads_x) vdw_all = 0.0 i = 0 while i < nbeads: j = 0 while j < i: dr = get_dr(aList_beads_x[i],aList_beads_y[i],aList_beads_z[i],\ aList_beads_x[j],aList_beads_y[j],aList_beads_z[j]) vdw = vdw_energy(econ12,econ6,e_width,dr,bead_sep3) vdw_all = vdw_all + vdw j = j + 1 i = i + 1 vdw_all = vdw_all/float(nbeads) return vdw_all; # Energy minimize def e_min(aList_beads_x,aList_beads_y,aList_beads_z,econ12,econ6,bead_sep,bead_sep3,aList_symm,e_width): eps = bead_sep/(2.0**(1.0/6.0)) eps12 = eps**12 eps6 = eps**6 step_max = bead_sep scale = 0.0 icount = -1 nbeads = len(aList_beads_x) num_ops = len(aList_symm) num_cycles = 51 i = 0 while i < num_cycles: icount = icount + 1 aList_beads_x_new = [] aList_beads_y_new = [] aList_beads_z_new = [] num_forces = 0.0 sum_forces = 0.0 sum_forces_scale = 0.0 k = 0 while k < nbeads - num_ops: xold = float(aList_beads_x[k]) yold = float(aList_beads_y[k]) zold = float(aList_beads_z[k]) fx = 0.0 fy = 0.0 fz = 0.0 j = 0 while j < nbeads: xj = aList_beads_x[j] yj = aList_beads_y[j] zj = aList_beads_z[j] dr = get_dr(xold,yold,zold,xj,yj,zj) # Truncate very steep dr = min(eps,dr) if dr < bead_sep3: dr_sq = dr*dr dr12 = dr_sq**6 dr6 = dr_sq**3 dx = xold - xj dy = yold - yj dz = zold - zj force = (1.0/dr_sq)*(eps12/dr12 - 0.5*eps6/dr6) fx = fx + force*dx fy = fy + force*dy fz = fz + force*dz sum_forces_scale = sum_forces_scale + abs(fx) + abs(fy) + abs(fz) j = j + 1 # xstep = scale*fx ystep = scale*fy zstep = scale*fz if xstep > 0.0: xstep = min(xstep,step_max) else: xstep = max(xstep,-step_max) if ystep > 0.0: ystep = min(ystep,step_max) else: ystep = max(ystep,-step_max) if zstep > 0.0: zstep = min(zstep,step_max) else: zstep = max(zstep,-step_max) xtest = xold + xstep ytest = yold + ystep ztest = zold + zstep aList_beads_x_new.append(xtest) aList_beads_y_new.append(ytest) aList_beads_z_new.append(ztest) # Apply shifs to symm positions l = 0 while l < num_ops: aList_s = aList_symm[l] m11 = float(aList_s[0]) m12 = float(aList_s[1]) m21 = float(aList_s[2]) m22 = float(aList_s[3]) xs = m11*xtest + m12*ytest ys = m21*xtest + m22*ytest zs = ztest aList_beads_x_new.append(xs) aList_beads_y_new.append(ys) aList_beads_z_new.append(zs) l = l + 1 # k = k + num_ops + 1 # Apply shifted positions after first cycle if i > 0: m = 0 while m < nbeads: aList_beads_x[m] = aList_beads_x_new[m] aList_beads_y[m] = aList_beads_y_new[m] aList_beads_z[m] = aList_beads_z_new[m] m = m + 1 # mean_force = (num_ops+1)*sum_forces_scale/(nbeads*3.0) scale = bead_sep/mean_force vdw_all = get_total_energy(aList_beads_x,aList_beads_y,aList_beads_z,econ12,econ6,bead_sep3,e_width) if icount == 0: aString = 'Emin cycle: ' + str(i) + ' Energy: ' + str('%4.2f'%(vdw_all)) print (aString) icount = -10 if vdw_all < 0.0: i = num_cycles i = i + 1 return; # Set up symmetry operators for rotational symmetry def make_symm(aList_symm,num_symm): angle_step = 360.0/float(num_symm) i = 1 while i < num_symm: theta = float(i) * angle_step theta = math.radians(theta) cos_theta = math.cos(theta) sin_theta = math.sin(theta) aList_s = [cos_theta,sin_theta,-sin_theta,cos_theta] aList_symm.append(aList_s) i = i + 1 return aList_symm; # Set up a shift vector in P(r) for a change in bead position def pr_shift_atom(aList_r,aList_pr_model_test2,x1,y1,z1,\ aList_beads_x,aList_beads_y,aList_beads_z,hist_grid,ii): num_hist = len(aList_r) max_dr = (float(num_hist)-1.0)*hist_grid num_beads = len(aList_beads_x) i = 0 while i < num_hist: aList_pr_model_test2[i] = 0.0 i = i + 1 i = 0 while i < num_beads: if i != ii: x2 = float(aList_beads_x[i]) y2 = float(aList_beads_y[i]) z2 = float(aList_beads_z[i]) dr = get_dr(x1,y1,z1,x2,y2,z2) dr = min(dr,max_dr) dr_grid = dr/hist_grid int_dr_grid = int(dr_grid) int_dr_grid = max(int_dr_grid,0) int_dr_grid = min(int_dr_grid,num_hist-2) ip_low = int_dr_grid ip_high = ip_low + 1 ip_high_frac = dr_grid - float(int_dr_grid) ip_low_frac = 1.0 - ip_high_frac aList_pr_model_test2[ip_low] = float(aList_pr_model_test2[ip_low]) + ip_low_frac aList_pr_model_test2[ip_high] = float(aList_pr_model_test2[ip_high]) + ip_high_frac i = i + 1 return; # Recenter set of beads to origin def recenter_pdb(aList_beads_x,aList_beads_y,aList_beads_z): nbeads = len(aList_beads_x) xsum = 0.0 ysum = 0.0 zsum = 0.0 i = 0 while i < nbeads: xsum = xsum + float(aList_beads_x[i]) ysum = ysum + float(aList_beads_y[i]) zsum = zsum + float(aList_beads_z[i]) i = i + 1 xmean = xsum/float(nbeads) ymean = ysum/float(nbeads) zmean = zsum/float(nbeads) i = 0 while i < nbeads: aList_beads_x[i] = float(aList_beads_x[i]) - xmean aList_beads_y[i] = float(aList_beads_y[i]) - ymean aList_beads_z[i] = float(aList_beads_z[i]) - zmean i = i + 1 return;