potential function, bug fixing, errorbar for the traning set

python-module/ML_of_Forces/DataModel/ForceConfs.py

@@ -56,14 +56,14 @@ class ForceConfs(DataModel):
         with open(filename, 'r') as jsonfile:
             data = json.load(jsonfile)['data']
 
-        self.confs = [timestep['confs'] for timestep in data]
-        self.force = [timestep['force'] for timestep in data]
+        confs = [np.array(timestep['confs']) for timestep in data]
+        force = [np.array(timestep['force']) for timestep in data]
 
-        self.confs = np.asarray(self.confs)
-        self.force = np.asarray(self.force)
+        self.confs = np.asarray(confs)
+        self.force = np.asarray(force)
 
 
-    def subsampling(self, train: int = 0, test: int = 0, seed: int = 0, random: bool = True) -> None:
+    def subsampling(self, train: int = 1, test: int = 1, seed: int = 0, random: bool = True) -> None:
 
         self.nsample = train + test
 

python-module/ML_of_Forces/GP_custom_gen.py

@@ -1045,9 +1045,9 @@ class GaussianProcess3:
             K = diag + off_diag + off_diag.T
             np.set_printoptions(precision = 3)
             self.K = K
-            print("Is K symmetric? ", (K == K.T).all())
+            # print("Is K symmetric? ", (K == K.T).all())
             eigs = np.linalg.eigvalsh(K)
-            print("Is K positive definite?  ", (eigs> 0).all())
+            # print("Is K positive definite?  ", (eigs> 0).all())
 
             return K
 
@@ -1211,7 +1211,7 @@ class GaussianProcess3:
                         kv = K2_ten[i]
                         var[i] = np.diag(mat_kernel_func(X[i], X[i]))
                         var[i] += nugget
-                        var[i] = np.diag(np.einsum('Uab, UbDc, Ddc -> ad', kv, self.inv_ten, kv)) 
+                        var[i] -= np.diag(np.einsum('Uab, UbDc, Ddc -> ad', kv, self.inv_ten, kv)) 
 
                 return pred, var
         else:

python-module/ML_of_Forces/Potential_draw.py

+import numpy as np
+import matplotlib.pyplot as plt
+import random
+import os.path
+import sys
+import datetime
+sys.path.insert(0, '../ML_of_Forces/3D')
+from GP_custom_gen import GaussianProcess3 as GPvec
+from ManyBodyExpansion import MBExp
+
+
+### Importing data from simulation ###
+
+InDir = "../OneDrive/ML_Data/3D/TB/Si/300K"
+
+forces = np.asarray(np.load(os.path.join(InDir,"forcest.npy")))
+confs = np.asarray(np.load(os.path.join(InDir,"confst.npy")))
+
+forces = np.reshape(forces, (len(forces)*len(forces[0]), 3))
+lenc = len(forces)
+print("Database length is: ", lenc)
+
+### Subsampling from big database ###
+
+ncal = 100
+
+ntot = ncal
+
+ind = np.arange(lenc)
+
+ind_ntot = np.random.choice(ind, size=ntot, replace=False)
+ind_ncal = ind_ntot[0:ncal]
+ind_ntest = ind_ntot[ncal:ntot]
+
+print("Ntot Database length is: " ,len(ind_ntot))
+
+### Spliting database in training/testing sets ###
+
+tr_confs = confs[ind_ncal]
+tr_forces = forces[ind_ncal]
+
+### Radial/angular distribution ###
+
+distances = []
+for d in np.arange(len(tr_confs)):
+
+  dist = np.linalg.norm(tr_confs[d], axis = 1)
+  dist = list(dist)
+  distances = distances + dist
+d_min = np.min(distances)
+print(d_min)
+
+
+
+### TRAINING ###
+t0_train = datetime.datetime.now()
+# m_theta0 = [None],
+
+gp = GPvec( ker=[ 'id'], fvecs =['cov_sim'] ,nugget = 1e-10, theta0=np.array([None]), m_theta0 = [1.], sig =.25, bounds = ((0.1,10.),),
+                      optimizer= None , calc_error = True, eval_grad = False)
+
+gp.fit(tr_confs, tr_forces)
+
+tf_train = datetime.datetime.now()
+print("Training computational time is", (tf_train-t0_train).total_seconds())
+
+
+### REMAPPING ###
+
+t0_remap = datetime.datetime.now()
+
+
+d0, df = 1.7, 5.
+Delta_d = 0.05
+
+xgrid = np.linspace(d0, df, (df-d0)/Delta_d)
+confs = np.zeros(((df-d0)/Delta_d, 1, 3))
+confs[:, 0, 0] = xgrid
+pred_forces, variances = gp.predict(confs)
+
+xforces = pred_forces[:, 0]
+xen = np.cumsum(xforces*Delta_d)
+
+pred_errors = 2*np.sqrt(np.sum(variances, axis= 1))
+en_err = np.cumsum(pred_errors*Delta_d)
+
+
+### PLOTS ###
+
+fig, ax1 = plt.subplots()
+
+ax1.plot(xgrid, xen )
+ax1.plot(xgrid, xen + pred_errors)
+ax1.plot(xgrid, xen - pred_errors)
+ 
+ax1.set_ylabel('Potential (eV)')
+ 
+ax1.set_xlabel('Distance (A)')
+
+ax2 = ax1.twinx()
+ax2.hist(distances, 40, alpha = 0.2, normed =1, label = "pdf")
+
+ax2.set_ylabel('p(d)')
+    
+plt.show()   

python-module/main.py

@@ -6,7 +6,7 @@ from ML_of_Forces.ManyBodyExpansion import MBExp
 
 
 data = ForceConfs()
-data.read_json_data(filename='../data-mlforce/Si_TB_300K.json')
+data.read_json_data(filename='../data-mlforce/Si_TB_1000K.json')
 
 train_confs, train_force, test_confs, test_force = data.subsampling(40, 10, random=True)
 
@@ -22,6 +22,8 @@ gp = GPvec(
     eval_grad=False,
 )
 
+print(train_confs.shape)
+
 
 gp.fit(train_confs, train_force)
 pred_force, var_force = gp.predict(test_confs)
@@ -42,8 +44,8 @@ pred_force, var_force = gp.predict(test_confs)
 
 # print(train_confs)
 
-fit = MBExp(gp)
-
-fit.PW_S_fit(0.1,4)
-
-fit.pair_potential_scalars(0.1, 4, Delta_d=0.1)
+# fit = MBExp(gp)
+#
+# fit.PW_S_fit(0.1,4)
+#
+# fit.pair_potential_scalars(0.1, 4, Delta_d=0.1)