Only in pyoz-0.2.2: _features.txt diff -u pyoz-0.2/_notes.txt pyoz-0.2.2/_notes.txt --- pyoz-0.2/_notes.txt 2009-05-14 11:08:33.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: 14. May 2009 -version: 0.2 -second release +released: 24. July 2009 +version: 0.2.2 check the website http://pyoz.vrbka.net @@ -13,27 +12,24 @@ - theory -changes from the previous release (0.1) - - simplified output - - improved convergence handling - - Newton-Raphson/conjugate gradients algorithm (for HNC) - - simplified handling of concentrations/densities (normalization now done within FT) - - further code cleanup, speed optimization - - fixed output of the pair potential (overflow for large numbers) - - new option - writing U(r) + gamma - total interaction potential +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.4 +pyoz.py ........................................... 0.1.6 main program pyoz_closure.py ................................... 0.1.1 definition of closures -pyoz_const.py ..................................... 0.1.3 +pyoz_const.py ..................................... 0.1.5 definition of constants -pyoz_dft.py ....................................... 0.1.3 +pyoz_dft.py ....................................... 0.1.5 definition of fourier transforms -pyoz_input.py ..................................... 0.1.5 +pyoz_input.py ..................................... 0.1.7 handling of input pyoz_misc.py ...................................... 0.1.3 miscellanous routines @@ -45,3 +41,4 @@ evaluation of thermodynamic properties pyoz_solver.py .................................... 0.1.1 optimized solvers for special cases + diff -u pyoz-0.2/pyoz_const.py pyoz-0.2.2/pyoz_const.py --- pyoz-0.2/pyoz_const.py 2009-05-08 16:58:32.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.3 +revision of this file: 0.1.5 """ from math import sqrt, pi from scipy import inf -pyoz_version = "0.2" +pyoz_version = "0.2.2" class const: """ @@ -152,7 +152,11 @@ nr_convergence_factor = 0.1 # max number of iterations in the inner NR/CG cycle - nr_max_iter = 50 + nr_max_iter = 10 + + # use non-converged dgam as increment if nr_max_iter has been reached? + # if false, the Picard iteration will be used instead + nr_noconv_incr = False # ************************************************************************************************ @@ -225,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 @@ -254,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/pyoz_dft.py pyoz-0.2.2/pyoz_dft.py --- pyoz-0.2/pyoz_dft.py 2009-05-08 23:28:29.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.3 + 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,19 +119,22 @@ # ********************************************************************************************** # 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 -_ft_3d_inverse = 1.0 +#_ft_3d_forward = 1/(2.0 * pi)**3 +#_ft_3d_inverse = 1.0 # unitary #_ft_3d_forward = sqrt(1/(2.0 * pi)**3) #_ft_3d_inverse = sqrt(1/(2.0 * pi)**3) # normalization within the inverse transformation -#_ft_3d_forward = 1.0 -#_ft_3d_inverse = 1/(2.0 * pi)**3 +_ft_3d_forward = 1.0 +_ft_3d_inverse = 1/(2.0 * pi)**3 # 3D FT -> Fourier-Bessel transform # both 4pi factors within the forward transformation @@ -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/pyoz_input.py pyoz-0.2.2/pyoz_input.py --- pyoz-0.2/pyoz_input.py 2009-05-14 15:26:37.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.5 +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 @@ -271,6 +292,7 @@ nr_mix_param = ctrl_def.nr_mix_param nr_convergence_factor = ctrl_def.nr_convergence_factor nr_max_iter = ctrl_def.nr_max_iter + nr_noconv_incr = ctrl_def.nr_noconv_incr # parse the inputs and (possibly) overwrite the input data try: @@ -303,61 +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)) @@ -371,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 @@ -390,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 @@ -403,15 +435,20 @@ ctrl['nr_mix_param'] = nr_mix_param ctrl['nr_convergence_factor'] = nr_convergence_factor ctrl['nr_max_iter'] = nr_max_iter + ctrl['nr_noconv_incr'] = nr_noconv_incr if (ctrl['do_nr']): print("\n\tNewton-Raphson iteration method will be used") print("\tNR max iterations\t\t%lu" % ctrl['nr_max_iter']) - print("\tNR initial mix parameter\t%f" % ctrl['nr_mix_param']) - print("\tNR relative convergence factor\t%f" % ctrl['nr_convergence_factor']) - print("\tfallback (Picard) mix parameter\t%f" % ctrl['mix_param']) + print("\tNR initial mix parameter\t%e" % ctrl['nr_mix_param']) + print("\tNR relative convergence factor\t%e" % ctrl['nr_convergence_factor']) + if (not nr_noconv_incr): + print("\tPicard will be used if NR fails") + print("\tfallback Picard mix parameter\t%e" % ctrl['mix_param']) + else: + print("\tnon-converged increment will be used if NR fails") else: print("\n\tPicard iteration method will be used") - print("\tPicard mix parameter\t\t%f" % ctrl['mix_param']) + print("\tPicard mix parameter\t\t%e" % ctrl['mix_param']) # use graphics? ctrl['do_graphics'] = do_graphics @@ -479,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)): @@ -497,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 @@ -561,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, ... @@ -827,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: @@ -838,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 @@ -896,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 @@ -1009,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: @@ -1020,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 @@ -1244,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: @@ -1255,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 @@ -1454,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: @@ -1465,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 @@ -1574,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 @@ -1609,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 @@ -1665,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 @@ -1687,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 @@ -1702,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: @@ -1716,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' @@ -1782,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 @@ -1824,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)): @@ -1838,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 @@ -1858,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)) @@ -1910,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/pyoz.py pyoz-0.2.2/pyoz.py --- pyoz-0.2/pyoz.py 2009-05-14 11:07:20.000000000 +0200 +++ pyoz-0.2.2/pyoz.py 2009-07-24 09:39:20.000000000 +0200 @@ -9,18 +9,9 @@ """ 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.4 +revision of this file: 0.1.6 -supported closure relations: - hnc - HyperNetted Chain - py - Percus-Yevick - -suported interaction potentials: - hard spheres -planned support - coulomb potential - van der Waals potential - external PMFs +check the file _features.txt and website http://pyoz.vrbka.net for further information in the whole program (except where stated otherwise), the numpy/scipy arrays are used. as a result, all mathematical operations are done element-wise. this applies also for multiplication @@ -131,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'])) @@ -173,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'])) @@ -305,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") @@ -329,7 +320,7 @@ # test for convergence and write the gamma if everything is OK norm_dsqn = convergence_dsqn(ctrl, syst, G_o_ij, G_n_ij) time_end = time() - print("%f sec - DSQN %e -" % ((time_end - time_beg), norm_dsqn)), + print("%f sec - DSQN %.3e -" % ((time_end - time_beg), norm_dsqn)), if (norm_dsqn > ctrl['max_dsqn'] or (not isfinite(norm_dsqn))): print("\nDSQN too large, calculation is probably diverging") print("check inputs and outputs and/or increase the value of max_dsqn (%e at the moment)" % ctrl['max_dsqn']) @@ -389,19 +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 = E_ij + 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() @@ -411,7 +415,6 @@ # we make one half-iteration (number 0) and then carry on with full cycles until convergence while (not nr_converged and nr_niter <= ctrl['nr_max_iter']): - #print("main\t%4u " % niter), print(" nr/cg\t :%-4u" % (nr_niter)), # show "progress bar" when output is redirected - this time with plus sign @@ -436,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']): @@ -457,15 +468,15 @@ # do the calculation Rcur = B - AX - #print("remainder %e" % abs(dotproduct(ctrl, syst, r, Rcur, Rcur))), # check for convergence here - if converged, abandon the cycle! # we check how far is Rcur from zero (Z is zero array with the same dimensions as Rcur) nr_norm_dsqn = convergence_dsqn(ctrl, syst, Rcur, Z) + #nr_norm_dsqn = convergence_dsqn(ctrl, syst, B, AX) nr_time_end = time() # convergence is tested relatively to the DSQN of the 'outer' cycle - print("%f sec - rel. DSQN %e -" % ((nr_time_end - nr_time_beg), nr_norm_dsqn/norm_dsqn)), + print("%f sec - rel. DSQN %.3e -" % ((nr_time_end - nr_time_beg), nr_norm_dsqn/norm_dsqn)), if (nr_norm_dsqn > ctrl['max_dsqn'] or (not isfinite(nr_norm_dsqn))): print("\n\tDSQN too large, calculation is probably diverging") break @@ -476,7 +487,7 @@ else: # the conjugate gradients algorithm needed - NR has not converged print("not converged") - + nr_time_beg = time() nr_niter += 1 @@ -484,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] @@ -499,16 +515,16 @@ 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 Lcur = abs(dotproduct(ctrl, syst, r, Rcur, Rcur))/abs(dotproduct(ctrl, syst, r, AtR, AtR)) - + # calculate Wnew = W(n+1) (except for first iteration) if (nr_niter != 1): # do the full calculation, in the first iteration the value is pre-set to 1.0 - Wpartial = 1.0 - Lcur * abs(dotproduct(ctrl, syst, r, Rcur, Rcur))/Lold * Wcur * abs(dotproduct(ctrl, syst, r, Rold, Rold)) + Wpartial = 1.0 - Lcur * abs(dotproduct(ctrl, syst, r, Rcur, Rcur))/(Lold * Wcur * abs(dotproduct(ctrl, syst, r, Rold, Rold))) Wnew = 1.0 / Wpartial # end calculation of Wnew @@ -517,21 +533,26 @@ # test if (AX,R) = (X,AtR) - is the adjoint OK? is giving me zero (always)! #if (nr_niter == 1): - # test1 = dotproduct(ctrl, syst, r, AX, Rcur) - # test2 = dotproduct(ctrl, syst, r, Xcur, AtR) - # diff = abs(test2-test1) - # print("\tadjoint test - difference %e" % abs(test2-test1)) + #test1 = dotproduct(ctrl, syst, r, AX, Rcur) + #test2 = dotproduct(ctrl, syst, r, Xcur, AtR) + #diff = abs(test2-test1) + #print("\tadjoint test - difference %e" % abs(test2-test1)) # end if (nr_norm_dsqn <= nr_convergence_crit) - handling of the else-branch (not converged) # end while (not nr-converged and nr_niter < ctrl['nr_max_iter']) # in case convergence was not reached, do Picard if (not nr_converged): - print("\tcouldn't converge NR/CG cycle, using Picard iteration instead") - G_r_ij = (1.0 - ctrl['mix_param']) * G_o_ij + ctrl['mix_param'] * G_n_ij + print("\tcouldn't converge NR/CG cycle,"), + if (not ctrl['nr_noconv_incr']): + print("using Picard iteration instead") + G_r_ij = (1.0 - ctrl['mix_param']) * G_o_ij + ctrl['mix_param'] * G_n_ij + else: + print("using non-converged increment") + G_r_ij = G_o_ij + Xnew else: - G_r_ij = G_o_ij + Xcur - #G_r_ij = G_o_ij + Xnew + #G_r_ij = G_o_ij + Xcur + G_r_ij = G_o_ij + Xnew # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # end else (Newton-Raphson method) @@ -566,7 +587,7 @@ # end while (not converged) - print("\niteration process completed after iteration %u" % niter) + print("\niteration process completed in iteration %u" % niter) if (converged): print("\tcalculation converged") else: @@ -574,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) @@ -605,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 @@ -664,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/tests and pyoz-0.2.2/tests