Basic Tutorial
# Load in the relevant packages
from ase.io import read,write
from ase.build import bulk
import numpy as np
import ase.db as db
from ase.visualize import view
from ase.calculators.lj import LennardJones
import matplotlib.pyplot as plt
Fit_on_forces=True #Enable/disable option for fitting CCS potential to atomic forces.
Generate the CCS_fit input file; structures.json
The next commands fetch the training set data from the ASE database containing the LJ structures and the corresponding energies. The reference data is called LJ.db since the reference energies and forces are obtained from a Lennard Jones potential.
from ccs_fit.scripts.ccs_fetch import ccs_fetch
ccs_fetch(mode="CCS",DFT_DB="LJ.db",include_forces=Fit_on_forces)
Fetching data: 0%| | 0/40 [00:00<?, ?it/s]
Fetching data: 18%|█▊ | 7/40 [00:00<00:00, 65.69it/s]
Fetching data: 35%|███▌ | 14/40 [00:00<00:00, 48.30it/s]
Fetching data: 50%|█████ | 20/40 [00:00<00:00, 42.40it/s]
Fetching data: 68%|██████▊ | 27/40 [00:00<00:00, 50.06it/s]
Fetching data: 82%|████████▎ | 33/40 [00:00<00:00, 42.96it/s]
Fetching data: 100%|██████████| 40/40 [00:00<00:00, 49.67it/s]
Fetching data: 100%|██████████| 40/40 [00:00<00:00, 48.25it/s]
Fit training data to Curvature Constrained Splines
Finally, the splines are fitted to the target defined in the structures.json input file. The splines can be restricted to be fully repulsive (“rep”), or have a turning point/switch (“sw”), which is defined by the “Swtype” key. A more comprehensive guide on the fitting options can be found in Advanced_Tutorials/Search_Mode.
### Generate input.json file
import json
input={
"General": {
"interface": "CCS",
"merging" : "True"
},
"Twobody": {
"Ne-Ne": {
"Rcut": 5.0,
"Resolution": 0.01,
"Swtype": "sw",
"const_type" : "Mono"
}
}
}
#SAVE TO FILE
with open('CCS_input.json', 'w') as f:
json.dump(input, f, indent=8)
#RUN FIT
from ccs_fit import ccs_fit
ccs_fit("CCS_input.json")
Generating one-body information from training-set.
Added elements: ['Ne']
Applying monotonic constraints for pair: Ne-Ne
Merging intervals for pair Ne-Ne; number of intervals reduced from 412 to 411.
Finding optimum switch: 0%| | 0/411 [00:00<?, ?it/s]
Finding optimum switch: 0%| | 1/411 [00:00<01:31, 4.50it/s]
Finding optimum switch: 0%| | 2/411 [00:00<01:27, 4.68it/s]
Finding optimum switch: 1%| | 3/411 [00:00<01:26, 4.72it/s]
Finding optimum switch: 1%| | 4/411 [00:00<01:27, 4.64it/s]
Finding optimum switch: 1%| | 5/411 [00:01<01:27, 4.66it/s]
Finding optimum switch: 1%|▏ | 6/411 [00:01<01:24, 4.78it/s]
Finding optimum switch: 2%|▏ | 7/411 [00:01<01:21, 4.93it/s]
Finding optimum switch: 2%|▏ | 8/411 [00:01<01:20, 5.04it/s]
Finding optimum switch: 2%|▏ | 9/411 [00:01<01:17, 5.19it/s]
Finding optimum switch: 2%|▏ | 10/411 [00:02<01:16, 5.24it/s]
Finding optimum switch: 3%|▎ | 11/411 [00:02<01:14, 5.34it/s]
Finding optimum switch: 3%|▎ | 12/411 [00:02<01:12, 5.50it/s]
Finding optimum switch: 3%|▎ | 13/411 [00:02<01:07, 5.92it/s]
Finding optimum switch: 3%|▎ | 14/411 [00:02<01:05, 6.03it/s]
Finding optimum switch: 4%|▎ | 15/411 [00:02<01:05, 6.09it/s]
Finding optimum switch: 4%|▍ | 16/411 [00:02<01:05, 6.00it/s]
Finding optimum switch: 4%|▍ | 17/411 [00:03<01:05, 6.02it/s]
Finding optimum switch: 4%|▍ | 18/411 [00:03<01:04, 6.05it/s]
Finding optimum switch: 5%|▍ | 19/411 [00:03<01:02, 6.29it/s]
Finding optimum switch: 5%|▍ | 20/411 [00:03<01:00, 6.44it/s]
Finding optimum switch: 5%|▌ | 21/411 [00:03<00:58, 6.70it/s]
Finding optimum switch: 5%|▌ | 22/411 [00:03<00:56, 6.90it/s]
Finding optimum switch: 6%|▌ | 23/411 [00:04<00:54, 7.09it/s]
Finding optimum switch: 6%|▌ | 24/411 [00:04<00:53, 7.28it/s]
Finding optimum switch: 6%|▌ | 25/411 [00:04<00:53, 7.26it/s]
Finding optimum switch: 6%|▋ | 26/411 [00:04<00:54, 7.09it/s]
Finding optimum switch: 7%|▋ | 27/411 [00:04<00:55, 6.91it/s]
Finding optimum switch: 7%|▋ | 28/411 [00:04<00:56, 6.78it/s]
Finding optimum switch: 7%|▋ | 29/411 [00:04<00:58, 6.56it/s]
Finding optimum switch: 7%|▋ | 30/411 [00:05<01:01, 6.22it/s]
Finding optimum switch: 8%|▊ | 31/411 [00:05<01:01, 6.18it/s]
Finding optimum switch: 8%|▊ | 32/411 [00:05<01:00, 6.21it/s]
Finding optimum switch: 8%|▊ | 33/411 [00:05<00:59, 6.37it/s]
Finding optimum switch: 8%|▊ | 34/411 [00:05<00:58, 6.41it/s]
Finding optimum switch: 9%|▊ | 35/411 [00:05<01:01, 6.13it/s]
Finding optimum switch: 9%|▉ | 36/411 [00:06<01:05, 5.72it/s]
Finding optimum switch: 9%|▉ | 37/411 [00:06<01:06, 5.67it/s]
Finding optimum switch: 9%|▉ | 38/411 [00:06<01:05, 5.69it/s]
Finding optimum switch: 9%|▉ | 39/411 [00:06<01:07, 5.52it/s]
Finding optimum switch: 10%|▉ | 40/411 [00:06<01:07, 5.51it/s]
Finding optimum switch: 10%|▉ | 41/411 [00:06<01:06, 5.55it/s]
Finding optimum switch: 10%|█ | 42/411 [00:07<01:04, 5.74it/s]
Finding optimum switch: 10%|█ | 43/411 [00:07<01:03, 5.80it/s]
Finding optimum switch: 11%|█ | 44/411 [00:07<01:02, 5.84it/s]
Finding optimum switch: 11%|█ | 45/411 [00:07<01:02, 5.89it/s]
Finding optimum switch: 11%|█ | 46/411 [00:07<01:01, 5.96it/s]
Finding optimum switch: 11%|█▏ | 47/411 [00:07<01:00, 6.05it/s]
Finding optimum switch: 12%|█▏ | 48/411 [00:08<00:59, 6.07it/s]
Finding optimum switch: 12%|█▏ | 49/411 [00:08<00:58, 6.14it/s]
Finding optimum switch: 12%|█▏ | 50/411 [00:08<00:58, 6.20it/s]
Finding optimum switch: 12%|█▏ | 51/411 [00:08<00:58, 6.15it/s]
Finding optimum switch: 13%|█▎ | 52/411 [00:08<00:57, 6.20it/s]
Finding optimum switch: 13%|█▎ | 53/411 [00:08<00:56, 6.35it/s]
Finding optimum switch: 13%|█▎ | 54/411 [00:09<00:57, 6.26it/s]
Finding optimum switch: 13%|█▎ | 55/411 [00:09<00:56, 6.26it/s]
Finding optimum switch: 14%|█▎ | 56/411 [00:09<00:56, 6.33it/s]
Finding optimum switch: 14%|█▍ | 57/411 [00:09<00:56, 6.27it/s]
Finding optimum switch: 14%|█▍ | 58/411 [00:09<00:56, 6.28it/s]
Finding optimum switch: 14%|█▍ | 59/411 [00:09<00:55, 6.35it/s]
Finding optimum switch: 15%|█▍ | 60/411 [00:10<00:56, 6.25it/s]
Finding optimum switch: 15%|█▍ | 61/411 [00:10<00:55, 6.27it/s]
Finding optimum switch: 15%|█▌ | 62/411 [00:10<00:54, 6.40it/s]
Finding optimum switch: 15%|█▌ | 63/411 [00:10<00:54, 6.37it/s]
Finding optimum switch: 16%|█▌ | 64/411 [00:10<00:55, 6.30it/s]
Finding optimum switch: 16%|█▌ | 65/411 [00:10<00:54, 6.30it/s]
Finding optimum switch: 16%|█▌ | 66/411 [00:10<00:54, 6.34it/s]
Finding optimum switch: 16%|█▋ | 67/411 [00:11<00:54, 6.30it/s]
Finding optimum switch: 17%|█▋ | 68/411 [00:11<00:53, 6.36it/s]
Finding optimum switch: 17%|█▋ | 69/411 [00:11<00:53, 6.35it/s]
Finding optimum switch: 17%|█▋ | 70/411 [00:11<00:53, 6.34it/s]
Finding optimum switch: 17%|█▋ | 71/411 [00:11<00:53, 6.39it/s]
Finding optimum switch: 18%|█▊ | 72/411 [00:11<00:53, 6.34it/s]
Finding optimum switch: 18%|█▊ | 73/411 [00:12<00:54, 6.15it/s]
Finding optimum switch: 18%|█▊ | 74/411 [00:12<00:54, 6.20it/s]
Finding optimum switch: 18%|█▊ | 75/411 [00:12<00:53, 6.28it/s]
Finding optimum switch: 18%|█▊ | 76/411 [00:12<00:53, 6.25it/s]
Finding optimum switch: 19%|█▊ | 77/411 [00:12<00:52, 6.39it/s]
Finding optimum switch: 19%|█▉ | 78/411 [00:12<00:50, 6.55it/s]
Finding optimum switch: 19%|█▉ | 79/411 [00:13<00:49, 6.68it/s]
Finding optimum switch: 19%|█▉ | 80/411 [00:13<00:49, 6.65it/s]
Finding optimum switch: 20%|█▉ | 81/411 [00:13<00:48, 6.74it/s]
Finding optimum switch: 20%|█▉ | 82/411 [00:13<00:47, 6.88it/s]
Finding optimum switch: 20%|██ | 83/411 [00:13<00:47, 6.91it/s]
Finding optimum switch: 20%|██ | 84/411 [00:13<00:46, 7.00it/s]
Finding optimum switch: 21%|██ | 85/411 [00:13<00:46, 7.06it/s]
Finding optimum switch: 21%|██ | 86/411 [00:14<00:46, 7.04it/s]
Finding optimum switch: 21%|██ | 87/411 [00:14<00:46, 6.91it/s]
Finding optimum switch: 21%|██▏ | 88/411 [00:14<00:47, 6.85it/s]
Finding optimum switch: 22%|██▏ | 89/411 [00:14<00:47, 6.82it/s]
Finding optimum switch: 22%|██▏ | 90/411 [00:14<00:47, 6.77it/s]
Finding optimum switch: 22%|██▏ | 91/411 [00:14<00:47, 6.76it/s]
Finding optimum switch: 22%|██▏ | 92/411 [00:14<00:47, 6.75it/s]
Finding optimum switch: 23%|██▎ | 93/411 [00:15<00:46, 6.81it/s]
Finding optimum switch: 23%|██▎ | 94/411 [00:15<00:46, 6.75it/s]
Finding optimum switch: 23%|██▎ | 95/411 [00:15<00:46, 6.81it/s]
Finding optimum switch: 23%|██▎ | 95/411 [00:15<00:51, 6.18it/s]
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[4], line 4
1 #RUN FIT
2 from ccs_fit import ccs_fit
----> 4 ccs_fit("CCS_input.json")
File ~/work/CCS/CCS/src/ccs_fit/fitting/main.py:411, in twp_fit(filename)
399 nn = Objective(
400 atom_pairs,
401 atom_onebodies,
(...)
407 force_ewald=ewald_forces,
408 )
410 # Solve QP problem
--> 411 predicted_energies, mse, xx_unfolded = nn.solution()
File ~/work/CCS/CCS/src/ccs_fit/fitting/objective.py:169, in Objective.solution(self)
167 hh = np.zeros(gg.shape[0])
168 bb = np.zeros(aa.shape[0])
--> 169 sol = self.solver(pp, qq, matrix(gg), matrix(hh),
170 matrix(aa), matrix(bb))
171 obj.append(float(self.eval_obj(sol["x"])))
173 obj = np.asarray(obj)
File ~/work/CCS/CCS/src/ccs_fit/fitting/objective.py:300, in Objective.solver(pp, qq, gg, hh, aa, bb, maxiter, tol)
298 sol = solvers.qp(pp, qq, gg, hh, aa, bb)
299 else:
--> 300 sol = solvers.qp(pp, qq, gg, hh)
302 return sol
File ~/.cache/pypoetry/virtualenvs/ccs-fit-oj5MS1KF-py3.9/lib/python3.9/site-packages/cvxopt/coneprog.py:4485, in qp(P, q, G, h, A, b, solver, kktsolver, initvals, **kwargs)
4475 pinfres, dinfres = None, None
4477 return {'status': status, 'x': x, 's': s, 'y': y, 'z': z,
4478 'primal objective': pcost, 'dual objective': dcost,
4479 'gap': gap, 'relative gap': relgap,
(...)
4482 'residual as primal infeasibility certificate': pinfres,
4483 'residual as dual infeasibility certificate': dinfres}
-> 4485 return coneqp(P, q, G, h, None, A, b, initvals, kktsolver = kktsolver, options = options)
File ~/.cache/pypoetry/virtualenvs/ccs-fit-oj5MS1KF-py3.9/lib/python3.9/site-packages/cvxopt/coneprog.py:2174, in coneqp(P, q, G, h, dims, A, b, initvals, kktsolver, xnewcopy, xdot, xaxpy, xscal, ynewcopy, ydot, yaxpy, yscal, **kwargs)
2172 f0 = 0.5 * (xdot(x, rx) + xdot(x, q))
2173 fA(y, rx, beta = 1.0, trans = 'T')
-> 2174 fG(z, rx, beta = 1.0, trans = 'T')
2175 resx = math.sqrt(xdot(rx, rx))
2177 # ry = A*x - b
File ~/.cache/pypoetry/virtualenvs/ccs-fit-oj5MS1KF-py3.9/lib/python3.9/site-packages/cvxopt/coneprog.py:1897, in coneqp.<locals>.fG(x, y, trans, alpha, beta)
1896 def fG(x, y, trans = 'N', alpha = 1.0, beta = 0.0):
-> 1897 misc.sgemv(G, x, y, dims, trans = trans, alpha = alpha,
1898 beta = beta)
File ~/.cache/pypoetry/virtualenvs/ccs-fit-oj5MS1KF-py3.9/lib/python3.9/site-packages/cvxopt/misc.py:830, in sgemv(A, x, y, dims, trans, alpha, beta, n, offsetA, offsetx, offsety)
828 if n is None: n = A.size[1]
829 if trans == 'T' and alpha: trisc(x, dims, offsetx)
--> 830 base.gemv(A, x, y, trans = trans, alpha = alpha, beta = beta, m = m,
831 n = n, offsetA = offsetA, offsetx = offsetx, offsety = offsety)
832 if trans == 'T' and alpha: triusc(x, dims, offsetx)
KeyboardInterrupt:
Validate your potential
Make sure your potential (at least) reproduce the data points in your training-set. Performing further tests on strucutres not included in the training set is recomended but not included in the tutorial.
from ccs_fit.scripts.ccs_validate import ccs_validate
ccs_validate(mode="CCS",CCS_params="CCS_params.json",DFT_DB="LJ.db")
0%| | 0/40 [00:00<?, ?it/s]/home/thism292/anaconda3/envs/CCS_fit/lib/python3.9/site-packages/ase/utils/__init__.py:62: FutureWarning: Please use atoms.cell.rank instead
warnings.warn(warning)
100%|██████████| 40/40 [00:01<00:00, 29.38it/s]
with open("CCS_params.json", "r") as f:
CCS_params = json.load(f)
with open("structures.json", "r") as f:
training_set = json.load(f)
r=np.array(CCS_params["Two_body"]["Ne-Ne"]["r"])
e=CCS_params["Two_body"]["Ne-Ne"]["spl_a"]
e_LJ= 4 * ((1 / r) ** 12 - (1 / r) ** 6)
plt.xlim(0.5,3)
plt.ylim(-1.5,1)
plt.xlabel('Distance (Å)')
plt.ylabel('Potential (eV)')
plt.plot(r,e_LJ,color='black',label="Ref. Lennard-Jones potential")
plt.plot(r,e,'--',color='red',label="Fitted potential")
plt.legend()
plt.show()
err=np.loadtxt("CCS_validate.dat")
plt.xlabel('Reference energy (eV)')
plt.ylabel('Validation energy (eV)')
plt.plot( [min(err[:,0]),max(err[:,0])],[min(err[:,0]),max(err[:,0])],'--',color='black' )
plt.scatter(err[:,0],err[:,1],facecolors='none', edgecolors='red')
plt.show()
plt.xlabel('Reference energy (eV)')
plt.ylabel('Error in fit (eV)')
plt.scatter(err[:,0],err[:,2],facecolors='none', edgecolors='red')
plt.show()
try:
err_F=np.loadtxt("CCS_error_forces.out")
plt.xlabel('Reference force (eV/Å)')
plt.ylabel('Fitted force (eV/Å)')
plt.plot( [min(err_F[:,0]),max(err_F[:,0])],[min(err_F[:,0]),max(err_F[:,0])],'--',color='black')
plt.scatter(err_F[:,0],err_F[:,1],facecolors='none', edgecolors='red',alpha=0.1 )
plt.show()
except:
pass
d=[]
for t in training_set["energies"]:
d.extend(training_set["energies"][t]["Ne-Ne"])
plt.hist(d)
plt.show()
from ccs_fit.scripts.ccs_export_FF import write_FF
write_FF("CCS_params.json")
Cleaning up
import glob
import os
def rm(file):
try:
os.remove(file)
except OSError:
pass
list=glob.glob("CALCULATED_DATA/*")
for file in list:
rm(file)
list=glob.glob("CCS_*")
for file in list:
rm(file)
rm("structures.json")
rm("file_list")
rm("ccs.spl")