Changeset 4101 for trunk/G2shapes.py
- Timestamp:
- Aug 19, 2019 2:28:38 PM (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/G2shapes.py
r4099 r4101 31 31 import io as StringIO 32 32 import numpy as np 33 nxs = np.newaxis 33 34 34 def G2shapes(Profile,ProfDict,Limits,data ):35 def G2shapes(Profile,ProfDict,Limits,data,dlg=None): 35 36 36 37 ########## FUNCTIONS ######################## … … 396 397 for XYZi in XYZ: 397 398 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))) 399 400 if count>1: 400 401 aList_box_x_all.append(XYZi[0]) … … 617 618 return vdw; 618 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 619 629 # Set a random distribution of beads in a box with maximum extent dmax 620 630 … … 858 868 vdw_all = 0.0 859 869 870 860 871 i = 0 861 872 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 870 890 871 891 vdw_all = vdw_all/float(nbeads) … … 906 926 zold = float(aList_beads_z[k]) 907 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 908 948 fx = 0.0 909 949 fy = 0.0 910 950 fz = 0.0 951 911 952 j = 0 912 953 while j < nbeads: … … 1045 1086 # num_beads = len(aList_beads_x) 1046 1087 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 1049 1095 XYZ = np.array([aList_beads_x,aList_beads_y,aList_beads_z]).T 1050 1096 xyz = np.array([x1,y1,z1]) … … 1062 1108 aList_pr_model_test2[ip_high] += ip_high_fracs[ip_high] 1063 1109 1064 # i = 0 1065 # while i < num_hist: 1066 # aList_pr_model_test2[i] = 0.0 1067 # i = i + 1 1068 # 1110 1069 1111 # i = 0 1070 1112 # while i < num_beads: … … 1672 1714 # Get initial contributions to P(r) 1673 1715 1716 aList_pr_model_test2 = num_hist*[0.0,] 1674 1717 pr_shift_atom(aList_pr_model_test2,xold,yold,zold,aList_beads_x,\ 1675 1718 aList_beads_y,aList_beads_z,hist_grid,ii) … … 1713 1756 # Get shifted contributions to P(r) 1714 1757 1758 aList_pr_model_test2 = num_hist*[0.0,] 1715 1759 pr_shift_atom(aList_pr_model_test2,xtest,ytest,ztest,aList_beads_x,\ 1716 1760 aList_beads_y,aList_beads_z,hist_grid,ii) … … 1806 1850 success_rate = success*float(num_symm)/float(nbeads) 1807 1851 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') 1813 1860 1814 1861 # Debug statitics. Weight of 10 gives about 1.0
Note: See TracChangeset
for help on using the changeset viewer.