| 1 | //
 | 
|---|
| 2 | // tformv3.cc
 | 
|---|
| 3 | //
 | 
|---|
| 4 | // Copyright (C) 2001 Edward Valeev
 | 
|---|
| 5 | //
 | 
|---|
| 6 | // Author: Edward Valeev <edward.valeev@chemistry.gatech.edu>
 | 
|---|
| 7 | // Maintainer: EV
 | 
|---|
| 8 | //
 | 
|---|
| 9 | // This file is part of the SC Toolkit.
 | 
|---|
| 10 | //
 | 
|---|
| 11 | // The SC Toolkit is free software; you can redistribute it and/or modify
 | 
|---|
| 12 | // it under the terms of the GNU Library General Public License as published by
 | 
|---|
| 13 | // the Free Software Foundation; either version 2, or (at your option)
 | 
|---|
| 14 | // any later version.
 | 
|---|
| 15 | //
 | 
|---|
| 16 | // The SC Toolkit is distributed in the hope that it will be useful,
 | 
|---|
| 17 | // but WITHOUT ANY WARRANTY; without even the implied warranty of
 | 
|---|
| 18 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 | 
|---|
| 19 | // GNU Library General Public License for more details.
 | 
|---|
| 20 | //
 | 
|---|
| 21 | // You should have received a copy of the GNU Library General Public License
 | 
|---|
| 22 | // along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
 | 
|---|
| 23 | // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
 | 
|---|
| 24 | //
 | 
|---|
| 25 | // The U.S. Government is granted a limited license as per AL 91-7.
 | 
|---|
| 26 | //
 | 
|---|
| 27 | 
 | 
|---|
| 28 | #if defined(__GNUC__)
 | 
|---|
| 29 | #pragma implementation
 | 
|---|
| 30 | #endif
 | 
|---|
| 31 | 
 | 
|---|
| 32 | #include <stdlib.h>
 | 
|---|
| 33 | #include <string.h>
 | 
|---|
| 34 | #include <math.h>
 | 
|---|
| 35 | 
 | 
|---|
| 36 | #include <util/misc/formio.h>
 | 
|---|
| 37 | #include <chemistry/qc/basis/integral.h>
 | 
|---|
| 38 | #include <chemistry/qc/cints/macros.h>
 | 
|---|
| 39 | #include <chemistry/qc/cints/tform.h>
 | 
|---|
| 40 | #include <chemistry/qc/cints/int1e.h>
 | 
|---|
| 41 | #include <chemistry/qc/cints/int2e.h>
 | 
|---|
| 42 | 
 | 
|---|
| 43 | using namespace std;
 | 
|---|
| 44 | using namespace sc;
 | 
|---|
| 45 | 
 | 
|---|
| 46 | inline int max(int a,int b) { return (a > b) ? a : b;}
 | 
|---|
| 47 | inline void fail()
 | 
|---|
| 48 | {
 | 
|---|
| 49 |   ExEnv::errn() << scprintf("failing module:\n%s",__FILE__) << endl;
 | 
|---|
| 50 |   abort();
 | 
|---|
| 51 | }
 | 
|---|
| 52 | 
 | 
|---|
| 53 | static void transform1e_1(SphericalTransformIter&, double*, double*, int);
 | 
|---|
| 54 | static void transform1e_2(SphericalTransformIter&, double*, double*, int, int);
 | 
|---|
| 55 | static void transform1e_vec_2(const int, SphericalTransformIter&, double*, double*, int, int);
 | 
|---|
| 56 | static void transform2e_1(SphericalTransformIter&, double*, double*, int);
 | 
|---|
| 57 | static void transform2e_2(SphericalTransformIter&, double*, double*, int, int, int);
 | 
|---|
| 58 | static void transform2e_3(SphericalTransformIter&, double*, double*, int, int, int);
 | 
|---|
| 59 | static void transform2e_4(SphericalTransformIter&, double*, double*, int, int);
 | 
