#!python -u import numpy as np from susmost import make_triangular_Oxy_plane_cell_neg, make_single_point_zmatrix, normal_lattice_task, transferm, transferm, SiteStateType, register_properties,make_tensor, solve_TM from math import sinh, exp, sqrt,atan2, pi import os from time import sleep def main(): #constants INF_E = 1E6 # prohibitive energy level k_B = 8.3144598E-3 # J/mol a = 1.0 # lattice parameter r_cut_off= 8.0 + .0001 #parameters N = 6 T = 400.0 mol_int_w1 = -12.5 mol_int_w2 = 5. mol_int_w3 = -5. cell,atoms = make_triangular_Oxy_plane_cell_neg(a) register_properties('coverage_one', 'coverage_two') #other beta = 1./(k_B*T) mus = [] covs = [] FORCE_GEN_TM = (os.getenv('FORCE_GEN_TM', 'False')!='False') int_matrix = {('one', 'two', 'right'): mol_int_w2, \ ('two', 'one', 'right'): mol_int_w1, \ ('one', 'two', 'left'): mol_int_w1, \ ('two', 'one', 'left'): mol_int_w2, \ ('one', 'two', 'upright'): mol_int_w1, \ ('two', 'one', 'upright'): mol_int_w2, \ ('one', 'two', 'downright'): mol_int_w1, \ ('two', 'one', 'downright'): mol_int_w2, \ ('one', 'two', 'downleft'): mol_int_w2, \ ('two', 'one', 'downleft'): mol_int_w1, \ ('one', 'two', 'upleft'): mol_int_w2, \ ('two', 'one', 'upleft'): mol_int_w1} def trimez_interactions(cc1, cc2): if cc1.lctype == 'Empty' or cc2.lctype == 'Empty': return None c1 = cc1.coords[0] c2 = cc2.coords[0] v = c2 - c1 r = np.linalg.norm(v) # distance if r > 1.001: return None if r < 0.999: return INF_E # forbidden distance if r < 1.001: if (cc1.lctype == cc2.lctype): return mol_int_w3 d = atan2(v[0],v[1]) if abs(d-pi/2) < 0.001: d ='right' elif abs(d+pi/2) < 0.001: d = 'left' elif 0 < d < pi/2: d = 'upright' elif pi/2 < d < pi: d = 'downright' elif -pi < d < -pi/2: d = 'downleft' elif -pi/2 < d < 0: d = 'upleft' return int_matrix[cc1.lctype,cc2.lctype,d] return None interaction = lambda cc1,cc2: trimez_interactions(cc1,cc2) one_state = SiteStateType('one', make_single_point_zmatrix(), -1., ['N'], coverage=1., coverage_one=1., coverage_two=0.) two_state = SiteStateType('two', make_single_point_zmatrix(), -1., ['C'], coverage=1., coverage_one=0., coverage_two=1.) empty_state = SiteStateType('Empty',make_single_point_zmatrix(), 0.0, ['H'], coverage=0., coverage_one=0., coverage_two=0.) lt = normal_lattice_task([one_state, two_state, empty_state], cell, atoms, interaction, r_cut_off, INF_E) assert lt.states_count == 3, lt for mu in np.arange(0., 50.+0.0001, 5.0 ): # [0 ... 50] one_state.props['ads_energy'].value = mu two_state.props['ads_energy'].value = mu if 1: Wg = make_tensor(lt, beta, force_generic=True) Ws = make_tensor(lt, beta, force_generic=False) TMg = transferm.make_TM_operator(Wg, 2) TMs = transferm.make_TM_operator(Ws, 2) x = np.arange(TMg.nN) assert abs(np.linalg.norm(TMg * x) - np.linalg.norm(TMs * x))/np.linalg.norm(TMs * x) < 1E-6 Wg = transferm.transpose(Wg) Ws = transferm.transpose(Ws) TMg = transferm.make_TM_operator(Wg, 2) TMs = transferm.make_TM_operator(Ws, 2) assert abs(np.linalg.norm(TMg * x) - np.linalg.norm(TMs * x))/np.linalg.norm(TMs * x) < 1E-6 W = make_tensor(lt, beta, force_generic=FORCE_GEN_TM) avg_cov = transferm.average_props(lt.states, N, 'coverage', beta) avg_cov_one = transferm.average_props(lt.states, N, 'coverage_one', beta) avg_cov_two = transferm.average_props(lt.states, N, 'coverage_two', beta) tm_sol = solve_TM(W,N) cov = sum(tm_sol.probs*np.array(avg_cov)) cov_one = sum(tm_sol.probs*np.array(avg_cov_one)) cov_two = sum(tm_sol.probs*np.array(avg_cov_two)) covs.append([cov, cov_one, cov_two]) mus.append(mu) print("mu ==", mu, "theta == ",cov, tm_sol.lr, tm_sol.ll, W[0]) print("covs", covs) etalon_covs = np.array([[0.99999619209074875, 0.49999809604538503, 0.49999809604536094], [0.99998276411443687, 0.49999138205722088, 0.49999138205721649], [0.99992001197541414, 0.49996000598770768, 0.49996000598770779], [0.99643341031567645, 0.49821670515784622, 0.49821670515783323], [0.6683049639253047, 0.3341524819626509, 0.33415248196265152], [0.66702973921804887, 0.33351486960902371, 0.33351486960902393], [0.66674654351560281, 0.3333732717578014, 0.33337327175780113], [0.66668050403392864, 0.33334025201696382, 0.33334025201696388], [1.2374283145107218e-05, 6.187141572553636e-06, 6.1871415725536343e-06], [2.6792031805800048e-06, 1.339601590289999e-06, 1.3396015902899996e-06], [5.9230868457531227e-07, 2.9615434228765688e-07, 2.9615434228765677e-07]]) print('max error = ', np.max(abs(covs-etalon_covs))) print('FORCE_GEN_TM=',FORCE_GEN_TM) assert np.allclose(covs, etalon_covs, atol=1e-2), np.sort(abs(covs-etalon_covs)) print("TEST OK") if __name__ == '__main__': main()