00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "QMCThreeBodyJastrow.h"
00014 #include "MathFunctions.h"
00015
00016 QMCThreeBodyJastrow::QMCThreeBodyJastrow()
00017 {
00018 }
00019
00020 QMCThreeBodyJastrow::~QMCThreeBodyJastrow()
00021 {
00022 }
00023
00024 void QMCThreeBodyJastrow::initialize(QMCInput * input)
00025 {
00026 }
00027
00028 double QMCThreeBodyJastrow::get_p3_xxa_ln(int ai)
00029 {
00030 return p3_xxa(ai);
00031 }
00032
00033 Array2D<double> * QMCThreeBodyJastrow::get_p2_xa_ln(int ai)
00034 {
00035 return &p2_xa(ai);
00036 }
00037
00038 double QMCThreeBodyJastrow::get_p_a_ln(int ai)
00039 {
00040 return p_a(ai);
00041 }
00042
00043 void QMCThreeBodyJastrow::evaluate(QMCJastrowParameters & JP,
00044 QMCWalkerData * wData, Array2D<double> & X)
00045 {
00046 wd = wData;
00047 if(wd->whichE == -1)
00048 updateAll(JP,X);
00049 else
00050 updateOne(JP,X);
00051 wd = 0;
00052 }
00053
00054 void QMCThreeBodyJastrow::updateOne(QMCJastrowParameters & JP, Array2D<double> & X)
00055 {
00056 int nalpha = globalInput.WF.getNumberElectrons(true);
00057 int nbeta = globalInput.WF.getNumberElectrons(false);
00058 int E = wd->whichE;
00059 bool isAlpha = E < nalpha ? true : false;
00060
00061 for(int E2=0; E2<wd->UijA.dim1(); E2++)
00062 {
00063 if(E == E2) continue;
00064 int EH = max(E,E2);
00065 int EL = min(E,E2);
00066
00067 wd->U -= wd->UijA(EH,EL);
00068 wd->UijA(EH,EL) = 0.0;
00069
00070 wd->U_xx -= wd->UijA_xx(EH,EL);
00071 wd->UijA_xx(EH,EL) = 0.0;
00072
00073 for(int xyz=0; xyz<3; xyz++)
00074 {
00075 wd->U_x(EH,xyz) -= (wd->UijA_x1(xyz))(EH,EL);
00076 (wd->UijA_x1(xyz))(EH,EL) = 0.0;
00077 wd->U_x(EL,xyz) -= (wd->UijA_x2(xyz))(EH,EL);
00078 (wd->UijA_x2(xyz))(EH,EL) = 0.0;
00079 }
00080 }
00081
00082 Array1D<string> * NucleiTypes = JP.getNucleiTypes();
00083
00084 EupEdnNuclear = JP.getElectronUpElectronDownNuclearParameters();
00085 for(int nuc=0; nuc<EupEdnNuclear->dim1(); nuc++)
00086 (*EupEdnNuclear)(nuc).zeroOutDerivatives();
00087
00088 EupEupNuclear = JP.getElectronUpElectronUpNuclearParameters();
00089 for(int nuc=0; nuc<EupEupNuclear->dim1(); nuc++)
00090 (*EupEupNuclear)(nuc).zeroOutDerivatives();
00091
00092 EdnEdnNuclear = JP.getElectronDownElectronDownNuclearParameters();
00093 for(int nuc=0; nuc<EdnEdnNuclear->dim1(); nuc++)
00094 (*EdnEdnNuclear)(nuc).zeroOutDerivatives();
00095
00096 Array1D<double> xyz(3);
00097 for(int Nuclei=0; Nuclei<globalInput.Molecule.getNumberAtoms(); Nuclei++)
00098 {
00099
00100 int NucleiType = -1;
00101 for( int i=0; i<NucleiTypes->dim1(); i++ )
00102 if( globalInput.Molecule.Atom_Labels(Nuclei) == (*NucleiTypes)(i) )
00103 {
00104 NucleiType = i;
00105 break;
00106 }
00107 bool used;
00108
00109
00110 if (nalpha > 0 && nbeta > 0)
00111 if(isAlpha)
00112 {
00113 for (int EB=nalpha; EB<X.dim1(); EB++)
00114 {
00115 for(int i=0; i<3; i++) xyz(i) = wd->riA_uvec(EB,Nuclei,i);
00116 used = (*EupEdnNuclear)(NucleiType).getThreeBodyCorrelationFunction()->setElectron(true, xyz, wd->riA(EB,Nuclei));
00117 if(!used) continue;
00118 collectForPair(EB, E, Nuclei,& (*EupEdnNuclear)(NucleiType));
00119 }
00120 } else {
00121 for (int EA=0; EA<nalpha; EA++)
00122 {
00123 for(int i=0; i<3; i++) xyz(i) = wd->riA_uvec(E,Nuclei,i);
00124 used = (*EupEdnNuclear)(NucleiType).getThreeBodyCorrelationFunction()->setElectron(true, xyz, wd->riA(E,Nuclei));
00125 if(!used) continue;
00126 collectForPair(E, EA, Nuclei, & (*EupEdnNuclear)(NucleiType));
00127 }
00128 }
00129
00130
00131
00132 if (nalpha > 1)
00133 if(isAlpha)
00134 {
00135 for(int EA=0; EA<E; EA++)
00136 {
00137 for(int i=0; i<3; i++) xyz(i) = wd->riA_uvec(E,Nuclei,i);
00138 used = (*EupEupNuclear)(NucleiType).getThreeBodyCorrelationFunction()->setElectron(true, xyz, wd->riA(E,Nuclei));
00139 if(!used) continue;
00140 collectForPair(E, EA, Nuclei, & (*EupEupNuclear)(NucleiType));
00141 }
00142 for(int EA=E+1; EA<nalpha; EA++)
00143 {
00144 for(int i=0; i<3; i++) xyz(i) = wd->riA_uvec(EA,Nuclei,i);
00145 used = (*EupEupNuclear)(NucleiType).getThreeBodyCorrelationFunction()->setElectron(true, xyz, wd->riA(EA,Nuclei));
00146 if(!used) continue;
00147 collectForPair(EA, E, Nuclei, & (*EupEupNuclear)(NucleiType));
00148 }
00149 }
00150
00151
00152 if (nbeta > 1)
00153 if(!isAlpha)
00154 {
00155 for(int EB=nalpha; EB<E; EB++)
00156 {
00157 for(int i=0; i<3; i++) xyz(i) = wd->riA_uvec(E,Nuclei,i);
00158 used = (*EdnEdnNuclear)(NucleiType).getThreeBodyCorrelationFunction()->setElectron(true, xyz, wd->riA(E,Nuclei));
00159 if(!used) continue;
00160 collectForPair(E, EB, Nuclei, & (*EdnEdnNuclear)(NucleiType));
00161 }
00162 for(int EB=E+1; EB<X.dim1(); EB++)
00163 {
00164 for(int i=0; i<3; i++) xyz(i) = wd->riA_uvec(EB,Nuclei,i);
00165 used = (*EdnEdnNuclear)(NucleiType).getThreeBodyCorrelationFunction()->setElectron(true, xyz, wd->riA(EB,Nuclei));
00166 if(!used) continue;
00167 collectForPair(EB, E, Nuclei, & (*EdnEdnNuclear)(NucleiType));
00168 }
00169 }
00170 }
00171
00172 for(int E2=0; E2<wd->UijA.dim1(); E2++)
00173 {
00174 if(E == E2) continue;
00175 int EH = max(E,E2);
00176 int EL = min(E,E2);
00177
00178 wd->U += wd->UijA(EH,EL);
00179 wd->U_xx += wd->UijA_xx(EH,EL);
00180 for(int xyz=0; xyz<3; xyz++)
00181 {
00182 wd->U_x(EH,xyz) += (wd->UijA_x1(xyz))(EH,EL);
00183 wd->U_x(EL,xyz) += (wd->UijA_x2(xyz))(EH,EL);
00184 }
00185 }
00186
00187 xyz.deallocate();
00188 }
00189
00190 double QMCThreeBodyJastrow::jastrowOnGrid(QMCJastrowParameters & JP,
00191 int E,
00192 Array2D<double> & R,
00193 Array2D<double> & grid,
00194 Array1D<double> & integrand)
00195 {
00196 int nalpha = globalInput.WF.getNumberElectrons(true);
00197 int nbeta = globalInput.WF.getNumberElectrons(false);
00198 bool isAlpha = E < nalpha ? true : false;
00199
00200 double denom = 0.0;
00201 Array1D<double> sumU(integrand.dim1());
00202 sumU = 0.0;
00203
00204 Array1D<string> * NucleiTypes = JP.getNucleiTypes();
00205
00206 EupEdnNuclear = JP.getElectronUpElectronDownNuclearParameters();
00207 EupEupNuclear = JP.getElectronUpElectronUpNuclearParameters();
00208 EdnEdnNuclear = JP.getElectronDownElectronDownNuclearParameters();
00209
00210 for(int Nuclei=0; Nuclei<globalInput.Molecule.getNumberAtoms(); Nuclei++)
00211 {
00212
00213 int NucleiType = -1;
00214 for( int i=0; i<NucleiTypes->dim1(); i++ )
00215 if( globalInput.Molecule.Atom_Labels(Nuclei) == (*NucleiTypes)(i) )
00216 {
00217 NucleiType = i;
00218 break;
00219 }
00220
00221 QMCThreeBodyCorrelationFunction *U_Function;
00222
00223 double r1 = MathFunctions::rij(R,globalInput.Molecule.Atom_Positions,E,Nuclei);
00224 double r1g[grid.dim1()];
00225 for(int gp=0; gp<grid.dim1(); gp++)
00226 r1g[gp] = MathFunctions::rij(grid,globalInput.Molecule.Atom_Positions,gp,Nuclei);
00227
00228
00229 if((*EupEdnNuclear)(NucleiType).isUsed()){
00230 U_Function = (*EupEdnNuclear)(NucleiType).getThreeBodyCorrelationFunction();
00231 if (nalpha > 0 && nbeta > 0)
00232 if(isAlpha){
00233 for (int EB=nalpha; EB<R.dim1(); EB++){
00234 double r2 = MathFunctions::rij(R,globalInput.Molecule.Atom_Positions,EB,Nuclei);
00235 denom += U_Function->getFunctionValue(MathFunctions::rij(R,EB,E),r1,r2);
00236 for(int gp=0; gp<grid.dim1(); gp++)
00237 sumU(gp) += U_Function->getFunctionValue(MathFunctions::rij(grid,R,gp,EB),r1g[gp],r2);
00238 }
00239 } else {
00240 for (int EA=0; EA<nalpha; EA++){
00241 double r2 = MathFunctions::rij(R,globalInput.Molecule.Atom_Positions,EA,Nuclei);
00242 denom += U_Function->getFunctionValue(MathFunctions::rij(R,EA,E),r1,r2);
00243 for(int gp=0; gp<grid.dim1(); gp++)
00244 sumU(gp) += U_Function->getFunctionValue(MathFunctions::rij(grid,R,gp,EA),r1g[gp],r2);
00245 }
00246 }
00247 }
00248
00249
00250 if((*EupEupNuclear)(NucleiType).isUsed()){
00251 U_Function = (*EupEupNuclear)(NucleiType).getThreeBodyCorrelationFunction();
00252 if (nalpha > 1)
00253 if(isAlpha)
00254 {
00255 for(int EA=0; EA<nalpha; EA++)
00256 {
00257 if(EA == E) continue;
00258 double r2 = MathFunctions::rij(R,globalInput.Molecule.Atom_Positions,EA,Nuclei);
00259 denom += U_Function->getFunctionValue(MathFunctions::rij(R,EA,E),r1,r2);
00260 for(int gp=0; gp<grid.dim1(); gp++)
00261 sumU(gp) += U_Function->getFunctionValue(MathFunctions::rij(grid,R,gp,EA),r1g[gp],r2);
00262 }
00263 }
00264 }
00265
00266
00267 if((*EdnEdnNuclear)(NucleiType).isUsed()){
00268 U_Function = (*EdnEdnNuclear)(NucleiType).getThreeBodyCorrelationFunction();
00269 if (nbeta > 1)
00270 if(!isAlpha)
00271 {
00272 for(int EB=nalpha; EB<R.dim1(); EB++)
00273 {
00274 if(EB == E) continue;
00275 double r2 = MathFunctions::rij(R,globalInput.Molecule.Atom_Positions,EB,Nuclei);
00276 denom += U_Function->getFunctionValue(MathFunctions::rij(R,EB,E),r1,r2);
00277 for(int gp=0; gp<grid.dim1(); gp++)
00278 sumU(gp) += U_Function->getFunctionValue(MathFunctions::rij(grid,R,gp,EB),r1g[gp],r2);
00279 }
00280 }
00281 }
00282 }
00283
00284 for(int gp=0; gp<grid.dim1(); gp++)
00285 integrand(gp) *= exp(sumU(gp));
00286 return exp(denom);
00287 }
00288
00289 void QMCThreeBodyJastrow::updateAll(QMCJastrowParameters & JP, Array2D<double> & X)
00290 {
00291 int nalpha = globalInput.WF.getNumberElectrons(true);
00292 int nbeta = globalInput.WF.getNumberElectrons(false);
00293
00294 for(int E1=0; E1<wd->UijA.dim1(); E1++)
00295 for(int E2=0; E2<E1; E2++)
00296 {
00297 wd->U -= wd->UijA(E1,E2);
00298 wd->U_xx -= wd->UijA_xx(E1,E2);
00299 for(int xyz=0; xyz<3; xyz++)
00300 {
00301 wd->U_x(E1,xyz) -= (wd->UijA_x1(xyz))(E1,E2);
00302 wd->U_x(E2,xyz) -= (wd->UijA_x2(xyz))(E1,E2);
00303 }
00304 }
00305
00306 wd->UijA = 0;
00307 wd->UijA_xx = 0;
00308 for (int i=0; i<3; i++)
00309 {
00310 wd->UijA_x1(i) = 0.0;
00311 wd->UijA_x2(i) = 0.0;
00312 }
00313
00314
00315
00316 Array1D<string> * NucleiTypes = JP.getNucleiTypes();
00317
00318 EupEdnNuclear = JP.getElectronUpElectronDownNuclearParameters();
00319 for(int nuc=0; nuc<EupEdnNuclear->dim1(); nuc++)
00320 (*EupEdnNuclear)(nuc).zeroOutDerivatives();
00321
00322 EupEupNuclear = JP.getElectronUpElectronUpNuclearParameters();
00323 for(int nuc=0; nuc<EupEupNuclear->dim1(); nuc++)
00324 (*EupEupNuclear)(nuc).zeroOutDerivatives();
00325
00326 EdnEdnNuclear = JP.getElectronDownElectronDownNuclearParameters();
00327 for(int nuc=0; nuc<EdnEdnNuclear->dim1(); nuc++)
00328 (*EdnEdnNuclear)(nuc).zeroOutDerivatives();
00329
00330 Array1D<double> xyz(3);
00331 for(int Nuclei=0; Nuclei<globalInput.Molecule.getNumberAtoms(); Nuclei++)
00332 {
00333
00334 int NucleiType = -1;
00335 for( int i=0; i<NucleiTypes->dim1(); i++ )
00336 if( globalInput.Molecule.Atom_Labels(Nuclei) == (*NucleiTypes)(i) )
00337 {
00338 NucleiType = i;
00339 break;
00340 }
00341
00342
00343 if (nalpha > 0 && nbeta > 0)
00344 for (int E2=nalpha; E2<nalpha+nbeta; E2++)
00345 {
00346 bool used;
00347 for(int i=0; i<3; i++)
00348 xyz(i) = wd->riA_uvec(E2,Nuclei,i);
00349 used = (*EupEdnNuclear)(NucleiType).getThreeBodyCorrelationFunction()->setElectron(true, xyz, wd->riA(E2,Nuclei));
00350 if(!used) continue;
00351
00352 for(int E1=0; E1<nalpha; E1++)
00353 collectForPair(E2, E1, Nuclei,
00354 & (*EupEdnNuclear)(NucleiType));
00355 }
00356
00357
00358
00359 if (nalpha > 1)
00360 for(int E1=1; E1<nalpha; E1++)
00361 {
00362 bool used;
00363 for(int i=0; i<3; i++)
00364 xyz(i) = wd->riA_uvec(E1,Nuclei,i);
00365 used = (*EupEupNuclear)(NucleiType).getThreeBodyCorrelationFunction()->setElectron(true, xyz, wd->riA(E1,Nuclei));
00366 if(!used) continue;
00367
00368 for(int E2=0; E2<E1; E2++)
00369 collectForPair(E1, E2, Nuclei,
00370 & (*EupEupNuclear)(NucleiType));
00371 }
00372
00373
00374 if (nbeta > 1)
00375 for(int E1=nalpha+1; E1<X.dim1(); E1++)
00376 {
00377 bool used;
00378 for(int i=0; i<3; i++)
00379 xyz(i) = wd->riA_uvec(E1,Nuclei,i);
00380 used = (*EdnEdnNuclear)(NucleiType).getThreeBodyCorrelationFunction()->setElectron(true, xyz, wd->riA(E1,Nuclei));
00381 if(!used) continue;
00382
00383 for(int E2=nalpha; E2<E1; E2++)
00384 collectForPair(E1, E2, Nuclei,
00385 & (*EdnEdnNuclear)(NucleiType));
00386 }
00387 }
00388
00389 for(int E1=0; E1<wd->UijA.dim1(); E1++)
00390 for(int E2=0; E2<E1; E2++)
00391 {
00392 wd->U += wd->UijA(E1,E2);
00393 wd->U_xx += wd->UijA_xx(E1,E2);
00394 for(int xyz=0; xyz<3; xyz++)
00395 {
00396 wd->U_x(E1,xyz) += (wd->UijA_x1(xyz))(E1,E2);
00397 wd->U_x(E2,xyz) += (wd->UijA_x2(xyz))(E1,E2);
00398 }
00399 }
00400
00401 xyz.deallocate();
00402 packageDerivatives();
00403 }
00404
00405 inline void QMCThreeBodyJastrow::collectForPair(int E1,
00406 int E2,
00407 int Nuclei,
00408 QMCThreeBodyCorrelationFunctionParameters * paramset)
00409 {
00410 if(!paramset->isUsed()) return;
00411
00412 QMCThreeBodyCorrelationFunction *U_Function = paramset->getThreeBodyCorrelationFunction();
00413
00414 Array1D<double> xyz2(3);
00415 Array1D<double> xyz12(3);
00416
00417 for(int i=0; i<3; i++)
00418 {
00419 xyz2(i) = wd->riA_uvec(E2,Nuclei,i);
00420 xyz12(i) = wd->rij_uvec(E1, E2,i);
00421 }
00422
00423 bool used;
00424 used = U_Function->setElectron(false,xyz2, wd->riA(E2,Nuclei));
00425 if(!used) return;
00426
00427 U_Function->evaluate(xyz12, wd->rij(E1,E2));
00428
00429 if(IeeeMath::isNaN( U_Function->getLaplacianValue() ) && false ){
00430 printf("(%2i,%2i) sum %20.10e cur %20.10e psi %20.10e r12 %20.10e r1 %20.10e r2 %20.10e \n",E1,E2,
00431 wd->UijA_xx(E1,E2),
00432 U_Function->getLaplacianValue(),
00433 U_Function->getFunctionValue(),
00434 wd->rij(E1,E2),wd->riA(E1,Nuclei),wd->riA(E2,Nuclei)
00435 );
00436 U_Function->print(cerr);
00437
00438 }
00439
00440 wd->UijA(E1,E2) += U_Function->getFunctionValue();
00441 wd->UijA_xx(E1,E2) += U_Function->getLaplacianValue();
00442
00443 Array1D<double> * grad1 = U_Function->getElectron1Gradient();
00444 Array1D<double> * grad2 = U_Function->getElectron2Gradient();
00445 for (int i=0; i<3; i++)
00446 {
00447 (wd->UijA_x1(i))(E1,E2) += (*grad1)(i);
00448 (wd->UijA_x2(i))(E1,E2) += (*grad2)(i);
00449 }
00450
00451 if(globalInput.flags.calculate_Derivatives == 1 &&
00452 globalInput.flags.optimize_NEE_Jastrows == 1)
00453 {
00454 for(int ai=0; ai<paramset->getNumberOfTotalParameters()+1; ai++)
00455 {
00456 paramset->pt_a(ai) += U_Function->get_p_a(ai);
00457 paramset->pt3_xxa(ai) += U_Function->get_p3_xxa(ai);
00458
00459 for(int i=0; i<3; i++)
00460 {
00461 (paramset->pt2_xa(ai))(E1,i) += U_Function->get_p2_xa(true,i,ai);
00462 (paramset->pt2_xa(ai))(E2,i) += U_Function->get_p2_xa(false,i,ai);
00463 }
00464 }
00465 }
00466 }
00467
00468 void QMCThreeBodyJastrow::packageDerivatives()
00469 {
00470 if(globalInput.flags.calculate_Derivatives != 1 ||
00471 globalInput.flags.optimize_NEE_Jastrows == 0)
00472 return;
00473
00474 int numNEE =
00475 globalInput.JP.getNumberNEupEdnParameters() +
00476 globalInput.JP.getNumberNEupEupParameters() +
00477 globalInput.JP.getNumberNEdnEdnParameters();
00478
00479 p_a.allocate(numNEE);
00480 p2_xa.allocate(numNEE);
00481 p3_xxa.allocate(numNEE);
00482 for(int i=0; i<p2_xa.dim1(); i++)
00483 {
00484 p2_xa(i).allocate(globalInput.WF.getNumberElectrons(),3);
00485 p2_xa(i) = 0.0;
00486 }
00487 p_a = 0.0;
00488 p3_xxa = 0.0;
00489
00490 int numNuc = EupEdnNuclear->dim1();
00491 int ai = 0;
00492 for(int nuc=0; nuc<numNuc; nuc++)
00493 {
00494 QMCThreeBodyCorrelationFunctionParameters * accum = & (*EupEdnNuclear)(nuc);
00495
00496 switch(globalInput.flags.link_NEE_Jastrows)
00497 {
00498 case 0:
00499 accum->totalDerivativesToFree(ai,p_a,p2_xa,p3_xxa);
00500 accum = & (*EupEupNuclear)(nuc);
00501 accum->totalDerivativesToFree(ai,p_a,p2_xa,p3_xxa);
00502 accum = & (*EdnEdnNuclear)(nuc);
00503 accum->totalDerivativesToFree(ai,p_a,p2_xa,p3_xxa);
00504 break;
00505
00506 case 1:
00507 accum->totalDerivativesToFree(ai,p_a,p2_xa,p3_xxa);
00508 accum = & (*EupEupNuclear)(nuc);
00509 accum->totalDerivativesAccumulate((*EdnEdnNuclear)(nuc));
00510 accum->totalDerivativesToFree(ai,p_a,p2_xa,p3_xxa);
00511 break;
00512
00513 case 2:
00514 accum->totalDerivativesAccumulate((*EupEupNuclear)(nuc));
00515 accum->totalDerivativesAccumulate((*EdnEdnNuclear)(nuc));
00516 accum->totalDerivativesToFree(ai,p_a,p2_xa,p3_xxa);
00517 break;
00518 }
00519 }
00520 }
00521