|---|
| 60 | 
 | 
|---|
| 61 | void Int1eCints::transform_contrquartets_(double * source_ints_buf, double *target_ints_buf)
 | 
|---|
| 62 | {
 | 
|---|
| 63 |   double *source1, *target1;
 | 
|---|
| 64 |   double *source2, *target2;
 | 
|---|
| 65 |   double *source = source_ints_buf;
 | 
|---|
| 66 |   double *target = target_ints_buf;
 | 
|---|
| 67 |   double *tmpbuf = tformbuf_;
 | 
|---|
| 68 | 
 | 
|---|
| 69 |   for (int gc1=0; gc1<int_shell1_->ncontraction(); gc1++) {
 | 
|---|
| 70 |     int am1 = int_shell1_->am(gc1);
 | 
|---|
| 71 |     int is_pure1 = int_shell1_->is_pure(gc1) ? 1 : 0;
 | 
|---|
| 72 |     int ncart1 = int_shell1_->ncartesian(gc1);
 | 
|---|
| 73 |     int nbf1 = int_shell1_->nfunction(gc1);
 | 
|---|
| 74 | 
 | 
|---|
| 75 |     int target_bf2_offset = 0;
 | 
|---|
| 76 |     for (int gc2=0; gc2<int_shell2_->ncontraction(); gc2++) {
 | 
|---|
| 77 |       int am2 = int_shell2_->am(gc2);
 | 
|---|
| 78 |       int is_pure2 = int_shell2_->is_pure(gc2) ? 1 : 0;
 | 
|---|
| 79 |       int ncart2 = int_shell2_->ncartesian(gc2);
 | 
|---|
| 80 |       int nbf2 = int_shell2_->nfunction(gc2);
 | 
|---|
| 81 |       
 | 
|---|
| 82 |       int transform_index = 2*is_pure1 + is_pure2;
 | 
|---|
| 83 |       switch (transform_index) {
 | 
|---|
| 84 |       case 0:
 | 
|---|
| 85 |         break;
 | 
|---|
| 86 | 
 | 
|---|
| 87 |       case 1:
 | 
|---|
| 88 |         source2 = source;
 | 
|---|
| 89 |         target2 = target;
 | 
|---|
| 90 |         break;
 | 
|---|
| 91 |         
 | 
|---|
| 92 |       case 2:
 | 
|---|
| 93 |         source1 = source;
 | 
|---|
| 94 |         target1 = target;
 | 
|---|
| 95 |         break;
 | 
|---|
| 96 | 
 | 
|---|
| 97 |       case 3:
 | 
|---|
| 98 |         source2 = source;
 | 
|---|
| 99 |         target2 = tmpbuf;
 | 
|---|
| 100 |         source1 = tmpbuf;
 | 
|---|
| 101 |         target1 = target;
 | 
|---|
| 102 |         break;
 | 
|---|
| 103 |       }
 | 
|---|
| 104 | 
 | 
|---|
| 105 |       if (is_pure2) {
 | 
|---|
| 106 |         SphericalTransformIter stiter(integral_->spherical_transform(am2));
 | 
|---|
| 107 |         transform1e_2(stiter,source2, target2, ncart1,ncart2);
 | 
|---|
| 108 |       }
 | 
|---|
| 109 |       if (is_pure1) {
 | 
|---|
| 110 |         SphericalTransformIter stiter(integral_->spherical_transform(am1));
 | 
|---|
| 111 |         transform1e_1(stiter,source1, target1, nbf2);
 | 
|---|
| 112 |       }
 | 
|---|
| 113 |       
 | 
|---|
| 114 |       source += (ncart1*ncart2);
 | 
|---|
| 115 |       target += (nbf1*nbf2);
 | 
|---|
| 116 |     }
 | 
|---|
| 117 |   }
 | 
|---|
| 118 | }
 | 
