From f81adf1e08c290bf2c0b7fb128e3d2aacc591432 Mon Sep 17 00:00:00 2001 From: "andrea.grisafi" Date: Thu, 18 Feb 2021 13:08:08 +0100 Subject: [PATCH 1/2] Centers sorting option included --- soapfast/get_power_spectrum.py | 56 ++++++++++++++----------- soapfast/utils/LODE/PS_utils.py | 12 +++--- soapfast/utils/LODE/direct_ewald.py | 4 +- soapfast/utils/LODE/direct_potential.py | 4 +- soapfast/utils/LODE/fourier_ewald.py | 2 +- soapfast/utils/LODE/parsing.py | 4 +- soapfast/utils/parsing.py | 4 +- 7 files changed, 48 insertions(+), 38 deletions(-) diff --git a/soapfast/get_power_spectrum.py b/soapfast/get_power_spectrum.py index 0f29108..49d0ce6 100755 --- a/soapfast/get_power_spectrum.py +++ b/soapfast/get_power_spectrum.py @@ -13,7 +13,7 @@ ############################################################################################################################### -def get_power_spectrum(lam,frames,nmax=8,lmax=6,rc=4.0,sg=0.3,ncut=-1,cw=1.0,periodic=False,outfile='',subset=['NO',None],initial=-1,sparsefile='',sparse_options=[''],cen=[],spec=[],atomic=[False,None],all_radial=[1.0,0.0,0.0],useall=False,xyz_slice=[],verbose=True,get_imag=False,norm=True,electro=False,sigewald=1.0,single_radial=False,radsize=50,lebsize=146,average=False): +def get_power_spectrum(lam,frames,nmax=8,lmax=6,rc=4.0,sg=0.3,ncut=-1,cw=1.0,periodic=False,outfile='',subset=['NO',None],initial=-1,sparsefile='',sparse_options=[''],cen=[],spec=[],atomic=[False,None],all_radial=[1.0,0.0,0.0],useall=False,xyz_slice=[],verbose=True,get_imag=False,norm=True,electro=False,sigewald=1.0,single_radial=False,radsize=50,lebsize=146,average=False,censort=False): # If we want a slice of the coordinates, do this BEFORE anything else. if (len(xyz_slice)!=0): @@ -230,26 +230,32 @@ def get_power_spectrum(lam,frames,nmax=8,lmax=6,rc=4.0,sg=0.3,ncut=-1,cw=1.0,per if not average: - # Reorder the power spectrum so that the ordering of the atoms matches their positions in the frame. - for i in range(npoints): - if (lam==0): - ps_row = np.zeros((len(power[0]),len(power[0,0])),dtype=float) - else: - ps_row = np.zeros((len(power[0]),len(power[0,0]),len(power[0,0,0])),dtype=float) - # Reorder this row - reorder_list = [[] for cen in centers] - atnum = frames[i].get_atomic_numbers() - for j in range(len(atnum)): - atom = atnum[j] - if (atom in centers): - place = list(centers).index(atom) - reorder_list[place].append(j) - reorder_list = sum(reorder_list,[]) - # The reordering list tells us where each column of the current power spectrum should go. - for j in range(len(reorder_list)): - ps_row[reorder_list[j]] = power[i,j] - # Insert the reordered row back into the power spectrum - power[i] = ps_row + + if not censort: + + # Reorder the power spectrum so that the ordering of the atoms matches their positions in the frame. + for i in range(npoints): + if (lam==0): + ps_row = np.zeros((len(power[0]),len(power[0,0])),dtype=float) + else: + ps_row = np.zeros((len(power[0]),len(power[0,0]),len(power[0,0,0])),dtype=float) + # Reorder this row + reorder_list = [[] for cen in centers] + atnum = frames[i].get_atomic_numbers() + for j in range(len(atnum)): + atom = atnum[j] + if (atom in centers): + place = list(centers).index(atom) + reorder_list[place].append(j) + reorder_list = sum(reorder_list,[]) + # The reordering list tells us where each column of the current power spectrum should go. + for j in range(len(reorder_list)): + ps_row[reorder_list[j]] = power[i,j] + # Insert the reordered row back into the power spectrum + power[i] = ps_row + else: + + print("Power spectrum entries sorted according to atomic centers list.") # Print power spectrum, if we have not asked for only a sample to be taken (we assume that taking a subset is just for the purpose of generating a sparsification) if (subset[0] == 'NO'): @@ -308,12 +314,12 @@ def main(): pname = os.path.dirname(os.path.realpath(__file__)) if os.path.isfile(pname + "/utils/LODE/gvectors.so"): args = parse.add_command_line_arguments_PS("Calculate power spectrum") - [nmax,lmax,rcut,sig,cen,spec,cweight,lam,periodic,ncut,sparsefile,frames,subset,sparse_options,outfile,initial,atomic,all_radial,useall,xyz_slice,get_imag,nonorm,electro,sigewald,single_radial,radsize,lebsize,average] = parse.set_variable_values_PS(args) - get_power_spectrum(lam,frames,nmax=nmax,lmax=lmax,rc=rcut,sg=sig,ncut=ncut,periodic=periodic,outfile=outfile,cw=cweight,subset=subset,initial=initial,sparsefile=sparsefile,sparse_options=sparse_options,cen=cen,spec=spec,atomic=atomic,all_radial=all_radial,useall=useall,xyz_slice=xyz_slice,get_imag=get_imag,norm=(not nonorm),electro=electro,sigewald=sigewald,single_radial=single_radial,radsize=radsize,lebsize=lebsize,average=average) + [nmax,lmax,rcut,sig,cen,spec,cweight,lam,periodic,ncut,sparsefile,frames,subset,sparse_options,outfile,initial,atomic,all_radial,useall,xyz_slice,get_imag,nonorm,electro,sigewald,single_radial,radsize,lebsize,average,censort] = parse.set_variable_values_PS(args) + get_power_spectrum(lam,frames,nmax=nmax,lmax=lmax,rc=rcut,sg=sig,ncut=ncut,periodic=periodic,outfile=outfile,cw=cweight,subset=subset,initial=initial,sparsefile=sparsefile,sparse_options=sparse_options,cen=cen,spec=spec,atomic=atomic,all_radial=all_radial,useall=useall,xyz_slice=xyz_slice,get_imag=get_imag,norm=(not nonorm),electro=electro,sigewald=sigewald,single_radial=single_radial,radsize=radsize,lebsize=lebsize,average=average,censort=censort) else: args = parse.add_command_line_arguments_PS("Calculate power spectrum") - [nmax,lmax,rcut,sig,cen,spec,cweight,lam,periodic,ncut,sparsefile,frames,subset,sparse_options,outfile,initial,atomic,all_radial,useall,xyz_slice,get_imag,nonorm] = parse.set_variable_values_PS(args) - get_power_spectrum(lam,frames,nmax=nmax,lmax=lmax,rc=rcut,sg=sig,ncut=ncut,periodic=periodic,outfile=outfile,cw=cweight,subset=subset,initial=initial,sparsefile=sparsefile,sparse_options=sparse_options,cen=cen,spec=spec,atomic=atomic,all_radial=all_radial,useall=useall,xyz_slice=xyz_slice,get_imag=get_imag,norm=(not nonorm)) + [nmax,lmax,rcut,sig,cen,spec,cweight,lam,periodic,ncut,sparsefile,frames,subset,sparse_options,outfile,initial,atomic,all_radial,useall,xyz_slice,get_imag,nonorm,censort] = parse.set_variable_values_PS(args) + get_power_spectrum(lam,frames,nmax=nmax,lmax=lmax,rc=rcut,sg=sig,ncut=ncut,periodic=periodic,outfile=outfile,cw=cweight,subset=subset,initial=initial,sparsefile=sparsefile,sparse_options=sparse_options,cen=cen,spec=spec,atomic=atomic,all_radial=all_radial,useall=useall,xyz_slice=xyz_slice,get_imag=get_imag,norm=(not nonorm),censort=censort) if __name__=="__main__": from utils import regression_utils diff --git a/soapfast/utils/LODE/PS_utils.py b/soapfast/utils/LODE/PS_utils.py index af0d559..c15310a 100644 --- a/soapfast/utils/LODE/PS_utils.py +++ b/soapfast/utils/LODE/PS_utils.py @@ -161,9 +161,9 @@ def compute_power_spectrum(nat,nneighmax,natmax,lam,lmax,npoints,nspecies,nnmax, else: - print("") - print("Computing potential projection for frame", i+1) - print("---------------------------------------------") + # print("") + # print("Computing potential projection for frame", i+1) + # print("---------------------------------------------") if np.sum(cell[i])==0: @@ -307,9 +307,9 @@ def compute_power_spectrum(nat,nneighmax,natmax,lam,lmax,npoints,nspecies,nnmax, else: - print("") - print("Computing potential projection for frame", i+1) - print("---------------------------------------------") + # print("") + # print("Computing potential projection for frame", i+1) + # print("---------------------------------------------") if np.sum(cell[i])==0: diff --git a/soapfast/utils/LODE/direct_ewald.py b/soapfast/utils/LODE/direct_ewald.py index afaa454..faf9c03 100644 --- a/soapfast/utils/LODE/direct_ewald.py +++ b/soapfast/utils/LODE/direct_ewald.py @@ -60,8 +60,8 @@ def direct_ewald(sigewald,cell,nat,nnmax,nspecies,lmax,centers,all_species,nneig omega_near = nearfield_ewald.nearfield(nat,nspecies,nmax,lmax,lebsize*radsize,nneigh_near,nnmax,alpha,coordx_near,spherical_grid,orthoradial,harmonics,integration_weights,sigewald) omega_near = np.transpose(omega_near,(4,3,2,1,0)) - print("direct space potential computed in", time.time()-start, "seconds") - print("") + #print("direct space potential computed in", time.time()-start, "seconds") + #print("") return omega_near diff --git a/soapfast/utils/LODE/direct_potential.py b/soapfast/utils/LODE/direct_potential.py index e2a9336..62dc3cd 100644 --- a/soapfast/utils/LODE/direct_potential.py +++ b/soapfast/utils/LODE/direct_potential.py @@ -100,8 +100,8 @@ def direct_potential(nat,nnmax,nspecies,lmax,centers,all_species,nneighmax,atom_ #print "-----------------------------------------" - print("Direct space potential computed in", time.time()-start, "seconds") - print("") + #print("Direct space potential computed in", time.time()-start, "seconds") + #print("") return omega_near diff --git a/soapfast/utils/LODE/fourier_ewald.py b/soapfast/utils/LODE/fourier_ewald.py index d6ae107..8b1dd5a 100644 --- a/soapfast/utils/LODE/fourier_ewald.py +++ b/soapfast/utils/LODE/fourier_ewald.py @@ -57,7 +57,7 @@ def fourier_ewald(sigewald,nat,nnmax,nspecies,lmax,centers,all_species,nneighmax # print "G-contraction :", time.time()-start4, "seconds" # print "-----------------------------------------" - print("reciprocal space potential computed in", time.time()-start, "seconds") + #print("reciprocal space potential computed in", time.time()-start, "seconds") omega *= 16.0*np.pi**2/volume diff --git a/soapfast/utils/LODE/parsing.py b/soapfast/utils/LODE/parsing.py index 9f3d830..682db1a 100644 --- a/soapfast/utils/LODE/parsing.py +++ b/soapfast/utils/LODE/parsing.py @@ -42,6 +42,7 @@ def add_command_line_arguments_PS(parsetext): parser.add_argument("-rad", "--radsize", type=int, default=50, help="Dimension of the Gauss-Legendre grid needed for the numerical radial integration of the direct-space potential") parser.add_argument("-leb", "--lebsize", type=int, default=146, help="Dimension of the Lebedev grid needed for the numerical angular integration of the direct-space potential. Choose among [6, 14, 26, 38, 50, 74, 86, 110, 146, 170, 194, 230, 266, 302, 350, 434, 590, 770, 974, 1202, 1454, 1730, 2030 (army grade), 2354, 2702, 3074, 3470, 3890, 4334, 4802, 5294, 5810]") parser.add_argument("-ave", "--average", action='store_true', help="Structural averaged features?") + parser.add_argument("-cs", "--centersorting", action='store_true', help="Sort power spectrum entries according to atomic centers list.") args = parser.parse_args() return args @@ -55,6 +56,7 @@ def set_variable_values_PS(args): srad = args.sradial ele = args.electro ave = args.average + censort = args.centersorting sew = args.sigewald # Gaussian width nmax = args.nmax # number of radial functions lmax = args.lmax # number of angular functions @@ -132,6 +134,6 @@ def set_variable_values_PS(args): sys.exit(0) xyz_slice = [args.slice[0],args.slice[1]] - return [nmax,lmax,rc,sg,cen,spec,cw,lam,periodic,ncut,sparsefile,frames,subset,sparse_options,outfile,args.initial,atomic,all_radial,not args.uselist,xyz_slice,args.imag,args.nonorm,ele,sew,srad,rad,leb,ave] + return [nmax,lmax,rc,sg,cen,spec,cw,lam,periodic,ncut,sparsefile,frames,subset,sparse_options,outfile,args.initial,atomic,all_radial,not args.uselist,xyz_slice,args.imag,args.nonorm,ele,sew,srad,rad,leb,ave,censort] ######################################################################### diff --git a/soapfast/utils/parsing.py b/soapfast/utils/parsing.py index bbd6902..ca90cf0 100644 --- a/soapfast/utils/parsing.py +++ b/soapfast/utils/parsing.py @@ -177,6 +177,7 @@ def add_command_line_arguments_PS(parsetext): parser.add_argument("-sl", "--slice", type=int, default=-1, nargs='+', help="Choose a slice of the input frames to calculate the power spectrum") parser.add_argument("-im", "--imag", action='store_true', help="Get imaginary power spectrum for building SO(3) kernel") parser.add_argument("-nn", "--nonorm", action='store_true', help="Do not normalize power spectrum") + parser.add_argument("-cs", "--centersorting", action='store_true', help="Sort power spectrum entries according to atomic centers list.") args = parser.parse_args() return args @@ -205,6 +206,7 @@ def set_variable_values_PS(args): fname = args.fname frames = read(fname,':') outfile = args.outfile + censort = args.centersorting sparse_options = [sparsefile] if sparsefile != '': @@ -261,7 +263,7 @@ def set_variable_values_PS(args): sys.exit(0) xyz_slice = [args.slice[0],args.slice[1]] - return [nmax,lmax,rc,sg,cen,spec,cw,lam,periodic,ncut,sparsefile,frames,subset,sparse_options,outfile,args.initial,atomic,all_radial,not args.uselist,xyz_slice,args.imag,args.nonorm] + return [nmax,lmax,rc,sg,cen,spec,cw,lam,periodic,ncut,sparsefile,frames,subset,sparse_options,outfile,args.initial,atomic,all_radial,not args.uselist,xyz_slice,args.imag,args.nonorm,censort] ######################################################################### From fe7b8cb4c83f076a284eeb7e22049d4bac63b52c Mon Sep 17 00:00:00 2001 From: "andrea.grisafi" Date: Thu, 18 Feb 2021 13:11:19 +0100 Subject: [PATCH 2/2] readme update --- README.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/README.rst b/README.rst index 8a9133b..aa55a30 100644 --- a/README.rst +++ b/README.rst @@ -1,7 +1,7 @@ -SA-GPR -====== +TENSOAP +======= -This repository contains a Python code for carrying out Symmetry-Adapted Gaussian Process Regression (SA-GPR) for the machine-learning of tensors. For more information, see: +This repository contains a Python code for carrying out Symmetry-Adapted Gaussian Process Regression (SA-GPR) for the machine-learning of tensors, using SOAP and LODE features. For more information, see: 1. Andrea Grisafi, David M. Wilkins, Gabor Csányi, Michele Ceriotti, "Symmetry-Adapted Machine Learning for Tensorial Properties of Atomistic Systems", Phys. Rev. Lett. 120, 036002 (2018) @@ -19,7 +19,7 @@ This repository contains a Python code for carrying out Symmetry-Adapted Gaussia Versions ======== -The current version of SOAPFAST (v3.0.1) is written in python 3. It has the same functionality as the previous release (v2.3), which is written in python 2. +The current version of TENSOAP (v3.0.1) is written in python 3. It has the same functionality as the previous release (v2.3), which is written in python 2. Requirements ============