#!python -u import numpy as np from math import log from susmost import make_square_plane_cell, make_single_point_zmatrix, \ normal_lattice_task, tensors, SiteStateType, nearest_int, solve_TM import os def main(): #constants INF_E = 1E6 # prohibitive energy level k_B = 1.0 a = 1.0 # lattice parameter r_cut_off= 8.0 + .0001 precission = 6 #parameters J = 1.0 dop_int = {('up','up'):-J/2,('down','down'):-J/2,('up','down'):+J/2,('down','up'):+J/2,} cell,atoms = make_square_plane_cell(a) up_state = SiteStateType('up', make_single_point_zmatrix(), 0.0, ['N']) down_state = SiteStateType('down',make_single_point_zmatrix(), 0.0, ['O']) interaction = lambda cc1,cc2: nearest_int(cc1,cc2,dop_int=dop_int) N = 8 #other Tcrit_exact = (2./log(1. + 2.**0.5)) * J / k_B omegas = [] M = [] Ts = [] theor_M = [] FORCE_GEN_TM = os.getenv('FORCE_GEN_TM', 'False' ) != 'False' lt = normal_lattice_task( [up_state, down_state], cell, atoms, interaction, r_cut_off, INF_E) assert lt.states_count == 2, lt pair_energies = tensors.calc_pair_energies(lt) for T in np.arange(2.4, 2.1, -0.01): beta = 1./(k_B*T) W = tensors.make_tensor_from_pair_energies(pair_energies, lt.states, beta, force_generic=FORCE_GEN_TM) tm_sol = solve_TM(W, N, symmetric=True) omega_tm = -log(tm_sol.lr)/beta omegas.append(omega_tm) Ts.append(T) M.append(tm_sol.pair_corr[0,0] + tm_sol.pair_corr[1,1]) theor_M.append( (Tcrit_exact - T)/Tcrit_exact ) #print 'T = ', T #print 'pair_corr:',tm_sol.pair_corr etalon_omegas = [-17.312521443256731, -17.282591191166375, -17.253036452687468, -17.223862874013804, -17.195075947016967, -17.166680987051809, -17.138683110388651, -17.111087211446048, -17.083897940014982, -17.057119678680518, -17.030756520659441, -17.004812248282295, -16.979290312354077, -16.954193812630269, -16.929525479642365, -16.905287658100335, -16.881482292087973, -16.858110912250329, -16.835174625151872, -16.812674104958543, -16.790609587567811, -16.768980867278863, -16.747787296060036, -16.72702778543453, -16.706700810968666, -16.686804419309617, -16.667336237684665, -16.64829348574002, -16.629672989566714, -16.611471197734275] etalon_M = [0.81793751768174694, 0.8207545076915479, 0.82360030653939487, 0.82647334727555066, 0.82937191710095082, 0.83229415799364692, 0.83523806866248895, 0.83820150791229009, 0.84118219948692741, 0.84417773843556954, 0.84718559902275259, 0.8502031441757536, 0.85322763643323307, 0.85625625032813102, 0.85928608610618151, 0.86231418465010612, 0.8653375434496019, 0.868353133429699, 0.87135791642596128, 0.8743488630753129, 0.87732297087680033, 0.88027728216802892, 0.88320890176075872, 0.88611501398343562, 0.88899289888919508, 0.89183994740478623, 0.89465367521838512, 0.89743173523159359, 0.90017192843209459, 0.9028722130773984] omega = np.array(omegas[:]) d_omega = omega[2:] - omega[:-2] d2_omega = d_omega[2:] - d_omega[:-2] print('\n'.join(map(str, list(zip(Ts, omegas, M, theor_M)) ))) Tc_idx = np.argmin(d2_omega) Tcrit = (Ts[2:-2])[Tc_idx] print("Tcrit = ", Tcrit, Tc_idx) print("Tcrit_exact = ", Tcrit_exact) print('FORCE_GEN_TM=',FORCE_GEN_TM) assert abs(Tcrit - 2.21) < 1e-10 assert np.allclose(omegas, etalon_omegas, atol=1e-9), (np.sort(abs(omegas - etalon_omegas)), omegas) assert np.allclose(M, etalon_M, atol=1e-9), (np.sort(abs(M - etalon_M)), M) print("TEST OK") if __name__ == '__main__': main()