#!python -u import numpy as np from susmost import make_single_point_zmatrix, make_square_plane_cell,\ normal_lattice_task, transferm, register_property, SiteStateType, make_tensor, sitestate, make_tensor, solve_TM from math import sinh, exp, sqrt import os from time import sleep 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 = 400.0 mol_int_w1 = -12.5 mol_int_w2 = 5. cell,atoms = make_square_plane_cell(a) register_property('coverage_rightleft','==', None, sitestate.join_cover) register_property('coverage_updown','==', None, sitestate.join_cover) #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 = cc1.coords[0] c2 = cc2.coords[0] r = np.linalg.norm(c1 - c2) # distance if r < 0.999: return INF_E # forbidden distance if r < 1.001: if abs(c1[0]-c2[0]) < 0.01: #vertical if cc1.lctype == 'updown' and cc2.lctype == 'updown': return mol_int_w1 else: return mol_int_w2 if abs(c1[1]-c2[1]) < 0.01: #horizontal if cc1.lctype == 'rightleft' and cc2.lctype == 'rightleft': return mol_int_w1 else: return mol_int_w2 assert(False) #imposible situation return None interaction = lambda cc1,cc2: porph_interactions(cc1,cc2) rightleft_state = SiteStateType('rightleft', make_single_point_zmatrix(), -1., ['N'], coverage=1.0, coverage_rightleft=1.0, coverage_updown=0.0) updown_state = SiteStateType('updown', make_single_point_zmatrix(), -1., ['C'], coverage=1.0, coverage_rightleft=0.0, coverage_updown=1.0) empty_state = SiteStateType('Empty',make_single_point_zmatrix(), 0.0, ['H'], coverage=0.0, coverage_rightleft=0.0, coverage_updown=0.0) lt = normal_lattice_task([rightleft_state, updown_state, empty_state], cell, atoms, interaction, r_cut_off, INF_E) for mu in np.arange(-5., 35.+0.0001, 4.0 ): # [-5 ... 35] rightleft_state.props['ads_energy'].value = mu updown_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) avg_cov_rightleft = transferm.average_props(lt.states, N, 'coverage_rightleft', beta) avg_cov_updown = transferm.average_props(lt.states, N, 'coverage_updown', beta) tm_sol = solve_TM(W,N) cov = sum(tm_sol.probs*np.array(avg_cov)) cov_rightleft = sum(tm_sol.probs*np.array(avg_cov_rightleft)) cov_updown = sum(tm_sol.probs*np.array(avg_cov_updown)) covs.append([cov,cov_rightleft,cov_updown]) mus.append(mu) print("mu ==", mu, "theta == ",cov) print("covs", covs) etalon_covs = np.array([[0.99997022769428412, 3.8788036971901128e-09, 0.99997022381547973], [0.99985167717251566, 8.4086284721576516e-08, 0.99985159308622595], [0.97393288643593434, 2.4995574587642413e-05, 0.97390789086134344], [0.54161048563509739, 0.00013659232290069842, 0.5414738933122043], [0.50121987965939618, 4.5059570531890215e-05, 0.50117482008886383], [0.50001327393077732, 6.1471516378165347e-05, 0.49995180241439741], [0.49882667172458006, 5.8926548791756338e-06, 0.49882077906969874], [0.45837522562477745, 4.0315271671395967e-05, 0.45833491035310436], [0.026851860460743324, 0.0012525782414728624, 0.025599282219270339], [0.00027591221106598859, 0.00012808176145751477, 0.00014783044960847382], [5.9512999668565377e-05, 2.9749222196030969e-05, 2.9763777472534411e-05]]) print('max error = ', np.max(abs(covs-etalon_covs))) print('FORCE_GEN_TM=',FORCE_GEN_TM) assert np.allclose(covs, etalon_covs, atol=1e-5) print("TEST OK") if __name__ == '__main__': main()