|---|
| 119 | 
 | 
|---|
| 120 | void Int1eCints::transform_contrquartets_vec_(const int ntypes, double * source_ints_buf, double *target_ints_buf)
 | 
|---|
| 121 | {
 | 
|---|
| 122 |   double *source1, *target1;
 | 
|---|
| 123 |   double *source2, *target2;
 | 
|---|
| 124 |   double *source = source_ints_buf;
 | 
|---|
| 125 |   double *target = target_ints_buf;
 | 
|---|
| 126 |   double *tmpbuf = tformbuf_;
 | 
|---|
| 127 | 
 | 
|---|
| 128 |   for (int gc1=0; gc1<int_shell1_->ncontraction(); gc1++) {
 | 
|---|
| 129 |     int am1 = int_shell1_->am(gc1);
 | 
|---|
| 130 |     int is_pure1 = int_shell1_->is_pure(gc1) ? 1 : 0;
 | 
|---|
| 131 |     int ncart1 = int_shell1_->ncartesian(gc1);
 | 
|---|
| 132 |     int nbf1 = int_shell1_->nfunction(gc1);
 | 
|---|
| 133 | 
 | 
|---|
| 134 |     int target_bf2_offset = 0;
 | 
|---|
| 135 |     for (int gc2=0; gc2<int_shell2_->ncontraction(); gc2++) {
 | 
|---|
| 136 |       int am2 = int_shell2_->am(gc2);
 | 
|---|
| 137 |       int is_pure2 = int_shell2_->is_pure(gc2) ? 1 : 0;
 | 
|---|
| 138 |       int ncart2 = int_shell2_->ncartesian(gc2);
 | 
|---|
| 139 |       int nbf2 = int_shell2_->nfunction(gc2);
 | 
|---|
| 140 |       
 | 
|---|
| 141 |       int transform_index = 2*is_pure1 + is_pure2;
 | 
|---|
| 142 |       switch (transform_index) {
 | 
|---|
| 143 |       case 0:
 | 
|---|
| 144 |         break;
 | 
|---|
| 145 | 
 | 
|---|
| 146 |       case 1:
 | 
|---|
| 147 |         source2 = source;
 | 
|---|
| 148 |         target2 = target;
 | 
|---|
| 149 |         break;
 | 
|---|
| 150 |         
 | 
|---|
| 151 |       case 2:
 | 
|---|
| 152 |         source1 = source;
 | 
|---|
| 153 |         target1 = target;
 | 
|---|
| 154 |         break;
 | 
|---|
| 155 | 
 | 
|---|
| 156 |       case 3:
 | 
|---|
| 157 |         source2 = source;
 | 
|---|
| 158 |         target2 = tmpbuf;
 | 
|---|
| 159 |         source1 = tmpbuf;
 | 
|---|
| 160 |         target1 = target;
 | 
|---|
| 161 |         break;
 | 
|---|
| 162 |       }
 | 
|---|
| 163 | 
 | 
|---|
| 164 |       if (is_pure2) {
 | 
|---|
| 165 |         SphericalTransformIter stiter(integral_->spherical_transform(am2));
 | 
|---|
| 166 |         transform1e_vec_2(ntypes, stiter, source2, target2, ncart1,ncart2);
 | 
|---|
| 167 |       }
 | 
|---|
| 168 |       if (is_pure1) {
 | 
|---|
| 169 |         SphericalTransformIter stiter(integral_->spherical_transform(am1));
 | 
|---|
| 170 |         transform1e_1(stiter, source1, target1, nbf2*ntypes);
 | 
|---|
| 171 |       }
 | 
|---|
| 172 |       
 | 
|---|
| 173 |       source += (ntypes*ncart1*ncart2);
 | 
|---|
| 174 |       target += (ntypes*nbf1*nbf2);
 | 
|---|
| 175 |     }
 | 
|---|
| 176 |   }
 | 
|---|
| 177 | }
 | 
|---|
| 178 | 
 | 
|---|
| 179 | void Int2eCints::transform_contrquartets_(double * source_ints_buf, double *target_ints_buf)
 | 
|---|
| 180 | {
 | 
|---|
| 181 |   double *source1, *target1;
 | 
|---|
| 182 |   double *source2, *target2;
 | 
|---|
| 183 |   double *source3, *target3;
 | 
|---|
| 184 |   double *source4, *target4;
 | 
|---|
| 185 |   double *source = source_ints_buf;
 | 
|---|
| 186 |   double *target = target_ints_buf;
 | 
|---|
| 187 |   double *tmpbuf = tformbuf_;
 | 
|---|
| 188 | 
 | 
|---|
| 189 |   for (int gc1=0; gc1<int_shell1_->ncontraction(); gc1++) {
 | 
|---|
| 190 |     int am1 = int_shell1_->am(gc1);
 | 
|---|
| 191 |     int is_pure1 = int_shell1_->is_pure(gc1) ? 1 : 0;
 | 
|---|
| 192 |     int ncart1 = int_shell1_->ncartesian(gc1);
 | 
|---|
| 193 |     int nbf1 = int_shell1_->nfunction(gc1);
 | 
|---|
| 194 | 
 | 
|---|
| 195 |     int target_bf2_offset = 0;
 | 
|---|
| 196 |     for (int gc2=0; gc2<int_shell2_->ncontraction(); gc2++) {
 | 
|---|
| 197 |       int am2 = int_shell2_->am(gc2);
 | 
|---|
| 198 |       int is_pure2 = int_shell2_->is_pure(gc2) ? 1 : 0;
 | 
|---|
| 199 |       int ncart2 = int_shell2_->ncartesian(gc2);
 | 
|---|
| 200 |       int nbf2 = int_shell2_->nfunction(gc2);
 | 
|---|
| 201 |       
 | 
|---|
| 202 |       int target_bf3_offset = 0;
 | 
|---|
| 203 |       for (int gc3=0; gc3<int_shell3_->ncontraction(); gc3++) {
 | 
|---|
| 204 |         int am3 = int_shell3_->am(gc3);
 | 
|---|
| 205 |         int is_pure3 = int_shell3_->is_pure(gc3) ? 1 : 0;
 | 
|---|
| 206 |         int ncart3 = int_shell3_->ncartesian(gc3);
 | 
|---|
| 207 |         int nbf3 = int_shell3_->nfunction(gc3);
 | 
|---|
| 208 |         
 | 
|---|
| 209 |         int target_bf4_offset = 0;
 | 
|---|
| 210 |         for (int gc4=0; gc4<int_shell4_->ncontraction(); gc4++) {
 | 
|---|
| 211 |           int am4 = int_shell4_->am(gc4);
 | 
|---|
| 212 |           int is_pure4 = int_shell4_->is_pure(gc4) ? 1 : 0;
 | 
|---|
| 213 |           int ncart4 = int_shell4_->ncartesian(gc4);
 | 
|---|
| 214 |           int nbf4 = int_shell4_->nfunction(gc4);
 | 
|---|
| 215 | 
 | 
|---|
| 216 |           int transform_index = 8*is_pure1 + 4*is_pure2 + 2*is_pure3 + is_pure4;
 | 
|---|
| 217 |           switch (transform_index) {
 | 
|---|
| 218 |             case 0:
 | 
|---|
| 219 |               break;
 | 
|---|
| 220 | 
 | 
|---|
| 221 |             case 1:
 | 
|---|
| 222 |               source4 = source;
 | 
|---|
| 223 |               target4 = target;
 | 
|---|
| 224 |               break;
 | 
|---|
| 225 | 
 | 
|---|
| 226 |             case 2:
 | 
|---|
| 227 |               source3 = source;
 | 
|---|
| 228 |               target3 = target;
 | 
|---|
| 229 |               break;
 | 
|---|
| 230 | 
 | 
|---|
| 231 |             case 3:
 | 
|---|
| 232 |               source4 = source;
 | 
|---|
| 233 |               target4 = tmpbuf;
 | 
|---|
| 234 |               source3 = tmpbuf;
 | 
|---|
| 235 |               target3 = target;
 | 
|---|
| 236 |               break;
 | 
|---|
| 237 |               
 | 
|---|
| 238 |             case 4:
 | 
|---|
| 239 |               source2 = source;
 | 
|---|
| 240 |               target2 = target;
 | 
|---|
| 241 |               break;
 | 
|---|
| 242 | 
 | 
|---|
| 243 |             case 5:
 | 
|---|
| 244 |               source4 = source;
 | 
|---|
| 245 |               target4 = tmpbuf;
 | 
|---|
| 246 |               source2 = tmpbuf;
 | 
|---|
| 247 |               target2 = target;
 | 
|---|
| 248 |               break;
 | 
|---|
| 249 | 
 | 
|---|
| 250 |             case 6:
 | 
|---|
| 251 |               source3 = source;
 | 
|---|
| 252 |               target3 = tmpbuf;
 | 
|---|
| 253 |               source2 = tmpbuf;
 | 
|---|
| 254 |               target2 = target;
 | 
|---|
| 255 |               break;
 | 
|---|
| 256 |             
 | 
|---|
| 257 |             case 7:
 | 
|---|
| 258 |               source4 = source;
 | 
|---|
| 259 |               target4 = tmpbuf;
 | 
|---|
| 260 |               source3 = tmpbuf;
 | 
|---|
| 261 |               target3 = source;
 | 
|---|
| 262 |               source2 = source;
 | 
|---|
| 263 |               target2 = target;
 | 
|---|
| 264 |               break;
 | 
|---|
| 265 | 
 | 
|---|
| 266 |             case 8:
 | 
|---|
| 267 |               source1 = source;
 | 
|---|
| 268 |               target1 = target;
 | 
|---|
| 269 |               break;
 | 
|---|
| 270 |               
 | 
|---|
| 271 |             case 9:
 | 
|---|
| 272 |               source4 = source;
 | 
|---|
| 273 |               target4 = tmpbuf;
 | 
|---|
| 274 |               source1 = tmpbuf;
 | 
|---|
| 275 |               target1 = target;
 | 
|---|
| 276 |               break;
 | 
|---|
| 277 | 
 | 
|---|
| 278 |             case 10:
 | 
|---|
| 279 |               source3 = source;
 | 
|---|
| 280 |               target3 = tmpbuf;
 | 
|---|
| 281 |               source1 = tmpbuf;
 | 
|---|
| 282 |               target1 = target;
 | 
|---|
| 283 |               break;
 | 
|---|
| 284 | 
 | 
|---|
| 285 |             case 11:
 | 
|---|
| 286 |               source4 = source;
 | 
|---|
| 287 |               target4 = tmpbuf;
 | 
|---|
| 288 |               source3 = tmpbuf;
 | 
|---|
| 289 |               target3 = source;
 | 
|---|
| 290 |               source1 = source;
 | 
|---|
| 291 |               target1 = target;
 | 
|---|
| 292 |               break;
 | 
|---|
| 293 |               
 | 
|---|
| 294 |             case 12:
 | 
|---|
| 295 |               source2 = source;
 | 
|---|
| 296 |               target2 = tmpbuf;
 | 
|---|
| 297 |               source1 = tmpbuf;
 | 
|---|
| 298 |               target1 = target;
 | 
|---|
| 299 |               break;
 | 
|---|
| 300 | 
 | 
|---|
| 301 |             case 13:
 | 
|---|
| 302 |               source4 = source;
 | 
|---|
| 303 |               target4 = tmpbuf;
 | 
|---|
| 304 |               source2 = tmpbuf;
 | 
|---|
| 305 |               target2 = source;
 | 
|---|
| 306 |               source1 = source;
 | 
|---|
| 307 |               target1 = target;
 | 
|---|
| 308 |               break;
 | 
|---|
| 309 | 
 | 
|---|
| 310 |             case 14:
 | 
|---|
| 311 |               source3 = source;
 | 
|---|
| 312 |               target3 = tmpbuf;
 | 
|---|
| 313 |               source2 = tmpbuf;
 | 
|---|
| 314 |               target2 = source;
 | 
|---|
| 315 |               source1 = source;
 | 
|---|
| 316 |               target1 = target;
 | 
|---|
| 317 |               break;
 | 
|---|
| 318 |             
 | 
|---|
| 319 |             case 15:
 | 
|---|
| 320 |               source4 = source;
 | 
|---|
| 321 |               target4 = tmpbuf;
 | 
|---|
| 322 |               source3 = tmpbuf;
 | 
|---|
| 323 |               target3 = source;
 | 
|---|
| 324 |               source2 = source;
 | 
|---|
| 325 |               target2 = tmpbuf;
 | 
|---|
| 326 |               source1 = tmpbuf;
 | 
|---|
| 327 |               target1 = target;
 | 
|---|
| 328 |               break;
 | 
|---|
| 329 |           }
 | 
|---|
| 330 | 
 | 
|---|
| 331 |           if (is_pure4) {
 | 
|---|
| 332 |             SphericalTransformIter stiter(integral_->spherical_transform(am4));
 | 
|---|
| 333 |             transform2e_4(stiter, source4, target4, ncart1*ncart2*ncart3,ncart4);
 | 
|---|
| 334 |           }
 | 
|---|
| 335 |           if (is_pure3) {
 | 
|---|
| 336 |             SphericalTransformIter stiter(integral_->spherical_transform(am3));
 | 
|---|
| 337 |             transform2e_3(stiter,source3, target3, ncart1*ncart2,ncart3,nbf4);
 | 
|---|
| 338 |           }
 | 
|---|
| 339 |           if (is_pure2) {
 | 
|---|
| 340 |             SphericalTransformIter stiter(integral_->spherical_transform(am2));
 | 
|---|
| 341 |             transform2e_2(stiter,source2, target2, ncart1,ncart2,nbf3*nbf4);
 | 
|---|
| 342 |           }
 | 
|---|
| 343 |           if (is_pure1) {
 | 
|---|
| 344 |             SphericalTransformIter stiter(integral_->spherical_transform(am1));
 | 
|---|
| 345 |             transform2e_1(stiter,source1, target1, nbf2*nbf3*nbf4);
 | 
|---|
| 346 |           }
 | 
|---|
| 347 |           
 | 
|---|
| 348 |           source += (ncart1*ncart2*ncart3*ncart4);
 | 
|---|
| 349 |           target += (nbf1*nbf2*nbf3*nbf4);
 | 
|---|
| 350 |         }
 | 
|---|
| 351 |       }
 | 
|---|
| 352 |     }
 | 
|---|
| 353 |   }
 | 
|---|
| 354 | }
 | 
|---|
| 355 | 
 | 
|---|
| 356 | 
 | 
|---|
| 357 | static void transform1e_1(SphericalTransformIter& sti, double *s, double *t, int nl)
 | 
|---|
| 358 | {
 | 
|---|
| 359 |   memset(t,0,sti.n()*nl*sizeof(double));
 | 
|---|
| 360 | 
 | 
|---|
| 361 |   for (sti.begin(); sti.ready(); sti.next()) {
 | 
|---|
| 362 |     double *sptr = s + sti.cartindex()*nl;
 | 
|---|
| 363 |     double *tptr = t + sti.pureindex()*nl;
 | 
|---|
| 364 |     double coef = sti.coef();
 | 
|---|
| 365 |     for(int l=0; l<nl; l++)
 | 
|---|
| 366 |       *(tptr++) += coef * *(sptr++);
 | 
|---|
| 367 |   }
 | 
|---|
| 368 | }
 | 
|---|
| 369 | 
 | 
|---|
| 370 | static void transform1e_2(SphericalTransformIter& sti, double *s, double *t, int nk, int nl)
 | 
|---|
| 371 | {
 | 
|---|
| 372 |   const int sl = nl;
 | 
|---|
| 373 |   const int tl = sti.n();
 | 
|---|
| 374 | 
 | 
|---|
| 375 |   memset(t,0,nk*tl*sizeof(double));
 | 
|---|
| 376 | 
 | 
|---|
| 377 |   for (sti.begin(); sti.ready(); sti.next()) {
 | 
|---|
| 378 |     double *sptr = s + sti.cartindex();
 | 
|---|
| 379 |     double *tptr = t + sti.pureindex();
 | 
|---|
| 380 |     double coef = sti.coef();
 | 
|---|
| 381 |     for(int k=0; k<nk; k++,sptr+=sl,tptr+=tl) {
 | 
|---|
| 382 |         *(tptr) += coef * *(sptr);
 | 
|---|
| 383 |     }
 | 
|---|
| 384 |   }
 | 
|---|
| 385 | }
 | 
