From 95cf4023dd5e9662be09aa020882bd2ab9f434d4 Mon Sep 17 00:00:00 2001 From: Bert de Jong Date: Thu, 29 Nov 2001 22:49:23 +0000 Subject: [PATCH] Python program that can be used to generate Coulomb fitting basis sets for atoms, and for XF and XH molecules. --- contrib/python/coulombfitting.XF.mol.nw | 188 ++++++++++++++++++++++++ contrib/python/coulombfitting.XH.mol.nw | 188 ++++++++++++++++++++++++ contrib/python/coulombfitting.atom.nw | 184 +++++++++++++++++++++++ 3 files changed, 560 insertions(+) create mode 100644 contrib/python/coulombfitting.XF.mol.nw create mode 100644 contrib/python/coulombfitting.XH.mol.nw create mode 100644 contrib/python/coulombfitting.atom.nw diff --git a/contrib/python/coulombfitting.XF.mol.nw b/contrib/python/coulombfitting.XF.mol.nw new file mode 100644 index 0000000000..92549a5e0a --- /dev/null +++ b/contrib/python/coulombfitting.XF.mol.nw @@ -0,0 +1,188 @@ +start coulombfitting + +# This python program generates Coulomb fitting basis sets +# for a XF molecule, with X being the atom of choice, +# using a uhf scf or dft wave function, both +# with Hartree-Fock exchange +# +# Optimize an even tempered like set of primitive functions +# for an molecule and return list of exponents and contraction +# coefficients. The molecule is atom + H at 1.5 Angstrom + +# exp[i+1] = exp[i] * beta * (1 + gamma * (i^2/(n+1)^2) +# +# with i = 0,...,n-2 + +# Input block after the python call: +# run_type = "scf" or "dft" +# atom = which atom to optimize for +# charge = molecular charge +# nopen = nr of open shells (hence 2S) +# basis = name of basis in nwchem library +# ecp = name of ecp in nwchem library +# l = maximum angular momentum +# n = list with nr of exponents per l-value +# freeze = l-values that should not be optimized +# 0 = optimize +# 1 = not optimize +# exp0 = list with exp[0] a for expansion per l-value, see formula above +# beta = list with beta a for expansion per l-value, see formula above +# gamma = list with gamma a for expansion per l-value, see formula above + +set tolguess 1e-10 + +set lindep:n_dep 0 + +set int:acc_std 1e-25 + +print none + +python + from mathutil import * + + run_type = "scf" + atom = "Be" + charge = 0 + nopen = 0 + basis = "Ahlrichs_pVDZ" + ecp = "" + l = 2 + n = [9,5,2] + freeze = [0,1,1] + exp0 = [0.25,0.40,0.60] + beta = [1.44,1.44,1.44] + gamma = [1.00,1.00,1.00] + + def nr_lvalue_list(k): + nr = -1 + for i in range(0,k+1): + if (n[i]>0 and freeze[i]<1): + nr = nr + 1 + return (nr) + + def put_exp0_beta_gamma(): + z = zerovector(3*(nr_lvalue_list(l)+1)) + for i in range(0,l+1): + if (n[i]>0 and freeze[i]<1): + nr = nr_lvalue_list(i) + z[0+3*nr] = sqrt(exp0[i] - 0.000001) + z[1+3*nr] = sqrt(beta[i] - 0.000001) + z[2+3*nr] = sqrt(gamma[i] - 0.000001) + return (z) + + def get_exp0_beta_gamma(z,k): + if (freeze[k]>0): + exp0_n = exp0[k] + beta_n = beta[k] + gamma_n = gamma[k] + else: + nr = nr_lvalue_list(k) + exp0_n = z[0+3*nr]*z[0+3*nr]+0.000001 + beta_n = z[1+3*nr]*z[1+3*nr]+0.000001 + gamma_n = z[2+3*nr]*z[2+3*nr]+0.000001 + return (exp0_n, beta_n, gamma_n) + + def make_exp(z): + exponents = zeromatrix(max(n),l+1) + for j in range(0,l+1): + (exp0_n, beta_n, gamma_n) = get_exp0_beta_gamma(z,j) + exponents[0][j] = exp0_n + for i in range(0,n[j]-1): + ee = exponents[i][j] * beta_n + ee = ee * (1 + gamma_n * i * i / ((n[j] + 1) * (n[j] + 1)) ) + exponents[i+1][j] = ee + return exponents + + def type_setup(type): + if (type == "scf"): + string = "scf; ri-scf coulomb; uhf; vectors input start.movecs output final.movecs;" + string = string + ('nopen %i; maxiter 1; thresh 1e-8; end\n\n') % (nopen) + else: + string = "dft; xc hfexch 1.0; grid xfine; tolerances tight; odft;" + string = string + "vectors input start.movecs output final.movecs;" + string = string + ('mult %i; iterations 1; convergence energy 1e-8; end\n\n') % (nopen+1) + return string + + def ref_energy(): + if (run_type == "scf"): + return task_coulomb_ref('scf') + else: + return task_coulomb_ref('dft') + + def energy(z): + function_type = ["s", "p", "d", "f", "g", "h", "i"] + exponents = make_exp(z) + basis = 'basis "riscf basis" spherical noprint;' + for j in range(0,l+1): + for i in range(0,n[j]): + basis = basis + ('%s %s ; %f 1;' % (atom, function_type[j], exponents[i][j])) + basis = basis + "F library Ahlrichs_Coulomb_Fitting;" + if (run_type == "scf"): + basis = basis + "end\n\n" + type_setup('scf') + input_parse(basis) + return task_coulomb('scf') + else: + basis = basis + "end\n\n" + type_setup('dft') + input_parse(basis) + return task_coulomb('dft') + + def printexp(z): + function_type = ["s", "p", "d", "f", "g", "h", "i"] + exponents = make_exp(z) + print ' Exponents:' + for j in range(0,l+1): + print "%s - functions" % (function_type[j]) + for i in range(0,n[j]): + print " %14.8f" % exponents[i][j] + print "" + print " " + + # Setup list of variables to be optimized + + z = zerovector(3*(nr_lvalue_list(l)+1)) + z = put_exp0_beta_gamma() + + # Place geometry, basis, and (possibly) ecp into input / rtdb + + calc_setup = "geometry;" + calc_setup = calc_setup + ('%s 0 0 0;') % (atom) + calc_setup = calc_setup + "F 0 0 1.3; end\n\n" + calc_setup = calc_setup + ('basis spherical noprint; %s library %s;') % (atom, basis) + calc_setup = calc_setup + "F library Ahlrichs_VTZ; end\n\n" + if (ecp != ""): + calc_setup = calc_setup + ('ecp spherical noprint; %s library %s; end\n\n') % (atom,ecp) + calc_setup = calc_setup + ('charge %i \n\n') % (charge) + input_parse(calc_setup) + + # Calculate reference energy without Coulomb fitting + # + # First setup for scf or dft + + if (run_type == "scf"): + string = "scf; uhf; vectors output start.movecs;" + string = string + ('nopen %i; maxiter 100; thresh 1e-8; end\n\n') % (nopen) + else: + string = "dft; xc hfexch 1.0; grid xfine; tolerances tight; odft;" + string = string + ('mult %i; vectors output start.movecs;') % (nopen+1) + string = string + "iterations 100; convergence energy 1e-8; end\n\n" + input_parse(string) + reference = ref_energy() + + # Optimize fitting basis set + + (value,z) = quasinr(energy, z, 1e-5, 1e-12, printexp) + + # Print the final results + + print '\n\n Results\n' + print ' l n exp0 beta gamma ' + print ' --- --- --------- --------- --------- ' + for j in range(0,l+1): + (exp0_n, beta_n, gamma_n) = get_exp0_beta_gamma(z,j) + print '%4i %4i %12.6f %10.6f %10.6f' % (j, n[j], exp0_n, beta_n, gamma_n) + print '\n Final energy difference = %12.6f' % value + print '\n Energy without Coulomb fitting = %12.6f' % reference + printexp(z) +end + +task python diff --git a/contrib/python/coulombfitting.XH.mol.nw b/contrib/python/coulombfitting.XH.mol.nw new file mode 100644 index 0000000000..d831a91c07 --- /dev/null +++ b/contrib/python/coulombfitting.XH.mol.nw @@ -0,0 +1,188 @@ +start coulombfitting + +# This python program generates Coulomb fitting basis sets +# for a XH molecule, with X being the atom of choice, +# using a uhf scf or dft wave function, both +# with Hartree-Fock exchange +# +# Optimize an even tempered like set of primitive functions +# for an molecule and return list of exponents and contraction +# coefficients. The molecule is atom + H at 1.5 Angstrom + +# exp[i+1] = exp[i] * beta * (1 + gamma * (i^2/(n+1)^2) +# +# with i = 0,...,n-2 + +# Input block after the python call: +# run_type = "scf" or "dft" +# atom = which atom to optimize for +# charge = molecular charge +# nopen = nr of open shells (hence 2S) +# basis = name of basis in nwchem library +# ecp = name of ecp in nwchem library +# l = maximum angular momentum +# n = list with nr of exponents per l-value +# freeze = l-values that should not be optimized +# 0 = optimize +# 1 = not optimize +# exp0 = list with exp[0] a for expansion per l-value, see formula above +# beta = list with beta a for expansion per l-value, see formula above +# gamma = list with gamma a for expansion per l-value, see formula above + +set tolguess 1e-10 + +set lindep:n_dep 0 + +set int:acc_std 1e-25 + +print none + +python + from mathutil import * + + run_type = "scf" + atom = "Be" + charge = 0 + nopen = 0 + basis = "Ahlrichs_pVDZ" + ecp = "" + l = 2 + n = [9,5,2] + freeze = [0,1,1] + exp0 = [0.25,0.40,0.60] + beta = [1.44,1.44,1.44] + gamma = [1.00,1.00,1.00] + + def nr_lvalue_list(k): + nr = -1 + for i in range(0,k+1): + if (n[i]>0 and freeze[i]<1): + nr = nr + 1 + return (nr) + + def put_exp0_beta_gamma(): + z = zerovector(3*(nr_lvalue_list(l)+1)) + for i in range(0,l+1): + if (n[i]>0 and freeze[i]<1): + nr = nr_lvalue_list(i) + z[0+3*nr] = sqrt(exp0[i] - 0.000001) + z[1+3*nr] = sqrt(beta[i] - 0.000001) + z[2+3*nr] = sqrt(gamma[i] - 0.000001) + return (z) + + def get_exp0_beta_gamma(z,k): + if (freeze[k]>0): + exp0_n = exp0[k] + beta_n = beta[k] + gamma_n = gamma[k] + else: + nr = nr_lvalue_list(k) + exp0_n = z[0+3*nr]*z[0+3*nr]+0.000001 + beta_n = z[1+3*nr]*z[1+3*nr]+0.000001 + gamma_n = z[2+3*nr]*z[2+3*nr]+0.000001 + return (exp0_n, beta_n, gamma_n) + + def make_exp(z): + exponents = zeromatrix(max(n),l+1) + for j in range(0,l+1): + (exp0_n, beta_n, gamma_n) = get_exp0_beta_gamma(z,j) + exponents[0][j] = exp0_n + for i in range(0,n[j]-1): + ee = exponents[i][j] * beta_n + ee = ee * (1 + gamma_n * i * i / ((n[j] + 1) * (n[j] + 1)) ) + exponents[i+1][j] = ee + return exponents + + def type_setup(type): + if (type == "scf"): + string = "scf; ri-scf coulomb; uhf; vectors input start.movecs output final.movecs;" + string = string + ('nopen %i; maxiter 1; thresh 1e-8; end\n\n') % (nopen) + else: + string = "dft; xc hfexch 1.0; grid xfine; tolerances tight; odft;" + string = string + "vectors input start.movecs output final.movecs;" + string = string + ('mult %i; iterations 1; convergence energy 1e-8; end\n\n') % (nopen+1) + return string + + def ref_energy(): + if (run_type == "scf"): + return task_coulomb_ref('scf') + else: + return task_coulomb_ref('dft') + + def energy(z): + function_type = ["s", "p", "d", "f", "g", "h", "i"] + exponents = make_exp(z) + basis = 'basis "riscf basis" spherical noprint;' + for j in range(0,l+1): + for i in range(0,n[j]): + basis = basis + ('%s %s ; %f 1;' % (atom, function_type[j], exponents[i][j])) + basis = basis + "H library Ahlrichs_Coulomb_Fitting;" + if (run_type == "scf"): + basis = basis + "end\n\n" + type_setup('scf') + input_parse(basis) + return task_coulomb('scf') + else: + basis = basis + "end\n\n" + type_setup('dft') + input_parse(basis) + return task_coulomb('dft') + + def printexp(z): + function_type = ["s", "p", "d", "f", "g", "h", "i"] + exponents = make_exp(z) + print ' Exponents:' + for j in range(0,l+1): + print "%s - functions" % (function_type[j]) + for i in range(0,n[j]): + print " %14.8f" % exponents[i][j] + print "" + print " " + + # Setup list of variables to be optimized + + z = zerovector(3*(nr_lvalue_list(l)+1)) + z = put_exp0_beta_gamma() + + # Place geometry, basis, and (possibly) ecp into input / rtdb + + calc_setup = "geometry;" + calc_setup = calc_setup + ('%s 0 0 0;') % (atom) + calc_setup = calc_setup + "H 0 0 1.0; end\n\n" + calc_setup = calc_setup + ('basis spherical noprint; %s library %s;') % (atom, basis) + calc_setup = calc_setup + "H library Ahlrichs_VTZ; end\n\n" + if (ecp != ""): + calc_setup = calc_setup + ('ecp spherical noprint; %s library %s; end\n\n') % (atom,ecp) + calc_setup = calc_setup + ('charge %i \n\n') % (charge) + input_parse(calc_setup) + + # Calculate reference energy without Coulomb fitting + # + # First setup for scf or dft + + if (run_type == "scf"): + string = "scf; uhf; vectors output start.movecs;" + string = string + ('nopen %i; maxiter 100; thresh 1e-8; end\n\n') % (nopen) + else: + string = "dft; xc hfexch 1.0; grid xfine; tolerances tight; odft;" + string = string + ('mult %i; vectors output start.movecs;') % (nopen+1) + string = string + "iterations 100; convergence energy 1e-8; end\n\n" + input_parse(string) + reference = ref_energy() + + # Optimize fitting basis set + + (value,z) = quasinr(energy, z, 1e-5, 1e-12, printexp) + + # Print the final results + + print '\n\n Results\n' + print ' l n exp0 beta gamma ' + print ' --- --- --------- --------- --------- ' + for j in range(0,l+1): + (exp0_n, beta_n, gamma_n) = get_exp0_beta_gamma(z,j) + print '%4i %4i %12.6f %10.6f %10.6f' % (j, n[j], exp0_n, beta_n, gamma_n) + print '\n Final energy difference = %12.6f' % value + print '\n Energy without Coulomb fitting = %12.6f' % reference + printexp(z) +end + +task python diff --git a/contrib/python/coulombfitting.atom.nw b/contrib/python/coulombfitting.atom.nw new file mode 100644 index 0000000000..0cbe413f98 --- /dev/null +++ b/contrib/python/coulombfitting.atom.nw @@ -0,0 +1,184 @@ +start coulombfitting + +# This python program generates Coulomb fitting basis sets +# for an atom using a uhf scf or dft wave function, both +# with Hartree-Fock exchange +# +# Optimize an even tempered like set of primitive functions +# for an atom and return list of exponents and contraction +# coefficients + +# exp[i+1] = exp[i] * beta * (1 + gamma * (i^2/(n+1)^2) +# +# with i = 0,...,n-2 + +# Input block after the python call: +# run_type = "scf" or "dft" +# atom = which atom to optimize for +# charge = atomic charge +# nopen = nr of open shells (hence 2S) +# basis = name of basis in nwchem library +# ecp = name of ecp in nwchem library +# l = maximum angular momentum +# n = list with nr of exponents per l-value +# freeze = l-values that should not be optimized +# 0 = optimize +# 1 = not optimize +# exp0 = list with exp[0] a for expansion per l-value, see formula above +# beta = list with beta a for expansion per l-value, see formula above +# gamma = list with gamma a for expansion per l-value, see formula above + +set tolguess 1e-10 + +set lindep:n_dep 0 + +set int:acc_std 1e-25 + +print none + +python + from mathutil import * + + run_type = "scf" + atom = "Be" + charge = 0 + nopen = 0 + basis = "Ahlrichs_pVDZ" + ecp = "" + l = 2 + n = [9,5,2] + freeze = [0,1,1] + exp0 = [0.25,0.40,0.60] + beta = [1.44,1.44,1.44] + gamma = [1.00,1.00,1.00] + + def nr_lvalue_list(k): + nr = -1 + for i in range(0,k+1): + if (n[i]>0 and freeze[i]<1): + nr = nr + 1 + return (nr) + + def put_exp0_beta_gamma(): + z = zerovector(3*(nr_lvalue_list(l)+1)) + for i in range(0,l+1): + if (n[i]>0 and freeze[i]<1): + nr = nr_lvalue_list(i) + z[0+3*nr] = sqrt(exp0[i] - 0.000001) + z[1+3*nr] = sqrt(beta[i] - 0.000001) + z[2+3*nr] = sqrt(gamma[i] - 0.000001) + return (z) + + def get_exp0_beta_gamma(z,k): + if (freeze[k]>0): + exp0_n = exp0[k] + beta_n = beta[k] + gamma_n = gamma[k] + else: + nr = nr_lvalue_list(k) + exp0_n = z[0+3*nr]*z[0+3*nr]+0.000001 + beta_n = z[1+3*nr]*z[1+3*nr]+0.000001 + gamma_n = z[2+3*nr]*z[2+3*nr]+0.000001 + return (exp0_n, beta_n, gamma_n) + + def make_exp(z): + exponents = zeromatrix(max(n),l+1) + for j in range(0,l+1): + (exp0_n, beta_n, gamma_n) = get_exp0_beta_gamma(z,j) + exponents[0][j] = exp0_n + for i in range(0,n[j]-1): + ee = exponents[i][j] * beta_n + ee = ee * (1 + gamma_n * i * i / ((n[j] + 1) * (n[j] + 1)) ) + exponents[i+1][j] = ee + return exponents + + def type_setup(type): + if (type == "scf"): + string = "scf; ri-scf coulomb; uhf; vectors input start.movecs output final.movecs;" + string = string + ('nopen %i; maxiter 1; thresh 1e-8; end\n\n') % (nopen) + else: + string = "dft; xc hfexch 1.0; grid xfine; tolerances tight; odft;" + string = string + "vectors input start.movecs output final.movecs;" + string = string + ('mult %i; iterations 1; convergence energy 1e-8; end\n\n') % (nopen+1) + return string + + def ref_energy(): + if (run_type == "scf"): + return task_coulomb_ref('scf') + else: + return task_coulomb_ref('dft') + + def energy(z): + function_type = ["s", "p", "d", "f", "g", "h", "i"] + exponents = make_exp(z) + basis = 'basis "riscf basis" spherical noprint;' + for j in range(0,l+1): + for i in range(0,n[j]): + basis = basis + ('%s %s ; %f 1;' % (atom, function_type[j], exponents[i][j])) + if (run_type == "scf"): + basis = basis + "end\n\n" + type_setup('scf') + input_parse(basis) + return task_coulomb('scf') + else: + basis = basis + "end\n\n" + type_setup('dft') + input_parse(basis) + return task_coulomb('dft') + + def printexp(z): + function_type = ["s", "p", "d", "f", "g", "h", "i"] + exponents = make_exp(z) + print ' Exponents:' + for j in range(0,l+1): + print "%s - functions" % (function_type[j]) + for i in range(0,n[j]): + print " %14.8f" % exponents[i][j] + print "" + print " " + + # Setup list of variables to be optimized + + z = zerovector(3*(nr_lvalue_list(l)+1)) + z = put_exp0_beta_gamma() + + # Place geometry, basis, and (possibly) ecp into input / rtdb + + calc_setup = "geometry;" + calc_setup = calc_setup + ('%s 0 0 0; end\n\n') % (atom) + calc_setup = calc_setup + ('basis spherical noprint; %s library %s; end\n\n') % (atom, basis) + if (ecp != ""): + calc_setup = calc_setup + ('ecp spherical noprint; %s library %s; end\n\n') % (atom,ecp) + calc_setup = calc_setup + ('charge %i \n\n') % (charge) + input_parse(calc_setup) + + # Calculate reference energy without Coulomb fitting + # + # First setup for scf or dft + + if (run_type == "scf"): + string = "scf; uhf; vectors output start.movecs;" + string = string + ('nopen %i; maxiter 100; thresh 1e-8; end\n\n') % (nopen) + else: + string = "dft; xc hfexch 1.0; grid xfine; tolerances tight; odft;" + string = string + ('mult %i; vectors output start.movecs;') % (nopen+1) + string = string + "iterations 100; convergence energy 1e-8; end\n\n" + input_parse(string) + reference = ref_energy() + + # Optimize fitting basis set + + (value,z) = quasinr(energy, z, 1e-5, 1e-12, printexp) + + # Print the final results + + print '\n\n Results\n' + print ' l n exp0 beta gamma ' + print ' --- --- --------- --------- --------- ' + for j in range(0,l+1): + (exp0_n, beta_n, gamma_n) = get_exp0_beta_gamma(z,j) + print '%4i %4i %12.6f %10.6f %10.6f' % (j, n[j], exp0_n, beta_n, gamma_n) + print '\n Final energy difference = %12.6f' % value + print '\n Energy without Coulomb fitting = %12.6f' % reference + printexp(z) +end + +task python