#!python -u import numpy as np from math import exp import os from susmost import make_single_point_zmatrix, make_square_plane_cell, \ normal_lattice_task, transferm, SiteStateType, nearest_int def main(): #constants INF_E = 1E6 # prohibitive energy level k_B = 8.3144598E-3 # kJ/mol a = 1.0 # lattice parameter r_cut_off= 8.0 + .0001 precission = 6 #parameters N = 4 T = 100.0 cell,atoms = make_square_plane_cell(a) interaction = lambda cc1,cc2: nearest_int(cc1,cc2) #other beta = 1./(k_B*T) mus = [] covs = [] FORCE_GEN_TM = (os.getenv('FORCE_GEN_TM', 'False')!='False') mono_state = SiteStateType('mono', make_single_point_zmatrix(), -1., ['N'], coverage=1.) empty_state = SiteStateType('Empty',make_single_point_zmatrix(), 0.0, ['H'], coverage=0.) lt = normal_lattice_task([mono_state, empty_state],cell, atoms, interaction, r_cut_off, INF_E,) assert lt.states_count == 1, lt for mu in np.arange(-3., 3.+0.0001, 0.1 ): mono_state.props['ads_energy'].value = mu cov = lt.states[0].get_prop('coverage', beta) covs.append(cov) mus.append(mu) print("mu ==", mu, "theta == ",cov, exp(-mu/(k_B*T))/(1.+exp(-mu/(k_B*T)))) print('mu theta exact_theta') covs = np.array(covs).round(6) print(list(covs)) etalon_covs = [] for mu,c in zip(mus, covs): theor_cov = exp(-mu/(k_B*T))/(1.+exp(-mu/(k_B*T))) etalon_covs.append(theor_cov) print(mu, c, theor_cov) print('max error = ', np.max(abs(covs-etalon_covs))) print('FORCE_GEN_TM=',FORCE_GEN_TM) assert np.allclose(covs, etalon_covs, atol=1e-6), np.sort(abs(covs-etalon_covs)) print("TEST OK") if __name__ == '__main__': main()