|---|
| 386 | 
 | 
|---|
| 387 | static void transform1e_vec_2(const int ntypes, SphericalTransformIter& sti, double *s, double *t, int nk, int nl)
 | 
|---|
| 388 | {
 | 
|---|
| 389 |   const int sl = nl*ntypes;
 | 
|---|
| 390 |   const int tl = sti.n()*ntypes;
 | 
|---|
| 391 | 
 | 
|---|
| 392 |   memset(t,0,nk*tl*sizeof(double));
 | 
|---|
| 393 | 
 | 
|---|
| 394 |   for (sti.begin(); sti.ready(); sti.next()) {
 | 
|---|
| 395 |     double *sptr = s + sti.cartindex()*ntypes;
 | 
|---|
| 396 |     double *tptr = t + sti.pureindex()*ntypes;
 | 
|---|
| 397 |     double coef = sti.coef();
 | 
|---|
| 398 |     for(int k=0; k<nk; k++,sptr+=sl,tptr+=tl) {
 | 
|---|
| 399 |       for(int type=0; type<ntypes; type++)
 | 
|---|
| 400 |         tptr[type] += coef * sptr[type];
 | 
|---|
| 401 |     }
 | 
|---|
| 402 |   }
 | 
|---|
| 403 | }
 | 
|---|
| 404 | 
 | 
