Notebook created

data/error_estimates/error_estimates_base.py

+from numpy import zeros, int32, unique, genfromtxt, sort, argsort, arange, \
+        append, array
+from ase import Atom
+from ase.db import connect
+
+def string2symbols(s):
+    """Convert string to list of chemical symbols."""
+    n = len(s)
+
+    if n == 0:
+        return []
+    
+    c = s[0]
+    
+    if c.isdigit():
+        i = 1
+        while i < n and s[i].isdigit():
+            i += 1
+        return int(s[:i]) * string2symbols(s[i:])
+
+    if c == '(':
+        p = 0
+        for i, c in enumerate(s):
+            if c == '(':
+                p += 1
+            elif c == ')':
+                p -= 1
+                if p == 0:
+                    break
+        j = i + 1
+        while j < n and s[j].isdigit():
+            j += 1
+        if j > i + 1:
+            m = int(s[i + 1:j])
+        else:
+            m = 1
+        return m * string2symbols(s[1:i]) + string2symbols(s[j:])
+
+    if c.isupper():
+        i = 1
+        if 1 < n and s[1].islower():
+            i += 1
+        j = i
+        while j < n and s[j].isdigit():
+            j += 1
+        if j > i:
+            m = int(s[i:j])
+        else:
+            m = 1
+        return m * [s[:i]] + string2symbols(s[j:])
+    else:
+        raise ValueError
+
+def get_data(con, name_dict, Z, code, category, keys, recommended_VASP,
+             name_dict_monos_gpaw, name_dict_bins_gpaw):
+    """Function for obtaining the data from the databases for each code
+    Parameters:
+
+    con: The database
+
+    name_dict: dict
+        The dictionary for el. solid/binaries
+
+    Z: array like
+        atomic numbers for el. solids/binaries
+
+    code: str
+        DFT-code
+
+    category: str
+        code specific category string
+
+    keys: code specific parameters
+  """
+    N_monos = len(name_dict)
+    data = zeros(N_monos)
+    if code == 'VASP':
+        rows = con.select(selection=[('category=' + category),
+                                     ('precision=' + keys[0]),
+                                     ('k_point_density=' + str(keys[1])),
+                                     ('functional=' + keys[2])])
+        # For VASP we have to select the recommended Pseudo-Potential:
+        for row in rows:
+            if category.find('binaries') == 0:
+                A = unique(row.numbers)[0] - 1
+                B = unique(row.numbers)[1] - 1
+                rec_POT = ('POTCAR' +
+                           recommended_VASP[A, 1].replace('POTCAR', '') +
+                           recommended_VASP[B, 1].replace('POTCAR', ''))
+                data_POT = row.potcar.replace('_h', '')
+                rec_POT = rec_POT.replace('_3', '')
+                rec_POT = rec_POT.replace('_2', '')
+            else:
+                rec_POT = recommended_VASP[Z[name_dict[row.name]]-1, 1]
+                data_POT = row.potcar
+            if data_POT == rec_POT:
+                data[name_dict[row.name]] = row.total_energy
+
+    elif code == 'FHI-aims':
+        rows = con.select(selection=[('category=' + category),
+                                     ('basis_set=' + keys[0]),
+                                     ('k_point_density=' + str(keys[1])),
+                                     ('functional=' + keys[2]),
+                                     ('tiers=' + keys[3]),
+                                     ('relativistic_treatment=' + keys[4])])
+        for row in rows:
+            data[name_dict[row.name]] = row.total_energy
+
+    elif code == 'exciting':
+        rows = con.select(selection=[('category=' + category),
+                                     ('total_precision','=',float(keys[0])),
+                                     ('k_point_density','=',float(keys[1])),
+                                     ('xc_functional=' + keys[2])])
+        for row in rows:
+            data[name_dict[row.name]] = row.total_energy
+
+    elif code == 'GPAW':
+        rows = con.select(selection=[('category=' + category),
+                                     ('ecut=' + str(keys[0])),
+                                     ('k_point_density=' + str(keys[1]))])
+        # Mapping from formula to system name
+        for row in rows:
+            if category.find('binaries') == 0:
+                data[name_dict[name_dict_bins_gpaw[row.formula]]] = row.energy
+            else:
+                data[name_dict[name_dict_monos_gpaw[row.formula]]] = row.energy
+    return data
+
+def get_rows(con, name_dict, Z, code, category, keys, recommended_VASP,
+             name_dict_monos_gpaw, name_dict_bins_gpaw):
+    """Function for obtaining the rows from the databases for each code
+    Parameters:
+
+    con: The database
+
+    name_dict: dict
+        The dictionary for el. solid/binaries
+
+    Z: array like
+        atomic numbers for el. solids/binaries
+
+    code: str
+        DFT-code
+
+    category: str
+        code specific category string
+
+    keys: code specific parameters
+  """
+    if code == 'VASP':
+        rows = con.select(selection=[('category=' + category),
+                                     ('precision=' + keys[0]),
+                                     ('k_point_density=' + str(keys[1])),
+                                     ('functional=' + keys[2])])
+
+    elif code == 'FHI-aims':
+        rows = con.select(selection=[('category=' + category),
+                                     ('basis_set=' + keys[0]),
+                                     ('k_point_density=' + str(keys[1])),
+                                     ('functional=' + keys[2]),
+                                     ('tiers=' + keys[3]),
+                                     ('relativistic_treatment=' + keys[4])])
+
+    elif code == 'exciting':
+        rows = con.select(selection=[('category=' + category),
+                                     ('total_precision','=',float(keys[0])),
+                                     ('k_point_density','=',float(keys[1])),
+                                     ('xc_functional=' + keys[2])])
+
+    elif code == 'GPAW':
+        rows = con.select(selection=[('category=' + category),
+                                     ('ecut=' + str(keys[0])),
+                                     ('k_point_density=' + str(keys[1]))])
+    return rows
+
+
+def get_keys(code, prec, kpt, xc, tiers, rel):
+    if code == 'VASP':
+        keys = [prec, str(kpt), xc, tiers]
+        ref_keys = [xc]
+    elif code == 'FHI-aims':
+        keys = [prec, str(kpt), xc, tiers, rel]
+        ref_keys = [xc, rel]
+    elif code == 'exciting':
+        keys = [prec, str(kpt), xc]
+        ref_keys = [xc]
+    elif code == 'GPAW':
+        keys = [prec, str(kpt)]
+        ref_keys = [xc]
+    return keys, ref_keys
+
+
+def get_xy(Z, N, data_ref, data):
+    plot_x = sort(Z)
+    plot_y = ((data_ref-data)/N)[argsort(Z)]
+    nonzero=((data!=0)*(data_ref!=0))[argsort(Z)]
+    plot_x = plot_x[nonzero]
+    plot_y = plot_y[nonzero]
+    #print(plot_x[abs(plot_y)>0.001])
+    return plot_x, abs(plot_y)
+
+
+def get_xy_predict(N, data_ref, data):
+    plot_x = ((data) / N)
+    plot_y = data_ref / N
+    nonzero = (plot_y != 0) * (plot_x != 0)
+    plot_x = plot_x[nonzero]
+    plot_y = plot_y[nonzero]
+    #print(plot_x[abs(plot_y)>0.001])
+    return abs(plot_x), abs(plot_y)
+
+
+def get_binary_error_from_solids(error, binaries_to_monos_min, N_bins_min,
+                                 binaries_to_monos_max, N_bins_max, pred):
+    error_A = error[binaries_to_monos_min]
+    error_B = error[binaries_to_monos_max]
+    N_A = N_bins_min
+    N_B = N_bins_max
+    N_AB = N_A + N_B
+
+    if pred == '1':
+        pred_error=(error_A * N_A + error_B * N_B)
+        rel_pred_error = pred_error / N_AB
+    elif pred == '2':
+        pred_error = (error_A * N_A + error_B * N_B)
+        rel_pred_error = pred_error / N_AB
+    return pred_error #, rel_pred_error
+
+
+def do_plot(fig, ax, Z, N, data_ref, data, lab):
+    plot_x = sort(Z)
+    plot_y = ((data_ref - data) / N)[argsort(Z)]
+    nonzero = plot_y != 0
+    plot_x = plot_x[nonzero]
+    plot_y = plot_y[nonzero]
+    #print(plot_x[abs(plot_y)>0.001])
+    ax.semilogy(plot_x, abs(plot_y), 'o', label=lab)
+    ax.legend(numpoints=1, loc=4, fontsize=8)
+    fig.canvas.draw_idle()
+
+
+def do_plot_predict(fig, ax, N, data_ref, data, lab):
+    plot_x = ((data) / N)
+    plot_y = data_ref / N
+    nonzero = (plot_y != 0) * (plot_x != 0)
+    plot_x = plot_x[nonzero]
+    plot_y = plot_y[nonzero]
+    #print(plot_x[abs(plot_y)>0.001])
+    ax.loglog(abs(plot_x), abs(plot_y), 'o', label=lab)
+    ax.plot(ax.get_xlim(), ax.get_xlim(), '-k')
+    ax.legend(numpoints=1, loc=4, fontsize=8)
+    fig.canvas.draw_idle()
+
+
+def get_xy(Z, N, data_ref, data):
+    plot_x = sort(Z)
+    plot_y = ((data_ref-data)/N)[argsort(Z)]
+    nonzero = plot_y != 0
+    plot_x = plot_x#[nonzero]
+    plot_y = plot_y#[nonzero]
+    #print(plot_x[abs(plot_y)>0.001])
+    return plot_x, abs(plot_y)
+
+
+def get_xy_predict(N, data_ref, data):
+    plot_x = ((data) / N)
+    plot_y = data_ref / N
+    nonzero = (plot_y != 0) * (plot_x != 0)
+    plot_x = plot_x[nonzero]
+    plot_y = plot_y[nonzero]
+    #print(plot_x[abs(plot_y)>0.001])
+    return abs(plot_x), abs(plot_y)
+
+def get_xy_Ecoh(Z, data_monos,data_bins,zeroinds,binaries_to_monos_min,binaries_to_monos_max,N_bins_min,N_bins_max):
+    error_A = data_monos[binaries_to_monos_min]
+    error_B = data_monos[binaries_to_monos_max]
+    N_A = N_bins_min
+    N_B = N_bins_max
+    N_AB = N_A + N_B
+    plot_x=Z[zeroinds]
+    plot_y=((data_bins)[zeroinds]-(error_A * N_A + error_B * N_B)[zeroinds])/N_AB[zeroinds]
+    return plot_x, plot_y
+  
+def get_mono_ind(formula,name_dict_monos):
+    symbols = string2symbols(formula) 
+    ind = zeros(0,dtype='int32')
+    N_f = len(symbols)
+    notinset=False
+    for s in symbols:
+        at = Atom(s)
+        if at.number < 10:
+            try:
+                ind = append(ind, name_dict_monos['0'+str(at.number)+'_'+s])
+            except KeyError:
+                print("Element "+s+" is not in the set of elementary solids.")
+                notinset=True
+                break
+        else:
+            try:
+                ind = append(ind, name_dict_monos[str(at.number)+'_'+s])
+            except KeyError:
+                print("Element "+s+" is not in the set of elementary solids.")
+                notinset=True
+                break
+    return N_f, ind, notinset
+
+def intermediate_delta_energy_exciting(elementname, intermediate_prec, con, name_dict_monos_exciting, keys, data_ref, N_mono):
+    import numpy as np
+    import sys
+    import matplotlib.pyplot as mpl
+    import scipy
+    import scipy.stats
+
+    available_precs=np.asarray([30,40,50,60,70,80]) # this list needs to be sorted
+    natoms=N_mono[name_dict_monos_exciting[elementname]]
+    total_energy_ref=data_ref['exciting','monomers',keys[2]][name_dict_monos_exciting[elementname]]
+    
+    if intermediate_prec < available_precs[0]:
+        #sys.exit('intermediate_delta_energy_exciting: precision of one element in a binary is below '+str(available_precs[0]))
+        number_of_points=2
+        x__considered_precs=available_precs[0:number_of_points]
+    elif intermediate_prec > available_precs[-1]:
+        #print('CAUTION: extrapolation! Precision of one element in a binary is above '+str(available_precs[-1]))
+        number_of_points=3
+        number_of_points=-1*number_of_points
+        x__considered_precs=available_precs[number_of_points:]
+    else:
+        for icount in range(1,len(available_precs)):
+            if intermediate_prec < available_precs[icount]:
+                x__considered_precs=available_precs[[icount-1,icount]]
+                break
+                
+    y__log_error_per_atom=[]
+    #for prec in x__considered_precs:
+    rows = con.select(selection=[('category=monomers'),
+                                        ('total_precision','<=',x__considered_precs[-1]),
+                                        ('total_precision','>=',x__considered_precs[0]),
+                                        ('k_point_density','=',float(keys[1])),
+                                        ('xc_functional=' + keys[2]),
+                                        ('element1='+elementname)])
+    for row in rows:
+        total_energy=row.total_energy        
+        y__log_error_per_atom.append(np.log10((total_energy - total_energy_ref)/natoms))
+
+    linregr_results = scipy.stats.linregress(x__considered_precs,y__log_error_per_atom)
+    p = np.asarray([linregr_results[0], linregr_results[1]])
+    delta_energy_per_atom = np.power(10,np.polyval(p,intermediate_prec))
+    # want to plot results of linregress?
+    if False:
+        # get all monomers data
+        monomers_delta_energy_per_atom=[]
+        rows = con.select(selection=[('category=monomers'),
+                                        ('total_precision','<=',available_precs[-1]),
+                                        ('total_precision','>=',available_precs[0]),                                        
+                                        ('k_point_density','=',float(keys[1])),
+                                        ('xc_functional=' + keys[2]),
+                                        ('element1='+elementname)])
+        for row in rows:
+            total_energy=row.total_energy        
+            monomers_delta_energy_per_atom.append((total_energy - total_energy_ref)/natoms)
+        # plot
+        fig3 = figure(3,(10,10))
+        rect = fig3.patch
+        rect.set_facecolor('white')
+        mpl.semilogy(available_precs, monomers_delta_energy_per_atom,'ko-', ms=7.0, markeredgewidth=1, linewidth=1.5, label="data")
+        mpl.semilogy(x__considered_precs, np.power(10,y__log_error_per_atom),'go', ms=11.0, markeredgewidth=1, label="accounted data in linregress")
+        mpl.semilogy(intermediate_prec, delta_energy_per_atom,'rx', ms=11.0, markeredgewidth=3, label="result from linregress")
+        axes = mpl.gca()
+        axes.grid(True)
+        mpl.xlabel("Precision %", fontsize = "x-large")
+        mpl.ylabel("$\Delta E$ per atom [eV]", fontsize = "x-large")
+        mpl.legend(loc="lower left", frameon=True, fontsize = "x-large", numpoints=1)
+        mpl.title(elementname)
+        mpl.show()
+    return delta_energy_per_atom

