1 ////////////////////////////////////////////////////////////////////////////////////// 2 // This file is distributed under the University of Illinois/NCSA Open Source License. 3 // See LICENSE file in top directory for details. 4 // 5 // Copyright (c) 2016 Jeongnim Kim and QMCPACK developers. 6 // 7 // File developed by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign 8 // Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign 9 // Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign 10 // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory 11 // Ye Luo, yeluo@anl.gov, Argonne National Laboratory 12 // 13 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign 14 ////////////////////////////////////////////////////////////////////////////////////// 15 16 17 #ifndef QMCPLUSPLUS_POLYNOMIAL3D_FUNCTOR_H 18 #define QMCPLUSPLUS_POLYNOMIAL3D_FUNCTOR_H 19 #include "Numerics/OptimizableFunctorBase.h" 20 #include "Numerics/HDFNumericAttrib.h" 21 #include "Utilities/ProgressReportEngine.h" 22 #include "OhmmsData/AttributeSet.h" 23 #include "Numerics/LinearFit.h" 24 #include "Numerics/DeterminantOperators.h" 25 #include <cstdio> 26 #include <algorithm> 27 28 namespace qmcplusplus 29 { 30 struct PolynomialFunctor3D : public OptimizableFunctorBase 31 { 32 typedef real_type value_type; 33 int N_eI, N_ee; 34 Array<real_type, 3> gamma; 35 // Permutation vector, used when we need to pivot 36 // columns 37 std::vector<int> GammaPerm; 38 39 Array<int, 3> index; 40 std::vector<bool> IndepVar; 41 std::vector<real_type> GammaVec, dval_Vec; 42 std::vector<TinyVector<real_type, 3>> dgrad_Vec; 43 std::vector<Tensor<real_type, 3>> dhess_Vec; 44 int NumConstraints, NumGamma; 45 Matrix<real_type> ConstraintMatrix; 46 std::vector<real_type> Parameters, d_valsFD; 47 std::vector<TinyVector<real_type, 3>> d_gradsFD; 48 std::vector<Tensor<real_type, 3>> d_hessFD; 49 std::vector<std::string> ParameterNames; 50 std::string iSpecies, eSpecies1, eSpecies2; 51 int ResetCount; 52 // Order of continuity 53 const int C; 54 real_type scale; 55 bool notOpt; 56 57 ///constructor 58 PolynomialFunctor3D(real_type ee_cusp = 0.0, real_type eI_cusp = 0.0) 59 : N_eI(0), N_ee(0), ResetCount(0), C(3), scale(1.0), notOpt(false) 60 { 61 if (std::abs(ee_cusp) > 0.0 || std::abs(eI_cusp) > 0.0) 62 { 63 app_error() << "PolynomialFunctor3D does not support nonzero cusp.\n"; 64 abort(); 65 } 66 cutoff_radius = 0.0; 67 } 68 makeClonePolynomialFunctor3D69 OptimizableFunctorBase* makeClone() const { return new PolynomialFunctor3D(*this); } 70 resizePolynomialFunctor3D71 void resize(int neI, int nee) 72 { 73 N_eI = neI; 74 N_ee = nee; 75 const double L = 0.5 * cutoff_radius; 76 gamma.resize(N_eI + 1, N_eI + 1, N_ee + 1); 77 index.resize(N_eI + 1, N_eI + 1, N_ee + 1); 78 NumGamma = ((N_eI + 1) * (N_eI + 2) / 2 * (N_ee + 1)); 79 NumConstraints = (2 * N_eI + 1) + (N_eI + N_ee + 1); 80 int numParams = NumGamma - NumConstraints; 81 Parameters.resize(numParams); 82 d_valsFD.resize(numParams); 83 d_gradsFD.resize(numParams); 84 d_hessFD.resize(numParams); 85 GammaVec.resize(NumGamma); 86 dval_Vec.resize(NumGamma); 87 dgrad_Vec.resize(NumGamma); 88 dhess_Vec.resize(NumGamma); 89 ConstraintMatrix.resize(NumConstraints, NumGamma); 90 // Assign indices 91 int num = 0; 92 for (int m = 0; m <= N_eI; m++) 93 for (int l = m; l <= N_eI; l++) 94 for (int n = 0; n <= N_ee; n++) 95 index(l, m, n) = index(m, l, n) = num++; 96 assert(num == NumGamma); 97 // std::cerr << "NumGamma = " << NumGamma << std::endl; 98 // Fill up contraint matrix 99 // For 3 constraints and 2 parameters, we would have 100 // |A00 A01 A02 A03 A04| |g0| |0 | 101 // |A11 A11 A12 A13 A14| |g1| |0 | 102 // |A22 A21 A22 A23 A24| |g2| = |0 | 103 // | 0 0 0 1 0 | |g3| |p0| 104 // | 0 0 0 0 1 | |g4| |p1| 105 ConstraintMatrix = 0.0; 106 // std::cerr << "ConstraintMatrix.size = " << ConstraintMatrix.size(0) 107 // << " by " << ConstraintMatrix.size(1) << std::endl; 108 // std::cerr << "index.size() = (" << index.size(0) << ", " 109 // << index.size(1) << ", " << index.size(2) << ").\n"; 110 int k; 111 // e-e no-cusp constraint 112 for (k = 0; k <= 2 * N_eI; k++) 113 { 114 for (int m = 0; m <= k; m++) 115 { 116 int l = k - m; 117 if (l <= N_eI && m <= N_eI) 118 { 119 int i = index(l, m, 1); 120 if (l > m) 121 ConstraintMatrix(k, i) = 2.0; 122 else if (l == m) 123 ConstraintMatrix(k, i) = 1.0; 124 } 125 } 126 } 127 // e-I no-cusp constraint 128 for (int kp = 0; kp <= N_eI + N_ee; kp++) 129 { 130 if (kp <= N_ee) 131 { 132 ConstraintMatrix(k + kp, index(0, 0, kp)) = (real_type)C; 133 ConstraintMatrix(k + kp, index(0, 1, kp)) = -L; 134 } 135 for (int l = 1; l <= kp; l++) 136 { 137 int n = kp - l; 138 if (n >= 0 && n <= N_ee && l <= N_eI) 139 { 140 ConstraintMatrix(k + kp, index(l, 0, n)) = (real_type)C; 141 ConstraintMatrix(k + kp, index(l, 1, n)) = -L; 142 } 143 } 144 } 145 // { 146 // fprintf (stderr, "Constraint matrix:\n"); 147 // for (int i=0; i<NumConstraints; i++) { 148 // for (int j=0; j<NumGamma; j++) 149 // fprintf (stderr, "%5.2f ", ConstraintMatrix(i,j)); 150 // fprintf(stderr, "\n"); 151 // } 152 // } 153 // Now, row-reduce constraint matrix 154 GammaPerm.resize(NumGamma); 155 IndepVar.resize(NumGamma, false); 156 // Set identity permutation 157 for (int i = 0; i < NumGamma; i++) 158 GammaPerm[i] = i; 159 int col = -1; 160 for (int row = 0; row < NumConstraints; row++) 161 { 162 int max_loc; 163 real_type max_abs; 164 do 165 { 166 col++; 167 max_loc = row; 168 max_abs = std::abs(ConstraintMatrix(row, col)); 169 for (int ri = row + 1; ri < NumConstraints; ri++) 170 { 171 real_type abs_val = std::abs(ConstraintMatrix(ri, col)); 172 if (abs_val > max_abs) 173 { 174 max_loc = ri; 175 max_abs = abs_val; 176 } 177 } 178 if (max_abs < 1.0e-6) 179 IndepVar[col] = true; 180 } while (max_abs < 1.0e-6); 181 ConstraintMatrix.swap_rows(row, max_loc); 182 real_type lead_inv = 1.0 / ConstraintMatrix(row, col); 183 for (int c = 0; c < NumGamma; c++) 184 ConstraintMatrix(row, c) *= lead_inv; 185 // Now, eliminate column entries 186 for (int ri = 0; ri < NumConstraints; ri++) 187 { 188 if (ri != row) 189 { 190 real_type val = ConstraintMatrix(ri, col); 191 for (int c = 0; c < NumGamma; c++) 192 ConstraintMatrix(ri, c) -= val * ConstraintMatrix(row, c); 193 } 194 } 195 } 196 for (int c = col + 1; c < NumGamma; c++) 197 IndepVar[c] = true; 198 // fprintf (stderr, "Reduced Constraint matrix:\n"); 199 // for (int i=0; i<NumConstraints; i++) { 200 // for (int j=0; j<NumGamma; j++) 201 // fprintf (stderr, "%5.2f ", ConstraintMatrix(i,j)); 202 // fprintf(stderr, "\n"); 203 // } 204 // fprintf (stderr, "Independent vars = \n"); 205 // for (int i=0; i<NumGamma; i++) 206 // if (IndepVar[i]) 207 // fprintf (stderr, "%d ", i); 208 // fprintf (stderr, "\n"); 209 // fprintf (stderr, "Inverse matrix:\n"); 210 // // Now, invert constraint matrix 211 // Invert(ConstraintMatrix.data(), NumGamma, NumGamma); 212 // for (int i=0; i<NumGamma; i++) { 213 // for (int j=0; j<NumGamma; j++) 214 // fprintf (stderr, "%5.2f ", ConstraintMatrix(i,j)); 215 // fprintf(stderr, "\n"); 216 // } 217 } 218 resetPolynomialFunctor3D219 void reset() 220 { 221 resize(N_eI, N_ee); 222 reset_gamma(); 223 } 224 reset_gammaPolynomialFunctor3D225 void reset_gamma() 226 { 227 // fprintf (stderr, "Parameters:\n"); 228 // for (int i=0; i<Parameters.size(); i++) 229 // fprintf (stderr, " %16.10e\n", Parameters[i]); 230 const double L = 0.5 * cutoff_radius; 231 std::fill(GammaVec.begin(), GammaVec.end(), 0.0); 232 // First, set all independent variables 233 int var = 0; 234 for (int i = 0; i < NumGamma; i++) 235 if (IndepVar[i]) 236 GammaVec[i] = scale * Parameters[var++]; 237 assert(var == Parameters.size()); 238 // Now, set dependent variables 239 var = 0; 240 // std::cerr << "NumConstraints = " << NumConstraints << std::endl; 241 for (int i = 0; i < NumGamma; i++) 242 if (!IndepVar[i]) 243 { 244 // fprintf (stderr, "constraintMatrix(%d,%d) = %1.10f\n", 245 // var, i, ConstraintMatrix(var,i)); 246 assert(std::abs(ConstraintMatrix(var, i) - 1.0) < 1.0e-6); 247 for (int j = 0; j < NumGamma; j++) 248 if (i != j) 249 GammaVec[i] -= ConstraintMatrix(var, j) * GammaVec[j]; 250 var++; 251 } 252 int num = 0; 253 for (int m = 0; m <= N_eI; m++) 254 for (int l = m; l <= N_eI; l++) 255 for (int n = 0; n <= N_ee; n++) 256 // gamma(m,l,n) = gamma(l,m,n) = unpermuted[num++]; 257 gamma(m, l, n) = gamma(l, m, n) = GammaVec[num++]; 258 // Now check that constraints have been satisfied 259 // e-e constraints 260 for (int k = 0; k <= 2 * N_eI; k++) 261 { 262 real_type sum = 0.0; 263 for (int m = 0; m <= k; m++) 264 { 265 int l = k - m; 266 if (l <= N_eI && m <= N_eI) 267 { 268 int i = index(l, m, 1); 269 if (l > m) 270 sum += 2.0 * GammaVec[i]; 271 else if (l == m) 272 sum += GammaVec[i]; 273 } 274 } 275 if (std::abs(sum) > 1.0e-9) 276 std::cerr << "error in k = " << k << " sum = " << sum << std::endl; 277 } 278 for (int k = 0; k <= 2 * N_eI; k++) 279 { 280 real_type sum = 0.0; 281 for (int l = 0; l <= k; l++) 282 { 283 int m = k - l; 284 if (m <= N_eI && l <= N_eI) 285 { 286 // fprintf (stderr, "k = %d gamma(%d, %d, 1) = %1.8f\n", k, l, m, 287 // gamma(l,m,1)); 288 sum += gamma(l, m, 1); 289 } 290 } 291 if (std::abs(sum) > 1.0e-6) 292 { 293 app_error() << "e-e constraint not satisfied in PolynomialFunctor3D: k=" << k << " sum=" << sum << std::endl; 294 abort(); 295 } 296 } 297 // e-I constraints 298 for (int k = 0; k <= N_eI + N_ee; k++) 299 { 300 real_type sum = 0.0; 301 for (int m = 0; m <= k; m++) 302 { 303 int n = k - m; 304 if (m <= N_eI && n <= N_ee) 305 { 306 sum += (real_type)C * gamma(0, m, n) - L * gamma(1, m, n); 307 // fprintf (stderr, 308 // "k = %d gamma(0,%d,%d) = %1.8f gamma(1,%d,%d)=%1.8f\n", 309 // k, m, n, gamma(0,m,n), m, n, gamma(1,m,n)); 310 } 311 } 312 if (std::abs(sum) > 1.0e-6) 313 { 314 app_error() << "e-I constraint not satisfied in PolynomialFunctor3D: k=" << k << " sum=" << sum << std::endl; 315 abort(); 316 } 317 } 318 } 319 evaluatePolynomialFunctor3D320 inline real_type evaluate(real_type r_12, real_type r_1I, real_type r_2I) const 321 { 322 constexpr real_type czero(0); 323 constexpr real_type cone(1); 324 constexpr real_type chalf(0.5); 325 326 const real_type L = chalf * cutoff_radius; 327 if (r_1I >= L || r_2I >= L) 328 return czero; 329 real_type val = czero; 330 real_type r2l(cone); 331 for (int l = 0; l <= N_eI; l++) 332 { 333 real_type r2m(r2l); 334 for (int m = 0; m <= N_eI; m++) 335 { 336 real_type r2n(r2m); 337 for (int n = 0; n <= N_ee; n++) 338 { 339 val += gamma(l, m, n) * r2n; 340 r2n *= r_12; 341 } 342 r2m *= r_2I; 343 } 344 r2l *= r_1I; 345 } 346 for (int i = 0; i < C; i++) 347 val *= (r_1I - L) * (r_2I - L); 348 return val; 349 } 350 351 // assume r_1I < L && r_2I < L, compression and screening is handled outside evaluateVPolynomialFunctor3D352 inline real_type evaluateV(int Nptcl, 353 const real_type* restrict r_12_array, 354 const real_type* restrict r_1I_array, 355 const real_type* restrict r_2I_array) const 356 { 357 constexpr real_type czero(0); 358 constexpr real_type cone(1); 359 constexpr real_type chalf(0.5); 360 361 const real_type L = chalf * cutoff_radius; 362 real_type val_tot = czero; 363 364 #pragma omp simd aligned(r_12_array, r_1I_array, r_2I_array: QMC_SIMD_ALIGNMENT) reduction(+ : val_tot) 365 for (int ptcl = 0; ptcl < Nptcl; ptcl++) 366 { 367 const real_type r_12 = r_12_array[ptcl]; 368 const real_type r_1I = r_1I_array[ptcl]; 369 const real_type r_2I = r_2I_array[ptcl]; 370 real_type val = czero; 371 real_type r2l(cone); 372 for (int l = 0; l <= N_eI; l++) 373 { 374 real_type r2m(r2l); 375 for (int m = 0; m <= N_eI; m++) 376 { 377 real_type r2n(r2m); 378 for (int n = 0; n <= N_ee; n++) 379 { 380 val += gamma(l, m, n) * r2n; 381 r2n *= r_12; 382 } 383 r2m *= r_2I; 384 } 385 r2l *= r_1I; 386 } 387 const real_type both_minus_L = (r_2I - L) * (r_1I - L); 388 for (int i = 0; i < C; i++) 389 val *= both_minus_L; 390 val_tot += val; 391 } 392 393 return val_tot; 394 } 395 evaluatePolynomialFunctor3D396 inline real_type evaluate(real_type r_12, 397 real_type r_1I, 398 real_type r_2I, 399 TinyVector<real_type, 3>& grad, 400 Tensor<real_type, 3>& hess) const 401 { 402 constexpr real_type czero(0); 403 constexpr real_type cone(1); 404 constexpr real_type chalf(0.5); 405 constexpr real_type ctwo(2); 406 407 const real_type L = chalf * cutoff_radius; 408 if (r_1I >= L || r_2I >= L) 409 { 410 grad = czero; 411 hess = czero; 412 return czero; 413 } 414 real_type val = czero; 415 grad = czero; 416 hess = czero; 417 real_type r2l(cone), r2l_1(czero), r2l_2(czero), lf(czero); 418 for (int l = 0; l <= N_eI; l++) 419 { 420 real_type r2m(cone), r2m_1(czero), r2m_2(czero), mf(czero); 421 for (int m = 0; m <= N_eI; m++) 422 { 423 real_type r2n(cone), r2n_1(czero), r2n_2(czero), nf(czero); 424 for (int n = 0; n <= N_ee; n++) 425 { 426 const real_type g = gamma(l, m, n); 427 const real_type g00x = g * r2l * r2m; 428 const real_type g10x = g * r2l_1 * r2m; 429 const real_type g01x = g * r2l * r2m_1; 430 const real_type gxx0 = g * r2n; 431 432 val += g00x * r2n; 433 grad[0] += g00x * r2n_1; 434 grad[1] += g10x * r2n; 435 grad[2] += g01x * r2n; 436 hess(0, 0) += g00x * r2n_2; 437 hess(0, 1) += g10x * r2n_1; 438 hess(0, 2) += g01x * r2n_1; 439 hess(1, 1) += gxx0 * r2l_2 * r2m; 440 hess(1, 2) += gxx0 * r2l_1 * r2m_1; 441 hess(2, 2) += gxx0 * r2l * r2m_2; 442 nf += cone; 443 r2n_2 = r2n_1 * nf; 444 r2n_1 = r2n * nf; 445 r2n *= r_12; 446 } 447 mf += cone; 448 r2m_2 = r2m_1 * mf; 449 r2m_1 = r2m * mf; 450 r2m *= r_2I; 451 } 452 lf += cone; 453 r2l_2 = r2l_1 * lf; 454 r2l_1 = r2l * lf; 455 r2l *= r_1I; 456 } 457 458 const real_type r_2I_minus_L = r_2I - L; 459 const real_type r_1I_minus_L = r_1I - L; 460 const real_type both_minus_L = r_2I_minus_L * r_1I_minus_L; 461 for (int i = 0; i < C; i++) 462 { 463 hess(0, 0) = both_minus_L * hess(0, 0); 464 hess(0, 1) = both_minus_L * hess(0, 1) + r_2I_minus_L * grad[0]; 465 hess(0, 2) = both_minus_L * hess(0, 2) + r_1I_minus_L * grad[0]; 466 hess(1, 1) = both_minus_L * hess(1, 1) + ctwo * r_2I_minus_L * grad[1]; 467 hess(1, 2) = both_minus_L * hess(1, 2) + r_1I_minus_L * grad[1] + r_2I_minus_L * grad[2] + val; 468 hess(2, 2) = both_minus_L * hess(2, 2) + ctwo * r_1I_minus_L * grad[2]; 469 grad[0] = both_minus_L * grad[0]; 470 grad[1] = both_minus_L * grad[1] + r_2I_minus_L * val; 471 grad[2] = both_minus_L * grad[2] + r_1I_minus_L * val; 472 val *= both_minus_L; 473 } 474 hess(1, 0) = hess(0, 1); 475 hess(2, 0) = hess(0, 2); 476 hess(2, 1) = hess(1, 2); 477 return val; 478 } 479 480 // assume r_1I < L && r_2I < L, compression and screening is handled outside evaluateVGLPolynomialFunctor3D481 inline void evaluateVGL(int Nptcl, 482 const real_type* restrict r_12_array, 483 const real_type* restrict r_1I_array, 484 const real_type* restrict r_2I_array, 485 real_type* restrict val_array, 486 real_type* restrict grad0_array, 487 real_type* restrict grad1_array, 488 real_type* restrict grad2_array, 489 real_type* restrict hess00_array, 490 real_type* restrict hess11_array, 491 real_type* restrict hess22_array, 492 real_type* restrict hess01_array, 493 real_type* restrict hess02_array) const 494 { 495 constexpr real_type czero(0); 496 constexpr real_type cone(1); 497 constexpr real_type chalf(0.5); 498 constexpr real_type ctwo(2); 499 500 const real_type L = chalf * cutoff_radius; 501 #pragma omp simd aligned(r_12_array, \ 502 r_1I_array, \ 503 r_2I_array, \ 504 val_array, \ 505 grad0_array, \ 506 grad1_array, \ 507 grad2_array, \ 508 hess00_array, \ 509 hess11_array, \ 510 hess22_array, \ 511 hess01_array, \ 512 hess02_array: QMC_SIMD_ALIGNMENT) 513 for (int ptcl = 0; ptcl < Nptcl; ptcl++) 514 { 515 const real_type r_12 = r_12_array[ptcl]; 516 const real_type r_1I = r_1I_array[ptcl]; 517 const real_type r_2I = r_2I_array[ptcl]; 518 519 real_type val(czero); 520 real_type grad0(czero); 521 real_type grad1(czero); 522 real_type grad2(czero); 523 real_type hess00(czero); 524 real_type hess11(czero); 525 real_type hess22(czero); 526 real_type hess01(czero); 527 real_type hess02(czero); 528 529 real_type r2l(cone), r2l_1(czero), r2l_2(czero), lf(czero); 530 for (int l = 0; l <= N_eI; l++) 531 { 532 real_type r2m(cone), r2m_1(czero), r2m_2(czero), mf(czero); 533 for (int m = 0; m <= N_eI; m++) 534 { 535 real_type r2n(cone), r2n_1(czero), r2n_2(czero), nf(czero); 536 for (int n = 0; n <= N_ee; n++) 537 { 538 const real_type g = gamma(l, m, n); 539 const real_type g00x = g * r2l * r2m; 540 const real_type g10x = g * r2l_1 * r2m; 541 const real_type g01x = g * r2l * r2m_1; 542 const real_type gxx0 = g * r2n; 543 544 val += g00x * r2n; 545 grad0 += g00x * r2n_1; 546 grad1 += g10x * r2n; 547 grad2 += g01x * r2n; 548 hess00 += g00x * r2n_2; 549 hess01 += g10x * r2n_1; 550 hess02 += g01x * r2n_1; 551 hess11 += gxx0 * r2l_2 * r2m; 552 hess22 += gxx0 * r2l * r2m_2; 553 nf += cone; 554 r2n_2 = r2n_1 * nf; 555 r2n_1 = r2n * nf; 556 r2n *= r_12; 557 } 558 mf += cone; 559 r2m_2 = r2m_1 * mf; 560 r2m_1 = r2m * mf; 561 r2m *= r_2I; 562 } 563 lf += cone; 564 r2l_2 = r2l_1 * lf; 565 r2l_1 = r2l * lf; 566 r2l *= r_1I; 567 } 568 569 const real_type r_2I_minus_L = r_2I - L; 570 const real_type r_1I_minus_L = r_1I - L; 571 const real_type both_minus_L = r_2I_minus_L * r_1I_minus_L; 572 for (int i = 0; i < C; i++) 573 { 574 hess00 = both_minus_L * hess00; 575 hess01 = both_minus_L * hess01 + r_2I_minus_L * grad0; 576 hess02 = both_minus_L * hess02 + r_1I_minus_L * grad0; 577 hess11 = both_minus_L * hess11 + ctwo * r_2I_minus_L * grad1; 578 hess22 = both_minus_L * hess22 + ctwo * r_1I_minus_L * grad2; 579 grad0 = both_minus_L * grad0; 580 grad1 = both_minus_L * grad1 + r_2I_minus_L * val; 581 grad2 = both_minus_L * grad2 + r_1I_minus_L * val; 582 val *= both_minus_L; 583 } 584 585 val_array[ptcl] = val; 586 grad0_array[ptcl] = grad0 / r_12; 587 grad1_array[ptcl] = grad1 / r_1I; 588 grad2_array[ptcl] = grad2 / r_2I; 589 hess00_array[ptcl] = hess00; 590 hess11_array[ptcl] = hess11; 591 hess22_array[ptcl] = hess22; 592 hess01_array[ptcl] = hess01 / (r_12 * r_1I); 593 hess02_array[ptcl] = hess02 / (r_12 * r_2I); 594 } 595 } 596 597 evaluatePolynomialFunctor3D598 inline real_type evaluate(const real_type r_12, 599 const real_type r_1I, 600 const real_type r_2I, 601 TinyVector<real_type, 3>& grad, 602 Tensor<real_type, 3>& hess, 603 TinyVector<Tensor<real_type, 3>, 3>& d3) 604 { 605 grad = 0.0; 606 hess = 0.0; 607 d3 = grad; 608 const real_type L = 0.5 * cutoff_radius; 609 if (r_1I >= L || r_2I >= L) 610 return 0.0; 611 real_type val = 0.0; 612 real_type r2l(1.0), r2l_1(0.0), r2l_2(0.0), r2l_3(0.0), lf(0.0); 613 for (int l = 0; l <= N_eI; l++) 614 { 615 real_type r2m(1.0), r2m_1(0.0), r2m_2(0.0), r2m_3(0.0), mf(0.0); 616 for (int m = 0; m <= N_eI; m++) 617 { 618 real_type r2n(1.0), r2n_1(0.0), r2n_2(0.0), r2n_3(0.0), nf(0.0); 619 for (int n = 0; n <= N_ee; n++) 620 { 621 real_type g = gamma(l, m, n); 622 val += g * r2l * r2m * r2n; 623 grad[0] += nf * g * r2l * r2m * r2n_1; 624 grad[1] += lf * g * r2l_1 * r2m * r2n; 625 grad[2] += mf * g * r2l * r2m_1 * r2n; 626 hess(0, 0) += nf * (nf - 1.0) * g * r2l * r2m * r2n_2; 627 hess(0, 1) += nf * lf * g * r2l_1 * r2m * r2n_1; 628 hess(0, 2) += nf * mf * g * r2l * r2m_1 * r2n_1; 629 hess(1, 1) += lf * (lf - 1.0) * g * r2l_2 * r2m * r2n; 630 hess(1, 2) += lf * mf * g * r2l_1 * r2m_1 * r2n; 631 hess(2, 2) += mf * (mf - 1.0) * g * r2l * r2m_2 * r2n; 632 d3[0](0, 0) += nf * (nf - 1.0) * (nf - 2.0) * g * r2l * r2m * r2n_3; 633 d3[0](0, 1) += nf * (nf - 1.0) * lf * g * r2l_1 * r2m * r2n_2; 634 d3[0](0, 2) += nf * (nf - 1.0) * mf * g * r2l * r2m_1 * r2n_2; 635 d3[0](1, 1) += nf * lf * (lf - 1.0) * g * r2l_2 * r2m * r2n_1; 636 d3[0](1, 2) += nf * lf * mf * g * r2l_1 * r2m_1 * r2n_1; 637 d3[0](2, 2) += nf * mf * (mf - 1.0) * g * r2l * r2m_2 * r2n_1; 638 d3[1](1, 1) += lf * (lf - 1.0) * (lf - 2.0) * g * r2l_3 * r2m * r2n; 639 d3[1](1, 2) += lf * (lf - 1.0) * mf * g * r2l_2 * r2m_1 * r2n; 640 d3[1](2, 2) += lf * mf * (mf - 1.0) * g * r2l_1 * r2m_2 * r2n; 641 d3[2](2, 2) += mf * (mf - 1.0) * (mf - 2.0) * g * r2l * r2m_3 * r2n; 642 r2n_3 = r2n_2; 643 r2n_2 = r2n_1; 644 r2n_1 = r2n; 645 r2n *= r_12; 646 nf += 1.0; 647 } 648 r2m_3 = r2m_2; 649 r2m_2 = r2m_1; 650 r2m_1 = r2m; 651 r2m *= r_2I; 652 mf += 1.0; 653 } 654 r2l_3 = r2l_2; 655 r2l_2 = r2l_1; 656 r2l_1 = r2l; 657 r2l *= r_1I; 658 lf += 1.0; 659 } 660 for (int i = 0; i < C; i++) 661 { 662 d3[0](0, 0) = (r_1I - L) * (r_2I - L) * d3[0](0, 0); 663 d3[0](0, 1) = (r_1I - L) * (r_2I - L) * d3[0](0, 1) + (r_2I - L) * hess(0, 0); 664 d3[0](0, 2) = (r_1I - L) * (r_2I - L) * d3[0](0, 2) + (r_1I - L) * hess(0, 0); 665 d3[0](1, 1) = (r_1I - L) * (r_2I - L) * d3[0](1, 1) + 2.0 * (r_2I - L) * hess(0, 1); 666 d3[0](1, 2) = (r_1I - L) * (r_2I - L) * d3[0](1, 2) + (r_1I - L) * hess(0, 1) + (r_2I - L) * hess(0, 2) + grad[0]; 667 d3[0](2, 2) = (r_1I - L) * (r_2I - L) * d3[0](2, 2) + 2.0 * (r_1I - L) * hess(0, 2); 668 d3[1](1, 1) = (r_1I - L) * (r_2I - L) * d3[1](1, 1) + 3.0 * (r_2I - L) * hess(1, 1); 669 d3[1](1, 2) = (r_1I - L) * (r_2I - L) * d3[1](1, 2) + 2.0 * (r_2I - L) * hess(1, 2) + 2.0 * grad[1] + 670 (r_1I - L) * hess(1, 1); 671 d3[1](2, 2) = (r_1I - L) * (r_2I - L) * d3[1](2, 2) + 2.0 * (r_1I - L) * hess(1, 2) + 2.0 * grad[2] + 672 (r_2I - L) * hess(2, 2); 673 d3[2](2, 2) = (r_1I - L) * (r_2I - L) * d3[2](2, 2) + 3.0 * (r_1I - L) * hess(2, 2); 674 hess(0, 0) = (r_1I - L) * (r_2I - L) * hess(0, 0); 675 hess(0, 1) = (r_1I - L) * (r_2I - L) * hess(0, 1) + (r_2I - L) * grad[0]; 676 hess(0, 2) = (r_1I - L) * (r_2I - L) * hess(0, 2) + (r_1I - L) * grad[0]; 677 hess(1, 1) = (r_1I - L) * (r_2I - L) * hess(1, 1) + 2.0 * (r_2I - L) * grad[1]; 678 hess(1, 2) = (r_1I - L) * (r_2I - L) * hess(1, 2) + (r_1I - L) * grad[1] + (r_2I - L) * grad[2] + val; 679 hess(2, 2) = (r_1I - L) * (r_2I - L) * hess(2, 2) + 2.0 * (r_1I - L) * grad[2]; 680 grad[0] = (r_1I - L) * (r_2I - L) * grad[0]; 681 grad[1] = (r_1I - L) * (r_2I - L) * grad[1] + (r_2I - L) * val; 682 grad[2] = (r_1I - L) * (r_2I - L) * grad[2] + (r_1I - L) * val; 683 val *= (r_1I - L) * (r_2I - L); 684 } 685 hess(1, 0) = hess(0, 1); 686 hess(2, 0) = hess(0, 2); 687 hess(2, 1) = hess(1, 2); 688 d3[0](1, 0) = d3[0](0, 1); 689 d3[0](2, 0) = d3[0](0, 2); 690 d3[0](2, 1) = d3[0](1, 2); 691 d3[1](0, 0) = d3[0](1, 1); 692 d3[1](0, 1) = d3[0](0, 1); 693 d3[1](0, 2) = d3[0](1, 2); 694 d3[1](1, 0) = d3[0](0, 1); 695 d3[1](2, 0) = d3[0](1, 2); 696 d3[1](2, 1) = d3[1](1, 2); 697 d3[2](0, 0) = d3[0](0, 2); 698 d3[2](0, 1) = d3[0](1, 2); 699 d3[2](0, 2) = d3[0](2, 2); 700 d3[2](1, 0) = d3[0](1, 2); 701 d3[2](1, 1) = d3[1](1, 2); 702 d3[2](1, 2) = d3[1](2, 2); 703 d3[2](2, 0) = d3[0](2, 2); 704 d3[2](2, 1) = d3[1](2, 2); 705 return val; 706 } 707 708 evaluatePolynomialFunctor3D709 inline real_type evaluate(const real_type r, const real_type rinv) { return 0.0; } 710 711 evaluateDerivativesFDPolynomialFunctor3D712 inline bool evaluateDerivativesFD(const real_type r_12, 713 const real_type r_1I, 714 const real_type r_2I, 715 std::vector<double>& d_vals, 716 std::vector<TinyVector<real_type, 3>>& d_grads, 717 std::vector<Tensor<real_type, 3>>& d_hess) 718 { 719 const real_type eps = 1.0e-6; 720 assert(d_vals.size() == Parameters.size()); 721 assert(d_grads.size() == Parameters.size()); 722 assert(d_hess.size() == Parameters.size()); 723 for (int ip = 0; ip < Parameters.size(); ip++) 724 { 725 real_type v_plus, v_minus; 726 TinyVector<real_type, 3> g_plus, g_minus; 727 Tensor<real_type, 3> h_plus, h_minus; 728 real_type save_p = Parameters[ip]; 729 Parameters[ip] = save_p + eps; 730 reset_gamma(); 731 v_plus = evaluate(r_12, r_1I, r_2I, g_plus, h_plus); 732 Parameters[ip] = save_p - eps; 733 reset_gamma(); 734 v_minus = evaluate(r_12, r_1I, r_2I, g_minus, h_minus); 735 Parameters[ip] = save_p; 736 reset_gamma(); 737 real_type dp_inv = 0.5 / eps; 738 d_vals[ip] = dp_inv * (v_plus - v_minus); 739 d_grads[ip] = dp_inv * (g_plus - g_minus); 740 d_hess[ip] = dp_inv * (h_plus - h_minus); 741 } 742 return true; 743 } 744 745 evaluateDerivativesPolynomialFunctor3D746 inline bool evaluateDerivatives(const real_type r_12, 747 const real_type r_1I, 748 const real_type r_2I, 749 std::vector<real_type>& d_vals, 750 std::vector<TinyVector<real_type, 3>>& d_grads, 751 std::vector<Tensor<real_type, 3>>& d_hess) 752 { 753 const real_type L = 0.5 * cutoff_radius; 754 if (r_1I >= L || r_2I >= L) 755 return false; 756 757 constexpr real_type czero(0); 758 constexpr real_type cone(1); 759 constexpr real_type ctwo(2); 760 761 real_type dval_dgamma; 762 TinyVector<real_type, 3> dgrad_dgamma; 763 Tensor<real_type, 3> dhess_dgamma; 764 765 for (int i = 0; i < dval_Vec.size(); i++) 766 { 767 dval_Vec[i] = czero; 768 dgrad_Vec[i] = czero; 769 dhess_Vec[i] = czero; 770 } 771 772 const real_type r_2I_minus_L = r_2I - L; 773 const real_type r_1I_minus_L = r_1I - L; 774 const real_type both_minus_L = r_2I_minus_L * r_1I_minus_L; 775 776 real_type r2l(cone), r2l_1(czero), r2l_2(czero), lf(czero); 777 for (int l = 0; l <= N_eI; l++) 778 { 779 real_type r2m(cone), r2m_1(czero), r2m_2(czero), mf(czero); 780 for (int m = 0; m <= N_eI; m++) 781 { 782 int num; 783 if (m > l) 784 num = ((2 * N_eI - l + 3) * l / 2 + m - l) * (N_ee + 1); 785 else 786 num = ((2 * N_eI - m + 3) * m / 2 + l - m) * (N_ee + 1); 787 real_type r2n(cone), r2n_1(czero), r2n_2(czero), nf(czero); 788 for (int n = 0; n <= N_ee; n++, num++) 789 { 790 dval_dgamma = r2l * r2m * r2n; 791 dgrad_dgamma[0] = r2l * r2m * r2n_1; 792 dgrad_dgamma[1] = r2l_1 * r2m * r2n; 793 dgrad_dgamma[2] = r2l * r2m_1 * r2n; 794 dhess_dgamma(0, 0) = r2l * r2m * r2n_2; 795 dhess_dgamma(0, 1) = r2l_1 * r2m * r2n_1; 796 dhess_dgamma(0, 2) = r2l * r2m_1 * r2n_1; 797 dhess_dgamma(1, 1) = r2l_2 * r2m * r2n; 798 dhess_dgamma(1, 2) = r2l_1 * r2m_1 * r2n; 799 dhess_dgamma(2, 2) = r2l * r2m_2 * r2n; 800 801 for (int i = 0; i < C; i++) 802 { 803 dhess_dgamma(0, 0) = both_minus_L * dhess_dgamma(0, 0); 804 dhess_dgamma(0, 1) = both_minus_L * dhess_dgamma(0, 1) + r_2I_minus_L * dgrad_dgamma[0]; 805 dhess_dgamma(0, 2) = both_minus_L * dhess_dgamma(0, 2) + r_1I_minus_L * dgrad_dgamma[0]; 806 dhess_dgamma(1, 1) = both_minus_L * dhess_dgamma(1, 1) + ctwo * r_2I_minus_L * dgrad_dgamma[1]; 807 dhess_dgamma(1, 2) = both_minus_L * dhess_dgamma(1, 2) + r_1I_minus_L * dgrad_dgamma[1] + 808 r_2I_minus_L * dgrad_dgamma[2] + dval_dgamma; 809 dhess_dgamma(2, 2) = both_minus_L * dhess_dgamma(2, 2) + ctwo * r_1I_minus_L * dgrad_dgamma[2]; 810 dgrad_dgamma[0] = both_minus_L * dgrad_dgamma[0]; 811 dgrad_dgamma[1] = both_minus_L * dgrad_dgamma[1] + r_2I_minus_L * dval_dgamma; 812 dgrad_dgamma[2] = both_minus_L * dgrad_dgamma[2] + r_1I_minus_L * dval_dgamma; 813 dval_dgamma *= both_minus_L; 814 } 815 816 // Now, pack into vectors 817 dval_Vec[num] += scale * dval_dgamma; 818 for (int i = 0; i < 3; i++) 819 { 820 dgrad_Vec[num][i] += scale * dgrad_dgamma[i]; 821 dhess_Vec[num](i, i) += scale * dhess_dgamma(i, i); 822 for (int j = i + 1; j < 3; j++) 823 { 824 dhess_Vec[num](i, j) += scale * dhess_dgamma(i, j); 825 dhess_Vec[num](j, i) = dhess_Vec[num](i, j); 826 } 827 } 828 829 nf += cone; 830 r2n_2 = r2n_1 * nf; 831 r2n_1 = r2n * nf; 832 r2n *= r_12; 833 } 834 mf += cone; 835 r2m_2 = r2m_1 * mf; 836 r2m_1 = r2m * mf; 837 r2m *= r_2I; 838 } 839 lf += cone; 840 r2l_2 = r2l_1 * lf; 841 r2l_1 = r2l * lf; 842 r2l *= r_1I; 843 } 844 // for (int i=0; i<dval_Vec.size(); i++) 845 // fprintf (stderr, "dval_Vec[%d] = %12.6e\n", i, dval_Vec[i]); 846 /////////////////////////////////////////// 847 // Now, compensate for constraint matrix // 848 /////////////////////////////////////////// 849 std::fill(d_vals.begin(), d_vals.end(), 0.0); 850 int var = 0; 851 for (int i = 0; i < NumGamma; i++) 852 if (IndepVar[i]) 853 { 854 d_vals[var] = dval_Vec[i]; 855 d_grads[var] = dgrad_Vec[i]; 856 d_hess[var] = dhess_Vec[i]; 857 var++; 858 } 859 int constraint = 0; 860 for (int i = 0; i < NumGamma; i++) 861 { 862 if (!IndepVar[i]) 863 { 864 int indep_var = 0; 865 for (int j = 0; j < NumGamma; j++) 866 if (IndepVar[j]) 867 { 868 d_vals[indep_var] -= ConstraintMatrix(constraint, j) * dval_Vec[i]; 869 d_grads[indep_var] -= ConstraintMatrix(constraint, j) * dgrad_Vec[i]; 870 d_hess[indep_var] -= ConstraintMatrix(constraint, j) * dhess_Vec[i]; 871 indep_var++; 872 } 873 else if (i != j) 874 assert(std::abs(ConstraintMatrix(constraint, j)) < 1.0e-10); 875 constraint++; 876 } 877 } 878 return true; 879 #ifdef DEBUG_DERIVS 880 evaluateDerivativesFD(r_12, r_1I, r_2I, d_valsFD, d_gradsFD, d_hessFD); 881 fprintf(stderr, "Param Analytic Finite diffference\n"); 882 for (int ip = 0; ip < Parameters.size(); ip++) 883 fprintf(stderr, " %3d %12.6e %12.6e\n", ip, d_vals[ip], d_valsFD[ip]); 884 fprintf(stderr, "Param Analytic Finite diffference\n"); 885 for (int ip = 0; ip < Parameters.size(); ip++) 886 fprintf(stderr, 887 " %3d %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e\n", 888 ip, 889 d_grads[ip][0], 890 d_gradsFD[ip][0], 891 d_grads[ip][1], 892 d_gradsFD[ip][1], 893 d_grads[ip][2], 894 d_gradsFD[ip][2]); 895 fprintf(stderr, "Param Analytic Finite diffference\n"); 896 for (int ip = 0; ip < Parameters.size(); ip++) 897 for (int dim = 0; dim < 3; dim++) 898 fprintf(stderr, 899 " %3d %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e\n", 900 ip, 901 d_hess[ip](0, dim), 902 d_hessFD[ip](0, dim), 903 d_hess[ip](1, dim), 904 d_hessFD[ip](1, dim), 905 d_hess[ip](2, dim), 906 d_hessFD[ip](2, dim)); 907 #endif 908 } 909 fPolynomialFunctor3D910 inline real_type f(real_type r) { return 0.0; } dfPolynomialFunctor3D911 inline real_type df(real_type r) { return 0.0; } 912 putPolynomialFunctor3D913 bool put(xmlNodePtr cur) 914 { 915 ReportEngine PRE("PolynomialFunctor3D", "put(xmlNodePtr)"); 916 // //CuspValue = -1.0e10; 917 // NumParams_eI = NumParams_ee = 0; 918 cutoff_radius = 0.0; 919 OhmmsAttributeSet rAttrib; 920 rAttrib.add(N_ee, "esize"); 921 rAttrib.add(N_eI, "isize"); 922 rAttrib.add(cutoff_radius, "rcut"); 923 rAttrib.put(cur); 924 if (N_eI == 0) 925 PRE.error("You must specify a positive number for \"isize\"", true); 926 if (N_ee == 0) 927 PRE.error("You must specify a positive number for \"esize\"", true); 928 app_summary() << " Ion: " << iSpecies << " electron-electron: " << eSpecies1 << " - " << eSpecies2 << std::endl; 929 app_summary() << " Number of parameters for e-e: " << N_ee << ", for e-I: " << N_eI << std::endl; 930 app_summary() << " Cutoff radius: " << cutoff_radius << std::endl; 931 app_summary() << std::endl; 932 resize(N_eI, N_ee); 933 // Now read coefficents 934 xmlNodePtr xmlCoefs = cur->xmlChildrenNode; 935 while (xmlCoefs != NULL) 936 { 937 std::string cname((const char*)xmlCoefs->name); 938 if (cname == "coefficients") 939 { 940 std::string type("0"), id("0"), opt("yes"); 941 OhmmsAttributeSet cAttrib; 942 cAttrib.add(id, "id"); 943 cAttrib.add(type, "type"); 944 cAttrib.add(opt, "optimize"); 945 cAttrib.put(xmlCoefs); 946 notOpt = (opt == "no"); 947 if (type != "Array") 948 { 949 PRE.error("Unknown correlation type " + type + " in PolynomialFunctor3D." + "Resetting to \"Array\""); 950 xmlNewProp(xmlCoefs, (const xmlChar*)"type", (const xmlChar*)"Array"); 951 } 952 std::vector<real_type> params; 953 putContent(params, xmlCoefs); 954 if (params.size() == Parameters.size()) 955 Parameters = params; 956 else 957 { 958 app_log() << "Expected " << Parameters.size() << " parameters," 959 << " but found " << params.size() << " in PolynomialFunctor3D.\n"; 960 if (params.size() != 0) 961 abort(); //you think you know what they should be but don't. 962 } 963 // Setup parameter names 964 for (int i = 0; i < Parameters.size(); i++) 965 { 966 std::stringstream sstr; 967 sstr << id << "_" << i; 968 ; 969 if (!notOpt) 970 myVars.insert(sstr.str(), Parameters[i], optimize::LOGLINEAR_P, true); 971 } 972 // for (int i=0; i< N_ee; i++) 973 // for (int j=0; j < N_eI; j++) 974 // for (int k=0; k<=j; k++) { 975 // std::stringstream sstr; 976 // sstr << id << "_" << i << "_" << j << "_" << k; 977 // myVars.insert(sstr.str(),Parameters[index],true); 978 // ParamArray(i,j,k) = ParamArray(i,k,j) = Parameters[index]; 979 // index++; 980 // } 981 if (!notOpt) 982 { 983 int left_pad_spaces = 6; 984 myVars.print(app_log(), left_pad_spaces, true); 985 app_log() << std::endl; 986 } 987 } 988 xmlCoefs = xmlCoefs->next; 989 } 990 reset_gamma(); 991 //print(); 992 return true; 993 } 994 resetParametersPolynomialFunctor3D995 void resetParameters(const opt_variables_type& active) 996 { 997 if (notOpt) 998 return; 999 for (int i = 0; i < Parameters.size(); ++i) 1000 { 1001 int loc = myVars.where(i); 1002 if (loc >= 0) { 1003 Parameters[i] = std::real( myVars[i] = active[loc] ); 1004 } 1005 } 1006 if (ResetCount++ == 100) 1007 { 1008 ResetCount = 0; 1009 //print(); 1010 } 1011 reset_gamma(); 1012 } 1013 checkInVariablesPolynomialFunctor3D1014 void checkInVariables(opt_variables_type& active) { active.insertFrom(myVars); } 1015 checkOutVariablesPolynomialFunctor3D1016 void checkOutVariables(const opt_variables_type& active) { myVars.getIndex(active); } 1017 printPolynomialFunctor3D1018 void print() 1019 { 1020 const int N = 100; 1021 std::string fname = iSpecies + ".J3.h5"; 1022 hid_t hid = H5Fcreate(fname.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); 1023 Array<real_type, 3> val(N, N, N); 1024 for (int i = 0; i < N; i++) 1025 { 1026 double r_12 = (real_type)i / (real_type)(N - 1); 1027 for (int j = 0; j < N; j++) 1028 { 1029 double r_1I = (real_type)j / (real_type)(N - 1) * 0.5 * cutoff_radius; 1030 for (int k = 0; k < N; k++) 1031 { 1032 double r_2I = (real_type)k / (real_type)(N - 1) * 0.5 * cutoff_radius; 1033 double rmin = std::abs(r_1I - r_2I); 1034 double rmax = std::abs(r_1I + r_2I); 1035 double r = rmin + r_12 * (rmax - rmin); 1036 val(i, j, k) = evaluate(r, r_1I, r_2I); 1037 } 1038 } 1039 } 1040 // HDFAttribIO<Array<real_type,3> > coefs_attrib (SplineCoefs); 1041 // HDFAttribIO<Array<real_type,3> > param_attrib (ParamArray); 1042 Array<double, 3> val_DP(N, N, N); 1043 val_DP = val; 1044 HDFAttribIO<Array<double, 3>> val_attrib(val_DP); 1045 val_attrib.write(hid, "val"); 1046 // coefs_attrib.write (hid, "coefs"); 1047 // param_attrib.write (hid, "params"); 1048 H5Fclose(hid); 1049 // std::string fname = (elementType != "") ? elementType : pairType; 1050 // fname = fname + ".dat"; 1051 // //cerr << "Writing " << fname << " file.\n"; 1052 // FILE *fout = fopen (fname.c_str(), "w"); 1053 // for (double r=0.0; r<cutoff_radius; r+=0.001) 1054 // fprintf (fout, "%8.3f %16.10f\n", r, evaluate(r)); 1055 // fclose(fout); 1056 } 1057 1058 printPolynomialFunctor3D1059 void print(std::ostream& os) 1060 { 1061 /* no longer correct. Ye Luo 1062 int n=100; 1063 real_type d=cutoff_radius/100.,r=0; 1064 real_type u,du,d2du; 1065 for(int i=0; i<n; ++i) 1066 { 1067 u=evaluate(r,du,d2du); 1068 os << std::setw(22) << r << std::setw(22) << u << std::setw(22) << du 1069 << std::setw(22) << d2du << std::endl; 1070 r+=d; 1071 } 1072 */ 1073 } 1074 getNumParametersPolynomialFunctor3D1075 inline int getNumParameters() { return Parameters.size(); } 1076 }; 1077 } // namespace qmcplusplus 1078 #endif 1079