|---|
| 405 | static void transform2e_1(SphericalTransformIter& sti, double *s, double *t, int njkl)
 | 
|---|
| 406 | {
 | 
|---|
| 407 |   memset(t,0,sti.n()*njkl*sizeof(double));
 | 
|---|
| 408 | 
 | 
|---|
| 409 |   for (sti.begin(); sti.ready(); sti.next()) {
 | 
|---|
| 410 |     double *sptr = s + sti.cartindex()*njkl;
 | 
|---|
| 411 |     double *tptr = t + sti.pureindex()*njkl;
 | 
|---|
| 412 |     double coef = sti.coef();
 | 
|---|
| 413 |     for(int jkl=0; jkl<njkl; jkl++)
 | 
|---|
| 414 |       *(tptr++) += coef * *(sptr++);
 | 
|---|
| 415 |   }
 | 
|---|
| 416 | }
 | 
|---|
| 417 | 
 | 
|---|
| 418 | static void transform2e_2(SphericalTransformIter& sti, double *s, double *t, int ni, int nj, int nkl)
 | 
|---|
| 419 | {
 | 
|---|
| 420 |   int sj = sti.n();
 | 
|---|
| 421 |   const int sjkl = nj*nkl;
 | 
|---|
| 422 |   const int tjkl = sj*nkl;
 | 
|---|
| 423 | 
 | 
|---|
| 424 |   memset(t,0,ni*tjkl*sizeof(double));
 | 
|---|
| 425 | 
 | 
|---|
| 426 |   for (sti.begin(); sti.ready(); sti.next()) {
 | 
|---|
| 427 |     double *sptr = s + sti.cartindex()*nkl;
 | 
|---|
| 428 |     double *tptr = t + sti.pureindex()*nkl;
 | 
|---|
| 429 |     double coef = sti.coef();
 | 
|---|
| 430 |     for(int i=0; i<ni; i++,sptr+=sjkl,tptr+=tjkl) {
 | 
|---|
| 431 |       for(int kl=0; kl<nkl; kl++)
 | 
|---|
| 432 |         tptr[kl] += coef * sptr[kl];
 | 
|---|
| 433 |     }
 | 
|---|
| 434 |   }
 | 
|---|
| 435 | }
 | 
