Only in pyoz-0.2.2: _features.txt diff -u pyoz-0.2.1/_notes.txt pyoz-0.2.2/_notes.txt --- pyoz-0.2.1/_notes.txt 2009-05-20 09:42:30.000000000 +0200 +++ pyoz-0.2.2/_notes.txt 2009-07-24 10:33:38.000000000 +0200 @@ -1,9 +1,8 @@ release notes for pyOZ ====================== -released: 20. May 2009 -version: 0.2.1 -bugfix to the version 0.2 +released: 24. July 2009 +version: 0.2.2 check the website http://pyoz.vrbka.net @@ -13,22 +12,24 @@ - theory -changes from the previous release (0.2) - - fixed bug in NR/CG handling causing very slow convergence - - new handling of non-converged NR/CG cycle (ctrl:nr_noconv_picard) +changes from the previous release (0.2.1) + - several bugfixes + - printout of direct correlation functions c(r) + - nr/cg temporarily disabled due to problems + - new file in distribution (_features.txt) revisions/descriptions of individual source files -pyoz.py ........................................... 0.1.5 +pyoz.py ........................................... 0.1.6 main program pyoz_closure.py ................................... 0.1.1 definition of closures -pyoz_const.py ..................................... 0.1.4 +pyoz_const.py ..................................... 0.1.5 definition of constants -pyoz_dft.py ....................................... 0.1.4 +pyoz_dft.py ....................................... 0.1.5 definition of fourier transforms -pyoz_input.py ..................................... 0.1.6 +pyoz_input.py ..................................... 0.1.7 handling of input pyoz_misc.py ...................................... 0.1.3 miscellanous routines diff -u pyoz-0.2.1/pyoz_const.py pyoz-0.2.2/pyoz_const.py --- pyoz-0.2.1/pyoz_const.py 2009-05-20 09:09:50.000000000 +0200 +++ pyoz-0.2.2/pyoz_const.py 2009-07-24 10:38:37.000000000 +0200 @@ -8,14 +8,14 @@ """ Module defining parameters and constants -revision of this file: 0.1.4 +revision of this file: 0.1.5 """ from math import sqrt, pi from scipy import inf -pyoz_version = "0.2.1" +pyoz_version = "0.2.2" class const: """ @@ -229,7 +229,8 @@ gr_suffix = '-gr.dat' ur_suffix = '-ur.dat' urtot_suffix = '-urtot.dat' - + cr_suffix = '-cr.dat' + # Gamma related stuff # write output? G_ij_write = False @@ -258,6 +259,11 @@ # filename Utot_ij_name = name + urtot_suffix + # c(r) related stuff + # write? + c_ij_write = False + # filename + c_ij_name = name + cr_suffix # ************************************************************************************************ diff -u pyoz-0.2.1/pyoz_dft.py pyoz-0.2.2/pyoz_dft.py --- pyoz-0.2.1/pyoz_dft.py 2009-05-20 08:55:31.000000000 +0200 +++ pyoz-0.2.2/pyoz_dft.py 2009-07-24 10:32:19.000000000 +0200 @@ -11,7 +11,7 @@ """ module for performing various types of discrete transforms - revision of this file: 0.1.4 + revision of this file: 0.1.5 for the 3D-Fourier transformation, there is a normalization factor of 1/(2*pi)^3 involved then, the forward FT can be normalized @@ -119,9 +119,12 @@ # ********************************************************************************************** # let's define the normalization constants as described above -# choose always only one combination # these are private +# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +# ! do not change the definitions since the program relies on the selection made here ! +# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + # 3D Fourier transform # normalization within the forward transformation #_ft_3d_forward = 1/(2.0 * pi)**3 @@ -180,6 +183,7 @@ # it's known that there is an additional constant factor involved, equal to # the reciprocal of the forward fourier transform normalization factor # we need to access this factor from the solver routine + # due to the chosen parameters, this is equal 1.0 at the moment ft_convolution_factor = 1.0/_ft_3d_forward def __init__(self, npoints, dr, dk, r, k): diff -u pyoz-0.2.1/pyoz_input.py pyoz-0.2.2/pyoz_input.py --- pyoz-0.2.1/pyoz_input.py 2009-05-20 09:10:47.000000000 +0200 +++ pyoz-0.2.2/pyoz_input.py 2009-07-24 09:17:45.000000000 +0200 @@ -8,7 +8,7 @@ """ Module with functions for input file parsing -revision of this file: 0.1.6 +revision of this file: 0.1.7 """ # import external modules @@ -30,6 +30,27 @@ from pyoz_closure import supported_closures # ********************************************************************************************** +# classes for exception handling - needed to disable deprecation warning for python > 2.4 + +class badValue(Exception): + def __init__(self, value): + self.value = value + def __str__(self): + return repr(self.value) + +class unkKeyword(Exception): + def __init__(self, value): + self.value = value + def __str__(self): + return repr(self.value) + +class pmfError(Exception): + def __init__(self, value): + self.value = value + def __str__(self): + return repr(self.value) + +# ********************************************************************************************** def parse_cmdline(argv): """usage information @@ -304,63 +325,66 @@ graphics_p_xmax = float((l.split())[1]) elif (match(r'npoints$', keyword)): npoints = int((l.split())[1]) - if (npoints % 2): - raise "badValue", (i, "npoints not a power of 2") - # the input should be a power of 2 # we need to then subtract one because of the particular fourier transform used # the last point will be added by the transform itself npoints -= 1 + if (npoints <= 0): + raise badValue((i, "cannot use zero/negative number of discretization points!")) elif (match(r'deltar$', keyword)): deltar = float((l.split())[1]) if (deltar <= 0.0): - raise "badValue", (i, "step size must be larger than 0") + raise badValue((i, "step size must be larger than 0")) elif (match(r'mix_param$', keyword)): mix_param = float((l.split())[1]) if ((mix_param <= 0) or (mix_param > 1)): - raise "badValue", (i, "mixing parameter not in allowed range 0 < mix_param ≤ 1") + raise badValue((i, "mixing parameter not in allowed range 0 < mix_param ≤ 1")) elif (match(r'conv_crit$', keyword)): conv_crit = float((l.split())[1]) if (conv_crit <= 0.0): - raise "badValue", (i, "convergence criterion has to be larger than 0") + raise badValue((i, "convergence criterion has to be larger than 0")) + if (conv_crit >= 1.0): + raise badValue((i, "convergence criterion has to be smaller than 1")) elif (match(r'max_iter$', keyword)): max_iter = int((l.split())[1]) if (max_iter <= 0): - raise "badValue", (i, "minimum number of 1 iterations has to be performed") + raise badValue((i, "minimum number of 1 iterations has to be performed")) elif (match(r'max_dsqn$', keyword)): max_dsqn = float((l.split())[1]) if (max_dsqn <= 0): - raise "badValue", (i, "only positive dsqn values are allowed") + raise badValue((i, "only positive dsqn values are allowed")) elif (match(r'use_nr$', keyword)): do_nr = True + print("NR/CG method has been disabled in this version due to severe problems!") + sys.exit(1) elif (match(r'nr_mix_param$', keyword)): nr_mix_param = float((l.split())[1]) if ((nr_mix_param < 0) or (nr_mix_param > 1)): - raise "badValue", (i, "NR mixing parameter not in allowed range 0 ≤ nr_mix_param ≤ 1") + raise badValue((i, "NR mixing parameter not in allowed range 0 ≤ nr_mix_param ≤ 1")) elif (match(r'nr_conv_crit$', keyword)): nr_convergence_factor = float((l.split())[1]) if (nr_convergence_factor <= 0.0): - raise "badValue", (i, "NR relative convergence criterion has to be larger than 0") + raise badValue((i, "NR relative convergence criterion has to be larger than 0")) elif (match(r'nr_max_iter$', keyword)): nr_max_iter = int((l.split())[1]) if (nr_max_iter <= 0): - raise "badValue", (i, "minimum number of 1 NR iterations has to be performed") + raise badValue((i, "minimum number of 1 NR iterations has to be performed")) elif (match(r'nr_noconv_incr$', keyword)): nr_noconv_incr = True else: - raise "UnkKeyword", (i) + raise unkKeyword((i, l)) # end for i in range(len(lines)): # end try except ValueError, msg: print("error processing line %lu (%s)" % (i + lineno_beg, l)) print(msg) sys.exit(2) - except "UnkKeyword", msg: - print("error processing line %lu (%s)" % (msg + lineno_beg, l)) + except unkKeyword, msg: + print("error processing line %lu (%s)" % (msg.value[0] + lineno_beg, l)) print("unknown keyword") sys.exit(2) - except "BadValue", msg: - print ("error processing line %lu (%s)" % (msg[0] + lineno_beg, l)) - print (msg[1]) + except badValue, msg: + print ("error processing line %lu (%s)" % (msg.value[0] + lineno_beg, l)) + print (msg.value[1]) sys.exit(2) except IndexError, msg: print("error processing line %lu (%s)" % (i + lineno_beg, l)) @@ -374,6 +398,11 @@ ctrl['npoints_exp'] = npoints_exp ctrl['npoints'] = npoints print("\tnumber of points\t\t%u" % (ctrl['npoints'] + 1)) + if (((npoints + 1) & npoints) != 0): + # logical AND + # the value was already decremented! + # the input should be a power of 2 + print("\t!!!!! you should use number of points which is a power of 2 !!!!!") # recalculate step in fourier space, and maximum distances in r and k max_r = deltar * npoints @@ -393,7 +422,7 @@ # mixing parameter ctrl['mix_param'] = mix_param # convergence criterion - ctrl['convergence_crit'] = convergence_crit + ctrl['convergence_crit'] = conv_crit # maximum value of DSQN (convergence criterion) ctrl['max_dsqn'] = max_dsqn # maximum number of iterations @@ -487,15 +516,15 @@ elif (match(r'temp$', keyword)): temp = float((l.split())[1]) if (temp < 0.0): - raise "badValue", (i, "negative temperature not allowed") + raise badValue((i, "negative temperature not allowed")) elif (match(r'epsilon_r$', keyword)): epsilon_r = float((l.split())[1]) if (epsilon_r <= 0.0): - raise "badValue", (i, "epsilon_r has to be larger than zero") + raise badValue((i, "epsilon_r has to be larger than zero")) elif (match(r'ncomp$', keyword)): ncomponents = int((l.split())[1]) if (ncomponents <= 0): - raise "badValue", (i, "number of components has to be larger than 0") + raise badValue((i, "number of components has to be larger than 0")) elif (match(r'closure$', keyword)): closure_name = (l.split())[1] elif (match(r'names$', keyword)): @@ -505,21 +534,21 @@ elif (match(r'conc_unit$', keyword)): conc_unit = (l.split())[1] else: - raise "UnkKeyword", (i) + raise unkKeyword((i, l)) # end for i in range(len(lines)): # end try except ValueError, msg: print("error processing line %lu (%s)" % (i + lineno_beg, l)) print(msg) sys.exit(2) - except "BadValue", msg: - print ("error processing line %lu (%s)" % (msg[0] + lineno_beg, l)) - print (msg[1]) - sys.exit(2) - except "UnkKeyword", msg: - print("error processing line %lu (%s)" % (i + lineno_beg, l)) + except unkKeyword, msg: + print("error processing line %lu (%s)" % (msg.value[0] + lineno_beg, l)) print("unknown keyword") sys.exit(2) + except badValue, msg: + print ("error processing line %lu (%s)" % (msg.value[0] + lineno_beg, l)) + print (msg.value[1]) + sys.exit(2) # end try/except block # set the temperature @@ -569,8 +598,9 @@ # newton-raphson can be used only with the HNC closure! if (ctrl['do_nr'] and syst['closure_name'] != 'hnc'): - print("\tillegal combination of Newton-Raphson and closure (only HNC supported!)") - sys.exit(1) + #print("\tillegal combination of Newton-Raphson and closure (only HNC supported!)") + print("\tsupport of Newton-Raphson with PY closure is highly experimental!") + #sys.exit(1) # names of the constituents # if not given, use A, B, ... @@ -835,7 +865,7 @@ hs_unit = (l.split())[1] hs_unit_orig = (orig_l.split())[1] else: - raise "UnkKeyword", (i) + raise unkKeyword((i,l)) # end if (match...) # end try except ValueError, msg: @@ -846,8 +876,8 @@ print("line %lu with input data too short: %s" % (i + lineno_pot, l)) print(msg) sys.exit(2) - except "UnkKeyword", msg: - print("error processing line %lu (%s)" % (i + lineno_pot, l)) + except unkKeyword, msg: + print("error processing line %lu (%s)" % (msg.value[0] + lineno_beg, l)) print("unknown keyword") sys.exit(2) # end try/except block @@ -904,6 +934,7 @@ #rest,result = modf(hs_ij[i,j] / ctrl['deltar']) if (hs_ij[i,j] != ctrl['deltar'] * (dis_index[i,j] + 1)): print("\tWARNING: hard sphere diameter (combination %u,%u)\n\t\t\t%f not an integer multiple of step size\n\t\t\t%f used instead!" % ((i + 1), (j + 1), hs_ij[i,j], ctrl['deltar'] * (dis_index[i,j] + 1))) + print("\tWARNING: detection not 100%, but functionality of program not affected!") #else: # print("good") # increase the combination number @@ -1017,7 +1048,7 @@ dip_alpha[:] = float_(array( (l.split())[ 1:syst['ncomponents']+1 ] )) dip_alpha_loaded = True else: - raise "UnkKeyword", (i) + raise unkKeyword((i,l)) # end if (match...) # end try except ValueError, msg: @@ -1028,8 +1059,8 @@ print("line %lu with input data too short: %s" % (i + lineno_pot, l)) print(msg) sys.exit(2) - except "UnkKeyword", msg: - print("error processing line %lu (%s)" % (i + lineno_pot, l)) + except unkKeyword, msg: + print("error processing line %lu (%s)" % (msg.value[0] + lineno_beg, l)) print("unknown keyword") sys.exit(2) # end try/except block @@ -1252,7 +1283,7 @@ comb_rule['epsilon'] = (l.split())[1] epsilon_rule_orig = (orig_l.split())[1] else: - raise "UnkKeyword", (i) + raise unkKeyword((i,l)) # end if (match...) # end try except ValueError, msg: @@ -1263,8 +1294,8 @@ print("line %lu with input values too short: %s" % (i + lineno_pot, l)) print(msg) sys.exit(2) - except "UnkKeyword", msg: - print("error processing line %lu (%s)" % (i + lineno_pot, l)) + except unkKeyword, msg: + print("error processing line %lu (%s)" % (msg.value[0] + lineno_beg, l)) print("unknown keyword") sys.exit(2) # end try/except block @@ -1462,7 +1493,7 @@ elif (match(r'pmf_after$', keyword)): pmf_after = float((l.split())[1]) else: - raise "UnkKeyword", (i) + raise unkKeyword((i,l)) # end if (match...) # end try except IndexError, msg: @@ -1473,8 +1504,8 @@ print("error converting numerical input data on line %lu (%s)" % (i + lineno_pot, l)) print(msg) sys.exit(2) - except "UnkKeyword", msg: - print("error processing line %lu (%s)" % (i + lineno_pot, l)) + except unkKeyword, msg: + print("error processing line %lu (%s)" % (msg.value[0] + lineno_beg, l)) print("unknown keyword") sys.exit(2) # end try/except block @@ -1582,7 +1613,7 @@ if (len(datasets) < 2): # we need at least 2 points - raise "pmfError", ("need at least 2 points, got %u " % len(datasets[0])) + raise pmfError(("need at least 2 points, got %u " % len(datasets[0]))) fr.close() # pre-allocate arrays # distance @@ -1617,7 +1648,7 @@ if (len(datasets[i]) < 2): # we need at least 2 points - raise "pmfError", ("need at least 2 points, got %u " % len(datasets[i])) + raise pmfError(("need at least 2 points, got %u " % len(datasets[i]))) fr.close() # modify the arrays to correspond to the length of the pmf # distance @@ -1673,7 +1704,7 @@ sline = line.split() if (len(sline) < num_val): # something is not correct - raise "pmfError", ("expected %u values, got %u: %s" % (num_val, len(sline), line)) + raise pmfError(("expected %u values, got %u: %s" % (num_val, len(sline), line))) # end if (len(sline) != num_val) elif (len(sline) > num_val): # skip extra data but issue a warning @@ -1695,7 +1726,7 @@ # the first r value has to be the smallest # and it must increase among the datasets if ((pmf_r[i][index] <= pmf_rmin[i]) or (pmf_r[i][index] <= pmf_rmax[i])): - raise "pmfError", ("datasets have to be sorted in increasing order: around line %s" % line) + raise pmfError(("datasets have to be sorted in increasing order: around line %s" % line)) else: pmf_rmax[i] = pmf_r[i][index] # end if (index == 0) - testing for correctness @@ -1710,8 +1741,8 @@ except IOError, errmsg: print("load of PMF failed - IO error: %s" % errmsg) sys.exit(2) - except "pmfError", errmsg: - print("load of PMF failed: %s" % errmsg) + except pmfError, msg: + print("load of PMF failed: %s" % msg.value) sys.exit(2) # end except pmfError - PMF-related problems except ValueError, errmsg: @@ -1724,11 +1755,6 @@ # print some statistics? # end try/except/else block -# for i in range(3): -# for index in range(len(pmf_r[i])): -# print("val\t\t\tind=%u\tr=%f\tpmf=%f" % (i, pmf_r[i][index]/pmfdist_conversion, pmf_ij[i][index]/pmf_conversion)) -# print("") - # store the values parm.append({}) parm[-1]['pot_type'] = 'pmf' @@ -1790,11 +1816,15 @@ # pair potential + Gamma Utot_ij_write = outp_def.Utot_ij_write Utot_ij_name = outp_def.Utot_ij_name + # direct correlatio function + c_ij_write = outp_def.c_ij_write + c_ij_name = outp_def.c_ij_name # apply correction according to the commandline if (cmdline['name'] != None): G_ij_name = cmdline['name'] + outp_def.gamma_suffix g_ij_name = cmdline['name'] + outp_def.gr_suffix + c_ij_name = cmdline['name'] + outp_def.cr_suffix U_ij_name = cmdline['name'] + outp_def.ur_suffix Utot_ij_name = cmdline['name'] + outp_def.urtot_suffix @@ -1832,7 +1862,7 @@ G_ij_freq = int((l.split())[1]) G_ij_write = True if (G_ij_freq < 0): - raise "badValue", (i, "the Gamma save frequency cannot be negative") + raise badValue((i, "the Gamma save frequency cannot be negative")) elif (match(r'gamma_binary$', keyword)): G_ij_binary = True elif (match(r'gamma_text$', keyword)): @@ -1846,10 +1876,6 @@ elif (match(r'gr_name$', keyword)): g_ij_write = True g_ij_name = (orig_l.split())[1] - #elif (match(r'gr_freq$', keyword)): - # g_ij_freq = int((l.split())[1]) - # if (g_ij_freq < 0): - # raise "badValue", (i, "the g(r) save frequency cannot be negative") # ************************************************************************** # u(r) - total potential @@ -1866,21 +1892,31 @@ Utot_ij_name = (orig_l.split())[1] # ************************************************************************** + # c(r) + + elif (match(r'cr_write$', keyword)): + c_ij_write = True + elif (match(r'cr_name$', keyword)): + c_ij_write = True + c_ij_name = (orig_l.split())[1] + + # ************************************************************************** + else: - raise "UnkKeyword", (i) + raise unkKeyword((i,l)) # end for i in range(len(lines)): # end try except ValueError, msg: print("error processing line %lu (%s)" % (i + lineno_beg, l)) print(msg) sys.exit(2) - except "UnkKeyword", msg: - print("error processing line %lu (%s)" % (msg + lineno_beg, l)) + except unkKeyword, msg: + print("error processing line %lu (%s)" % (msg.value[0] + lineno_beg, l)) print("unknown keyword") sys.exit(2) - except "BadValue", msg: - print ("error processing line %lu (%s)" % (msg[0] + lineno_beg, l)) - print (msg[1]) + except badValue, msg: + print ("error processing line %lu (%s)" % (msg.value[0] + lineno_beg, l)) + print (msg.value[1]) sys.exit(2) except IndexError, msg: print("error processing line %lu (%s)" % (i + lineno_beg, l)) @@ -1918,7 +1954,15 @@ print("\tg(r) will be saved to file %s" % (outp['g_ij_name'])) else: print("\tg(r) will not be written!") - + + # c(r) related stuff + outp['c_ij_write'] = c_ij_write + outp['c_ij_name'] = c_ij_name + if (outp['c_ij_write']): + print("\tc(r) will be saved to file %s" % (outp['c_ij_name'])) + #else: + # print("\tc(r) will not be written!") + # U(r) related stuff outp['U_ij_write'] = U_ij_write outp['U_ij_name'] = U_ij_name diff -u pyoz-0.2.1/pyoz.py pyoz-0.2.2/pyoz.py --- pyoz-0.2.1/pyoz.py 2009-05-20 11:36:40.000000000 +0200 +++ pyoz-0.2.2/pyoz.py 2009-07-24 09:39:20.000000000 +0200 @@ -9,7 +9,7 @@ """ Program for numerical solution of the Ornstein-Zernike equation version of the program is defined in pyoz_const.py -revision of this file: 0.1.5 +revision of this file: 0.1.6 check the file _features.txt and website http://pyoz.vrbka.net for further information @@ -122,38 +122,35 @@ # calculate the erf-correction contribution exp(U_erf_ij) # store all in a dictionary modMayerFunc modMayerFunc = potential.def_modMayerFunc(syst, U_ij, U_ij_individual, U_discontinuity, U_erf_ij['real']) + # store the mayer function itself for the purpose of CG procedure with PY closure + M = modMayerFunc['u_ij'] - 1.0 # allocate arrays with direct, total and pair corr. functions, Gamma function - # it would be possible to work with smaller number of arrays - # but it would certainly obscure the code - # also, no memory would be probably saved, since x += 1 creates a copy of - # x anyway # some arrays will emerge from arithmetic operations, - # but defining them here makes the code clearer + # showing them here makes the code clearer # real space: _r_ # Fourier space _f_ # direct correlation function without (c) and with (C) density factor applied # direct correlation function with Ng-correction (cs) applied # w/o density correction in real space; w/ density correction in Fourier space - c_r_ij = zeros((syst['ncomponents'], syst['ncomponents'], ctrl['npoints'])) - c_f_ij = zeros((syst['ncomponents'], syst['ncomponents'], ctrl['npoints'])) + #c_r_ij C_f_ij = zeros((syst['ncomponents'], syst['ncomponents'], ctrl['npoints'])) - cs_r_ij = zeros((syst['ncomponents'], syst['ncomponents'], ctrl['npoints'])) Cs_f_ij = zeros((syst['ncomponents'], syst['ncomponents'], ctrl['npoints'])) + #cs_r_ij # pair correlation function - g_r_ij = zeros((syst['ncomponents'], syst['ncomponents'], ctrl['npoints'])) + #g_r_ij # h in Fourier space without (h) and with (H) density factor applied - # h not needed - it is later created from H - H_f_ij = zeros((syst['ncomponents'], syst['ncomponents'], ctrl['npoints'])) - h_f_ij = zeros((syst['ncomponents'], syst['ncomponents'], ctrl['npoints'])) + #H_f_ij + # matrix of partial structure factors + #S # actual, old (o) and new (n) values for Gamma # these are short-ranged Gamma (see the code for more details) G_r_ij = zeros((syst['ncomponents'], syst['ncomponents'], ctrl['npoints'])) - G_o_ij = zeros((syst['ncomponents'], syst['ncomponents'], ctrl['npoints'])) + #G_o_ij G_n_ij = zeros((syst['ncomponents'], syst['ncomponents'], ctrl['npoints'])) G_f_ij = zeros((syst['ncomponents'], syst['ncomponents'], ctrl['npoints'])) + # identity matrix for the solver (dr copies!) - # ones matrix for the solver (dr copies!) e_ij = eye(syst['ncomponents']) #o_ij = ones((syst['ncomponents'],syst['ncomponents'])) E_ij = zeros((syst['ncomponents'], syst['ncomponents'], ctrl['npoints'])) @@ -164,7 +161,7 @@ # zero array for the newton-raphson Z = zeros((syst['ncomponents'], syst['ncomponents'], ctrl['npoints'])) # arrays for the Newton-Raphson procedure (allocated even if not used) - HXq = empty((syst['ncomponents'], syst['ncomponents'], ctrl['npoints'])) + CFXq = empty((syst['ncomponents'], syst['ncomponents'], ctrl['npoints'])) AX = empty((syst['ncomponents'], syst['ncomponents'], ctrl['npoints'])) AXq = empty((syst['ncomponents'], syst['ncomponents'], ctrl['npoints'])) Rq = empty((syst['ncomponents'], syst['ncomponents'], ctrl['npoints'])) @@ -296,15 +293,18 @@ # H - aCH = C # (E - aC)H = C # H = {E - aC}^-1 * C + # however, thanks to the normalization chosen so that for FT it is 1, we can write + # H = {E - C}^-1 * C + # E + H = {E - C}^-1 * (C + E - C) + # S = {E - C}^-1 * E # for some cases (1,2 components), own routines are used - H_f_ij = solver_function((E_ij - dft.ft_convolution_factor * C_f_ij), C_f_ij, ctrl['npoints']) - # equivalent formalism involves the matrix of structure factors S = E + H(k) where 1 is identity matrix and h(k) contains the density prefactors - # the equation to solve is then in the continuous case S = 1/(1-C) - # in general case, we have extra constant factors involved in the convolution, and therefore also in (1-aC) => the situation is not that simple - # the S-matrix will be needed later (nr procedure) + #H_f_ij = solver_function((E_ij - dft.ft_convolution_factor * C_f_ij), C_f_ij, ctrl['npoints']) + S = solver_function((E_ij - C_f_ij), E_ij, ctrl['npoints']) # convert H to short ranged Gamma G(k) = H(k) - Cs(k) - G_f_ij = H_f_ij - Cs_f_ij + #G_f_ij = H_f_ij - Cs_f_ij + # convert S to short ranged Gamma G(k) = S(k) - E - Cs(k) + G_f_ij = S - E_ij - Cs_f_ij # FT G_f_ij to G_r_ij #print("\tiFT Gamma_f -> Gamma_r") @@ -380,25 +380,32 @@ # we don't need Wcur here, but let's define it... Wcur = 1.0 - # we will also need matrices S and H for the operator A (defined in the nr-cycle) + # we will also need matrices S (both HNC and PY) and H (HNC) or M (PY - Mayer function) for the operator A (defined in the nr-cycle) # we use functions from the main cycle - # S = 1 + H(k) # it's not possible to have 3D-matrix => we do it with arrays here and convert to matrices where needed and appropriate # - # S is defined with H(k) calculated with forward FT having the normalization constant of 1.0 - # the 1/(2pi)**3 norm. is then done in the inverse transformation - # ft_convolution factor is reciprocal of the forward normalization used in this program -> we just need to multiply - # H(k) with it to get the needed result - # tested with 1.0|(2pi)**3 and (2pi)**3|1.0 pairs - S = E_ij + dft.ft_convolution_factor*H_f_ij # G = H - C => H = G + C = Gs + Ucorr + Cs - Ucorr # in the algorithm both G and cs are short-ranged => the corrections cancel and provide the correct H # we have to take the H using the old Gamma! i.e., taking g(r) provided by the closure shortly after the main iteration # cycle is started and subtracting 1 H = g_r_ij - 1.0 + # modMayerFunc is MayerFunc + 1.0 + # this has been done before, so commenting out + #M = modMayerFunc['u_ij'] - 1.0 # the arrays for the operations involving operators A and At were created during the initialization - + + # the code will operate with the matrix CF (closure factor), which is set according to closure + # to either H (total correlation function, for HNC) or to M (mayer function, for PY) + # check below for the algorithm details + if (syst['closure_name'] == 'hnc'): + CF = H.copy() + elif (syst['closure_name'] == 'py'): + CF = M + else: + sys.stderr.write("unsupported closure! \n") + sys.exit(1) + # timing purposes nr_time_beg = time() @@ -432,18 +439,26 @@ # for matricial relations (part of operators A and At) we need to do everything separately for each discretization step # perform the calculation of AX # this will be done in several steps since + # !!!!!!!!!!!!!!!!!!! in HNC !!!!!!!!!!!!!!!!!!! # AX = 1X - iFT ( S FT(HX) S - FT(HX)) - + # !!!!!!!!!!!!!!!!!!! in PY !!!!!!!!!!!!!!!!!!! + # AX = 1X - iFT ( S FT(MX) S - FT(MX)) + # where M is the Mayer function exp(-betaU) - 1 + + # the code will operate with the matrix CF (closure factor), which is set according to closure + # to either H (total correlation function, for HNC) or to M (mayer function, for PY) + # this was done outside of this cycle + # H.X is not matricial product! - HX = H*Xcur + CFX = CF*Xcur for i in range(syst['ncomponents']): for j in range(syst['ncomponents']): - HXq[i,j] = dft.dfbt(HX[i,j])[0] + CFXq[i,j] = dft.dfbt(CFX[i,j])[0] # matricial products here for dr in range(ctrl['npoints']): - AXq[:,:,dr] = mat(S[:,:,dr])*mat(HXq[:,:,dr])*mat(S[:,:,dr]) - mat(HXq[:,:,dr]) + AXq[:,:,dr] = mat(S[:,:,dr])*mat(CFXq[:,:,dr])*mat(S[:,:,dr]) - mat(CFXq[:,:,dr]) for i in range(syst['ncomponents']): for j in range(syst['ncomponents']): @@ -480,7 +495,12 @@ # perform the calculation of AtR # this will be done in several steps since and At is an adjoint of the operator A + # !!!!! in HNC !!!!! # AT R = 1R - FT (S iFT(R) S - iFT(R))H + # !!!!! in PY !!!!! + # AT R = 1R - FT (S iFT(R) S - iFT(R))M + + # the operator works with the matrix CF, set according to the used closure for i in range(syst['ncomponents']): for j in range(syst['ncomponents']): Rq[i,j] = dft.dfbt(Rcur[i,j])[0] @@ -495,7 +515,7 @@ SRSq[i,j] = dft.idfbt(SRS[i,j]) # not a matricial product! - AtR = Rcur - SRSq*H + AtR = Rcur - SRSq*CF # calculate Lcur = alpha(n) = (R(n),R(n))/(AtR(n),AtR(n)) # where (Y,Z) is inner product \sum_ij rho_i rho_j \int Y_ij Z_ij 4 \pi r^2 dr @@ -575,7 +595,9 @@ # do closure cs_r_ij, g_r_ij = syst['closure'](syst, r, modMayerFunc, U_discontinuity, G_r_ij) - + # and evaluate uncorrected c(r) as well + c_r_ij = cs_r_ij - U_erf_ij['real'] + # update the plot if requested if (ctrl['do_graphics']): pyoz_plot.plot_update(syst, const, U_r=None, U_erf=None, G_r=G_r_ij, c_r=cs_r_ij, g_r=g_r_ij, c_f=None) @@ -606,6 +628,27 @@ fw.write("\n") fw.close() + # save c_r_ij to file + if (outp['c_ij_write']): + print("\tdirect correlation function\t(%s)" % outp['c_ij_name']) + fw = open(outp['c_ij_name'], "wt") + fw.write("%8.3f" % 0.0) + for i in range(syst['ncomponents']): + for j in range(i, syst['ncomponents']): + fw.write("%10.5f" % 0.0) + fw.write("\n") + for dr in range(ctrl['npoints']): + fw.write("%8.3f" % r[dr]) + for i in range(syst['ncomponents']): + for j in range(i, syst['ncomponents']): + # careful, we have to write the complete c(r), i.e., we need to compensate for the + # Ng-correction! + # cs(r) = c(r) + Ucorr(r) => c(r) = cs(r) - Ucorr(r) + fw.write("%10.5f" % c_r_ij[i, j, dr]) + # end for i,j in range ncomponents... + fw.write("\n") + fw.close() + if (outp['G_ij_write']): print("\tGamma function\t\t\t(%s)" % outp['G_ij_name']) # store the Gamma function if required @@ -665,13 +708,14 @@ index = i if (index >= 0): print("\t\tfound, using short-ranged c(r)\n") - # firstly, we get rid of ng, (subtract U_erf_ij) converting cs_r_ij to c_r_ij + # ----- was done already! firstly, we get rid of ng, (subtract U_erf_ij) converting cs_r_ij to c_r_ij # then we add the coulomb (lr correction) as shown above, to get c_r_ij_sr - c_r_ij_sr = cs_r_ij -U_erf_ij['real'] + U_ij_individual[index] + #c_r_ij_sr = cs_r_ij -U_erf_ij['real'] + U_ij_individual[index] + c_r_ij_sr = c_r_ij + U_ij_individual[index] else: print("\t\tnot found, using original c(r)\n") # in this case, ng-correction is zero and cs_r_ij is already the function we need - c_r_ij_sr = cs_r_ij + c_r_ij_sr = c_r_ij # kirkwood-buff integrals properties.kirkwood_buff(ctrl, syst, r, g_r_ij) Common subdirectories: pyoz-0.2.1/tests and pyoz-0.2.2/tests