make_ff_water.py

# -*- coding: utf-8 -*-
"""
Created on Tue Jan 20 13:14:20 2026

@author: n.patsalidis
"""
import pickle 
import sys
from matplotlib import pyplot as plt
#sys.path.append('../../code refactoring/main_scripts')
import FF_Develop as ff
import pandas as pd
import numpy as np
ele ='Ag'
molecule ='water'
inputfile = f'{ele}{molecule}.in'
#inputfile='Results/EnergyForcesDEA_InductiveBias_0.2/runned.in'
with open(f'init_pickles/{ele}-{molecule}_init.pickle', 'rb') as f:
    data = pickle.load(f)
    data = data[data['sys_name'] == f'{ele}{molecule}']
    data_all = data.copy()
    
x = data_all['energy_error'].to_numpy()
y = data_all['gradient_error'].to_numpy()
nall = len(data_all) 
data_all = data_all [ x <0.1]
nec = len(data_all)
data_all = data_all[ data_all['scf_correction'] < 0.5]
data = data_all
nef = len(data_all)
print(f'data all = {nall}, after energy cleaning --> {nec}  after force cleaning {nef}')
#data = data_all.sample(n=300)
#data = data [ data['surface'] == '111' ]
#data = pd.concat( ( data_all.sample(n=100), data), ignore_index=True )
#data = data.sample(n=1000)
#data = data.sample(n=10)
 
setup = ff.Setup_Interfacial_Optimization(inputfile)


ff.al_help.make_absolute_Energy_to_interaction(data, setup)
#data = data [ data['Energy']<10]
#e = data['Energy'].to_numpy()
print('Number of data points: ',len(data))

#data = data.sample(n=400, weights = np.exp(-0.2*e ))
#data = pd.concat( ( data_opt, data), ignore_index=True )

data,  optimizer = ff.al_help.solve_model(data,setup)

fe, fu, filt = optimizer.get_Forces_and_ForceClass('all')
from matplotlib import pyplot as plt 

_ = plt.figure(figsize=(3.3,3.3), dpi=300)
x = list(data['at_type'])[0]
plt.title('Force error distributions per atom type')
for t,color in zip(['H','O'], ['#1b9e77','#d95f02'] ):
    filt_idx = t ==x
    f =  list(filt_idx)*len(data) 
    u = fe[f] -fu[f]
    print(f'{t} MAE ={ abs(u).mean() :5.4f} BIAS = {u.mean() :5.4f} STD = { u.std() :5.4f}, Min = {u.min() : 5.4f}, Max = {u.max() : 5.4f} ')
    plt.hist(u.flatten(), bins =50,  label=t, histtype='step', color=color, density=True)
plt.xlabel(r'$F_{True}-F_{predict}$ $(kcal/mol/\AA)$')
plt.legend(fontsize=7.4, frameon=False)
plt.savefig(f'{setup.runpath}/force_error_distribution.png', bbox_inches='tight')
plt.show()