|---|
| 436 | 
 | 
|---|
| 437 | static void transform2e_3(SphericalTransformIter& sti, double *s, double *t, int nij, int nk, int nl)
 | 
|---|
| 438 | {
 | 
|---|
| 439 |   int sk = sti.n();
 | 
|---|
| 440 |   const int skl = nk*nl;
 | 
|---|
| 441 |   const int tkl = sk*nl;
 | 
|---|
| 442 | 
 | 
|---|
| 443 |   memset(t,0,nij*tkl*sizeof(double));
 | 
|---|
| 444 | 
 | 
|---|
| 445 |   for (sti.begin(); sti.ready(); sti.next()) {
 | 
|---|
| 446 |     double *sptr = s + sti.cartindex()*nl;
 | 
|---|
| 447 |     double *tptr = t + sti.pureindex()*nl;
 | 
|---|
| 448 |     double coef = sti.coef();
 | 
|---|
| 449 |     for(int ij=0; ij<nij; ij++,sptr+=skl,tptr+=tkl) {
 | 
|---|
| 450 |       for(int l=0; l<nl; l++)
 | 
|---|
| 451 |         tptr[l] += coef * sptr[l];
 | 
|---|
| 452 |     }
 | 
|---|
| 453 |   }
 | 
|---|
| 454 | }
 | 
|---|
| 455 | 
 | 
|---|
| 456 | static void transform2e_4(SphericalTransformIter& sti, double *s, double *t, int nijk, int nl)
 | 
|---|
| 457 | {
 | 
|---|
| 458 |   const int sl = nl;
 | 
|---|
| 459 |   const int tl = sti.n();
 | 
|---|
| 460 | 
 | 
|---|
| 461 |   memset(t,0,nijk*tl*sizeof(double));
 | 
|---|
| 462 | 
 | 
|---|
| 463 |   for (sti.begin(); sti.ready(); sti.next()) {
 | 
|---|
| 464 |     double *sptr = s + sti.cartindex();
 | 
|---|
| 465 |     double *tptr = t + sti.pureindex();
 | 
|---|
| 466 |     double coef = sti.coef();
 | 
|---|
| 467 |     for(int ijk=0; ijk<nijk; ijk++,sptr+=sl,tptr+=tl) {
 | 
|---|
| 468 |         *(tptr) += coef * *(sptr);
 | 
|---|
| 469 |     }
 | 
|---|
| 470 |   }
 | 
|---|
| 471 | }
 | 
|---|
| 472 | 
 | 
|---|
| 473 | 
 | 
|---|
| 474 | /////////////////////////////////////////////////////////////////////////////
 | 
|---|
| 475 | 
 | 
|---|
| 476 | // Local Variables:
 | 
|---|
| 477 | // mode: c++
 | 
|---|
| 478 | // c-file-style: "CLJ-CONDENSED"
 | 
|---|
| 479 | // End:
 | 
|---|