diff -u pyoz-0.2.2/_features.txt pyoz-0.2.3/_features.txt --- pyoz-0.2.2/_features.txt 2009-07-24 10:33:57.000000000 +0200 +++ pyoz-0.2.3/_features.txt 2009-07-29 10:00:10.000000000 +0200 @@ -1,14 +1,14 @@ list of features of pyOZ ======================== -released: 24. July 2009 -version: 0.2.2 +released: 29. July 2009 +version: 0.2.3 * general features o bulk calculations for simple atomic/ionic systems (no atomic details in molecules and molecular ions) on McMillan-Mayer level o simple iteration scheme (direct Picard iteration) - o improved iteration scheme (Newton-Raphson/conjugate gradients) for HNC closure + o (not available at the moment) improved iteration scheme (Newton-Raphson/conjugate gradients) for HNC closure o graphics output using Matplotlib * supported closure relations o Hypernetted Chain (HNC) @@ -24,5 +24,10 @@ o isothermal compressibility o excess chemical potential o Kirkwood-Buff factor - + * output + o pair correlation functions g(r) + o direct correlation functions c(r) + o partial structure factors S(k) + o pair potentials U(r) + o total interaction potential - U(r) + indirect correlation function diff -u pyoz-0.2.2/_notes.txt pyoz-0.2.3/_notes.txt --- pyoz-0.2.2/_notes.txt 2009-07-24 10:33:38.000000000 +0200 +++ pyoz-0.2.3/_notes.txt 2009-07-29 10:00:21.000000000 +0200 @@ -1,8 +1,9 @@ release notes for pyOZ ====================== -released: 24. July 2009 -version: 0.2.2 +released: 29. July 2009 +version: 0.2.3 +... check the website http://pyoz.vrbka.net @@ -12,24 +13,24 @@ - theory -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) +changes from the previous release (0.2.2) + - improved readability of the main code + - bugfixes + - output of partial structure factors S(k) revisions/descriptions of individual source files + - changed files denoted with star (*) -pyoz.py ........................................... 0.1.6 +pyoz.py ........................................... 0.1.7* main program pyoz_closure.py ................................... 0.1.1 definition of closures -pyoz_const.py ..................................... 0.1.5 +pyoz_const.py ..................................... 0.1.6* definition of constants -pyoz_dft.py ....................................... 0.1.5 +pyoz_dft.py ....................................... 0.1.4 definition of fourier transforms -pyoz_input.py ..................................... 0.1.7 +pyoz_input.py ..................................... 0.1.8* handling of input pyoz_misc.py ...................................... 0.1.3 miscellanous routines @@ -39,6 +40,6 @@ definition and handling of pair potentials pyoz_property.py .................................. 0.1.5 evaluation of thermodynamic properties -pyoz_solver.py .................................... 0.1.1 +pyoz_solver.py .................................... 0.1.2* optimized solvers for special cases diff -u pyoz-0.2.2/pyoz_const.py pyoz-0.2.3/pyoz_const.py --- pyoz-0.2.2/pyoz_const.py 2009-07-24 10:38:37.000000000 +0200 +++ pyoz-0.2.3/pyoz_const.py 2009-07-29 10:03:35.000000000 +0200 @@ -8,14 +8,14 @@ """ Module defining parameters and constants -revision of this file: 0.1.5 +revision of this file: 0.1.6 """ from math import sqrt, pi from scipy import inf -pyoz_version = "0.2.2" +pyoz_version = "0.2.3" class const: """ @@ -230,6 +230,7 @@ ur_suffix = '-ur.dat' urtot_suffix = '-urtot.dat' cr_suffix = '-cr.dat' + S_suffix = '-sk.dat' # Gamma related stuff # write output? @@ -265,6 +266,12 @@ # filename c_ij_name = name + cr_suffix + # S(k) related stuff + # write? + S_ij_write = False + # filename + S_ij_name = name + S_suffix + # ************************************************************************************************ if __name__ == "__main__": diff -u pyoz-0.2.2/pyoz_dft.py pyoz-0.2.3/pyoz_dft.py --- pyoz-0.2.2/pyoz_dft.py 2009-07-24 10:32:19.000000000 +0200 +++ pyoz-0.2.3/pyoz_dft.py 2009-07-24 09:23:29.000000000 +0200 @@ -11,7 +11,7 @@ """ module for performing various types of discrete transforms - revision of this file: 0.1.5 + revision of this file: 0.1.4 for the 3D-Fourier transformation, there is a normalization factor of 1/(2*pi)^3 involved then, the forward FT can be normalized diff -u pyoz-0.2.2/pyoz_input.py pyoz-0.2.3/pyoz_input.py --- pyoz-0.2.2/pyoz_input.py 2009-07-24 09:17:45.000000000 +0200 +++ pyoz-0.2.3/pyoz_input.py 2009-08-04 16:51:51.000000000 +0200 @@ -8,7 +8,7 @@ """ Module with functions for input file parsing -revision of this file: 0.1.7 +revision of this file: 0.1.8 """ # import external modules @@ -1819,12 +1819,16 @@ # direct correlatio function c_ij_write = outp_def.c_ij_write c_ij_name = outp_def.c_ij_name + # structure factors + S_ij_write = outp_def.S_ij_write + S_ij_name = outp_def.S_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 + S_ij_name = cmdline['name'] + outp_def.S_suffix U_ij_name = cmdline['name'] + outp_def.ur_suffix Utot_ij_name = cmdline['name'] + outp_def.urtot_suffix @@ -1901,6 +1905,15 @@ c_ij_name = (orig_l.split())[1] # ************************************************************************** + # S(k) + + elif (match(r'sk_write$', keyword)): + S_ij_write = True + elif (match(r'sk_name$', keyword)): + S_ij_write = True + S_ij_name = (orig_l.split())[1] + + # ************************************************************************** else: raise unkKeyword((i,l)) @@ -1963,6 +1976,14 @@ #else: # print("\tc(r) will not be written!") + # S(k) related stuff + outp['S_ij_write'] = S_ij_write + outp['S_ij_name'] = S_ij_name + if (outp['S_ij_write']): + print("\tS(k) will be saved to file %s" % (outp['S_ij_name'])) + #else: + # print("\tS(k) 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.2/pyoz.py pyoz-0.2.3/pyoz.py --- pyoz-0.2.2/pyoz.py 2009-07-24 09:39:20.000000000 +0200 +++ pyoz-0.2.3/pyoz.py 2009-07-28 15:24:12.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.6 +revision of this file: 0.1.7 check the file _features.txt and website http://pyoz.vrbka.net for further information @@ -19,11 +19,7 @@ """ # python modules -# probably will have to be made faster by importing -# individual methods/attributes only import sys -# mathematical constants -from math import pi # timing functions from time import time #import scipy @@ -128,11 +124,10 @@ # allocate arrays with direct, total and pair corr. functions, Gamma function # some arrays will emerge from arithmetic operations, # showing them here makes the code clearer - # real space: _r_ - # Fourier space _f_ - # direct correlation function without (c) and with (C) density factor applied + # real space: _r_, Fourier space _f_ # direct correlation function with Ng-correction (cs) applied - # w/o density correction in real space; w/ density correction in Fourier space + # direct correlation function without (c) and with (C) density factor applied + # w/o density correction in real space; w/ density correction in Fourier space #c_r_ij C_f_ij = zeros((syst['ncomponents'], syst['ncomponents'], ctrl['npoints'])) Cs_f_ij = zeros((syst['ncomponents'], syst['ncomponents'], ctrl['npoints'])) @@ -267,7 +262,6 @@ # FT c_r_ij to c_f_ij # we are using Fourier-sine transform; there are some steps involved in between FBT (Bessel) and FST - # which lead to use of constants and extra functions # this will not be discussed here, check the documentation and pyoz_dft.py for further information # the whole program is using FTs normalized with the density prefactors # sqrt(rho_i * rho_j) in order to have dimensionless functions in k-space @@ -297,9 +291,10 @@ # 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']) - S = solver_function((E_ij - C_f_ij), E_ij, ctrl['npoints']) + H_f_ij = solver_function((E_ij - dft.ft_convolution_factor * C_f_ij), C_f_ij, ctrl['npoints']) + from math import pi + S = E_ij + H_f_ij + #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 @@ -307,7 +302,6 @@ G_f_ij = S - E_ij - Cs_f_ij # FT G_f_ij to G_r_ij - #print("\tiFT Gamma_f -> Gamma_r") for i in range(syst['ncomponents']): for j in range(syst['ncomponents']): # perform the inverse Fourier transform of the Gamma function @@ -382,10 +376,7 @@ # 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 - # it's not possible to have 3D-matrix => we do it with arrays here and convert to matrices where needed and appropriate - # - # 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 do it with arrays here and convert to matrices where needed and appropriate # 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 @@ -414,9 +405,7 @@ # 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(" nr/cg\t :%-4u" % (nr_niter)), - # show "progress bar" when output is redirected - this time with plus sign # normal iterations are done with "." as an indicator if (cmdline['output'] != None): @@ -529,15 +518,7 @@ # end calculation of Wnew # calculate X(n+1) = X(n-1) + W(n+1)(alpha(n)ATR(n) + X(n) + X(n-1) - Xnew = Xold + Wnew * (Lcur * AtR + Xcur - Xold) - - # 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)) - + Xnew = Xold + Wnew * (Lcur * AtR + Xcur - Xold) # 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']) @@ -556,11 +537,6 @@ # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # end else (Newton-Raphson method) - - # symmetrize the matrix - takes too long and shouldn't actually be needed! - #for dr in range(ctrl['npoints']): - # G_r_ij[:,:,dr] = (G_r_ij[:,:,dr] + G_r_ij[:,:,dr].transpose()) / 2 - # end else (calculation not converged) # output data @@ -604,11 +580,7 @@ # end if(do_graphics): print("\nsaving outputs") - # some error checking should be added here! for both g, G - # allow more functions to be output - # probably g(r), Gamma, PMF (U_ij, U_ij_sr)... - # new section to input file %output %end - + # some error checking should be added here! for both g, G try: # save g_r_ij to file if (outp['g_ij_write']): @@ -649,6 +621,27 @@ fw.write("\n") fw.close() + # save S to file + if (outp['S_ij_write']): + print("\tpartial structure factors\t(%s)" % outp['S_ij_name']) + fw = open(outp['S_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" % k[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" % S[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 diff -u pyoz-0.2.2/pyoz_solver.py pyoz-0.2.3/pyoz_solver.py --- pyoz-0.2.2/pyoz_solver.py 2009-04-28 15:30:45.000000000 +0200 +++ pyoz-0.2.3/pyoz_solver.py 2009-07-28 16:00:21.000000000 +0200 @@ -8,7 +8,7 @@ """ Module defining solver-related functions -revision of this file: 0.1.1 +revision of this file: 0.1.2 currently contained are solvers for 1,2,3,n components: solver_1 @@ -54,7 +54,8 @@ # do the calculation and return immediately # we do not need the np here - return((b/a)/dens_factor) + #return((b/a)/dens_factor) + return(b/a) # two components def solver_2(a, b, np):