1;;;; 2;;;; Copyright (C) 2013- Qiming Sun <osirpt.sun@gmail.com> 3;;;; Description: 4;;; dump to output file 5 6(load "utility.cl") 7(load "parser.cl") 8(load "derivator.cl") 9 10(defun gen-subscript (cells-streamer raw-script) 11 (labels ((gen-tex-iter (raw-script) 12 (cond ((null raw-script) raw-script) 13 ((vector? raw-script) 14 (list (gen-tex-iter (comp-x raw-script)) 15 (gen-tex-iter (comp-y raw-script)) 16 (gen-tex-iter (comp-z raw-script)))) 17 ((cells? raw-script) 18 (funcall cells-streamer raw-script)) 19 (t (mapcar cells-streamer raw-script))))) 20 (gen-tex-iter raw-script))) 21(defun flatten-raw-script (raw-script) 22 (let ((terms ())) 23 (gen-subscript (lambda (x) (push x terms)) raw-script) 24 (reverse terms))) 25 26(defun convert-from-n-sys (ls n) 27 (reduce (lambda (x y) (+ (* x n) y)) ls 28 :initial-value 0)) 29 30(defun xyz-to-ternary (xyzs) 31 (cond ((eql xyzs 'x) 0) 32 ((eql xyzs 'y) 1) 33 ((eql xyzs 'z) 2) 34 (t (error " unknown subscript ~a" xyzs)))) 35 36(defun ternary-subscript (ops) 37 "convert the polynomial xyz to the ternary" 38 (cond ((null ops) ops) 39 (t (convert-from-n-sys (mapcar #'xyz-to-ternary 40 (remove-if (lambda (x) (eql x 's)) 41 (scripts-of ops))) 42 3)))) 43(defun cell-converter (cell fout &optional (with-grids nil)) 44 (let ((fac (realpart (phase-of cell))) 45 (const@3 (ternary-subscript (consts-of cell))) 46 (op@3 (if with-grids 47 (format nil "ig+GRID_BLKSIZE*~a" (ternary-subscript (ops-of cell))) 48 (ternary-subscript (ops-of cell))))) 49 (cond ((equal fac 1) 50 (cond ((null const@3) 51 (if (null op@3) 52 (format fout " + s\[0\]" ) 53 (format fout " + s\[~a\]" op@3))) 54 ((null op@3) 55 (format fout " + c\[~a\]*s\[0\]" const@3)) 56 (t (format fout " + c\[~a\]*s\[~a\]" const@3 op@3)))) 57 ((equal fac -1) 58 (cond ((null const@3) 59 (if (null op@3) 60 (format fout " - s\[0\]" ) 61 (format fout " - s\[~a\]" op@3))) 62 ((null op@3) 63 (format fout " - c\[~a\]*s\[0\]" const@3)) 64 (t (format fout " - c\[~a\]*s\[~a\]" const@3 op@3)))) 65 ((< fac 0) 66 (cond ((null const@3) 67 (if (null op@3) 68 (format fout " ~a*s\[0\]" fac) 69 (format fout " ~a*s\[~a\]" fac op@3))) 70 ((null op@3) 71 (format fout " ~a*c\[~a\]*s\[0\]" fac const@3)) 72 (t (format fout " ~a*c\[~a\]*s\[~a\]" fac const@3 op@3)))) 73 (t 74 (cond ((null const@3) 75 (if (null op@3) 76 (format fout " + ~a*s\[0\]" fac) 77 (format fout " + ~a*s\[~a\]" fac op@3))) 78 ((null op@3) 79 (format fout " + ~a*c\[~a\]*s\[0\]" fac const@3)) 80 (t (format fout " + ~a*c\[~a\]*s\[~a\]" fac const@3 op@3))))))) 81 82(defun to-c-code-string (fout c-converter flat-script &optional (with-grids nil)) 83 (flet ((c-streamer (cs) 84 (with-output-to-string (tmpout) 85 (cond ((null cs) (format tmpout " 0")) 86 ((cell? cs) (funcall c-converter cs tmpout with-grids)) 87 (t (mapcar (lambda (c) (funcall c-converter c tmpout with-grids)) cs)))))) 88 (mapcar #'c-streamer flat-script))) 89 90(defun gen-c-block (fout flat-script &optional (with-grids nil)) 91 (let ((assemb (to-c-code-string fout #'cell-converter flat-script with-grids)) 92 (comp (length flat-script))) 93 (loop for s in assemb 94 for gid from 0 do 95 (if with-grids 96 (format fout "gout[ig+bgrids*(n*~a+~a)] =~a;~%" comp gid s) 97 (format fout "gout[n*~a+~a] =~a;~%" comp gid s))))) 98(defun gen-c-block+ (fout flat-script &optional (with-grids nil)) 99 (let ((assemb (to-c-code-string fout #'cell-converter flat-script with-grids)) 100 (comp (length flat-script))) 101 (loop for s in assemb 102 for gid from 0 do 103 (if with-grids 104 (format fout "gout[ig+bgrids*(n*~a+~a)] +=~a;~%" comp gid s) 105 (format fout "gout[n*~a+~a] +=~a;~%" comp gid s))))) 106(defun gen-c-block-with-empty (fout flat-script) 107 (format fout "if (gout_empty) {~%") 108 (gen-c-block fout flat-script) 109 (format fout "} else {~%") 110 (gen-c-block+ fout flat-script) 111 (format fout "}")) 112;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 113 114;;; effective keys are p,r,ri,... 115(defun effect-keys (ops) 116 (remove-if-not (lambda (x) 117 (or (member x *act-left-right*) 118 (member x *intvar*))) 119 ops)) 120(defun g?e-of (key) 121 (case key 122 ((p ip nabla px py pz p* ip* nabla* px* py* pz*) "D_") 123 ((r x y z) "R_") ; the vector origin is on the center of the basis it acts on 124 ((ri rj rk rl) "RC") ; the vector origin is R[ijkl] 125 ((xi xj xk xl) "RC") 126 ((yi yj yk yl) "RC") 127 ((zi zj zk zl) "RC") 128 ((r0 x0 y0 z0 g) "R0") ; R0 ~ the vector origin is (0,0,0) 129 ((rc xc yc zc) "RC") ; the vector origin is set in env[PTR_COMMON_ORIG] 130 ((nabla-rinv nabla-r12 breit-r1 breit-r2) "D_") 131 (otherwise (error "unknown key ~a~%" key)))) 132 133(defun dump-header (fout) 134 (format fout "/* 135 * Copyright (C) 2013- Qiming Sun <osirpt.sun@gmail.com> 136 * Description: code generated by gen-code.cl 137 */ 138#include <stdlib.h> 139#include <stdio.h> 140#include \"cint_bas.h\" 141#include \"cart2sph.h\" 142#include \"g1e.h\" 143#include \"g1e_grids.h\" 144#include \"g2e.h\" 145#include \"optimizer.h\" 146#include \"cint1e.h\" 147#include \"cint2e.h\" 148#include \"misc.h\" 149#include \"c2f.h\" 150")) 151 152(defun dump-declare-dri-for-rc (fout i-ops symb) 153 (when (intersection '(rc xc yc zc) i-ops) 154 (format fout "double dr~a[3];~%" symb) 155 (format fout "dr~a[0] = envs->r~a[0] - envs->env[PTR_COMMON_ORIG+0];~%" symb symb) 156 (format fout "dr~a[1] = envs->r~a[1] - envs->env[PTR_COMMON_ORIG+1];~%" symb symb) 157 (format fout "dr~a[2] = envs->r~a[2] - envs->env[PTR_COMMON_ORIG+2];~%" symb symb)) 158 (when (intersection '(ri xi yi zi) i-ops) 159 (if (intersection '(rc xc yc zc) i-ops) 160 (error "Cannot declare dri because rc and ri coexist")) 161 (format fout "double dr~a[3];~%" symb) 162 (format fout "dr~a[0] = envs->r~a[0] - envs->ri[0];~%" symb symb) 163 (format fout "dr~a[1] = envs->r~a[1] - envs->ri[1];~%" symb symb) 164 (format fout "dr~a[2] = envs->r~a[2] - envs->ri[2];~%" symb symb)) 165 (when (intersection '(rj xj yj zj) i-ops) 166 (if (intersection '(rc xc yc zc) i-ops) 167 (error "Cannot declare drj because rc and rj coexist")) 168 (format fout "double dr~a[3];~%" symb) 169 (format fout "dr~a[0] = envs->r~a[0] - envs->rj[0];~%" symb symb) 170 (format fout "dr~a[1] = envs->r~a[1] - envs->rj[1];~%" symb symb) 171 (format fout "dr~a[2] = envs->r~a[2] - envs->rj[2];~%" symb symb)) 172 (when (intersection '(rk xk yk zk) i-ops) 173 (if (intersection '(rc xc yc zc) i-ops) 174 (error "Cannot declare drk because rc and rk coexist")) 175 (format fout "double dr~a[3];~%" symb) 176 (format fout "dr~a[0] = envs->r~a[0] - envs->rk[0];~%" symb symb) 177 (format fout "dr~a[1] = envs->r~a[1] - envs->rk[1];~%" symb symb) 178 (format fout "dr~a[2] = envs->r~a[2] - envs->rk[2];~%" symb symb)) 179 (when (intersection '(rl xl yl zl) i-ops) 180 (if (intersection '(rc xc yc zc) i-ops) 181 (error "Cannot declare drl because rc and rl coexist")) 182 (format fout "double dr~a[3];~%" symb) 183 (format fout "dr~a[0] = envs->r~a[0] - envs->rl[0];~%" symb symb) 184 (format fout "dr~a[1] = envs->r~a[1] - envs->rl[1];~%" symb symb) 185 (format fout "dr~a[2] = envs->r~a[2] - envs->rl[2];~%" symb symb))) 186 187(defun dump-declare-giao-ij (fout bra ket) 188 (let ((n-giao (count 'g (append bra ket)))) 189 (when (> n-giao 0) 190 (format fout "double rirj[3], c[~a];~%" (expt 3 n-giao)) 191 (format fout "rirj[0] = envs->ri[0] - envs->rj[0];~%" ) 192 (format fout "rirj[1] = envs->ri[1] - envs->rj[1];~%" ) 193 (format fout "rirj[2] = envs->ri[2] - envs->rj[2];~%" ) 194 (loop 195 for i upto (1- (expt 3 n-giao)) do 196 (format fout "c[~a] = 1" i) 197 (loop 198 for j from (1- n-giao) downto 0 199 and res = i then (multiple-value-bind (int res) (floor res (expt 3 j)) 200 (format fout " * rirj[~a]" int) 201 res)) 202 (format fout ";~%"))))) 203 204; l-combo searches op_bit from left to right 205; o100 o010 o001|g...> => |g...> = o100 |g0..> 206; a operator can only be applied to the left of the existed ones 207; |g100,l> = o100 |g000,l+1> 208; |g101,l> = o100 |g001,l+1> 209; |g110,l> = o100 |g010,l+1> 210; |g111,l> = o100 |g011,l+1> 211; as a result, g* intermediates are generated from the previous one whose 212; id (in binary) can be obtained by removing the first bit 1 from current 213; id (in binary), see combo-bra function, eg 214; 000 g0 215; 001 g1 from g0 (000) 216; 010 g2 from g0 (000) 217; 011 g3 from g1 (001) 218; 100 g4 from g0 (000) 219; 101 g5 from g1 (001) 220; 110 g6 from g2 (010) 221; 111 g7 from g3 (011) 222; r-combo searches op_bit from right to left 223; o100 o010 o001|g...> => |g...> = o001 |g..0> 224; a operator can only be applied to the right of the exsited ones 225; |g100,l+2> = o100 |g000,l+3> 226; |g101,l > = o001 |g100,l+1> 227; |g110,l+1> = o010 |g100,l+2> 228; |g111,l > = o001 |g110,l+1> 229; as a result, g* intermediates are generated from the previous one whose 230; id (in binary) can be obtained by removing the last bit 1 from current 231; id (in binary), see combo-ket function, eg 232; 000 g0 233; 001 g1 from g0 (000) 234; 010 g2 from g0 (000) 235; 011 g3 from g2 (010) 236; 100 g4 from g0 (000) 237; 101 g5 from g4 (100) 238; 110 g6 from g4 (100) 239; 111 g7 from g6 (110) 240; [lr]-combo have no connection with <bra| or |ket> 241; def l_combinator(self, ops, ig, mask, template): 242(defun first-bit1 (n) 243 (loop 244 for i upto 31 245 thereis (if (zerop (ash n (- i))) (1- i)))) 246(defun last-bit1 (n) 247 ; how many 0s follow the last bit 1 248 (loop 249 for i upto 31 250 thereis (if (oddp (ash n (- i))) i))) 251(defun combo-bra (fout fmt ops-rev n-ops ig mask) 252 (let* ((right (first-bit1 (ash ig (- mask)))) 253 (left (- n-ops right 1)) 254 (ig0 (- ig (ash 1 (+ mask right)))) 255 (op (nth right ops-rev))) 256 (format fout fmt (g?e-of op) ig ig0 left))) 257(defun combo-ket (fout fmt ops-rev i-len ig mask) 258 (let* ((right (last-bit1 (ash ig (- mask)))) 259 (ig0 (- ig (ash 1 (+ mask right)))) 260 (op (nth right ops-rev))) 261 (format fout fmt (g?e-of op) ig ig0 i-len right))) 262(defun combo-opj (fout fmt-op fmt-j opj-rev i-len j-len ig mask) 263 (let ((right (last-bit1 (ash ig (- mask))))) 264 (if (< right j-len) ; does not reach op yet 265 (combo-ket fout fmt-j opj-rev i-len ig mask) 266 (let ((ig0 (- ig (ash 1 (+ mask right)))) 267 (op (nth right opj-rev))) 268 (if (member op *act-left-right*) 269 (format fout fmt-op 270 (g?e-of op) ig ig0 i-len right 271 (g?e-of op) (1+ ig) ig0 i-len right 272 ig (1+ ig)) 273 (if (and (intersection *act-left-right* opj-rev) 274 (< (1+ right) (length opj-rev))) ; ops have *act-left-right* but the rightmost op is not 275 (format fout fmt-j (g?e-of op) ig ig0 (1+ i-len) right) 276 (format fout fmt-j (g?e-of op) ig ig0 i-len right))))))) 277 278(defun power2-range (n &optional (shift 0)) 279 (range (+ shift (ash 1 n)) (+ shift (ash 1 (1+ n))))) 280(defun dump-combo-braket (fout fmt-i fmt-op fmt-j i-rev op-rev j-rev mask) 281 (let* ((i-len (length i-rev)) 282 (j-len (length j-rev)) 283 (op-len (length op-rev)) 284 (opj-rev (append j-rev op-rev))) 285 (loop 286 for right from mask to (+ mask j-len op-len -1) do 287 (loop 288 for ig in (power2-range right) do 289 (combo-opj fout fmt-op fmt-j opj-rev i-len j-len ig mask))) 290 (let ((shft (+ op-len j-len mask))) 291 (loop 292 for right from shft to (+ shft i-len -1) do 293 (loop 294 for ig in (power2-range right) do 295 (combo-bra fout fmt-i i-rev i-len ig shft)))))) 296 297(defun dec-to-ybin (n) 298 (parse-integer (substitute #\0 #\2 (write-to-string n :base 3)) 299 :radix 2)) 300(defun dec-to-zbin (n) 301 (parse-integer (substitute #\1 #\2 302 (substitute #\0 #\1 303 (write-to-string n :base 3))) 304 :radix 2)) 305(defun dump-s-1e (fout n) 306 (format fout "for (n = 0; n < nf; n++, idx+=3) { 307ix = idx[0]; 308iy = idx[1]; 309iz = idx[2];~%") 310 (loop 311 for i upto (1- (expt 3 n)) do 312 (let* ((ybin (dec-to-ybin i)) 313 (zbin (dec-to-zbin i)) 314 (xbin (- (ash 1 n) 1 ybin zbin))) 315 (format fout "s[~a] = g~a[ix] * g~a[iy] * g~a[iz];~%" 316 i xbin ybin zbin)))) 317 318(defun name-c2sor (sfx sp sf ts) 319 (cond ((eql sp 'spinor) 320 (if (eql sf 'sf) 321 (if (eql ts 'ts) 322 (format nil "&c2s_sf_~a" sfx) 323 (format nil "&c2s_sf_~ai" sfx)) 324 (if (eql ts 'ts) 325 (format nil "&c2s_si_~a" sfx) 326 (format nil "&c2s_si_~ai" sfx)))) 327 ((eql sp 'spheric) 328 (format nil "&c2s_sph_~a" sfx)) 329 (t (format nil "&c2s_cart_~a" sfx)))) 330 331;;; g-intermediates are g0, g1, g2, ... 332(defun num-g-intermediates (tot-bits op i-len j-len) 333 (if (and (intersection *act-left-right* op) 334 ; nabla-rinv is the last one in the operation list 335 (<= (+ i-len (length op) j-len) 1)) 336 (ash 1 tot-bits) 337 (1- (ash 1 tot-bits)))) 338 339;;; Output function signatures for C header file 340(defun write-functype-header (fout intname &optional doc) 341 (if doc 342 (format fout "~a~%" doc)) 343 (format fout "extern CINTOptimizerFunction ~a_optimizer; 344extern CINTIntegralFunction ~a_cart; 345extern CINTIntegralFunction ~a_sph; 346extern CINTIntegralFunction ~a_spinor;~%~%" intname intname intname intname)) 347 348;!!! Be very cautious of the reverse order on i-operators and k operators! 349;!!! When multiple tensor components (>= rank 2) provided by the operators 350;!!! on bra functions, the ordering of the multiple tensor components are 351;!!! also reversed in the generated integral code 352(defun gen-code-gout1e (fout intname raw-infix flat-script) 353 (destructuring-bind (op bra-i ket-j bra-k ket-l) 354 (split-int-expression raw-infix) 355 (let* ((i-rev (effect-keys bra-i)) ;<i| already in reverse order 356 (j-rev (reverse (effect-keys ket-j))) 357 (op-rev (reverse (effect-keys op))) 358 (i-len (length i-rev)) 359 (j-len (length j-rev)) 360 (op-len (length op-rev)) 361 (tot-bits (+ i-len j-len op-len)) 362 (goutinc (length flat-script))) 363 (format fout "void CINTgout1e_~a" intname) 364 (format fout "(double *gout, double *g, FINT *idx, CINTEnvVars *envs, FINT gout_empty) { 365FINT nf = envs->nf; 366FINT ix, iy, iz, n; 367double *g0 = g;~%") 368 (loop 369 for i in (range (num-g-intermediates tot-bits op i-len j-len)) do 370 (format fout "double *g~a = g~a + envs->g_size * 3;~%" (1+ i) i)) 371 (dump-declare-dri-for-rc fout bra-i "i") 372 (dump-declare-dri-for-rc fout (append op ket-j) "j") 373 (dump-declare-giao-ij fout bra-i (append op ket-j)) 374 (format fout "double s[~a];~%" (expt 3 tot-bits)) 375;;; generate g_(bin) 376;;; for the operators act on the |ket>, the reversed scan order and r_combinator 377;;; is required; for the operators acto on the <bra|, the normal scan order and 378 (let ((fmt-i "G1E_~aI(g~a, g~a, envs->i_l+~a, envs->j_l, 0);~%") 379 (fmt-op (mkstr "G1E_~aJ(g~a, g~a, envs->i_l+~d, envs->j_l+~a, 0); 380G1E_~aI(g~a, g~a, envs->i_l+~d, envs->j_l+~a, 0); 381for (ix = 0; ix < envs->g_size * 3; ix++) {g~a[ix] += g~a[ix];}~%")) 382 (fmt-j (mkstr "G1E_~aJ(g~a, g~a, envs->i_l+~d, envs->j_l+~a, 0);~%"))) 383 (dump-combo-braket fout fmt-i fmt-op fmt-j i-rev op-rev j-rev 0)) 384 (format fout "for (n = 0; n < nf; n++) { 385ix = idx[0+n*3]; 386iy = idx[1+n*3]; 387iz = idx[2+n*3];~%") 388 (dump-s-for-nroots fout tot-bits 1) 389 (gen-c-block-with-empty fout flat-script) 390 (format fout "}}~%") 391 goutinc))) 392 393(defun gen-code-gout1e-rinv (fout intname raw-infix flat-script) 394 (destructuring-bind (op bra-i ket-j bra-k ket-l) 395 (split-int-expression raw-infix) 396 (let* ((i-rev (effect-keys bra-i)) ;<i| already in reverse order 397 (j-rev (reverse (effect-keys ket-j))) 398 (op-rev (reverse (effect-keys op))) 399 (i-len (length i-rev)) 400 (j-len (length j-rev)) 401 (op-len (length op-rev)) 402 (tot-bits (+ i-len j-len op-len)) 403 (goutinc (length flat-script))) 404 (format fout "void CINTgout1e_~a" intname) 405 (format fout "(double *gout, double *g, FINT *idx, CINTEnvVars *envs, FINT gout_empty) { 406FINT nf = envs->nf; 407FINT nrys_roots = envs->nrys_roots; 408FINT ix, iy, iz, n, i; 409double *g0 = g;~%") 410 (loop 411 for i in (range (num-g-intermediates tot-bits op i-len j-len)) do 412 (format fout "double *g~a = g~a + envs->g_size * 3;~%" (1+ i) i)) 413 (dump-declare-dri-for-rc fout bra-i "i") 414 (dump-declare-dri-for-rc fout (append op ket-j) "j") 415 (dump-declare-giao-ij fout bra-i (append op ket-j)) 416;;; generate g_(bin) 417;;; for the operators act on the |ket>, the reversed scan order and r_combinator 418;;; is required; for the operators acto on the <bra|, the normal scan order and 419 (let ((fmt-i "G2E_~aI(g~a, g~a, envs->i_l+~a, envs->j_l, 0, 0);~%") 420 (fmt-op (mkstr "G2E_~aJ(g~a, g~a, envs->i_l+~d, envs->j_l+~a, 0, 0); 421G2E_~aI(g~a, g~a, envs->i_l+~d, envs->j_l+~a, 0, 0); 422for (ix = 0; ix < envs->g_size * 3; ix++) {g~a[ix] += g~a[ix];}~%")) 423 (fmt-j (mkstr "G2E_~aJ(g~a, g~a, envs->i_l+~d, envs->j_l+~a, 0, 0);~%"))) 424 (dump-combo-braket fout fmt-i fmt-op fmt-j i-rev op-rev j-rev 0)) 425 (dump-s-2e fout tot-bits 0) 426 (gen-c-block-with-empty fout flat-script) 427 (format fout "}}~%") 428 goutinc))) 429 430(defun gen-code-int1e (fout intname raw-infix) 431 (destructuring-bind (op bra-i ket-j bra-k ket-l) 432 (split-int-expression raw-infix) 433 (let* ((i-rev (effect-keys bra-i)) ;<i| already in reverse order 434 (j-rev (reverse (effect-keys ket-j))) 435 (op-rev (reverse (effect-keys op))) 436 (i-len (length i-rev)) 437 (j-len (length j-rev)) 438 (op-len (length op-rev)) 439 (tot-bits (+ i-len j-len op-len)) 440 (raw-script (eval-int raw-infix)) 441 (flat-script (flatten-raw-script (last1 raw-script))) 442 (ts (car raw-script)) 443 (sf (cadr raw-script)) 444 (goutinc (length flat-script)) 445 (e1comps (if (eql sf 'sf) 1 4)) 446 (tensors (if (eql sf 'sf) goutinc (/ goutinc 4))) 447 (int1e-type (cond ((member 'nuc raw-infix) 2) 448 ((or (member 'rinv raw-infix) 449 (member 'nabla-rinv raw-infix)) 1) 450 (t 0))) 451 (ngdef (with-output-to-string (tmpout) 452 (if (or (member 'nuc raw-infix) 453 (member 'rinv raw-infix) 454 (member 'nabla-rinv raw-infix)) 455 (if (intersection *act-left-right* op) 456 (format tmpout "FINT ng[] = {~d, ~d, 0, 0, ~d, ~d, 0, ~d};~%" 457 (1+ i-len) (+ op-len j-len) tot-bits e1comps tensors) 458 (format tmpout "FINT ng[] = {~d, ~d, 0, 0, ~d, ~d, 0, ~d};~%" 459 i-len (+ op-len j-len) tot-bits e1comps tensors)) 460 (format tmpout "FINT ng[] = {~d, ~d, 0, 0, ~d, ~d, 1, ~d};~%" 461 i-len (+ op-len j-len) tot-bits e1comps tensors)))) 462 (envs-common (with-output-to-string (tmpout) 463 (format tmpout ngdef) 464 (format tmpout "CINTEnvVars envs;~%") 465 (format tmpout "CINTinit_int1e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);~%") 466 (format tmpout "envs.f_gout = &CINTgout1e_~a;~%" intname) 467 (unless (eql (factor-of raw-infix) 1) 468 (format tmpout "envs.common_factor *= ~a;~%" (factor-of raw-infix)))))) 469 (write-functype-header 470 t intname (format nil "/* <~{~a ~}i|~{~a ~}|~{~a ~}j> */" bra-i op ket-j)) 471 (format fout "/* <~{~a ~}i|~{~a ~}|~{~a ~}j> */~%" bra-i op ket-j) 472 (cond ((or (member 'nuc raw-infix) 473 (member 'rinv raw-infix) 474 (member 'nabla-rinv raw-infix)) 475 (gen-code-gout1e-rinv fout intname raw-infix flat-script)) 476 (t (gen-code-gout1e fout intname raw-infix flat-script))) 477 (format fout "void ~a_optimizer(CINTOpt **opt, FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env) {~%" intname) 478 (format fout ngdef) 479 (format fout "CINTall_1e_optimizer(opt, ng, atm, natm, bas, nbas, env);~%}~%") 480;;; _cart 481 (format fout "CACHE_SIZE_T ~a_cart(double *out, FINT *dims, FINT *shls, 482FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {~%" intname) 483 (format fout envs-common) 484 (when (member 'g raw-infix) 485 (when (or (member 'g bra-i) (member 'g ket-j)) 486 (format fout "if (out != NULL && envs.shls[0] == envs.shls[1]) { 487FINT i, nout; 488FINT counts[4]; 489counts[0] = envs.nfi * envs.x_ctr[0]; 490counts[1] = envs.nfj * envs.x_ctr[1]; 491counts[2] = 1; 492counts[3] = 1; 493if (dims == NULL) { dims = counts; } 494nout = dims[0] * dims[1]; 495for (i = 0; i < envs.ncomp_e1 * envs.ncomp_tensor; i++) { 496c2s_dset0(out+nout*i, dims, counts); } 497return 0; }~%"))) 498 (format fout "return CINT1e_drv(out, dims, &envs, cache, &c2s_cart_1e, ~d); 499} // ~a_cart~%" int1e-type intname) 500;;; _sph 501 (format fout "CACHE_SIZE_T ~a_sph(double *out, FINT *dims, FINT *shls, 502FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {~%" intname) 503 (format fout envs-common) 504 (when (member 'g raw-infix) 505 (when (or (member 'g bra-i) (member 'g ket-j)) 506 (format fout "if (out != NULL && envs.shls[0] == envs.shls[1]) { 507FINT i, nout; 508FINT counts[4]; 509counts[0] = (envs.i_l*2+1) * envs.x_ctr[0]; 510counts[1] = (envs.j_l*2+1) * envs.x_ctr[1]; 511counts[2] = 1; 512counts[3] = 1; 513if (dims == NULL) { dims = counts; } 514nout = dims[0] * dims[1]; 515for (i = 0; i < envs.ncomp_e1 * envs.ncomp_tensor; i++) { 516c2s_dset0(out+nout*i, dims, counts); } 517return 0; }~%"))) 518 (format fout "return CINT1e_drv(out, dims, &envs, cache, &c2s_sph_1e, ~d); 519} // ~a_sph~%" int1e-type intname) 520;;; _spinor 521 (format fout "CACHE_SIZE_T ~a_spinor(double complex *out, FINT *dims, FINT *shls, 522FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {~%" intname) 523 (format fout envs-common) 524 (when (member 'g raw-infix) 525 (when (or (member 'g bra-i) (member 'g ket-j)) 526 (format fout "if (out != NULL && envs.shls[0] == envs.shls[1]) { 527FINT i, nout; 528FINT counts[4]; 529counts[0] = CINTcgto_spinor(envs.shls[0], envs.bas); 530counts[1] = CINTcgto_spinor(envs.shls[1], envs.bas); 531counts[2] = 1; 532counts[3] = 1; 533if (dims == NULL) { dims = counts; } 534nout = dims[0] * dims[1]; 535for (i = 0; i < envs.ncomp_tensor; i++) { 536c2s_zset0(out+nout*i, dims, counts); } 537return 0; }~%"))) 538 (format fout "return CINT1e_spinor_drv(out, dims, &envs, cache, ~a, ~d); 539} // ~a_spinor~%" (name-c2sor "1e" 'spinor sf ts) int1e-type intname))) 540;;; int1e -> cint1e 541 (format fout "ALL_CINT1E(~a)~%" intname) 542 (format fout "ALL_CINT1E_FORTRAN_(~a)~%" intname)) 543;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 544 545(defun dump-declare-giao-ijkl (fout opi opj opk opl) 546 (let ((n-ij (count 'g (append opi opj))) 547 (n-kl (count 'g (append opk opl)))) 548 (when (> n-ij 0) 549 (format fout "double rirj[3];~%") 550 (format fout "rirj[0] = envs->ri[0] - envs->rj[0];~%" ) 551 (format fout "rirj[1] = envs->ri[1] - envs->rj[1];~%" ) 552 (format fout "rirj[2] = envs->ri[2] - envs->rj[2];~%" )) 553 (when (> n-kl 0) 554 (format fout "double rkrl[3];~%") 555 (format fout "rkrl[0] = envs->rk[0] - envs->rl[0];~%" ) 556 (format fout "rkrl[1] = envs->rk[1] - envs->rl[1];~%" ) 557 (format fout "rkrl[2] = envs->rk[2] - envs->rl[2];~%" )) 558 (when (> (+ n-ij n-kl) 0) 559 (format fout "double c[~a];~%" (expt 3 (+ n-ij n-kl))) 560 (loop 561 for i upto (1- (expt 3 (+ n-ij n-kl))) do 562 (format fout "c[~a] = 1" i) 563 (loop 564 for j from (+ n-ij n-kl -1) downto n-kl 565 and res = i then (multiple-value-bind (int res) (floor res (expt 3 j)) 566 (format fout " * rirj[~a]" int) 567 res)) 568 (loop 569 for j from (1- n-kl) downto 0 570 and res = (nth-value 1 (floor i (expt 3 n-kl))) 571 then (multiple-value-bind (int res) (floor res (expt 3 j)) 572 (format fout " * rkrl[~a]" int) 573 res)) 574 (format fout ";~%"))))) 575 576(defun dump-s-for-nroots (fout tot-bits nroots) 577 (loop 578 for i upto (1- (expt 3 tot-bits)) do 579 (let* ((ybin (dec-to-ybin i)) 580 (zbin (dec-to-zbin i)) 581 (xbin (- (ash 1 tot-bits) 1 ybin zbin))) 582 (format fout "s[~a] = " i) 583 (loop 584 for k upto (1- nroots) do 585 (format fout "+ g~a[ix+~a]*g~a[iy+~a]*g~a[iz+~a]" 586 xbin k ybin k zbin k)) 587 (format fout ";~%")))) 588(defun dump-s-loop (fout tot-bits &optional (with-grids nil)) 589 (loop 590 for i upto (1- (expt 3 tot-bits)) do 591 (let* ((ybin (dec-to-ybin i)) 592 (zbin (dec-to-zbin i)) 593 (xbin (- (ash 1 tot-bits) 1 ybin zbin))) 594 (if with-grids 595 (format fout "s[ig+GRID_BLKSIZE*~a] += g~a[ix+ig+i*GRID_BLKSIZE] * g~a[iy+ig+i*GRID_BLKSIZE] * g~a[iz+ig+i*GRID_BLKSIZE];~%" 596 i xbin ybin zbin) 597 (format fout "s[~a] += g~a[ix+i] * g~a[iy+i] * g~a[iz+i];~%" 598 i xbin ybin zbin))))) ; end do i = 1, envs->nrys_roots 599(defun dump-s-2e (fout tot-bits &optional (deriv-max 2)) 600 (format fout "double s[~a];~%" (expt 3 tot-bits)) 601 (format fout "for (n = 0; n < nf; n++) { 602ix = idx[0+n*3]; 603iy = idx[1+n*3]; 604iz = idx[2+n*3];~%") 605 (if (< tot-bits deriv-max) 606 (progn 607 (format fout "switch (nrys_roots) {~%") 608 (loop 609 for i from 1 to 4 do 610 (format fout "case ~a:~%" i) 611 (dump-s-for-nroots fout tot-bits i) 612 (format fout "break;~%" )) 613 (format fout "default: 614for (i = 0; i < ~a; i++) { s[i] = 0; } 615for (i = 0; i < nrys_roots; i++) {~%" (expt 3 tot-bits)) 616 (dump-s-loop fout tot-bits) 617 (format fout "} break;}~%")) ; else 618 (progn 619 (format fout "for (i = 0; i < ~a; i++) { s[i] = 0; } 620for (i = 0; i < nrys_roots; i++) {~%" (expt 3 tot-bits)) 621 (dump-s-loop fout tot-bits) 622 (format fout "}~%")))) 623 624;;; generate function gout2e 625(defun gen-code-gout4c2e (fout intname raw-infix flat-script) 626 (destructuring-bind (op bra-i ket-j bra-k ket-l) 627 (split-int-expression raw-infix) 628 (let* ((i-rev (effect-keys bra-i)) 629 (j-rev (reverse (effect-keys ket-j))) 630 (k-rev (effect-keys bra-k)) 631 (l-rev (reverse (effect-keys ket-l))) 632 (op-rev (reverse (effect-keys op))) 633 (op-len (length op-rev)) 634 (i-len (length i-rev)) 635 (j-len (length j-rev)) 636 (k-len (length k-rev)) 637 (l-len (length l-rev)) 638 (tot-bits (+ i-len j-len op-len k-len l-len)) 639 (goutinc (length flat-script))) 640 (format fout "void CINTgout2e_~a(double *gout, 641double *g, FINT *idx, CINTEnvVars *envs, FINT gout_empty) {~%" intname) 642 (format fout "FINT nf = envs->nf; 643FINT nrys_roots = envs->nrys_roots; 644FINT ix, iy, iz, i, n; 645double *g0 = g;~%") 646 (loop 647 for i in (range (num-g-intermediates tot-bits op i-len j-len)) do 648 (format fout "double *g~a = g~a + envs->g_size * 3;~%" (1+ i) i)) 649 (dump-declare-dri-for-rc fout bra-i "i") 650 (dump-declare-dri-for-rc fout ket-j "j") 651 (dump-declare-dri-for-rc fout bra-k "k") 652 (dump-declare-dri-for-rc fout ket-l "l") 653 (dump-declare-giao-ijkl fout bra-i ket-j bra-k ket-l) 654 (if (intersection *act-left-right* op) 655 (let ((fmt-k (mkstr "G2E_~aK(g~a, g~a, envs->i_l+" (1+ i-len) ", envs->j_l+" (1+ j-len) 656 ", envs->k_l+~a, envs->l_l);~%")) 657 (fmt-op "") 658 (fmt-l (mkstr "G2E_~aL(g~a, g~a, envs->i_l+" (1+ i-len) ", envs->j_l+" (1+ j-len) 659 ", envs->k_l+~a, envs->l_l+~a);~%"))) 660 (dump-combo-braket fout fmt-k fmt-op fmt-l k-rev '() l-rev 0)) 661 (let ((fmt-k (mkstr "G2E_~aK(g~a, g~a, envs->i_l+" i-len ", envs->j_l+" (+ op-len j-len) 662 ", envs->k_l+~a, envs->l_l);~%")) 663 (fmt-op "") 664 (fmt-l (mkstr "G2E_~aL(g~a, g~a, envs->i_l+" i-len ", envs->j_l+" (+ op-len j-len) 665 ", envs->k_l+~a, envs->l_l+~a);~%"))) 666 (dump-combo-braket fout fmt-k fmt-op fmt-l k-rev '() l-rev 0))) 667 (let ((fmt-i "G2E_~aI(g~a, g~a, envs->i_l+~a, envs->j_l, envs->k_l, envs->l_l);~%") 668 (fmt-op (mkstr "G2E_~aJ(g~a, g~a, envs->i_l+~d, envs->j_l+~a, envs->k_l, envs->l_l); 669G2E_~aI(g~a, g~a, envs->i_l+~a, envs->j_l+~a, envs->k_l, envs->l_l); 670for (ix = 0; ix < envs->g_size * 3; ix++) {g~a[ix] += g~a[ix];}~%")) 671 (fmt-j (mkstr "G2E_~aJ(g~a, g~a, envs->i_l+~d, envs->j_l+~a, envs->k_l, envs->l_l);~%"))) 672 (dump-combo-braket fout fmt-i fmt-op fmt-j i-rev op-rev j-rev (+ k-len l-len))) 673 (dump-s-2e fout tot-bits) 674 (gen-c-block-with-empty fout flat-script) 675 (format fout "}}~%") 676 goutinc))) 677 678(defun gen-code-int4c2e (fout intname raw-infix) 679 (destructuring-bind (op bra-i ket-j bra-k ket-l) 680 (split-int-expression raw-infix) 681 (let* ((i-rev (effect-keys bra-i)) 682 (j-rev (reverse (effect-keys ket-j))) 683 (k-rev (effect-keys bra-k)) 684 (l-rev (reverse (effect-keys ket-l))) 685 (op-rev (reverse (effect-keys op))) 686 (op-len (length op-rev)) 687 (i-len (length i-rev)) 688 (j-len (length j-rev)) 689 (k-len (length k-rev)) 690 (l-len (length l-rev)) 691 (tot-bits (+ i-len j-len op-len k-len l-len)) 692 (raw-script (eval-int raw-infix)) 693 (ts1 (car raw-script)) 694 (sf1 (cadr raw-script)) 695 (ts2 (caddr raw-script)) 696 (sf2 (cadddr raw-script)) 697 (flat-script (flatten-raw-script (last1 raw-script))) 698 (goutinc (length flat-script)) 699 (e1comps (if (eql sf1 'sf) 1 4)) 700 (e2comps (if (eql sf2 'sf) 1 4)) 701 (tensors (cond ((and (eql sf1 'sf) (eql sf2 'sf)) goutinc) 702 ((and (eql sf1 'si) (eql sf2 'si)) (/ goutinc 16)) 703 (t (/ goutinc 4)))) 704 (ngdef (with-output-to-string (tmpout) 705 (if (intersection *act-left-right* op) 706 (format tmpout "FINT ng[] = {~d, ~d, ~d, ~d, ~d, ~d, ~d, ~d};~%" 707 (1+ i-len) (1+ j-len) k-len l-len tot-bits e1comps e2comps tensors) 708 (format tmpout "FINT ng[] = {~d, ~d, ~d, ~d, ~d, ~d, ~d, ~d};~%" 709 i-len (+ op-len j-len) k-len l-len tot-bits e1comps e2comps tensors)))) 710 (envs-common (with-output-to-string (tmpout) 711 (format tmpout ngdef) 712 (format tmpout "CINTEnvVars envs;~%") 713 (format tmpout "CINTinit_int2e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);~%") 714 (format tmpout "envs.f_gout = &CINTgout2e_~a;~%" intname) 715 (unless (eql (factor-of raw-infix) 1) 716 (format tmpout "envs.common_factor *= ~a;~%" (factor-of raw-infix)))))) 717 (write-functype-header 718 t intname (format nil "/* (~{~a ~}i ~{~a ~}j|~{~a ~}|~{~a ~}k ~{~a ~}l) */" 719 bra-i ket-j op bra-k ket-l)) 720 (format fout "/* <~{~a ~}k ~{~a ~}i|~{~a ~}|~{~a ~}j ~{~a ~}l> : i,j \\in electron 1; k,l \\in electron 2~%" 721 bra-k bra-i op ket-j ket-l) 722 (format fout " * = (~{~a ~}i ~{~a ~}j|~{~a ~}|~{~a ~}k ~{~a ~}l) */~%" 723 bra-i ket-j op bra-k ket-l) 724 (gen-code-gout4c2e fout intname raw-infix flat-script) 725 (format fout "void ~a_optimizer(CINTOpt **opt, FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env) {~%" intname) 726 (format fout ngdef) 727 (format fout "CINTall_2e_optimizer(opt, ng, atm, natm, bas, nbas, env);~%}~%") 728;;; _cart 729 (format fout "CACHE_SIZE_T ~a_cart(double *out, FINT *dims, FINT *shls, 730FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {~%" intname) 731 (format fout envs-common) 732 (when (and (eql sf1 'si) (eql ts1 'tas) 733 (eql sf2 'si) (eql ts2 'tas)) 734 (format fout "envs.common_factor *= -1;~%")) 735 (when (member 'g raw-infix) 736 (format fout "FINT i, nout; 737FINT counts[4];~%") 738 (when (or (member 'g bra-i) (member 'g ket-j)) 739 (format fout "if (out != NULL && envs.shls[0] == envs.shls[1]) { 740counts[0] = envs.nfi * envs.x_ctr[0]; 741counts[1] = envs.nfj * envs.x_ctr[1]; 742counts[2] = envs.nfk * envs.x_ctr[2]; 743counts[3] = envs.nfl * envs.x_ctr[3]; 744if (dims == NULL) { dims = counts; } 745nout = dims[0] * dims[1] * dims[2] * dims[3]; 746for (i = 0; i < envs.ncomp_e1 * envs.ncomp_e2 * envs.ncomp_tensor; i++) { 747c2s_dset0(out+nout*i, dims, counts); } 748return 0; }~%")) 749 (when (or (member 'g bra-k) (member 'g ket-l)) 750 (format fout "if (out != NULL && envs.shls[2] == envs.shls[3]) { 751counts[0] = envs.nfi * envs.x_ctr[0]; 752counts[1] = envs.nfj * envs.x_ctr[1]; 753counts[2] = envs.nfk * envs.x_ctr[2]; 754counts[3] = envs.nfl * envs.x_ctr[3]; 755if (dims == NULL) { dims = counts; } 756nout = dims[0] * dims[1] * dims[2] * dims[3]; 757for (i = 0; i < envs.ncomp_e1 * envs.ncomp_e2 * envs.ncomp_tensor; i++) { 758c2s_dset0(out+nout*i, dims, counts); } 759return 0; }~%"))) 760 (format fout "return CINT2e_drv(out, dims, &envs, opt, cache, &c2s_cart_2e1);~%} // ~a_cart~%" intname) 761;;; _sph 762 (format fout "CACHE_SIZE_T ~a_sph(double *out, FINT *dims, FINT *shls, 763FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {~%" intname) 764 (format fout envs-common) 765 (when (and (eql sf1 'si) (eql ts1 'tas) 766 (eql sf2 'si) (eql ts2 'tas)) 767 (format fout "envs.common_factor *= -1;~%")) 768 (when (member 'g raw-infix) 769 (format fout "FINT i, nout; 770FINT counts[4];~%") 771 (when (or (member 'g bra-i) (member 'g ket-j)) 772 (format fout "if (out != NULL && envs.shls[0] == envs.shls[1]) { 773counts[0] = (envs.i_l*2+1) * envs.x_ctr[0]; 774counts[1] = (envs.j_l*2+1) * envs.x_ctr[1]; 775counts[2] = (envs.k_l*2+1) * envs.x_ctr[2]; 776counts[3] = (envs.l_l*2+1) * envs.x_ctr[3]; 777if (dims == NULL) { dims = counts; } 778nout = dims[0] * dims[1] * dims[2] * dims[3]; 779for (i = 0; i < envs.ncomp_e1 * envs.ncomp_e2 * envs.ncomp_tensor; i++) { 780c2s_dset0(out+nout*i, dims, counts); } 781return 0; }~%")) 782 (when (or (member 'g bra-k) (member 'g ket-l)) 783 (format fout "if (out != NULL && envs.shls[2] == envs.shls[3]) { 784counts[0] = (envs.i_l*2+1) * envs.x_ctr[0]; 785counts[1] = (envs.j_l*2+1) * envs.x_ctr[1]; 786counts[2] = (envs.k_l*2+1) * envs.x_ctr[2]; 787counts[3] = (envs.l_l*2+1) * envs.x_ctr[3]; 788if (dims == NULL) { dims = counts; } 789nout = dims[0] * dims[1] * dims[2] * dims[3]; 790for (i = 0; i < envs.ncomp_e1 * envs.ncomp_e2 * envs.ncomp_tensor; i++) { 791c2s_dset0(out+nout*i, dims, counts); } 792return 0; }~%"))) 793 (format fout "return CINT2e_drv(out, dims, &envs, opt, cache, &c2s_sph_2e1);~%} // ~a_sph~%" intname) 794;;; _spinor 795 (format fout "CACHE_SIZE_T ~a_spinor(double complex *out, FINT *dims, FINT *shls, 796FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {~%" intname) 797 (format fout envs-common) 798 (when (member 'g raw-infix) 799 (format fout "FINT i, nout; 800FINT counts[4];~%") 801 (when (or (member 'g bra-i) (member 'g ket-j)) 802 (format fout "if (out != NULL && envs.shls[0] == envs.shls[1]) { 803counts[0] = CINTcgto_spinor(envs.shls[0], envs.bas); 804counts[1] = CINTcgto_spinor(envs.shls[1], envs.bas); 805counts[2] = CINTcgto_spinor(envs.shls[2], envs.bas); 806counts[3] = CINTcgto_spinor(envs.shls[3], envs.bas); 807if (dims == NULL) { dims = counts; } 808nout = dims[0] * dims[1] * dims[2] * dims[3]; 809for (i = 0; i < envs.ncomp_tensor; i++) { 810c2s_zset0(out+nout*i, dims, counts); } 811return 0; }~%")) 812 (when (or (member 'g bra-k) (member 'g ket-l)) 813 (format fout "if (out != NULL && envs.shls[2] == envs.shls[3]) { 814counts[0] = CINTcgto_spinor(envs.shls[0], envs.bas); 815counts[1] = CINTcgto_spinor(envs.shls[1], envs.bas); 816counts[2] = CINTcgto_spinor(envs.shls[2], envs.bas); 817counts[3] = CINTcgto_spinor(envs.shls[3], envs.bas); 818if (dims == NULL) { dims = counts; } 819nout = dims[0] * dims[1] * dims[2] * dims[3]; 820for (i = 0; i < envs.ncomp_tensor; i++) { 821c2s_zset0(out+nout*i, dims, counts); } 822return 0; }~%"))) 823 (format fout "return CINT2e_spinor_drv(out, dims, &envs, opt, cache, ~a, ~a); 824} // ~a_spinor~%" (name-c2sor "2e1" 'spinor sf1 ts1) (name-c2sor "2e2" 'spinor sf2 ts2) intname))) 825;;; int2e -> cint2e 826 (format fout "ALL_CINT(~a)~%" intname) 827 (format fout "ALL_CINT_FORTRAN_(~a)~%" intname)) 828 829(defun gen-code-gout3c2e (fout intname raw-infix flat-script) 830 (destructuring-bind (op bra-i ket-j bra-k ket-l) 831 (split-int-expression raw-infix) 832 (let* ((i-rev (effect-keys bra-i)) 833 (j-rev (reverse (effect-keys ket-j))) 834 (k-rev (effect-keys bra-k)) 835 (op-rev (reverse (effect-keys op))) 836 (op-len (length op-rev)) 837 (i-len (length i-rev)) 838 (j-len (length j-rev)) 839 (k-len (length k-rev)) 840 (tot-bits (+ i-len j-len op-len k-len)) 841 (goutinc (length flat-script))) 842 (format fout "void CINTgout2e_~a(double *gout, 843double *g, FINT *idx, CINTEnvVars *envs, FINT gout_empty) {~%" intname) 844 (format fout "FINT nf = envs->nf; 845FINT nrys_roots = envs->nrys_roots; 846FINT ix, iy, iz, i, n; 847double *g0 = g;~%") 848 (loop 849 for i in (range (num-g-intermediates tot-bits op i-len j-len)) do 850 (format fout "double *g~a = g~a + envs->g_size * 3;~%" (1+ i) i)) 851 (dump-declare-dri-for-rc fout bra-i "i") 852 (dump-declare-dri-for-rc fout ket-j "j") 853 (dump-declare-dri-for-rc fout bra-k "k") 854 (dump-declare-giao-ij fout bra-i ket-j) 855 (if (intersection *act-left-right* op) 856 (let ((fmt-k (mkstr "G2E_~aK(g~a, g~a, envs->i_l+" (1+ i-len) ", envs->j_l+" (1+ j-len) 857 ", envs->k_l+~a, 0);~%")) 858 (fmt-op "") 859 (fmt-l "")) 860 (dump-combo-braket fout fmt-k fmt-op fmt-l k-rev '() '() 0)) 861 (let ((fmt-k (mkstr "G2E_~aK(g~a, g~a, envs->i_l+" i-len ", envs->j_l+" (+ op-len j-len) 862 ", envs->k_l+~a, 0);~%")) 863 (fmt-op "") 864 (fmt-l "")) 865 (dump-combo-braket fout fmt-k fmt-op fmt-l k-rev '() '() 0))) 866 (let ((fmt-i "G2E_~aI(g~a, g~a, envs->i_l+~a, envs->j_l, envs->k_l, 0);~%") 867 (fmt-op (mkstr "G2E_~aJ(g~a, g~a, envs->i_l+~d, envs->j_l+~a, envs->k_l, 0); 868G2E_~aI(g~a, g~a, envs->i_l+~a, envs->j_l+~a, envs->k_l, 0); 869for (ix = 0; ix < envs->g_size * 3; ix++) {g~a[ix] += g~a[ix];}~%")) 870 (fmt-j (mkstr "G2E_~aJ(g~a, g~a, envs->i_l+~d, envs->j_l+~a, envs->k_l, 0);~%"))) 871 (dump-combo-braket fout fmt-i fmt-op fmt-j i-rev op-rev j-rev k-len)) 872 (dump-s-2e fout tot-bits) 873 (gen-c-block-with-empty fout flat-script) 874 (format fout "}}~%") 875 goutinc))) 876 877(defun gen-code-int3c2e (fout intname raw-infix) 878 (destructuring-bind (op bra-i ket-j bra-k ket-l) 879 (split-int-expression raw-infix) 880 (let* ((i-rev (effect-keys bra-i)) 881 (j-rev (reverse (effect-keys ket-j))) 882 (k-rev (effect-keys bra-k)) 883 (op-rev (reverse (effect-keys op))) 884 (i-len (length i-rev)) 885 (j-len (length j-rev)) 886 (k-len (length k-rev)) 887 (op-len (length op-rev)) 888 (tot-bits (+ i-len j-len op-len k-len)) 889 (raw-script (eval-int raw-infix)) 890 (ts1 (car raw-script)) 891 (sf1 (cadr raw-script)) 892 (ts2 (caddr raw-script)) 893 (sf2 (cadddr raw-script)) 894 (flat-script (flatten-raw-script (last1 raw-script))) 895 (goutinc (length flat-script)) 896 (e1comps (if (eql sf1 'sf) 1 4)) 897 (e2comps (if (eql sf2 'sf) 1 4)) 898 (tensors (cond ((and (eql sf1 'sf) (eql sf2 'sf)) goutinc) 899 ((and (eql sf1 'si) (eql sf2 'si)) (/ goutinc 16)) 900 (t (/ goutinc 4)))) 901 (ngdef (with-output-to-string (tmpout) 902 (if (intersection *act-left-right* op) 903 (format tmpout "FINT ng[] = {~d, ~d, ~d, 0, ~d, ~d, ~d, ~d};~%" 904 (1+ i-len) (1+ j-len) k-len tot-bits e1comps e2comps tensors) 905 (format tmpout "FINT ng[] = {~d, ~d, ~d, 0, ~d, ~d, ~d, ~d};~%" 906 i-len (+ op-len j-len) k-len tot-bits e1comps e2comps tensors)))) 907 (envs-common (with-output-to-string (tmpout) 908 (format tmpout ngdef) 909 (format tmpout "CINTEnvVars envs;~%") 910 (format tmpout "CINTinit_int3c2e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);~%") 911 (format tmpout "envs.f_gout = &CINTgout2e_~a;~%" intname) 912 (unless (eql (factor-of raw-infix) 1) 913 (format tmpout "envs.common_factor *= ~a;~%" (factor-of raw-infix)))))) 914 (write-functype-header 915 t intname (format nil "/* (~{~a ~}i ~{~a ~}j|~{~a ~}|~{~a ~}k) */" 916 bra-i ket-j op bra-k)) 917 (format fout "/* (~{~a ~}i ~{~a ~}j|~{~a ~}|~{~a ~}k) */~%" 918 bra-i ket-j op bra-k) 919 (gen-code-gout3c2e fout intname raw-infix flat-script) 920 (format fout "void ~a_optimizer(CINTOpt **opt, FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env) {~%" intname) 921 (format fout ngdef) 922 (format fout "CINTall_3c2e_optimizer(opt, ng, atm, natm, bas, nbas, env);~%}~%") 923;;; _cart 924 (format fout "CACHE_SIZE_T ~a_cart(double *out, FINT *dims, FINT *shls, 925FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {~%" intname) 926 (format fout envs-common) 927 (when (member 'g raw-infix) 928 (when (or (member 'g bra-i) (member 'g ket-j)) 929 (format fout "if (out != NULL && envs.shls[0] == envs.shls[1]) { 930FINT i, nout; 931FINT counts[4]; 932counts[0] = envs.nfi * envs.x_ctr[0]; 933counts[1] = envs.nfj * envs.x_ctr[1]; 934counts[2] = envs.nfk * envs.x_ctr[2]; 935counts[3] = 1; 936if (dims == NULL) { dims = counts; } 937nout = dims[0] * dims[1] * dims[2]; 938for (i = 0; i < envs.ncomp_e1 * envs.ncomp_e2 * envs.ncomp_tensor; i++) { 939c2s_dset0(out+nout*i, dims, counts); } 940return 0; }~%"))) 941 (format fout "return CINT3c2e_drv(out, dims, &envs, opt, cache, &c2s_cart_3c2e1, 0);~%} // ~a_cart~%" intname) 942;;; _sph 943 (format fout "CACHE_SIZE_T ~a_sph(double *out, FINT *dims, FINT *shls, 944FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {~%" intname) 945 (format fout envs-common) 946 (when (member 'g raw-infix) 947 (when (or (member 'g bra-i) (member 'g ket-j)) 948 (format fout "if (out != NULL && envs.shls[0] == envs.shls[1]) { 949FINT i, nout; 950FINT counts[4]; 951counts[0] = (envs.i_l*2+1) * envs.x_ctr[0]; 952counts[1] = (envs.j_l*2+1) * envs.x_ctr[1]; 953counts[2] = (envs.k_l*2+1) * envs.x_ctr[2]; 954counts[3] = 1; 955if (dims == NULL) { dims = counts; } 956nout = dims[0] * dims[1] * dims[2]; 957for (i = 0; i < envs.ncomp_e1 * envs.ncomp_e2 * envs.ncomp_tensor; i++) { 958c2s_dset0(out+nout*i, dims, counts); } 959return 0; }~%"))) 960 (format fout "return CINT3c2e_drv(out, dims, &envs, opt, cache, &c2s_sph_3c2e1, 0); 961} // ~a_sph~%" intname) 962;;; _spinor 963 (format fout "CACHE_SIZE_T ~a_spinor(double complex *out, FINT *dims, FINT *shls, 964FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {~%" intname) 965 (format fout envs-common) 966 (when (member 'g raw-infix) 967 (when (or (member 'g bra-i) (member 'g ket-j)) 968 (format fout "if (out != NULL && envs.shls[0] == envs.shls[1]) { 969FINT i, nout; 970FINT counts[4]; 971counts[0] = CINTcgto_spinor(envs.shls[0], envs.bas); 972counts[1] = CINTcgto_spinor(envs.shls[1], envs.bas); 973counts[2] = (envs.k_l*2+1) * envs.x_ctr[2]; 974counts[3] = 1; 975if (dims == NULL) { dims = counts; } 976nout = dims[0] * dims[1] * dims[2]; 977for (i = 0; i < envs.ncomp_tensor; i++) { 978c2s_zset0(out+nout*i, dims, counts); } 979return 0; }~%"))) 980 (format fout "return CINT3c2e_spinor_drv(out, dims, &envs, opt, cache, ~a, 0); 981} // ~a_spinor~%" (name-c2sor "3c2e1" 'spinor sf1 ts1) intname))) 982;;; int2e -> cint2e 983 (format fout "ALL_CINT(~a)~%" intname) 984 (format fout "ALL_CINT_FORTRAN_(~a)~%" intname)) 985 986(defun gen-code-gout2c2e (fout intname raw-infix flat-script) 987 (destructuring-bind (op bra-i ket-j bra-k ket-l) 988 (split-int-expression raw-infix) 989 (let* ((i-rev (effect-keys bra-i)) 990 (k-rev (effect-keys bra-k)) 991 (op-rev (reverse (effect-keys op))) 992 (op-len (length op-rev)) 993 (i-len (length i-rev)) 994 (k-len (length k-rev)) 995 (tot-bits (+ i-len op-len k-len)) 996 (goutinc (length flat-script))) 997 (format fout "static void CINTgout2e_~a(double *gout, 998double *g, FINT *idx, CINTEnvVars *envs, FINT gout_empty) {~%" intname) 999 (format fout "FINT nf = envs->nf; 1000FINT nrys_roots = envs->nrys_roots; 1001FINT ix, iy, iz, i, n; 1002double *g0 = g;~%") 1003 (loop 1004 for i in (range (num-g-intermediates tot-bits op i-len 0)) do 1005 (format fout "double *g~a = g~a + envs->g_size * 3;~%" (1+ i) i)) 1006 (dump-declare-dri-for-rc fout bra-i "i") 1007 (dump-declare-dri-for-rc fout bra-k "k") 1008 (if (intersection *act-left-right* op) 1009 (let ((fmt-k (mkstr "G2E_~aK(g~a, g~a, envs->i_l+" (1+ i-len) ", 0, envs->k_l+~a, 0);~%")) 1010 (fmt-op "") 1011 (fmt-l "")) 1012 (dump-combo-braket fout fmt-k fmt-op fmt-l k-rev '() '() 0)) 1013 (let ((fmt-k (mkstr "G2E_~aK(g~a, g~a, envs->i_l+" i-len ", 0, envs->k_l+~a, 0);~%")) 1014 (fmt-op "") 1015 (fmt-l "")) 1016 (dump-combo-braket fout fmt-k fmt-op fmt-l k-rev '() '() 0))) 1017 (let ((fmt-i "G2E_~aI(g~a, g~a, envs->i_l+~a, 0, envs->k_l, 0);~%") 1018 (fmt-op (mkstr "G2E_~aJ(g~a, g~a, envs->i_l+~d, 0+~a, envs->k_l, 0); 1019G2E_~aI(g~a, g~a, envs->i_l+~a, 0+~a, envs->k_l, 0); 1020for (ix = 0; ix < envs->g_size * 3; ix++) {g~a[ix] += g~a[ix];}~%")) 1021 (fmt-j (mkstr "G2E_~aJ(g~a, g~a, envs->i_l+~d, 0, envs->k_l, 0);~%"))) 1022 (dump-combo-braket fout fmt-i fmt-op fmt-j i-rev op-rev '() k-len)) 1023 (dump-s-2e fout tot-bits) 1024 (gen-c-block-with-empty fout flat-script) 1025 (format fout "}}~%") 1026 goutinc))) 1027 1028(defun gen-code-int2c2e (fout intname raw-infix) 1029 (destructuring-bind (op bra-i ket-j bra-k ket-l) 1030 (split-int-expression raw-infix) 1031 (let* ((i-rev (effect-keys bra-i)) 1032 (k-rev (effect-keys bra-k)) 1033 (op-rev (reverse (effect-keys op))) 1034 (i-len (length i-rev)) 1035 (k-len (length k-rev)) 1036 (op-len (length op-rev)) 1037 (tot-bits (+ i-len op-len k-len)) 1038 (raw-script (eval-int raw-infix)) 1039 (ts1 (car raw-script)) 1040 (sf1 (cadr raw-script)) 1041 (ts2 (caddr raw-script)) 1042 (sf2 (cadddr raw-script)) 1043 (flat-script (flatten-raw-script (last1 raw-script))) 1044 (goutinc (length flat-script)) 1045 (e1comps (if (eql sf1 'sf) 1 4)) 1046 (e2comps (if (eql sf2 'sf) 1 4)) 1047 (tensors (cond ((and (eql sf1 'sf) (eql sf2 'sf)) goutinc) 1048 ((and (eql sf1 'si) (eql sf2 'si)) (/ goutinc 16)) 1049 (t (/ goutinc 4)))) 1050 (ngdef (with-output-to-string (tmpout) 1051 (if (intersection *act-left-right* op) 1052 (format tmpout "FINT ng[] = {~d, 1, ~d, 0, ~d, ~d, ~d, ~d};~%" 1053 (1+ i-len) k-len tot-bits e1comps e2comps tensors) 1054 (format tmpout "FINT ng[] = {~d, ~d, ~d, 0, ~d, ~d, ~d, ~d};~%" 1055 i-len op-len k-len tot-bits e1comps e2comps tensors)))) 1056 (envs-common (with-output-to-string (tmpout) 1057 (format tmpout ngdef) 1058 (format tmpout "CINTEnvVars envs;~%") 1059 (format tmpout "CINTinit_int2c2e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);~%") 1060 (format tmpout "envs.f_gout = &CINTgout2e_~a;~%" intname) 1061 (unless (eql (factor-of raw-infix) 1) 1062 (format tmpout "envs.common_factor *= ~a;~%" (factor-of raw-infix)))))) 1063 (write-functype-header 1064 t intname (format nil "/* (~{~a ~}i |~{~a ~}|~{~a ~}j) */" 1065 bra-i op bra-k)) 1066 (format fout "/* (~{~a ~}i |~{~a ~}|~{~a ~}j) */~%" 1067 bra-i op bra-k) 1068 (gen-code-gout2c2e fout intname raw-infix flat-script) 1069 (format fout "void ~a_optimizer(CINTOpt **opt, FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env) {~%" intname) 1070 (format fout ngdef) 1071 (format fout "CINTall_2c2e_optimizer(opt, ng, atm, natm, bas, nbas, env);~%}~%") 1072;;; _cart 1073 (format fout "CACHE_SIZE_T ~a_cart(double *out, FINT *dims, FINT *shls, 1074FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {~%" intname) 1075 (format fout envs-common) 1076 (format fout "return CINT2c2e_drv(out, dims, &envs, opt, cache, &c2s_cart_1e);~%} // ~a_cart~%" intname) 1077;;; _sph 1078 (format fout "CACHE_SIZE_T ~a_sph(double *out, FINT *dims, FINT *shls, 1079FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {~%" intname) 1080 (format fout envs-common) 1081 (format fout "return CINT2c2e_drv(out, dims, &envs, opt, cache, &c2s_sph_1e);~%} // ~a_sph~%" intname) 1082;;; _spinor 1083 (format fout "CACHE_SIZE_T ~a_spinor(double complex *out, FINT *dims, FINT *shls, 1084FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {~%" intname) 1085; (format fout envs-common) 1086; ; c2s_ function is incorrect if e1/e2 are spin-included operators 1087; (format fout "return CINT2c2e_spinor_drv(out, dims, &envs, opt, cache, ~a); 1088;} // ~a_spinor~%" (name-c2sor "1e" 'spinor sf1 ts1) intname))) 1089 (format fout "fprintf(stderr, \"~a_spinor not implemented\n\"); 1090return 0; 1091}~%" (name-c2sor "1e" 'spinor sf1 ts1) intname))) 1092;;; int2e -> cint2e 1093 (format fout "ALL_CINT(~a)~%" intname) 1094 (format fout "ALL_CINT_FORTRAN_(~a)~%" intname)) 1095 1096 1097(defun gen-code-gout3c1e (fout intname raw-infix flat-script) 1098 (destructuring-bind (op bra-i ket-j bra-k ket-l) 1099 (split-int-expression raw-infix) 1100 (let* ((i-rev (effect-keys bra-i)) 1101 (j-rev (reverse (effect-keys ket-j))) 1102 (k-rev (effect-keys bra-k)) 1103 (op-rev (reverse (effect-keys op))) 1104 (op-len (length op-rev)) 1105 (i-len (length i-rev)) 1106 (j-len (length j-rev)) 1107 (k-len (length k-rev)) 1108 (tot-bits (+ i-len j-len op-len k-len)) 1109 (goutinc (length flat-script))) 1110 (format fout "static void CINTgout1e_~a(double *gout, 1111double *g, FINT *idx, CINTEnvVars *envs, FINT gout_empty) {~%" intname) 1112 (format fout "FINT nf = envs->nf; 1113FINT ix, iy, iz, n; 1114double *g0 = g;~%") 1115 (loop 1116 for i in (range (num-g-intermediates tot-bits op i-len j-len)) do 1117 (format fout "double *g~a = g~a + envs->g_size * 3;~%" (1+ i) i)) 1118 (dump-declare-dri-for-rc fout bra-i "i") 1119 (dump-declare-dri-for-rc fout (append op ket-j) "j") 1120 (dump-declare-dri-for-rc fout bra-k "k") 1121 (dump-declare-giao-ij fout bra-i (append op ket-j)) 1122 (format fout "double s[~a];~%" (expt 3 tot-bits)) 1123 (let ((fmt-k (mkstr "G1E_~aK(g~a, g~a, envs->i_l+" i-len ", envs->j_l+" (+ op-len j-len) 1124 ", envs->k_l+~a);~%")) 1125 (fmt-op "") 1126 (fmt-l "")) 1127 (dump-combo-braket fout fmt-k fmt-op fmt-l k-rev '() '() 0)) 1128 (let ((fmt-i "G1E_~aI(g~a, g~a, envs->i_l+~a, envs->j_l, envs->k_l);~%") 1129 (fmt-op (mkstr "G1E_~aJ(g~a, g~a, envs->i_l+~d, envs->j_l+~a, envs->k_l); 1130G1E_~aI(g~a, g~a, envs->i_l+~a, envs->j_l+~a, envs->k_l); 1131for (ix = 0; ix < envs->g_size * 3; ix++) {g~a[ix] += g~a[ix];}~%")) 1132 (fmt-j (mkstr "G1E_~aJ(g~a, g~a, envs->i_l+~d, envs->j_l+~a, envs->k_l);~%"))) 1133 (dump-combo-braket fout fmt-i fmt-op fmt-j i-rev op-rev j-rev k-len)) 1134 (format fout "for (n = 0; n < nf; n++) { 1135ix = idx[0+n*3]; 1136iy = idx[1+n*3]; 1137iz = idx[2+n*3];~%") 1138 (dump-s-for-nroots fout tot-bits 1) 1139 (gen-c-block-with-empty fout flat-script) 1140 (format fout "}}~%") 1141 goutinc))) 1142(defun gen-code-gout3c1e-nuc (fout intname raw-infix flat-script) 1143 (gen-code-gout3c1e fout intname raw-infix flat-script)) 1144(defun gen-code-gout3c1e-rinv (fout intname raw-infix flat-script) 1145 (gen-code-gout3c1e fout intname raw-infix flat-script)) 1146 1147(defun gen-code-int3c1e (fout intname raw-infix) 1148 (destructuring-bind (op bra-i ket-j bra-k ket-l) 1149 (split-int-expression raw-infix) 1150 (let* ((i-rev (effect-keys bra-i)) 1151 (j-rev (reverse (effect-keys ket-j))) 1152 (k-rev (effect-keys bra-k)) 1153 (op-rev (reverse (effect-keys op))) 1154 (i-len (length i-rev)) 1155 (j-len (length j-rev)) 1156 (k-len (length k-rev)) 1157 (op-len (length op-rev)) 1158 (tot-bits (+ i-len j-len op-len k-len)) 1159 (raw-script (eval-int raw-infix)) 1160 (ts (car raw-script)) 1161 (sf (cadr raw-script)) 1162 (flat-script (flatten-raw-script (last1 raw-script))) 1163 (goutinc (length flat-script)) 1164 (e1comps (if (eql sf 'sf) 1 4)) 1165 (tensors (if (eql sf 'sf) goutinc (/ goutinc 4))) 1166 (int1e-type (cond ((member 'nuc raw-infix) 2) 1167 ((member 'rinv raw-infix) 1) 1168 ((member 'nabla-rinv raw-infix) 1169 (error "nabla-rinv for 3c1e not support")) 1170 (t 0))) 1171 (ngdef (with-output-to-string (tmpout) 1172 (if (or (member 'nuc raw-infix) 1173 (member 'rinv raw-infix) 1174 (member 'nabla-rinv raw-infix)) 1175 (format tmpout "FINT ng[] = {~d, ~d, ~d, 0, ~d, ~d, 0, ~d};~%" 1176 i-len (+ op-len j-len) k-len tot-bits e1comps tensors) 1177 (format tmpout "FINT ng[] = {~d, ~d, ~d, 0, ~d, ~d, 1, ~d};~%" 1178 i-len (+ op-len j-len) k-len tot-bits e1comps tensors)))) 1179 (envs-common (with-output-to-string (tmpout) 1180 (format tmpout ngdef) 1181 (format tmpout "CINTEnvVars envs;~%") 1182 (format tmpout "CINTinit_int3c1e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);~%") 1183 (format tmpout "envs.f_gout = &CINTgout1e_~a;~%" intname) 1184 (unless (eql (factor-of raw-infix) 1) 1185 (format tmpout "envs.common_factor *= ~a;~%" (factor-of raw-infix)))))) 1186 (cond ((member 'nuc raw-infix) 1187 (gen-code-gout3c1e-nuc fout intname raw-infix flat-script)) 1188 ((or (member 'rinv raw-infix) 1189 (member 'nabla-rinv raw-infix)) 1190 (gen-code-gout3c1e-rinv fout intname raw-infix flat-script)) 1191 (t (gen-code-gout3c1e fout intname raw-infix flat-script))) 1192 (write-functype-header 1193 t intname (format nil "/* 3-center 1-electron integral <(~{~a ~}i) (~{~a ~}j) (~{~a ~}k)> */" 1194 bra-i ket-j bra-k)) 1195 (format fout "void ~a_optimizer(CINTOpt **opt, FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env) {~%" intname) 1196 (format fout ngdef) 1197 (format fout "CINTall_3c1e_optimizer(opt, ng, atm, natm, bas, nbas, env);~%}~%") 1198;;; _cart 1199 (format fout "CACHE_SIZE_T ~a_cart(double *out, FINT *dims, FINT *shls, 1200FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {~%" intname) 1201 (format fout envs-common) 1202 (when (member 'g raw-infix) 1203 (when (or (member 'g bra-i) (member 'g ket-j)) 1204 (format fout "if (out != NULL && envs.shls[0] == envs.shls[1]) { 1205FINT i, nout; 1206FINT counts[4]; 1207counts[0] = envs.nfi * envs.x_ctr[0]; 1208counts[1] = envs.nfj * envs.x_ctr[1]; 1209counts[2] = envs.nfk * envs.x_ctr[2]; 1210counts[3] = 1; 1211if (dims == NULL) { dims = counts; } 1212nout = dims[0] * dims[1] * dims[2]; 1213for (i = 0; i < envs.ncomp_e1 * envs.ncomp_e2 * envs.ncomp_tensor; i++) { 1214c2s_dset0(out+nout*i, dims, counts); } 1215return 0; }~%"))) 1216 (format fout "return CINT3c1e_drv(out, dims, &envs, opt, cache, &c2s_cart_3c1e, ~d, 0);~%} 1217// ~a_cart~%" int1e-type intname) 1218;;; _sph 1219 (format fout "CACHE_SIZE_T ~a_sph(double *out, FINT *dims, FINT *shls, 1220FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {~%" intname) 1221 (format fout envs-common) 1222 (when (member 'g raw-infix) 1223 (when (or (member 'g bra-i) (member 'g ket-j)) 1224 (format fout "if (out != NULL && envs.shls[0] == envs.shls[1]) { 1225FINT i, nout; 1226FINT counts[4]; 1227counts[0] = (envs.i_l*2+1) * envs.x_ctr[0]; 1228counts[1] = (envs.j_l*2+1) * envs.x_ctr[1]; 1229counts[2] = (envs.k_l*2+1) * envs.x_ctr[2]; 1230counts[3] = 1; 1231if (dims == NULL) { dims = counts; } 1232nout = dims[0] * dims[1] * dims[2]; 1233for (i = 0; i < envs.ncomp_e1 * envs.ncomp_e2 * envs.ncomp_tensor; i++) { 1234c2s_dset0(out+nout*i, dims, counts); } 1235return 0; }~%"))) 1236 (format fout "return CINT3c1e_drv(out, dims, &envs, opt, cache, &c2s_sph_3c1e, ~d, 0); 1237} // ~a_sph~%" int1e-type intname) 1238;;; _spinor 1239 (format fout "CACHE_SIZE_T ~a_spinor(double complex *out, FINT *dims, FINT *shls, 1240FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {~%" intname) 1241 (format fout envs-common) 1242 (when (member 'g raw-infix) 1243 (when (or (member 'g bra-i) (member 'g ket-j)) 1244 (format fout "if (out != NULL && envs.shls[0] == envs.shls[1]) { 1245FINT i, nout; 1246FINT counts[4]; 1247counts[0] = CINTcgto_spinor(envs.shls[0], envs.bas); 1248counts[1] = CINTcgto_spinor(envs.shls[1], envs.bas); 1249counts[2] = (envs.k_l*2+1) * envs.x_ctr[2]; 1250counts[3] = 1; 1251if (dims == NULL) { dims = counts; } 1252nout = dims[0] * dims[1] * dims[2]; 1253for (i = 0; i < envs.ncomp_tensor; i++) { 1254c2s_zset0(out+nout*i, dims, counts); } 1255return 0; }~%"))) 1256 (format fout "return CINT3c1e_spinor_drv(out, dims, &envs, opt, cache, ~a, ~d, 0); 1257} // ~a_spinor~%" (name-c2sor "3c2e1" 'spinor sf ts) int1e-type intname))) 1258 (format fout "ALL_CINT(~a)~%" intname) 1259 (format fout "ALL_CINT_FORTRAN_(~a)~%" intname)) 1260 1261;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1262 1263(defun gen-code-gout1e-grids (fout intname raw-infix flat-script) 1264 (destructuring-bind (op bra-i ket-j bra-k ket-l) 1265 (split-int-expression raw-infix) 1266 (let* ((i-rev (effect-keys bra-i)) ;<i| already in reverse order 1267 (j-rev (reverse (effect-keys ket-j))) 1268 (op-rev (reverse (effect-keys op))) 1269 (i-len (length i-rev)) 1270 (j-len (length j-rev)) 1271 (op-len (length op-rev)) 1272 (tot-bits (+ i-len j-len op-len)) 1273 (comp (expt 3 tot-bits)) 1274 (goutinc (length flat-script))) 1275 (format fout "void CINTgout1e_~a" intname) 1276 (format fout "(double *gout, double *g, FINT *idx, CINTEnvVars *envs, FINT gout_empty) { 1277FINT ngrids = envs->ngrids; 1278FINT bgrids = MIN(ngrids - envs->grids_offset, GRID_BLKSIZE); 1279FINT nrys_roots = envs->nrys_roots; 1280FINT nf = envs->nf; 1281FINT ix, iy, iz, n, i, ig; 1282double *g0 = g;~%") 1283 (loop 1284 for i in (range (num-g-intermediates tot-bits op i-len j-len)) do 1285 (format fout "double *g~a = g~a + envs->g_size * 3;~%" (1+ i) i)) 1286 (dump-declare-dri-for-rc fout bra-i "i") 1287 (dump-declare-dri-for-rc fout (append op ket-j) "j") 1288 (dump-declare-giao-ij fout bra-i (append op ket-j)) 1289 (format fout "double s[GRID_BLKSIZE * ~a];~%" comp) 1290;;; generate g_(bin) 1291;;; for the operators act on the |ket>, the reversed scan order and r_combinator 1292;;; is required; for the operators acto on the <bra|, the normal scan order and 1293 (let ((fmt-i "G1E_GRIDS_~aI(g~a, g~a, envs->i_l+~a, envs->j_l);~%") 1294 (fmt-op (mkstr "G1E_GRIDS_~aJ(g~a, g~a, envs->i_l+~d, envs->j_l+~a); 1295G1E_~aI(g~a, g~a, envs->i_l+~d, envs->j_l+~a, 0); 1296for (ix = 0; ix < envs->g_size * 3; ix++) {g~a[ix] += g~a[ix];}~%")) 1297 (fmt-j (mkstr "G1E_GRIDS_~aJ(g~a, g~a, envs->i_l+~d, envs->j_l+~a);~%"))) 1298 (dump-combo-braket fout fmt-i fmt-op fmt-j i-rev op-rev j-rev 0)) 1299 (format fout "for (n = 0; n < nf; n++) { 1300ix = idx[0+n*3]; 1301iy = idx[1+n*3]; 1302iz = idx[2+n*3];~%") 1303 (format fout "for (i = 0; i < ~a; i++) { 1304for (ig = 0; ig < bgrids; ig++) { s[ig+i*GRID_BLKSIZE] = 0; }} 1305for (i = 0; i < nrys_roots; i++) { 1306for (ig = 0; ig < bgrids; ig++) {~%" comp) 1307 (dump-s-loop fout tot-bits t) 1308 (format fout "}};~%") 1309 1310 ; (gen-c-block-with-empty fout flat-script) 1311 (format fout "if (gout_empty) { 1312for (ig = 0; ig < bgrids; ig++) {~%") 1313 (gen-c-block fout flat-script t) 1314 (format fout "}} else { 1315for (ig = 0; ig < bgrids; ig++) {~%") 1316 (gen-c-block+ fout flat-script t) 1317 (format fout "}}}}~%") 1318 goutinc))) 1319 1320(defun gen-code-int1e-grids (fout intname raw-infix) 1321 (destructuring-bind (op bra-i ket-j bra-k ket-l) 1322 (split-int-expression raw-infix) 1323 (let* ((i-rev (effect-keys bra-i)) ;<i| already in reverse order 1324 (j-rev (reverse (effect-keys ket-j))) 1325 (op-rev (reverse (effect-keys op))) 1326 (i-len (length i-rev)) 1327 (j-len (length j-rev)) 1328 (op-len (length op-rev)) 1329 (tot-bits (+ i-len j-len op-len)) 1330 (raw-script (eval-int raw-infix)) 1331 (flat-script (flatten-raw-script (last1 raw-script))) 1332 (ts (car raw-script)) 1333 (sf (cadr raw-script)) 1334 (goutinc (length flat-script)) 1335 (e1comps (if (eql sf 'sf) 1 4)) 1336 (tensors (if (eql sf 'sf) goutinc (/ goutinc 4))) 1337 (ngdef (with-output-to-string (tmpout) 1338 (if (intersection *act-left-right* op) 1339 (format tmpout "FINT ng[] = {~d, ~d, 0, 0, ~d, ~d, 0, ~d};~%" 1340 (1+ i-len) (+ op-len j-len) tot-bits e1comps tensors) 1341 (format tmpout "FINT ng[] = {~d, ~d, 0, 0, ~d, ~d, 0, ~d};~%" 1342 i-len (+ op-len j-len) tot-bits e1comps tensors)))) 1343 (envs-common (with-output-to-string (tmpout) 1344 (format tmpout ngdef) 1345 (format tmpout "CINTEnvVars envs;~%") 1346 (format tmpout "CINTinit_int1e_grids_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);~%") 1347 (format tmpout "envs.f_gout = &CINTgout1e_~a;~%" intname) 1348 (unless (eql (factor-of raw-infix) 1) 1349 (format tmpout "envs.common_factor *= ~a;~%" (factor-of raw-infix)))))) 1350 (write-functype-header 1351 t intname (format nil "/* <~{~a ~}i| 1/r_{grids} |~{~a ~}j> */" bra-i ket-j)) 1352 (format fout "/* <~{~a ~}i| 1/r_{grids} |~{~a ~}j> */~%" bra-i ket-j) 1353 (gen-code-gout1e-grids fout intname raw-infix flat-script) 1354 (format fout "void ~a_optimizer(CINTOpt **opt, FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env) {~%" intname) 1355 (format fout ngdef) 1356 (format fout "CINTall_1e_grids_optimizer(opt, ng, atm, natm, bas, nbas, env);~%}~%") 1357;;; _cart 1358 (format fout "CACHE_SIZE_T ~a_cart(double *out, FINT *dims, FINT *shls, 1359FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {~%" intname) 1360 (format fout envs-common) 1361 (when (member 'g raw-infix) 1362 (when (or (member 'g bra-i) (member 'g ket-j)) 1363 (format fout "if (out != NULL && envs.shls[0] == envs.shls[1]) { 1364FINT i, nout; 1365FINT counts[4]; 1366counts[0] = envs.ngrids; 1367counts[1] = envs.nfi * envs.x_ctr[0]; 1368counts[2] = envs.nfj * envs.x_ctr[1]; 1369counts[3] = 1; 1370if (dims == NULL) { dims = counts; } 1371nout = dims[0] * dims[1] * dims[2]; 1372for (i = 0; i < envs.ncomp_e1 * envs.ncomp_tensor; i++) { 1373c2s_dset0(out+nout*i, dims, counts); } 1374return 0; }~%"))) 1375 (format fout "return CINT1e_grids_drv(out, dims, &envs, cache, &c2s_cart_1e_grids); 1376} // ~a_cart~%" intname) 1377;;; _sph 1378 (format fout "CACHE_SIZE_T ~a_sph(double *out, FINT *dims, FINT *shls, 1379FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {~%" intname) 1380 (format fout envs-common) 1381 (when (member 'g raw-infix) 1382 (when (or (member 'g bra-i) (member 'g ket-j)) 1383 (format fout "if (out != NULL && envs.shls[0] == envs.shls[1]) { 1384FINT i, nout; 1385FINT counts[4]; 1386counts[0] = envs.ngrids; 1387counts[1] = (envs.i_l*2+1) * envs.x_ctr[0]; 1388counts[2] = (envs.j_l*2+1) * envs.x_ctr[1]; 1389counts[3] = 1; 1390if (dims == NULL) { dims = counts; } 1391nout = dims[0] * dims[1] * dims[2]; 1392for (i = 0; i < envs.ncomp_e1 * envs.ncomp_tensor; i++) { 1393c2s_dset0(out+nout*i, dims, counts); } 1394return 0; }~%"))) 1395 (format fout "return CINT1e_grids_drv(out, dims, &envs, cache, &c2s_sph_1e_grids); 1396} // ~a_sph~%" intname) 1397;;; _spinor 1398 (format fout "CACHE_SIZE_T ~a_spinor(double complex *out, FINT *dims, FINT *shls, 1399FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {~%" intname) 1400 (format fout envs-common) 1401 (when (member 'g raw-infix) 1402 (when (or (member 'g bra-i) (member 'g ket-j)) 1403 (format fout "if (out != NULL && envs.shls[0] == envs.shls[1]) { 1404FINT i, nout; 1405FINT counts[4]; 1406counts[0] = envs.ngrids; 1407counts[1] = CINTcgto_spinor(envs.shls[0], envs.bas); 1408counts[2] = CINTcgto_spinor(envs.shls[1], envs.bas); 1409counts[3] = 1; 1410if (dims == NULL) { dims = counts; } 1411nout = dims[0] * dims[1] * dims[2]; 1412for (i = 0; i < envs.ncomp_tensor; i++) { 1413c2s_zset0(out+nout*i, dims, counts); } 1414return 0; }~%"))) 1415 (format fout "return CINT1e_grids_spinor_drv(out, dims, &envs, cache, ~a); 1416} // ~a_spinor~%" (name-c2sor "1e_grids" 'spinor sf ts) intname))) 1417;;; int1e -> cint1e 1418 (format fout "ALL_CINT1E(~a)~%" intname) 1419 (format fout "ALL_CINT1E_FORTRAN_(~a)~%" intname)) 1420 1421;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1422 1423(defun gen-cint (filename &rest items) 1424 "sp can be one of 'spinor 'spheric 'cart" 1425 (with-open-file (fout (mkstr filename) 1426 :direction :output :if-exists :supersede) 1427 (dump-header fout) 1428 (flet ((gen-code (item) 1429 (let ((intname (mkstr (car item))) 1430 (raw-infix (cadr item))) 1431 (cond ((int3c1e? raw-infix) 1432 (gen-code-int3c1e fout intname raw-infix)) 1433 ((int1e-grids? raw-infix) 1434 (gen-code-int1e-grids fout intname raw-infix)) 1435 ((one-electron-int? raw-infix) 1436 (gen-code-int1e fout intname raw-infix)) 1437 ((int4c2e? raw-infix) 1438 (gen-code-int4c2e fout intname raw-infix)) 1439 ((int3c2e? raw-infix) 1440 (gen-code-int3c2e fout intname raw-infix)) 1441 ((int2c2e? raw-infix) 1442 (gen-code-int2c2e fout intname raw-infix)))))) 1443 (mapcar #'gen-code items)))) 1444 1445; gcl -load sigma.o -batch -eval "( .. )" 1446 1447;; vim: ft=lisp 1448