QMcBeaver: QMCJastrowElectronElectron.cpp Source File

00001 //            QMcBeaver
00002 //
00003 //         Constructed by
00004 //
00005 //     Michael Todd Feldmann
00006 //              and
00007 //   David Randall "Chip" Kent IV
00008 //
00009 // Copyright 2000.  All rights reserved.
00010 //
00011 // drkent@users.sourceforge.net mtfeldmann@users.sourceforge.net
00012 
00013 #include "QMCJastrowElectronElectron.h"
00014 #include "MathFunctions.h"
00015 
00016 QMCJastrowElectronElectron::QMCJastrowElectronElectron(){}
00017 
00018 QMCJastrowElectronElectron::~QMCJastrowElectronElectron()
00019 {
00020   for(int i=0; i<p2_xa.dim1(); i++)
00021     p2_xa(i).deallocate();
00022 
00023   p_a.deallocate();
00024   p2_xa.deallocate();
00025   p3_xxa.deallocate();
00026 }
00027 
00028 void QMCJastrowElectronElectron::initialize(QMCInput * input)
00029 {
00030   wd = 0;
00031 
00032   int numEE = globalInput.JP.getNumberEEParameters();
00033   p_a.allocate(numEE);
00034   p2_xa.allocate(numEE);
00035   p3_xxa.allocate(numEE);
00036 
00037   for(int i=0; i<p2_xa.dim1(); i++)
00038     p2_xa(i).allocate(globalInput.WF.getNumberElectrons(),3);
00039 }
00040 
00041 double QMCJastrowElectronElectron::get_p3_xxa_ln(int ai)
00042 {
00043   return p3_xxa(ai);
00044 }
00045 
00046 Array2D<double> * QMCJastrowElectronElectron::get_p2_xa_ln(int ai)
00047 {
00048   return &p2_xa(ai);
00049 }
00050 
00051 double QMCJastrowElectronElectron::get_p_a_ln(int ai)
00052 {
00053   return p_a(ai);
00054 }
00055 
00056 void QMCJastrowElectronElectron::evaluate(QMCJastrowParameters & JP,
00057                                           QMCWalkerData * wData,
00058                                           Array2D<double> & X)
00059 {
00060   wd = wData;
00061   if(wd->whichE != -1)
00062     updateOne(JP,X);
00063   else
00064     updateAll(JP,X);
00065   wd = 0;
00066 }
00067 
00068 void QMCJastrowElectronElectron::updateOne(QMCJastrowParameters & JP,
00069                                            Array2D<double> & X)
00070 {
00071   int E = wd->whichE;
00072 
00073   int nalpha = globalInput.WF.getNumberElectrons(true);
00074   int nbeta  = globalInput.WF.getNumberElectrons(false);
00075   bool isAlpha = E < nalpha ? true : false;
00076   
00077   // Get values from JP that will be needed during the calc
00078 
00079   QMCCorrelationFunctionParameters * EupEdn = 0;
00080 
00081   if(nalpha > 0 && nbeta > 0)
00082     EupEdn = JP.getElectronUpElectronDownParameters();
00083   
00084   QMCCorrelationFunctionParameters * EupEup = 0;
00085   
00086   if(nalpha > 1)
00087     EupEup = JP.getElectronUpElectronUpParameters();
00088   
00089   QMCCorrelationFunctionParameters * EdnEdn = 0;
00090   
00091   if(nbeta > 1 )
00092     EdnEdn = JP.getElectronDownElectronDownParameters();
00093 
00094   QMCCorrelationFunction *U_Function = 0;
00095 
00096   int index = 0;
00097   int numP  = 0;
00098 
00099   if(EupEdn != 0)
00100     {
00101       U_Function = EupEdn->getCorrelationFunction();
00102       numP = EupEdn->getTotalNumberOfParameters();
00103 
00104       if(isAlpha)
00105         {
00106           for(int EB=nalpha; EB<X.dim1(); EB++)
00107             collectForPair(EB,E,U_Function,X,index,numP);
00108         } else {
00109           for(int EA=0; EA<nalpha; EA++)
00110             collectForPair(E,EA,U_Function,X,index,numP);
00111         }
00112     }
00113   
00114   index += numP;
00115   
00116   if(EupEup != 0)
00117     {
00118       numP  = EupEup->getTotalNumberOfParameters();
00119       U_Function = EupEup->getCorrelationFunction();
00120 
00121       if(isAlpha)
00122         {
00123           for(int EA=0; EA<E; EA++)
00124             collectForPair(E,EA,U_Function,X,index,numP);
00125           for(int EA=E+1; EA<nalpha; EA++)
00126             collectForPair(EA,E,U_Function,X,index,numP);
00127         }
00128     } else {
00129       numP = 0;
00130     }
00131 
00132   if(globalInput.flags.link_Jastrow_parameters == 0)
00133     index += numP;
00134 
00135   if(EdnEdn != 0)
00136     {
00137       numP  = EdnEdn->getTotalNumberOfParameters();
00138       U_Function = EdnEdn->getCorrelationFunction();
00139 
00140       if(!isAlpha)
00141         {
00142           for(int EB=nalpha; EB<E; EB++)
00143             collectForPair(E,EB,U_Function,X,index,numP);      
00144           for(int EB=E+1; EB<X.dim1(); EB++)
00145             collectForPair(EB,E,U_Function,X,index,numP);
00146         }
00147     } else {
00148       numP  = 0;
00149     }
00150 }
00151 
00152 void QMCJastrowElectronElectron::updateAll(QMCJastrowParameters & JP,
00153                                            Array2D<double> & X)
00154 {
00155   // initialize the results
00156   p_a              = 0.0;
00157   p3_xxa           = 0.0;
00158   for(int ai=0; ai<p2_xa.dim1(); ai++)
00159     p2_xa(ai)      = 0.0;
00160 
00161   int nalpha = globalInput.WF.getNumberElectrons(true);
00162   int nbeta  = globalInput.WF.getNumberElectrons(false);
00163 
00164   // Get values from JP that will be needed during the calc
00165 
00166   QMCCorrelationFunctionParameters * EupEdn = 0;
00167 
00168   if(nalpha > 0 && nbeta > 0)
00169     EupEdn = JP.getElectronUpElectronDownParameters();
00170   
00171   QMCCorrelationFunctionParameters * EupEup = 0;
00172   
00173   if(nalpha > 1)
00174     EupEup = JP.getElectronUpElectronUpParameters();
00175   
00176   QMCCorrelationFunctionParameters * EdnEdn = 0;
00177   
00178   if(nbeta > 1 )
00179     EdnEdn = JP.getElectronDownElectronDownParameters();
00180 
00181   QMCCorrelationFunction *U_Function = 0;
00182 
00183   int index = 0;
00184   int numP  = 0;
00185 
00186   if(EupEdn != 0)
00187     {
00188       U_Function = EupEdn->getCorrelationFunction();
00189       numP = EupEdn->getTotalNumberOfParameters();
00190 
00191       for(int E1=0; E1<nalpha; E1++)
00192         for(int E2=nalpha; E2<X.dim1(); E2++)
00193           collectForPair(E2,E1,U_Function,X,index,numP);
00194     }
00195   
00196   index += numP;
00197   
00198   if(EupEup != 0)
00199     {
00200       numP  = EupEup->getTotalNumberOfParameters();
00201       U_Function = EupEup->getCorrelationFunction();
00202 
00203       for(int E1=0; E1<nalpha; E1++)
00204         for(int E2=0; E2<E1; E2++)
00205           collectForPair(E1,E2,U_Function,X,index,numP);
00206     } else {
00207       numP = 0;
00208     }
00209   
00210   if(globalInput.flags.link_Jastrow_parameters == 0)
00211     index += numP;
00212 
00213   if(EdnEdn != 0)
00214     {
00215       numP  = EdnEdn->getTotalNumberOfParameters();
00216       U_Function = EdnEdn->getCorrelationFunction();
00217 
00218       for(int E1=nalpha; E1<X.dim1(); E1++)
00219         for(int E2=nalpha; E2<E1; E2++)
00220           collectForPair(E1,E2,U_Function,X,index,numP);
00221     } else {
00222       numP  = 0;
00223     }
00224 }
00225 
00226 inline void QMCJastrowElectronElectron::collectForPair(int E1, 
00227                                                        int E2,
00228                                                        QMCCorrelationFunction *U_Function,
00229                                                        Array2D<double> & X,
00230                                                        int index, int numP)
00231 {
00232   double temp;
00233   double Uij, Uij_x, Uij_xx;
00234   double r = wd->rij(E1, E2);
00235 
00236   U_Function->evaluate(r);      
00237   Uij    = U_Function->getFunctionValue();
00238   Uij_x  = U_Function->getFirstDerivativeValue();
00239   Uij_xx = 2.0*(2.0/r * Uij_x + U_Function->getSecondDerivativeValue());
00240   
00241   //Subtract the old value
00242   wd->U -= wd->Uij(E1, E2);
00243   //Add the new value
00244   wd->U += Uij;
00245   //Save the new value
00246   wd->Uij(E1, E2)    = Uij;
00247   
00248   for(int i=0; i<3; i++)
00249     {
00250       temp = wd->Uij_x(E1, E2, i);
00251       wd->U_x(E1,i) -= temp;
00252       wd->U_x(E2,i) += temp;
00253       
00254       temp = Uij_x * wd->rij_uvec(E1, E2,i);
00255       wd->Uij_x(E1, E2, i) = temp;
00256       
00257       wd->U_x(E1,i) += temp;
00258       wd->U_x(E2,i) -= temp;
00259     }
00260   
00261   wd->U_xx -= wd->Uij_xx(E1, E2);
00262   wd->U_xx += Uij_xx;
00263   wd->Uij_xx(E1, E2) = Uij_xx;
00264 
00265   /*
00266     These are for calculating the derivative of the energy
00267     with respect to parameter ai.
00268   */
00269   if(globalInput.flags.calculate_Derivatives == 1 &&
00270      globalInput.flags.optimize_EE_Jastrows == 1)
00271     {
00272       double firstDeriv;
00273 
00274       for(int ai=0; ai<numP; ai++)
00275         {         
00276           p_a(ai+index)    += U_Function->get_p_a(ai);
00277           firstDeriv        = U_Function->get_p2_xa(ai);
00278           p3_xxa(ai+index) += 2.0*(2.0/r * firstDeriv +
00279                                    U_Function->get_p3_xxa(ai));
00280           
00281           for(int i=0; i<3; i++)
00282             {
00283               temp = firstDeriv * wd->rij_uvec(E1, E2,i);
00284               
00285               (p2_xa(ai+index))(E1,i) += temp;
00286               (p2_xa(ai+index))(E2,i) -= temp;
00287             }
00288         }
00289     }
00290 }
00291 
00292 double QMCJastrowElectronElectron::jastrowOnGrid(QMCJastrowParameters & JP,
00293                                                  int E,
00294                                                  Array2D<double> & R,
00295                                                  Array2D<double> & grid,
00296                                                  Array1D<double> & integrand)
00297 {
00298   int nalpha = globalInput.WF.getNumberElectrons(true);
00299   int nbeta  = globalInput.WF.getNumberElectrons(false);
00300   double denom = 0.0;
00301 
00302   bool isAlpha = E < nalpha ? true : false;
00303 
00304   Array1D<double> sumU(integrand.dim1());
00305   sumU = 0.0;
00306 
00307   if(integrand.dim1() != grid.dim1()){
00308     cout << "Integrand dimensions don't match grid!\n";
00309   }
00310 
00311   // Get values from JP that will be needed during the calc
00312 
00313   QMCCorrelationFunctionParameters * EupEdn = 0;
00314 
00315   if(nalpha > 0 && nbeta > 0)
00316     EupEdn = JP.getElectronUpElectronDownParameters();
00317   
00318   QMCCorrelationFunctionParameters * EupEup = 0;
00319   
00320   if(nalpha > 1)
00321     EupEup = JP.getElectronUpElectronUpParameters();
00322   
00323   QMCCorrelationFunctionParameters * EdnEdn = 0;
00324   
00325   if(nbeta > 1 )
00326     EdnEdn = JP.getElectronDownElectronDownParameters();
00327 
00328   QMCCorrelationFunction *U_Function = 0;
00329 
00330   if(EupEdn != 0)
00331     {
00332       U_Function = EupEdn->getCorrelationFunction();
00333 
00334       if(isAlpha)
00335         {
00336           for(int EB=nalpha; EB<R.dim1(); EB++)
00337             {
00338               denom += U_Function->getFunctionValue(MathFunctions::rij(R,EB,E));
00339               for(int gr=0; gr<grid.dim1(); gr++)               
00340                 sumU(gr) += U_Function->getFunctionValue(MathFunctions::rij(grid,R,gr,EB));
00341             }
00342         } else {
00343           for(int EA=0; EA<nalpha; EA++)
00344             {
00345               denom += U_Function->getFunctionValue(MathFunctions::rij(R,E,EA));
00346               for(int gr=0; gr<grid.dim1(); gr++)
00347                 sumU(gr) += U_Function->getFunctionValue(MathFunctions::rij(grid,R,gr,EA));
00348             }
00349         }
00350     }  
00351   
00352   if(EupEup != 0)
00353     {
00354       U_Function = EupEup->getCorrelationFunction();
00355 
00356       if(isAlpha)
00357         {
00358           for(int EA=0; EA<nalpha; EA++)
00359             {
00360               if(EA == E) continue;
00361               denom += U_Function->getFunctionValue(MathFunctions::rij(R,E,EA));          
00362               for(int gr=0; gr<grid.dim1(); gr++)
00363                 sumU(gr) += U_Function->getFunctionValue(MathFunctions::rij(grid,R,gr,EA));
00364             }
00365         }
00366     }
00367 
00368   if(EdnEdn != 0)
00369     {
00370       U_Function = EdnEdn->getCorrelationFunction();
00371 
00372       if(!isAlpha)
00373         {
00374           for(int EB=nalpha; EB<R.dim1(); EB++)
00375             {
00376               if(EB == E) continue;
00377               denom += U_Function->getFunctionValue(MathFunctions::rij(R,E,EB));
00378               for(int gr=0; gr<grid.dim1(); gr++)
00379                 sumU(gr) += U_Function->getFunctionValue(MathFunctions::rij(grid,R,gr,EB));         
00380             }            
00381         }
00382     }
00383 
00384   for(int gr=0; gr<grid.dim1(); gr++)
00385     integrand(gr) *= exp(sumU(gr));
00386   return exp(denom);
00387 }