FF-optimization example

[qzhu2017]

# disable the lammps output
# /Users/qzhu8/miniconda3/envs/ost/lib/python3.9/site-packages/lammps/pylammps.py

from pyocse.parameters import ForceFieldParameters
from time import time
import os

params = ForceFieldParameters(
    smiles=['CC(=O)OC1=CC=CC=C1C(=O)O'],
    f_coef=1.0,
    s_coef=1.0,
    style = 'openff',
    ref_evaluator='mace',
    ncpu=1)

p0, errors = params.load_parameters("dataset/parameters.xml")
ref_dics = params.load_references("dataset/references.xml")[:200]

os.makedirs("ASP", exist_ok=True)
os.chdir("ASP")

t0 = time()
e_offset, params_opt = params.optimize_offset(ref_dics, p0)
params.update_ff_parameters(params_opt)
print(params.get_objective(ref_dics, e_offset, obj="MSE"))
errs = params.plot_ff_results("performance_init.png", ref_dics, [params_opt])

# Stepwise optimization of each term
for data in [
    (["bond", "angle", "proper"], 100),
    (["proper", "vdW", "charge"], 100),
    (["bond", "angle", "proper", "vdW", "charge"], 100),
]:
    (terms, steps) = data
    opt_dict = params.get_opt_dict(terms, None, params_opt)
    x, fun, values, it = params.optimize_global(ref_dics,
                                                opt_dict,
                                                params_opt,
                                                obj="MSE",
                                                t0=0.1,
                                                steps=25)
    params_opt = params.set_sub_parameters(values, terms, params_opt)

    opt_dict = params.get_opt_dict(terms, None, params_opt)
    x, fun, values, it = params.optimize_local(ref_dics,
                                               opt_dict,
                                               params_opt,
                                               obj="MSE",
                                               steps=steps)
    params_opt = params.set_sub_parameters(values, terms, params_opt)
    _, params_opt = params.optimize_offset(ref_dics, params_opt)

    t = (time() - t0) / 60
    print(f"\nFF optimization {t:.2f} min ", data)

#print(params.get_objective(ref_dics, e_offset))
errs = params.plot_ff_results("performance_opt.png", ref_dics, [params_opt])
params.plot_ff_parameters("parameters.png", [params_opt])
params.export_parameters("parameters_opt.xml", params_opt, errs[0])

 Opt. Energy (124)
Unit: [eV/mole]
RMSE: 0.3873
R2:   0.5196

 Opt. Forces (7056)
Unit: [eV/A]
RMSE: 0.6579
R2:   0.4196

 Opt. Stress (618)
Unit: [GPa]
RMSE: 0.0101
R2:   0.8258

Min_values: -139.8203 -139.5780
Init obj 4605.335512772165
Step    0  0.10 4762.6432 4605.3355
Step   10  0.09 4083.9351 4083.9351
Step   20  0.08 3464.5167 3461.1658
Best results after 25 steps: 3461.1658
Init obj 3461.165812874149
Step   10 3443.0312
Step   20 3434.6277
Step   30 3417.9817
Step   40 3318.0512
Step   50 3245.5106
Step   60 3184.1237
Step   70 3161.7609
Step   80 3154.2597
Step   90 3138.6318
Final Obj 3080.241259386598
optimized offset -137.22517993563605

FF optimization 4.98 min  (['bond', 'angle', 'proper'], 100)
Init obj 3072.0159005550972
Step    0  0.10 3088.5674 3072.0159
Step   10  0.09 3061.3139 3034.0383
Step   20  0.08 3013.4101 3012.0551
Best results after 25 steps: 3009.9888
Init obj 3009.988757736076
Step   10 2989.2609
Step   20 2985.5334
Step   30 2976.2584
Step   40 2963.0684
Step   50 2947.4123
Step   60 2948.2801
Step   70 2934.2236
Step   80 2923.3754
Step   90 2916.9483
Final Obj 2898.0347214881026
optimized offset -137.86219364361216

FF optimization 9.36 min  (['proper', 'vdW', 'charge'], 100)
Init obj 2921.9803266137183
Step    0  0.10 3061.3015 2921.9803
Step   10  0.09 2898.9479 2878.0570
Step   20  0.08 2854.0136 2854.0136
Best results after 25 steps: 2854.0136
Init obj 2854.0136226021814
Step   10 3052.9837
Step   20 2894.6201
Step   30 2857.3746
Step   40 2833.1030
Step   50 2832.8452
Step   60 2828.5256
Step   70 2822.4829
Step   80 2815.7560
Step   90 2806.7617
Final Obj 2776.942565526142
optimized offset -138.3939336948552

FF optimization 15.26 min  (['bond', 'angle', 'proper', 'vdW', 'charge'], 100)
Number of reference structures 124
(0.48054640530065185, 0.6188915013619671, 0.8193665614865759)

 Opt. Energy (124)
Unit: [eV/mole]
RMSE: 0.3940
R2:   0.4805

 Opt. Forces (7056)
Unit: [eV/A]
RMSE: 0.4334
R2:   0.6189

 Opt. Stress (618)
Unit: [GPa]
RMSE: 0.0102
R2:   0.8194

[qzhu2017]

• Better objective function (encourage good fitting on the low-E region, weights, ....)