#!python -u import numpy as np from time import time import os from mpi4py import MPI from math import exp, sqrt,pi from susmost import make_single_point_zmatrix, \ make_square_plane_cell, normal_lattice_task, mc, savexyz, SiteStateType, register_properties from susmost.mc import save_cells_as_xyz from susmost.zmatrix import make_trimer_zmatrix from susmost import phash import random INF_E = 1E6 # prohibitive energy level a = 1.0 # lattice parameter r_cut_off= 6.0 + .0001 precission = 6 k_B = 8.3144598E-3 # kJ/(mol*K) w_HB_perpendicular = -5 w_HB_parallel = w_HB_perpendicular #w_HB_parallel = w_HB_perpendicular *2*sqrt(2)/3.0 w_coord_bond = -20 w_coord_bond_diag = w_coord_bond/sqrt(2.0) w_MM = INF_E w_MM_diag = 0 xyz_fn = 'cur_cells.xyz' cell,atoms = make_square_plane_cell(a) register_properties('tpa_count','fe_count') random.seed(12321) def precover_random(cells, m, cnts): m.recalc_curs() assert m.curE.sum() < INF_E*0.5 s = [] for i,n in cnts.items(): s += [i]*int(n) random.shuffle(s) all_indices = list(range(len(cells)))[:] random.shuffle(all_indices) for i in s: while True: j = random.randint(0, len(cells)-1) old = cells[j] if old not in list(cnts.keys()): cells[j] = i m.recalc_curs() if m.curE.sum() >= INF_E*0.5: # reject cells[j] = old # revert m.recalc_curs() else: break # accept return 0 # Interaction potential for TPA-Fe on Cu(100) system def interaction_function(cc1, cc2): if cc1.lctype == 'Empty' or cc2.lctype == 'Empty': return None # There is no interaction with empty site if cc1.lctype == 'TPA' and cc2.lctype == 'TPA': # If any distance between -COO groups of the neighboring molecules is less than 2a, energy is infinite if any([np.linalg.norm(c1-c2) < 1.99*a for c1 in cc1.coords[1:] for c2 in cc2.coords[1:]]): return INF_E v1 = cc1.coords[2] - cc1.coords[1] # Orientation vector of the first molecule v2 = cc2.coords[2] - cc2.coords[1] # Orientation vector of the second molecule cosa = np.dot(v1,v2) / (np.linalg.norm(v1)* np.linalg.norm(v2)) # Angle between the molecular axis if abs(cosa) > 0.999: # If the molecules are parallel r = np.linalg.norm(cc1.coords[0] - cc2.coords[0]) if r < 2.0*sqrt(2.0)*a*0.999: return INF_E elif r > 3.0*a: return None else: return w_HB_parallel/2.0 elif abs(cosa) < 0.001: # If the molecules are perpendicular r = np.linalg.norm(cc1.coords[0] - cc2.coords[0]) if r > sqrt(10)*a*0.99 and r < sqrt(10)*a*1.01: # Single hydrogen bond O(-)...H-C return w_HB_perpendicular/4.0 elif r > 3.0*a*0.99 and r < 3*a*1.01: # Double hydrogen bond O(-)...H-C return w_HB_perpendicular/2.0 elif r > 3*a*1.01: return None if cc1.lctype == 'Fe' and cc2.lctype == 'TPA': cc1,cc2 = cc2,cc1 if cc1.lctype == 'TPA' and cc2.lctype == 'Fe': # If the Fe-atom is located on the left/right side of the molecule, the energy is infinite # Or any distance between -COO group and Fe atom is less than sqrt(2)*a, then energy is infinite c2 = cc2.coords[0] r = np.linalg.norm(cc1.coords[0] - cc2.coords[0]) if r < sqrt(2.0)*a*1.001: #print 'TPA-Fe', r, 'INF_E' return INF_E elif any([np.linalg.norm(c1-c2) > r and np.linalg.norm(c1-c2) < a*sqrt(5)*1.01 for c1 in cc1.coords[1:]]): return INF_E elif any([np.linalg.norm(c1-c2) < a*1.01 for c1 in cc1.coords[1:]]): #print 'TPA-Fe', r, w_coord_bond/2.0 return w_coord_bond/2.0 elif any([np.linalg.norm(c1-c2) < sqrt(2)*a*1.01 and np.linalg.norm(c1-c2) > a*1.01 for c1 in cc1.coords[1:]]): #print 'TPA-Fe', r, w_coord_bond_diag/2.0 return w_coord_bond_diag/2.0 elif any([np.linalg.norm(c1 - c2) > sqrt(2)*a*1.001 for c1 in cc1.coords[1:]]): #print 'TPA-Fe', r, None return None assert False if cc1.lctype == 'Fe' and cc2.lctype == 'Fe': r = np.linalg.norm(cc1.coords[0] - cc2.coords[0]) if r < 0.99*a: return INF_E elif r < 1.01*a: return w_MM/2.0 elif r < sqrt(2)*a*1.01: return w_MM_diag/2.0 elif r > sqrt(2)*a*1.01: return None return None fe_state = SiteStateType('Fe', make_single_point_zmatrix(), -1., ['Fe'], coverage=1., tpa_count=0., fe_count=1.) # The first sigment of TPA is a central one # S - aromatic ring of the TPA molecule, N is a -COO group tpa_state = SiteStateType('TPA', make_trimer_zmatrix(a,a,pi), -2., ['S','N','N'], coverage=3., tpa_count=1., fe_count=0) empty_state = SiteStateType('Empty',make_single_point_zmatrix(), 0.0, ['H'], coverage=0., tpa_count=0., fe_count=0.) lt = normal_lattice_task([fe_state, tpa_state, empty_state], cell, atoms, interaction_function, r_cut_off, INF_E) assert lt.states_count == 4, lt lt = MPI.COMM_WORLD.bcast(lt, root=0) fe_state_index = lt.find_state_index('fe_count', 1) tpa_state_index = lt.find_state_index('tpa_count', 1) def calc(muTPA, muFe): res_calc = [] xyz_fn = "muFe_" + str(muFe) + ".xyz" lt.set_ads_energy('Fe', muFe) lt.set_ads_energy('TPA', muTPA) T_list = [300.0, 600.0, 900.0, 1200.0] m = mc.make_metropolis(lt, 30, T_list, k_B, param_names=['coverage','tpa_count','fe_count'], seed=12321) m.phys_diffusion = True m.E_inf = 1E6 m.PT_period = 100 # Each m.PT_period steps we try to switch the temperature m.diff_cui_min = 1 # Minimal m.diff_cui_max = 2 m.min_phys_diff_edge_idx = 1 m.max_phys_diff_edge_idx = 8 def callback(m): if MPI.COMM_WORLD.Get_rank() == 0: save_cells_as_xyz(xyz_fn, m) #pass print('m.param_limits', MPI.COMM_WORLD.Get_rank(), m.param_limits[0].enabled, m.param_limits[0].param_index, "[", m.param_limits[0].lower_bound, m.param_limits[0].upper_bound, "]", m.cur_params, m.param_limits[0].check(m.cur_params[m.param_limits[0].param_index] )) #return np.array([1.]) precov_Fe_cnt = m.cells_count * 0.01 precov_TPA_cnt = m.cells_count * 0.05 unplaced = precover_random(m.cells, m, {fe_state_index:precov_Fe_cnt, tpa_state_index:precov_TPA_cnt}) assert unplaced == 0, unplaced save_cells_as_xyz("precov_{}.xyz".format(MPI.COMM_WORLD.Get_rank()), m) m.param_limits = [mc.param_limit(2, precov_Fe_cnt*5.0-0.001, precov_Fe_cnt*30+0.001, False)] log_periods_cnt = 100 # Amount of points for averaging # log_period_steps Amount of trial steps between adding the data to statistics # Total amount of trials equals to log_period_steps * log_periods_cnt stat_log = mc.run(m, log_periods_cnt, log_period_steps = 1000, log_callback=callback, relaxation_steps = 200000, traj_fns=None) xyz_str = save_cells_as_xyz("finish_{}.xyz".format(T_list[MPI.COMM_WORLD.Get_rank()]), m) xyz_hashes = MPI.COMM_WORLD.gather(phash(xyz_str),root=0) if MPI.COMM_WORLD.Get_rank() == 0: etalon_hashes = """c37973d7a5fcd779ba75d9c48b6407cd 6bee3fd86a9860f318fb2c3164c99d54 fcbfa9d328ea4583adfe06fa8b7cd739 a96379a33592db53ce426f2b13e36f38""".split('\n') print("xyz_hashes:\n" + "\n".join(xyz_hashes)) assert sorted(xyz_hashes) == sorted(etalon_hashes) # print "AVGE for T", T, (sum(stat_log_i)/len(stat_log_i))[1] # Print the energy for T, param_log in zip(T_list, m.full_params_log): stat = np.mean(param_log,axis=0)[[0,-2, 1, 2]] # cov E, TPA, Fe res_i = [T, muTPA, muFe, stat] res_calc.append(res_i) if MPI.COMM_WORLD.Get_rank() == 0: print('res_i',res_i) return res_calc #os.remove(xyz_fn) #res = calc(muTPA = 5., muFe = 10.) #print 'RES', res muTPA = -10.0 if __name__ == '__main__': # res = [] # if MPI.COMM_WORLD.Get_rank() == 0: # f = open('results.dat','w') # f.write('T muTPA muFe coverage E TPA_density Fe_density''\n') # f.close() # for muFe in np.arange(60., -60.-0.0001, -5. ): for muFe in [30.]: res_cur = calc(muTPA, muFe) cur_stat = 0 if MPI.COMM_WORLD.Get_rank() == 0: for i, (T, muTPA, muFe, stat) in enumerate(sorted(res_cur)): if abs(T - 300.) < 0.01: # s = ('{}\t'*7).format(T, muTPA, muFe, *stat) # f = open('results.dat','a') # f.write( s + '\n') # f.close() # print "RES:", s cur_stat = stat[2:] # TPA, Fe # print "cur_stat", cur_stat if muFe == 60.: etalon_stat = np.array([0.119938888889, 0.000380555555556]) elif muFe == 30.: etalon_stat = np.array([0.100244444444, 0.0604555555556]) elif muFe == -40.: etalon_stat = np.array([0., 0.5]) print('cur_stat = ', str(cur_stat)) print('etalon_stat = ', str(etalon_stat)) print('error value = ', abs(cur_stat-etalon_stat)) assert np.allclose(cur_stat, etalon_stat, atol=2e-2), abs(cur_stat-etalon_stat) if MPI.COMM_WORLD.Get_rank() == 0: print("TEST OK") # res += res_cur # print res_cov # print res_E