#!python -u import numpy as np from susmost import make_single_point_zmatrix, make_dimer_zmatrix, \ make_square_plane_cell, normal_lattice_task, transferm, SiteStateType, make_tensor, solve_TM from math import sinh, exp, sqrt import os 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 #parameters N = 6 T = 300.0 mol_int_w1 = -15. mol_int_w2 = 5. mol_int_w3 = 5. cell,atoms = make_square_plane_cell(a) interaction = lambda cc1,cc2: porph_interactions(cc1,cc2) #other beta = 1./(k_B*T) mus = [] covs = [] FORCE_GEN_TM = (os.getenv('FORCE_GEN_TM', 'False')!='False') def porph_interactions(cc1, cc2): if cc1.lctype == 'Empty' or cc2.lctype == 'Empty': return None c1_1 = cc1.coords[0] c1_2 = cc1.coords[1] c2_1 = cc2.coords[0] c2_2 = cc2.coords[1] first = sum(abs(c1_2-c2_1)) second = sum(abs(c2_2-c1_1)) r = np.linalg.norm(c1_1 - c2_1) # distance if r < 0.999: return INF_E # forbidden distance if r < 1.001: if (first < 0.01) and (second < 0.01): return mol_int_w1 elif (first < 0.01) or (second < 0.01): return mol_int_w3 else: return mol_int_w2 return None porph_state = SiteStateType('porph', make_dimer_zmatrix(a), -1., ['N','O'], coverage=1.) empty_state = SiteStateType('Empty',make_single_point_zmatrix(), 0.0, ['H'], coverage=0.) lt = normal_lattice_task([porph_state, empty_state],cell, atoms, interaction, r_cut_off, INF_E) assert lt.states_count == 5, lt for mu in np.concatenate((np.arange(-30., -20.+0.0001, 1.0 ), np.arange(-20., 20.+0.0001, 5.0 ))): # [-30 ... 20] porph_state.props['ads_energy'].value = mu W = make_tensor(lt, beta, force_generic=FORCE_GEN_TM) avg_cov = transferm.average_props(lt.states, N, 'coverage', beta) tm_sol = solve_TM(W,N, simple_iteration=True, k=2, precond_power = 2, ncv = 40) cov = sum(tm_sol.probs*np.array(avg_cov)) covs.append(cov) mus.append(mu) print("mu ==", mu, "theta == ",cov, tm_sol.lr, tm_sol.ll) print("covs", covs) etalon_covs = np.array([0.9984879232910383, 0.99686317316018114, 0.99371046019443809, 0.98815518414391024, 0.9795170646475182, 0.96758763918707336, 0.95244442413987973, 0.93425479219476137, 0.91336439324939711, 0.89038401237318232, 0.86610471681776846, 0.86610471681777124, 0.74811924800728968, 0.67363609127931368, 0.6432160769829115, 0.5093638238173207, 0.49970026682484309, 0.41999534473451194, 0.26416592940438305, 0.053050517482918592]) print('errs', abs(covs-etalon_covs)) print('max error = ', np.max(abs(covs-etalon_covs))) print('FORCE_GEN_TM=',FORCE_GEN_TM) scipy_version = list(map(int, transferm.scipy.__version__.split('.'))) print('transferm.scipy.__version__ == ', transferm.scipy.__version__, scipy_version) assert np.allclose(covs, etalon_covs, atol=1e-5) print("TEST OK") if __name__ == '__main__': main()