Sorted species

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
 ============

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

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:
 

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
 

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
 

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      
 

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]
 
 #########################################################################

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]
 
 #########################################################################