Changeset 4101


Ignore:
Timestamp:
Aug 19, 2019 2:28:38 PM (2 years ago)
Author:
vondreele
Message:

progress dilog for SHAPES & more speedup

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/G2shapes.py

    r4099 r4101  
    3131import io as StringIO
    3232import numpy as np
     33nxs = np.newaxis
    3334
    34 def G2shapes(Profile,ProfDict,Limits,data):   
     35def G2shapes(Profile,ProfDict,Limits,data,dlg=None):   
    3536
    3637########## FUNCTIONS ########################
     
    396397        for XYZi in XYZ:
    397398            dsq = np.sum((xyz-XYZi)**2,axis=1)
    398             count = int(np.sum(np.array([1 for dist in dsq if dist < search_sq])))
     399            count = int(np.sum(np.where(dsq<search_sq,1,0)))
    399400            if count>1:
    400401                aList_box_x_all.append(XYZi[0])
     
    617618        return vdw;
    618619   
     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       
    619629    # Set a random distribution of beads in a box with maximum extent dmax
    620630   
     
    858868        vdw_all = 0.0
    859869       
     870       
    860871        i = 0
    861872        while i < nbeads:
    862             j = 0
    863             while j < i:
    864                 dr = get_dr(aList_beads_x[i],aList_beads_y[i],aList_beads_z[i],\
    865                             aList_beads_x[j],aList_beads_y[j],aList_beads_z[j])
    866                 vdw = vdw_energy(econ12,econ6,e_width,dr,bead_sep3)           
    867                 vdw_all = vdw_all + vdw
    868                 j = j + 1
    869             i = i + 1
     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
    870890           
    871891        vdw_all = vdw_all/float(nbeads)
     
    906926                zold = float(aList_beads_z[k])
    907927   
     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               
    908948                fx = 0.0
    909949                fy = 0.0
    910950                fz = 0.0
     951               
    911952                j = 0
    912953                while j < nbeads:
     
    10451086#        num_beads = len(aList_beads_x)
    10461087       
    1047         aList_pr_model_test2 = num_hist*[0.0,]
    1048        
     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           
    10491095        XYZ = np.array([aList_beads_x,aList_beads_y,aList_beads_z]).T
    10501096        xyz = np.array([x1,y1,z1])
     
    10621108            aList_pr_model_test2[ip_high] += ip_high_fracs[ip_high]
    10631109   
    1064 #        i = 0
    1065 #        while i < num_hist:
    1066 #            aList_pr_model_test2[i] = 0.0
    1067 #            i = i + 1
    1068 #   
     1110   
    10691111#        i = 0
    10701112#        while i < num_beads:
     
    16721714                # Get initial contributions to P(r) 
    16731715   
     1716                aList_pr_model_test2 = num_hist*[0.0,]
    16741717                pr_shift_atom(aList_pr_model_test2,xold,yold,zold,aList_beads_x,\
    16751718                              aList_beads_y,aList_beads_z,hist_grid,ii)
     
    17131756                # Get shifted contributions to P(r) 
    17141757   
     1758                aList_pr_model_test2 = num_hist*[0.0,]
    17151759                pr_shift_atom(aList_pr_model_test2,xtest,ytest,ztest,aList_beads_x,\
    17161760                              aList_beads_y,aList_beads_z,hist_grid,ii)
     
    18061850            success_rate = success*float(num_symm)/float(nbeads)
    18071851            success_rate_all = success_rate_all + success_rate
    1808    
    1809             aString = 'Cycle ' + str(count_it+1) + ' Moves ' + str('%.2f'%(success_rate)) + \
    1810                       ' Possibles ' + str('%.2f'%(count_hist_yes)) + ' rms P(r) '+ str('%4.3f'%(hist_score)) + \
    1811                       ' Energy ' + str('%4.2f'%(vdw_all))
    1812             print (aString)
     1852
     1853            if not (count_it+1)%10:   
     1854                aString = 'Cycle ' + str(count_it+1) + ' Moves ' + str('%.2f'%(success_rate)) + \
     1855                          ' Possibles ' + str('%.2f'%(count_hist_yes)) + ' rms P(r) '+ str('%4.3f'%(hist_score)) + \
     1856                          ' Energy ' + str('%4.2f'%(vdw_all))
     1857                print (aString)
     1858            if dlg:
     1859                dlg.Update(count_it+1,newmsg='Cycle no.: '+str(count_it)+' of 160')
    18131860   
    18141861            # Debug statitics. Weight of 10 gives about 1.0
  • trunk/GSASIIpwdGUI.py

    r4097 r4101  
    54375437      doi: 10.1107/S1600576719009774''',
    54385438      caption='Program Shapes',style=wx.ICON_INFORMATION)
     5439            dlg = wx.ProgressDialog('Running SHAPES','Cycle no.: 0 of 160',161,
     5440                style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME)
     5441               
    54395442            data['Pair']['Result'] = []       #clear old results (if any) for now
    5440             data['Pair']['Result'] = G2shapes.G2shapes(Profile,ProfDict,Limits,data)
     5443            data['Pair']['Result'] = G2shapes.G2shapes(Profile,ProfDict,Limits,data,dlg)
    54415444            wx.CallAfter(UpdateModelsGrid,G2frame,data)
    54425445           
Note: See TracChangeset for help on using the changeset viewer.