diff -u pyoz-0.2/_notes.txt pyoz-0.2.1/_notes.txt --- pyoz-0.2/_notes.txt 2009-05-14 11:08:33.000000000 +0200 +++ pyoz-0.2.1/_notes.txt 2009-05-20 09:42:30.000000000 +0200 @@ -1,9 +1,9 @@ release notes for pyOZ ====================== -released: 14. May 2009 -version: 0.2 -second release +released: 20. May 2009 +version: 0.2.1 +bugfix to the version 0.2 check the website http://pyoz.vrbka.net @@ -13,27 +13,22 @@ - 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) + - fixed bug in NR/CG handling causing very slow convergence + - new handling of non-converged NR/CG cycle (ctrl:nr_noconv_picard) revisions/descriptions of individual source files -pyoz.py ........................................... 0.1.4 +pyoz.py ........................................... 0.1.5 main program pyoz_closure.py ................................... 0.1.1 definition of closures -pyoz_const.py ..................................... 0.1.3 +pyoz_const.py ..................................... 0.1.4 definition of constants -pyoz_dft.py ....................................... 0.1.3 +pyoz_dft.py ....................................... 0.1.4 definition of fourier transforms -pyoz_input.py ..................................... 0.1.5 +pyoz_input.py ..................................... 0.1.6 handling of input pyoz_misc.py ...................................... 0.1.3 miscellanous routines @@ -45,3 +40,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.1/pyoz_const.py --- pyoz-0.2/pyoz_const.py 2009-05-08 16:58:32.000000000 +0200 +++ pyoz-0.2.1/pyoz_const.py 2009-05-20 09:09:50.000000000 +0200 @@ -8,14 +8,14 @@ """ Module defining parameters and constants -revision of this file: 0.1.3 +revision of this file: 0.1.4 """ from math import sqrt, pi from scipy import inf -pyoz_version = "0.2" +pyoz_version = "0.2.1" 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 # ************************************************************************************************ diff -u pyoz-0.2/pyoz_dft.py pyoz-0.2.1/pyoz_dft.py --- pyoz-0.2/pyoz_dft.py 2009-05-08 23:28:29.000000000 +0200 +++ pyoz-0.2.1/pyoz_dft.py 2009-05-20 08:55:31.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.4 for the 3D-Fourier transformation, there is a normalization factor of 1/(2*pi)^3 involved then, the forward FT can be normalized @@ -124,14 +124,14 @@ # 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 diff -u pyoz-0.2/pyoz_input.py pyoz-0.2.1/pyoz_input.py --- pyoz-0.2/pyoz_input.py 2009-05-14 15:26:37.000000000 +0200 +++ pyoz-0.2.1/pyoz_input.py 2009-05-20 09:10:47.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.6 """ # import external modules @@ -271,6 +271,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: @@ -343,6 +344,8 @@ 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") + elif (match(r'nr_noconv_incr$', keyword)): + nr_noconv_incr = True else: raise "UnkKeyword", (i) # end for i in range(len(lines)): @@ -403,15 +406,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 diff -u pyoz-0.2/pyoz.py pyoz-0.2.1/pyoz.py --- pyoz-0.2/pyoz.py 2009-05-14 11:07:20.000000000 +0200 +++ pyoz-0.2.1/pyoz.py 2009-05-20 11:36:40.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.5 -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 @@ -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']) @@ -393,7 +384,13 @@ # 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 + # + # 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 @@ -411,7 +408,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 @@ -457,15 +453,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 +472,7 @@ else: # the conjugate gradients algorithm needed - NR has not converged print("not converged") - + nr_time_beg = time() nr_niter += 1 @@ -504,11 +500,11 @@ # 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 +513,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 +567,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: