QMcBeaver: QMCJastrowElectronNuclear.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 "QMCJastrowElectronNuclear.h"
00014 #include "MathFunctions.h"
00015 
00016 QMCJastrowElectronNuclear::QMCJastrowElectronNuclear(){}
00017 
00018 QMCJastrowElectronNuclear::~QMCJastrowElectronNuclear()
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 QMCJastrowElectronNuclear::initialize(QMCInput * input)
00029 {
00030   int numNE = globalInput.JP.getNumberNEParameters();
00031   p_a.allocate(numNE);
00032   p2_xa.allocate(numNE);
00033   p3_xxa.allocate(numNE);
00034 
00035   for(int i=0; i<p2_xa.dim1(); i++)
00036     p2_xa(i).allocate(globalInput.WF.getNumberElectrons(),3);
00037 }
00038 
00039 double QMCJastrowElectronNuclear::get_p3_xxa_ln(int ai)
00040 {
00041   return p3_xxa(ai);
00042 }
00043 
00044 Array2D<double> * QMCJastrowElectronNuclear::get_p2_xa_ln(int ai)
00045 {
00046   return &p2_xa(ai);
00047 }
00048 
00049 double QMCJastrowElectronNuclear::get_p_a_ln(int ai)
00050 {
00051   return p_a(ai);
00052 }
00053 
00054 void QMCJastrowElectronNuclear::evaluate(QMCJastrowParameters & JP,
00055                                          QMCWalkerData * wd,
00056                                          Array2D<double> & X)
00057 {
00058   // initialize the results
00059   double temp;
00060   double UiA, UiA_x, UiA_xx;
00061 
00062   p_a              = 0.0;
00063   p3_xxa           = 0.0;
00064   for(int ai=0; ai<p2_xa.dim1(); ai++)
00065     p2_xa(ai)      = 0.0;
00066 
00067   int nalpha = globalInput.WF.getNumberElectrons(true);
00068   int nbeta  = globalInput.WF.getNumberElectrons(false);
00069 
00070   // Get values from JP that will be needed during the calc
00071 
00072   Array1D<string> * NucleiTypes = JP.getNucleiTypes();
00073 
00074   Array1D<QMCCorrelationFunctionParameters> * EupNuclear = 0;
00075   if (nalpha > 0)
00076     EupNuclear = JP.getElectronUpNuclearParameters();
00077   
00078   Array1D<QMCCorrelationFunctionParameters> * EdnNuclear = 0;
00079   if (nbeta > 0)
00080     EdnNuclear = JP.getElectronDownNuclearParameters();
00081 
00082   // Loop over each atom calculating the e-n jastrow function
00083   double r;
00084   for(int Nuclei=0; Nuclei<globalInput.Molecule.getNumberAtoms(); Nuclei++)
00085     {
00086 
00087       /*
00088         For labeling the ai parameters, we need to make sure that the order
00089         we place the ai derivatives in our vectors lines up with how the ai
00090         are placed in the parameter vector in QMCJastrowParameters.
00091 
00092         In particular, the loops are inverted relative to what is done here.
00093         The parameters need to have the nuclei indices as the "fast loop".
00094         Eup comes first, then if the Jastrows are not linked, we have the Edn
00095         parameters.
00096 
00097         nuc_E??_shift will bring us to the right nucleus
00098         ai_shift will bring us to the right Jastrow.
00099         ai indexes the parameter within the Jastrow.
00100       */
00101       int nuc_Eup_shift = 0;
00102       int nuc_Edn_shift = 0;
00103       int numP;
00104 
00105       // Find the number of the current nucleus in the nuclei list
00106       int CurrentNucleiNumber = -1;
00107       for( int i=0; i<NucleiTypes->dim1(); i++ )
00108         if( globalInput.Molecule.Atom_Labels(Nuclei) == (*NucleiTypes)(i) )
00109           {
00110             CurrentNucleiNumber = i;
00111             break;
00112           } else {
00113             nuc_Eup_shift += (*EupNuclear)(i).getTotalNumberOfParameters();
00114             nuc_Edn_shift += (*EdnNuclear)(i).getTotalNumberOfParameters();
00115           }
00116 
00117       for(int Electron=0; Electron<X.dim1(); Electron++)
00118         {
00119           if(wd->whichE != Electron && wd->whichE != -1)
00120             continue;
00121 
00122           r = wd->riA(Electron,Nuclei);
00123 
00124           // Get the correct correlation function to use and evaluate it
00125           QMCCorrelationFunction *U_Function = 0;
00126 
00127           int ai_shift;
00128           if( Electron < nalpha )
00129             {
00130               U_Function = 
00131                 (*EupNuclear)(CurrentNucleiNumber).getCorrelationFunction();
00132 
00133               ai_shift = nuc_Eup_shift;
00134               numP = (*EupNuclear)(CurrentNucleiNumber).getTotalNumberOfParameters();
00135             }
00136           else
00137             {
00138               U_Function =
00139                 (*EdnNuclear)(CurrentNucleiNumber).getCorrelationFunction();
00140 
00141               ai_shift = nuc_Edn_shift;
00142               numP = (*EdnNuclear)(CurrentNucleiNumber).getTotalNumberOfParameters();
00143               if(globalInput.flags.link_Jastrow_parameters == 0)
00144                 ai_shift += globalInput.JP.getNumberNEupParameters();
00145             }
00146 
00147           U_Function->evaluate(r);
00148 
00149           UiA    = U_Function->getFunctionValue();
00150           UiA_x  = U_Function->getFirstDerivativeValue();
00151           UiA_xx = 2.0/r * UiA_x + U_Function->getSecondDerivativeValue();
00152 
00153           // Update the values being calculated ...
00154 
00155           wd->U -= wd->UiA(Electron, Nuclei);
00156           wd->U += UiA;
00157           wd->UiA(Electron, Nuclei) = UiA;
00158 
00159           wd->U_xx -= wd->UiA_xx(Electron,Nuclei);
00160           wd->U_xx += UiA_xx;
00161           wd->UiA_xx(Electron,Nuclei) = UiA_xx;
00162 
00163           for(int i=0; i<3; i++)
00164             {
00165               temp = wd->UiA_x(Electron,Nuclei,i);
00166               wd->U_x(Electron,i) -= temp;
00167               temp = UiA_x * wd->riA_uvec(Electron,Nuclei,i);
00168               wd->UiA_x(Electron,Nuclei,i) = temp;
00169               wd->U_x(Electron,i) += temp;
00170             }
00171 
00172           if(globalInput.flags.calculate_Derivatives == 1 &&
00173              globalInput.flags.optimize_EN_Jastrows == 1)
00174             {
00175               for(int ai=0; ai<numP; ai++)
00176                 {
00177                   p_a(ai+ai_shift)    += U_Function->get_p_a(ai);
00178                   UiA_x           = U_Function->get_p2_xa(ai);
00179                   p3_xxa(ai+ai_shift) += 2.0/r * UiA_x + U_Function->get_p3_xxa(ai);
00180                   
00181                   for(int i=0; i<3; i++)
00182                     (p2_xa(ai+ai_shift))(Electron,i) += UiA_x * wd->riA_uvec(Electron,Nuclei,i);
00183                 }
00184             }
00185         }
00186     }
00187 }
00188 
00189 double QMCJastrowElectronNuclear::jastrowOnGrid(QMCJastrowParameters & JP,
00190                                                 int Electron,
00191                                                 Array2D<double> & R,
00192                                                 Array2D<double> & grid,
00193                                                 Array1D<double> & integrand)
00194 {
00195   // initialize the results
00196   double denom = 0.0;
00197   double r;
00198   Array1D<double> sumU(integrand.dim1());
00199   sumU = 0.0;
00200 
00201   int nalpha = globalInput.WF.getNumberElectrons(true);
00202   int nbeta  = globalInput.WF.getNumberElectrons(false);
00203 
00204   // Get values from JP that will be needed during the calc
00205 
00206   Array1D<string> * NucleiTypes = JP.getNucleiTypes();
00207 
00208   Array1D<QMCCorrelationFunctionParameters> * EupNuclear = 0;
00209   if (nalpha > 0)
00210     EupNuclear = JP.getElectronUpNuclearParameters();
00211   
00212   Array1D<QMCCorrelationFunctionParameters> * EdnNuclear = 0;
00213   if (nbeta > 0)
00214     EdnNuclear = JP.getElectronDownNuclearParameters();
00215 
00216   // Loop over each atom calculating the e-n jastrow function
00217   for(int Nuclei=0; Nuclei<globalInput.Molecule.getNumberAtoms(); Nuclei++)
00218     {
00219       // Find the number of the current nucleus in the nuclei list
00220       int CurrentNucleiNumber = -1;
00221       for( int i=0; i<NucleiTypes->dim1(); i++ )
00222         if( globalInput.Molecule.Atom_Labels(Nuclei) == (*NucleiTypes)(i) )
00223           {
00224             CurrentNucleiNumber = i;
00225             break;
00226           }      
00227 
00228       // Get the correct correlation function to use and evaluate it
00229       QMCCorrelationFunction *U_Function = 0;
00230       
00231       if( Electron < nalpha )
00232         {
00233           U_Function = 
00234             (*EupNuclear)(CurrentNucleiNumber).getCorrelationFunction();
00235         }
00236       else
00237         {
00238           U_Function =
00239             (*EdnNuclear)(CurrentNucleiNumber).getCorrelationFunction();
00240         }
00241       
00242       r = MathFunctions::rij(R,globalInput.Molecule.Atom_Positions,
00243                              Electron,Nuclei);
00244       denom += U_Function->getFunctionValue(r);
00245       
00246       for(int gr=0; gr<grid.dim1(); gr++){
00247         r = MathFunctions::rij(grid,globalInput.Molecule.Atom_Positions,
00248                                gr,Nuclei);
00249         sumU(gr) += U_Function->getFunctionValue(r);
00250       }
00251     }
00252 
00253   for(int gr=0; gr<grid.dim1(); gr++)
00254     integrand(gr) *= exp(sumU(gr));
00255   return exp(denom);
00256 }
00257