data/error_estimates/recommended_POTCAR.dat

+H POTCAR 250  1
+He POTCAR 479  2
+Li POTCAR_sv 499  3
+Be POTCAR 248  2
+B POTCAR 319  3
+C POTCAR 400  4
+N POTCAR 400  5
+O POTCAR 400  6
+F POTCAR 400  7
+Ne POTCAR 344  8
+Na POTCAR_pv 260  7
+Mg POTCAR 200  2
+Al POTCAR 240  3
+Si POTCAR 245  4
+P POTCAR 255  5
+S POTCAR 259  6
+Cl POTCAR 262  7
+Ar POTCAR 266  8
+K POTCAR_sv 259  9
+Ca POTCAR_sv 267  10
+Sc POTCAR_sv 223  11
+Ti POTCAR_sv 275  12
+V POTCAR_sv 264  13
+Cr POTCAR_pv 266  12
+Mn POTCAR_pv 270  13
+Fe POTCAR 268  8
+Co POTCAR 268  9
+Ni POTCAR 270  10
+Cu POTCAR 295  11
+Zn POTCAR 277  12
+Ga POTCAR_d 283  13
+Ge POTCAR_d 310  14
+As POTCAR 209  5
+Se POTCAR 212  6
+Br POTCAR 216  7
+Kr POTCAR 185  8
+Rb POTCAR_sv 220  9
+Sr POTCAR_sv 229  10
+Y POTCAR_sv 203  11
+Zr POTCAR_sv 230  12
+Nb POTCAR_sv 293  13
+Mo POTCAR_sv 243  14
+Tc POTCAR_pv 264  13
+Ru POTCAR_pv 240  14
+Rh POTCAR_pv 247  15
+Pd POTCAR 251  10
+Ag POTCAR 250  11
+Cd POTCAR 274  12
+In POTCAR_d 239  13
+Sn POTCAR_d 241  14
+Sb POTCAR 172  5
+Te POTCAR 175  6
+I POTCAR 176  7
+Xe POTCAR 153  8
+Cs POTCAR_sv 220  9
+Ba POTCAR_sv 187  10
+La POTCAR 219  11
+Ce POTCAR 273  12
+Pr POTCAR_3 182  11
+Nd POTCAR_3 183  11
+Pm POTCAR_3 177  11
+Sm POTCAR_3 177  11
+Eu POTCAR_3 129  9
+Gd POTCAR_3 154  9
+Tb POTCAR_3 156  9
+Dy POTCAR_3 156  9
+Ho POTCAR_3 154  9
+Er POTCAR_3 155  9
+Tm POTCAR_3 149  9
+Yb POTCAR_2 113  8
+Lu POTCAR_3 155  9
+Hf POTCAR_pv 220  10
+Ta POTCAR_pv 224  11
+W POTCAR_sv 223  12
+Re POTCAR 226  7
+Os POTCAR 228  8
+Ir POTCAR 211  9
+Pt POTCAR 230  10
+Au POTCAR 230  11
+Hg POTCAR 233  12
+Tl POTCAR_d 237  13
+Pb POTCAR_d 238  14
+Bi POTCAR_d 243  15
+Po POTCAR_d 265  16
+At POTCAR_d 266  17
+Rn POTCAR 152  8
+Fr POTCAR_sv 215  9
+Ra POTCAR_sv 237  10
+Ac POTCAR 172  11
+Th POTCAR 247  12
+Pa POTCAR 252  13
+U POTCAR 253  14
+Np POTCAR 254  15
+Pu POTCAR 254  16
+Am POTCAR 256  17
+Cm POTCAR 258  18
\ No newline at end of file