1 2from collections import namedtuple, defaultdict 3from sympy import * 4 5# See the GaussianOrbitals notebook in the qmc_algorithms repo for more explanation, 6# especially about the normalization. 7 8def single_gaussian(): 9 r = Symbol('r',nonnegative=True) 10 alpha = Symbol('alpha') 11 12 alpha_val = 3.0 13 r_val = 1.2 14 print 'alpha = ',alpha_val 15 print 'r = ',r_val 16 17 phi = exp(-alpha*r**2) 18 print phi 19 print 'f == ',phi.subs(alpha, alpha_val).subs(r,r_val) 20 21 d_phi = diff(phi, r) 22 print 'symbolic df = ',d_phi 23 print 'df == ',d_phi.subs(alpha, alpha_val).subs(r,r_val) 24 25 dd_phi = diff(phi, r, 2) 26 print 'symbolic d2f = ',dd_phi 27 print 'd2f == ',dd_phi.subs(alpha, alpha_val).subs(r,r_val) 28 29 d3d_phi = diff(phi, r, 3) 30 print 'symbolic d3f = ',d3d_phi 31 print 'd3f == ',d3d_phi.subs(alpha, alpha_val).subs(r,r_val) 32 33 34CG_basis = namedtuple('CG_basis',['orbtype','nbasis','zeta','contraction_coeff']) 35 36# Read basis sets in Gamess format 37def read_gms_basis(fname): 38 element = '' 39 basis_sets = dict() 40 with open(fname,'r') as f: 41 lines = f.readlines() 42 idx = 0 43 while idx < len(lines): 44 line = lines[idx].strip() 45 if line.startswith('!') or len(line) == 0: 46 idx += 1 47 continue 48 if line.startswith('$DATA'): 49 idx += 1 50 element = lines[idx].strip() 51 basis_sets[element] = defaultdict(list) 52 while idx < len(lines): 53 idx += 1 54 if lines[idx].startswith('$END'): 55 break 56 orbtype, s_nval = lines[idx].strip().split() 57 nval = int(s_nval) 58 coef_list = [] 59 zeta_list = [] 60 for i in range(nval): 61 idx += 1 62 s_orbidx,s_zeta,s_c = lines[idx].strip().split() 63 orbidx = int(s_orbidx) 64 zeta = float(s_zeta) 65 c = float(s_c) 66 zeta_list.append(zeta) 67 coef_list.append(c) 68 assert(len(coef_list) == nval) 69 assert(len(zeta_list) == nval) 70 cg = CG_basis(orbtype, nval, zeta_list, coef_list) 71 basis_sets[element][orbtype].append(cg) 72 idx +=1 73 return basis_sets 74 75def create_full_gto_norm_symbolic(): 76 # full normalization 77 n1 = sympify('(2*alpha/pi)**(3/4)') 78 n2 = sympify('(8*alpha)**(i+j+k)') 79 80 i = Symbol('i') 81 j = Symbol('j') 82 k = Symbol('k') 83 n3 = factorial(i)*factorial(j)*factorial(k) 84 d1 = factorial(2*i)*factorial(2*j)*factorial(2*k) 85 norm = n1*sqrt(n2*n3/d1) 86 return norm 87 88def create_radial_gto_norm_symbolic(): 89 # just the radial part of the normalization 90 n1 = sympify('2*alpha**(3/4)*2**(3/4)* (1/pi)**(1/4)') 91 #n2 = sympify('(8*alpha)**(i+j+k)') 92 n2 = sympify('(8*alpha)**L') 93 #i = Symbol('i') 94 #j = Symbol('j') 95 #k = Symbol('k') 96 L = Symbol('L') 97 #n3 = factorial(i)*factorial(j)*factorial(k) 98 #d1 = factorial(2*i)*factorial(2*j)*factorial(2*k) 99 n3 = factorial(L) 100 d1 = factorial(2*L) 101 #d2 = 2*(i+j+k) + 1 102 103 d2 = 2*L + 1 104 norm = n1*sqrt(n2*n3/d1/d2) 105 return norm 106 107def create_angular_gto_norm_symbolic(): 108 return create_full_gto_norm_symbolic() / create_radial_gto_norm_symbolic() 109 110def create_gto_symbolic(): 111 pre = sympify('x**i * y**j * z**k * exp(-alpha *r**2)') 112 return pre 113 114 115def create_sym(ijk=[0,0,0]): 116 gto_sym = create_gto_symbolic() 117 gto = gto_sym.subs({Symbol('i'):ijk[0], Symbol('j'):ijk[1],Symbol('k'):ijk[2]}) 118 norm = create_radial_gto_norm_symbolic() 119 l_val = sum(ijk) 120 norm_s = norm.subs({Symbol('i'):ijk[0], Symbol('j'):ijk[1],Symbol('k'):ijk[2], Symbol('L'):l_val}) 121 print 'norm_s',norm_s 122 norm = lambdify(Symbol('alpha'), norm_s) 123 124 i = Symbol('i',integer=True) 125 c = IndexedBase('c') 126 alpha2 = IndexedBase('alpha') 127 norm2 = IndexedBase('N') 128 N_basis = Symbol('N_b', integer=True) 129 cg_sym = Sum(norm2[i]*c[i]*gto.subs(Symbol('alpha'),alpha2[i]),(i,1,N_basis)) 130 return cg_sym,norm,norm2,c,alpha2,N_basis 131 132def eval_sym(cg_sym, norm, norm2, N_basis, c, alpha2, h_basis): 133 cg_unroll = cg_sym.subs(N_basis, h_basis.nbasis).doit() 134 for i in range(h_basis.nbasis): 135 cc = h_basis.contraction_coeff[i] 136 cz = h_basis.zeta[i] 137 cg_unroll = cg_unroll.subs(c[i+1],cc).subs(alpha2[i+1],cz).subs(norm2[i+1],norm(cz)) 138 print cc,cz,norm(cz),'normL',norm(1.0) 139 return cg_unroll 140 141def compute_from_sym(h_basis, ijk=[0,0,0]): 142 143 cg_sym, norm, norm2, c, alpha2, N_basis = create_sym(ijk) 144 print 'norm',norm 145 cg = eval_sym(cg_sym, norm, norm2, N_basis, c, alpha2, h_basis) 146 r = Symbol('r') 147 x = Symbol('x') 148 # setting x to 1.0 is important to compute just the radial part 149 slist = {r:1.3, x:1.0} 150 151 print cg 152 print 'f = ',cg.subs(slist) 153 154 d_cg = diff(cg, r); 155 print 'df = ',d_cg.subs(slist) 156 157 dd_cg = diff(cg, r, 2); 158 print 'ddf = ',dd_cg.subs(slist) 159 160 d3_cg = diff(cg, r, 3); 161 print 'd3f = ',d3_cg.subs(slist) 162 163# generated from read_order.py 164def get_ijk(): 165 ijk = [] 166 # S 167 ijk.append( (0,0,0,"S") ) 168 # P 169 ijk.append( (1,0,0,"X") ) 170 ijk.append( (0,1,0,"Y") ) 171 ijk.append( (0,0,1,"Z") ) 172 # D 173 ijk.append( (2,0,0,"XX") ) 174 ijk.append( (0,2,0,"YY") ) 175 ijk.append( (0,0,2,"ZZ") ) 176 ijk.append( (1,1,0,"XY") ) 177 ijk.append( (1,0,1,"XZ") ) 178 ijk.append( (0,1,1,"YZ") ) 179 # F 180 ijk.append( (3,0,0,"XXX") ) 181 ijk.append( (0,3,0,"YYY") ) 182 ijk.append( (0,0,3,"ZZZ") ) 183 ijk.append( (2,1,0,"XXY") ) 184 ijk.append( (2,0,1,"XXZ") ) 185 ijk.append( (1,2,0,"YYX") ) 186 ijk.append( (0,2,1,"YYZ") ) 187 ijk.append( (1,0,2,"ZZX") ) 188 ijk.append( (0,1,2,"ZZY") ) 189 ijk.append( (1,1,1,"XYZ") ) 190 # G 191 ijk.append( (4,0,0,"XXXX") ) 192 ijk.append( (0,4,0,"YYYY") ) 193 ijk.append( (0,0,4,"ZZZZ") ) 194 ijk.append( (3,1,0,"XXXY") ) 195 ijk.append( (3,0,1,"XXXZ") ) 196 ijk.append( (1,3,0,"YYYX") ) 197 ijk.append( (0,3,1,"YYYZ") ) 198 ijk.append( (1,0,3,"ZZZX") ) 199 ijk.append( (0,1,3,"ZZZY") ) 200 ijk.append( (2,2,0,"XXYY") ) 201 ijk.append( (2,0,2,"XXZZ") ) 202 ijk.append( (0,2,2,"YYZZ") ) 203 ijk.append( (2,1,1,"XXYZ") ) 204 ijk.append( (1,2,1,"YYXZ") ) 205 ijk.append( (1,1,2,"ZZXY") ) 206 # H 207 ijk.append( (5,0,0,"XXXXX") ) 208 ijk.append( (0,5,0,"YYYYY") ) 209 ijk.append( (0,0,5,"ZZZZZ") ) 210 ijk.append( (4,1,0,"XXXXY") ) 211 ijk.append( (4,0,1,"XXXXZ") ) 212 ijk.append( (1,4,0,"YYYYX") ) 213 ijk.append( (0,4,1,"YYYYZ") ) 214 ijk.append( (1,0,4,"ZZZZX") ) 215 ijk.append( (0,1,4,"ZZZZY") ) 216 ijk.append( (3,2,0,"XXXYY") ) 217 ijk.append( (3,0,2,"XXXZZ") ) 218 ijk.append( (2,3,0,"YYYXX") ) 219 ijk.append( (0,3,2,"YYYZZ") ) 220 ijk.append( (2,0,3,"ZZZXX") ) 221 ijk.append( (0,2,3,"ZZZYY") ) 222 ijk.append( (3,1,1,"XXXYZ") ) 223 ijk.append( (1,3,1,"YYYXZ") ) 224 ijk.append( (1,1,3,"ZZZXY") ) 225 ijk.append( (2,2,1,"XXYYZ") ) 226 ijk.append( (2,1,2,"XXZZY") ) 227 ijk.append( (1,2,2,"YYZZX") ) 228 # I 229 ijk.append( (6,0,0,"X6") ) 230 ijk.append( (0,6,0,"Y6") ) 231 ijk.append( (0,0,6,"Z6") ) 232 ijk.append( (5,1,0,"X5Y") ) 233 ijk.append( (5,0,1,"X5Z") ) 234 ijk.append( (1,5,0,"Y5X") ) 235 ijk.append( (0,5,1,"Y5Z") ) 236 ijk.append( (1,0,5,"Z5X") ) 237 ijk.append( (0,1,5,"Z5Y") ) 238 ijk.append( (4,2,0,"X4Y2") ) 239 ijk.append( (4,0,2,"X4Z2") ) 240 ijk.append( (2,4,0,"Y4X2") ) 241 ijk.append( (0,4,2,"Y4Z2") ) 242 ijk.append( (2,0,4,"Z4X2") ) 243 ijk.append( (0,2,4,"Z4Y2") ) 244 ijk.append( (4,1,1,"X4YZ") ) 245 ijk.append( (1,4,1,"Y4XZ") ) 246 ijk.append( (1,1,4,"Z4XY") ) 247 ijk.append( (3,3,0,"X3Y3") ) 248 ijk.append( (3,0,3,"X3Z3") ) 249 ijk.append( (0,3,3,"Y3Z3") ) 250 ijk.append( (3,2,1,"X3Y2Z") ) 251 ijk.append( (3,1,2,"X3Z2Y") ) 252 ijk.append( (2,3,1,"Y3X2Z") ) 253 ijk.append( (1,3,2,"Y3Z2X") ) 254 ijk.append( (2,1,3,"Z3X2Y") ) 255 ijk.append( (1,2,3,"Z3Y2X") ) 256 ijk.append( (2,2,2,"X2Y2Z2") ) 257 258 return ijk 259 260 261# To create test cases for each routine 262# evaluate print_value = True, others False 263# evaluateAll print_value,print_grad,print_lap = True, others False 264# evaluateWithHessian print_value,print_grad,print_hess = True, others False 265# evaluateWithThirdDeriv print_value,print_grad,print_hess,print_ggg = True 266# evaluateThirdDerivOnly print_ggg = True, others False 267 268def compute_radial_values(p, print_value=False, print_grad=False, print_lap=False, 269 print_hess=False, print_ggg=False, using_soa=False): 270 271 out = '' 272 gto_s = create_gto_symbolic() 273 # just the 'angular' part 274 gto_s = gto_s.subs(Symbol('alpha'),0) 275 #print gto_s 276 277 norm_s = create_angular_gto_norm_symbolic() 278 #print 'norm',norm_s 279 280 281 ijk_s = get_ijk() 282 for idx, (i,j,k,s) in enumerate(ijk_s): 283 expr = gto_s * norm_s 284 l_val = i + j + k 285 slist = {Symbol('i'):i, Symbol('j'):j, Symbol('k'):k, Symbol('L'):l_val} 286 rlist = {Symbol('x'):p[0], Symbol('y'):p[1], Symbol('z'):p[2]} 287 val = expr.subs(slist).subs(rlist).evalf() 288 if print_value: 289 if using_soa: 290 out += " REQUIRE(XYZ[%d] == Approx(%.12g));\n"%(idx, val) 291 else: 292 out += " REQUIRE(ct.getYlm(%d) == Approx(%.12g));\n"%(idx, val) 293 dx = diff(expr, Symbol('x')) 294 dy = diff(expr, Symbol('y')) 295 dz = diff(expr, Symbol('z')) 296 dx_val = dx.subs(slist).subs(rlist).evalf() 297 dy_val = dy.subs(slist).subs(rlist).evalf() 298 dz_val = dz.subs(slist).subs(rlist).evalf() 299 if print_grad: 300 if using_soa: 301 out += " REQUIRE(gr0[%d] == Approx(%.12g));\n"%(idx, dx_val) 302 out += " REQUIRE(gr1[%d] == Approx(%.12g));\n"%(idx, dy_val) 303 out += " REQUIRE(gr2[%d] == Approx(%.12g));\n"%(idx, dz_val) 304 else: 305 out += " REQUIRE(ct.getGradYlm(%d)[0] == Approx(%.12g));\n"%(idx, dx_val) 306 out += " REQUIRE(ct.getGradYlm(%d)[1] == Approx(%.12g));\n"%(idx, dy_val) 307 out += " REQUIRE(ct.getGradYlm(%d)[2] == Approx(%.12g));\n"%(idx, dz_val) 308 309 lap = diff(expr, Symbol('x'), 2) + diff(expr, Symbol('y'), 2) + diff(expr, Symbol('z'), 2) 310 lap_val = lap.subs(slist).subs(rlist).evalf() 311 if print_lap: 312 if using_soa: 313 out += " REQUIRE(lap[%d] == Approx(%.12g));\n"%(idx, lap_val) 314 else: 315 out += " REQUIRE(ct.getLaplYlm(%d) == Approx(%.12g));\n"%(idx, lap_val) 316 317 318 if print_hess: 319 out += '\n' 320 axis_syms = [Symbol('x'), Symbol('y'), Symbol('z')] 321 for ii,si in enumerate(axis_syms): 322 for jj,sj in enumerate(axis_syms): 323 h_s = diff(diff(expr, si), sj) 324 hess_val = h_s.subs(slist).subs(rlist).evalf() 325 326 #print ii,jj,hess_val 327 #if hess_val != 0: 328 # print " REQUIRE(ct.getHessYlm(%d)(%d,%d) == Approx(%.12g));"%(idx, ii, jj, hess_val) 329 if using_soa: 330 if ii <= jj: 331 out += " REQUIRE(h%d%d[%d] == Approx(%.12g));\n"%(ii, jj, idx, hess_val) 332 else: 333 out += " REQUIRE(ct.getHessYlm(%d)(%d,%d) == Approx(%.12g));\n"%(idx, ii, jj, hess_val) 334 335 out += '\n' 336 337 if print_ggg: 338 out += '\n' 339 axis_syms = [Symbol('x'), Symbol('y'), Symbol('z')] 340 for ii,si in enumerate(axis_syms): 341 for jj,sj in enumerate(axis_syms): 342 for kk,sk in enumerate(axis_syms): 343 ggg_s = diff(diff(diff(expr, si), sj), sk) 344 ggg_val = ggg_s.subs(slist).subs(rlist).evalf() 345 346 out += " REQUIRE(ct.getGGGYlm(%d)[%d](%d,%d) == Approx(%.12g));\n"%(idx, ii, jj, kk, ggg_val) 347 348 out += '\n' 349 350 return out 351 352# A simple template replacement engine. 353# Template items to be replaced start on a line with '%'. 354def run_template(fname_in, fname_out, bodies): 355 out = '' 356 with open(fname_in, 'r') as f: 357 for line in f: 358 if line.startswith('%'): 359 key = line.strip()[1:] 360 if key in bodies: 361 line = bodies[key] 362 else: 363 print 'Error, template item not found, key:',key, ' line = ',line 364 out += line 365 366 with open(fname_out, 'w') as f: 367 f.write(out) 368 369 370def create_test_full_cartesian_tensor(): 371 p = [1.3,1.2,-0.5] 372 #compute_radial_values([1.3,1.2,-0.5]) 373 bodies = dict() 374 bodies['test_evaluate'] = compute_radial_values(p, print_value=True) 375 bodies['test_evaluate_all'] = compute_radial_values(p, print_value=True, 376 print_grad=True, 377 print_lap=True) 378 bodies['test_evaluate_with_hessian'] = compute_radial_values(p, print_value=True, 379 print_grad=True, 380 print_hess=True) 381 bodies['test_evaluate_with_third_deriv'] = compute_radial_values(p, print_value=True, 382 print_grad=True, 383 print_hess=True, 384 print_ggg=True) 385 bodies['test_evaluate_third_deriv_only'] = compute_radial_values(p, print_ggg=True) 386 387 fname_in = 'test_full_cartesian_tensor.cpp.in' 388 fname_out = 'test_full_cartesian_tensor.cpp' 389 390 run_template(fname_in, fname_out, bodies) 391 392def create_test_full_soa_cartesian_tensor(): 393 p = [1.3,1.2,-0.5] 394 #compute_radial_values([1.3,1.2,-0.5]) 395 bodies = dict() 396 bodies['test_evaluateV'] = compute_radial_values(p, print_value=True, using_soa=True) 397 bodies['test_evaluateVGL'] = compute_radial_values(p, print_value=True, 398 print_grad=True, 399 print_lap=True, 400 using_soa=True) 401 bodies['test_evaluateVGH'] = compute_radial_values(p, print_value=True, 402 print_grad=True, 403 print_hess=True, 404 using_soa=True) 405 406 fname_in = 'test_full_soa_cartesian_tensor.cpp.in' 407 fname_out = 'test_full_soa_cartesian_tensor.cpp' 408 409 run_template(fname_in, fname_out, bodies) 410 411 412 413 414if __name__ == '__main__': 415 # Generate data for unit tests below 416 417 # For test_gaussian_basis 418 419 # Data for 'Gaussian Combo' 420 #basis_sets = read_gms_basis('sto3g.txt') 421 #compute_from_sym(basis_sets['HYDROGEN']['S'][0]) 422 423 # Data for 'Gaussian Combo P' 424 #basis_sets_cc = read_gms_basis('cc-pVDZ.txt') 425 #compute_from_sym(basis_sets_cc['CARBON']['P'][0],[1,0,0]) 426 427 # Data for 'Gaussian Combo D' 428 #basis_sets_cc = read_gms_basis('cc-pVDZ.txt') 429 #compute_from_sym(basis_sets_cc['CARBON']['D'][0],[2,0,0]) 430 431 432 # Data for test_cartesian_tensor 433 #compute_radial_values([1.3,1.2,-0.5]) 434 435 # Create full test 436 create_test_full_cartesian_tensor() 437 create_test_full_soa_cartesian_tensor() 438 439