wavefunction.cpp

Go to the documentation of this file.

 
#include "wavefunction.h"
#include "path.h"
 
// ---------------------------------------------------------------------------
// ---------------------------------------------------------------------------
// WAVEFUNCTION BASE CLASS ---------------------------------------------------
// ---------------------------------------------------------------------------
// ---------------------------------------------------------------------------
 
/**************************************************************************/
WaveFunctionBase::WaveFunctionBase (const Path &_path, LookupTable &_lookup,
        string _name) :
    name(_name),
    path(_path),
    lookup(_lookup)
{
    // empty constructor
}
 
/**************************************************************************/
WaveFunctionBase::~WaveFunctionBase() {
    // empty destructor
}
 
// ---------------------------------------------------------------------------
// ---------------------------------------------------------------------------
// SECH WAVEFUNCTION CLASS ---------------------------------------------------
// ---------------------------------------------------------------------------
// ---------------------------------------------------------------------------
 
/**************************************************************************/
SechWaveFunction::SechWaveFunction(const Path &_path,LookupTable &_lookup, string _name) :
    WaveFunctionBase(_path,_lookup,_name)
{
    /* Set the parameter to its optimized value */
    a = sqrt(0.5*M_PI);
}
 
/**************************************************************************/
SechWaveFunction::~SechWaveFunction() {
    // empty destructor
}
 
/**************************************************************************/
double SechWaveFunction::PsiTrial(const int slice) {
 
    /* The cumulative value */
    double psiT = 1.0;
    beadLocator beadIndex;
    dVec pos;
    double r;
 
    for (int ptcl = 0; ptcl < path.numBeadsAtSlice(slice); ptcl++) {
        beadIndex = slice,ptcl;
        pos = path(beadIndex);
        r = sqrt(dot(pos,pos));
        psiT *= 1.0/cosh(a*r);
    }
 
    return psiT;
}
 
 
// ---------------------------------------------------------------------------
// ---------------------------------------------------------------------------
// JASTROW WAVEFUNCTION CLASS ---------------------------------------------------
// ---------------------------------------------------------------------------
// ---------------------------------------------------------------------------
 
/**************************************************************************/
JastrowWaveFunction::JastrowWaveFunction(const Path &_path,LookupTable &_lookup, string _name) :
WaveFunctionBase(_path,_lookup,_name)
{
    /* Set the parameter to its optimized value */
    alpha = 19.0;
    beta = 0.12;
    //beta = 3.07;
}
 
/**************************************************************************/
JastrowWaveFunction::~JastrowWaveFunction() {
    // empty destructor
}
 
/**************************************************************************/
double JastrowWaveFunction::PsiTrial(const double r) {
    
    double psiT = exp( -0.5*alpha/(1.0+beta*pow(r,5.0)));
    return psiT;
}
 
/**************************************************************************/
double JastrowWaveFunction::delPsiTrial(const double r) {
    
    double delpsiT = 2.5*alpha*beta*pow(r,4.0)/pow((1.0+beta*pow(r,5.0)),2.0);
    return delpsiT;
}
 
/**************************************************************************/
double JastrowWaveFunction::delSqPsiTrial(const double r) {
    
    double delSqPsiT = -1.25*alpha*beta*pow(r,3.0)*(-8.0+
                    beta*(4.0-5.0*alpha)*pow(r,5.0)+12.0*pow(beta,2.0)*pow(r,10.0) )/
                        ( pow((1.0+beta*pow(r,5.0)),4.0) );
    return delSqPsiT;
}
 
/**************************************************************************/
double JastrowWaveFunction::PsiTrial(const int slice) {
    
    /* The cumulative value */
    double psiT = 1.0;
    int numParticles = path.numBeadsAtSlice(slice);    
    dVec sep;                       // The spatial separation between beads.
    double r;                       // Distance between beads
    beadLocator bead1,bead2;
    bead1[0] = bead2[0] = slice;
        
    for (bead1[1] = 0; bead1[1] < numParticles; bead1[1]++) {
        /* The loop over all other particles, to find the total interaction
         * potential */
        for (bead2[1] = bead1[1]+1; bead2[1] < numParticles; bead2[1]++) {
            sep = path.getSeparation(bead2,bead1);
            r = sqrt(dot(sep,sep));
            psiT *= PsiTrial(r);
        } // bead2
        
    } // bead1
    
    return psiT;
}
 
/**************************************************************************/
double JastrowWaveFunction::gradSqPsiTrial(const int slice) {
    
    /* The cumulative value */
    double gradSqPsiT = 1.0;
    double delPsi12;
    int numParticles = path.numBeadsAtSlice(slice);
    dVec sep;                       // The spatial separation between beads.
    double r;                       // Distance between beads
    beadLocator bead1,bead2,bead3;
    bead1[0] = bead2[0] = bead3[0] = slice;
    
    for (bead1[1] = 0; bead1[1] < numParticles; bead1[1]++) {
        /* The loop over all other particles, to find the total interaction
         * potential */
        for (bead2[1] = bead1[1]+1; bead2[1] < numParticles; bead2[1]++) {
            sep = path.getSeparation(bead2,bead1);
            r = sqrt(dot(sep,sep));
            gradSqPsiT += 2.0*delSqPsiTrial(r);
            delPsi12 = delPsiTrial(r);
            for (bead3[1] = bead2[1]+1; bead3[1] < numParticles; bead3[1]++) {
                sep = path.getSeparation(bead3,bead1);
                r = sqrt(dot(sep,sep));
                gradSqPsiT += 2.0*delPsi12*delPsiTrial(r);
            }// bead 3
        } // bead2
        
    } // bead1
    
    return gradSqPsiT;
}
 
 
// ---------------------------------------------------------------------------
// ---------------------------------------------------------------------------
// LiebLiniger WAVEFUNCTION CLASS --------------------------------------------
// ---------------------------------------------------------------------------
// ---------------------------------------------------------------------------
 
/**************************************************************************/
LiebLinigerWaveFunction::LiebLinigerWaveFunction(const Path &_path,LookupTable &_lookup, string _name) :
WaveFunctionBase(_path,_lookup,_name)
{
    /* Set the parameter to its optimized value */
    R = constants()->R_LL_wfn();
    k = constants()->k_LL_wfn();
}
 
/**************************************************************************/
LiebLinigerWaveFunction::~LiebLinigerWaveFunction() {
    // empty destructor
}
 
/**************************************************************************/
double LiebLinigerWaveFunction::PsiTrial(const double r) {
    
    double psiT = 1.0;
    if(r<R)
        psiT = cos(k*(abs(r)-R));
    return psiT;
}
 
/**************************************************************************/
double LiebLinigerWaveFunction::delPsiTrial(const double r) {
    
    return 0.0;
}
 
/**************************************************************************/
double LiebLinigerWaveFunction::delSqPsiTrial(const double r) {
    
    return 0.0;
}
 
/**************************************************************************/
double LiebLinigerWaveFunction::PsiTrial(const int slice) {
    
    /* The cumulative value */
    double psiT = 1.0;
    Array <bool,1> doParticles(path.numBeadsAtSlice(slice));
    doParticles = true;
 
    dVec sep;                       // The spatial separation between beads.
    double r;                       // Distance between beads
    beadLocator bead1,bead2;
    bead1[0] = bead2[0] = slice;
    
    /* No cutoff */
//    for (bead1[1] = 0; bead1[1] < numParticles; bead1[1]++) {
//        /* The loop over all other particles, to find the total interaction
//         * potential */
//        for (bead2[1] = bead1[1]+1; bead2[1] < numParticles; bead2[1]++) {
//            sep = path.getSeparation(bead2,bead1);
//            r = sqrt(dot(sep,sep));
//            psiT *= PsiTrial(r);
//        } // bead2
//        
//  } // bead1
    
    /* Using cutoff */
    for (bead1[1] = 0; bead1[1] < path.numBeadsAtSlice(slice); bead1[1]++) {
        
        doParticles(bead1[1]) = false;
        
        /* Get the interaction list */
        lookup.updateInteractionList(path,bead1);
        
        /* Sum the interaction potential over all NN beads */
        for (int n = 0; n < lookup.numBeads; n++) {
            bead2 = lookup.beadList(n);
            if (doParticles(bead2[1])) {
                sep = path.getSeparation(bead2,bead1);
                r = sqrt(dot(sep,sep));
                psiT *= PsiTrial(r);
            }
        } // n
        
    } // bead1
    
    return psiT;
}
 
/**************************************************************************/
double LiebLinigerWaveFunction::PsiTrial(const beadLocator &bead1) {
    
    /* The cumulative value */
    double psiT = 1.0;
    
    dVec sep;                       // The spatial separation between beads.
    double r;                       // Distance between beads
    
    /* No cutoff */
    //    for (bead1[1] = 0; bead1[1] < numParticles; bead1[1]++) {
    //        /* The loop over all other particles, to find the total interaction
    //         * potential */
    //        for (bead2[1] = bead1[1]+1; bead2[1] < numParticles; bead2[1]++) {
    //            sep = path.getSeparation(bead2,bead1);
    //            r = sqrt(dot(sep,sep));
    //            psiT *= PsiTrial(r);
    //        } // bead2
    //
    //  } // bead1
    
    /* We only continue if bead1 is turned on */
    if (path.worm.beadOn(bead1)) {
        
        /* Fill up th nearest neighbor list */
        lookup.updateInteractionList(path,bead1);
        
        /* Sum the interaction potential over all NN beads */
        for (int n = 0; n < lookup.numBeads; n++) {
            
            sep = path.getSeparation(bead1,lookup.beadList(n));
            r = sqrt(dot(sep,sep));
            psiT *= PsiTrial(r);
        }
    }
    
    return psiT;
}
 
 
/**************************************************************************/
double LiebLinigerWaveFunction::gradSqPsiTrial(const int slice) {
    
    return 0.0;
}
 
 
// ---------------------------------------------------------------------------
// ---------------------------------------------------------------------------
// SUTHERLAND WAVEFUNCTION CLASS ---------------------------------------------
// ---------------------------------------------------------------------------
// ---------------------------------------------------------------------------
 
/**************************************************************************/
SutherlandWaveFunction::SutherlandWaveFunction(const Path &_path, LookupTable &_lookup, 
        double _lambda, string _name) :
WaveFunctionBase(_path,_lookup,_name)
{
    // The Sutherland model value of the interaction paramter \lambda
    lambda = _lambda;
 
    // pi/L
    pioL = M_PI/constants()->L();
}
 
/**************************************************************************/
SutherlandWaveFunction::~SutherlandWaveFunction() {
    // empty destructor
}
 
/**************************************************************************/
double SutherlandWaveFunction::PsiTrial(const int slice) {
    
    /* The cumulative value */
    double psiT = 1.0;
 
    /* The number of particles */
    int numParticles = path.numBeadsAtSlice(slice);    
 
    dVec sep;                       // The spatial separation between beads.
    double r;                       // Distance between beads
    beadLocator bead1,bead2;
    bead1[0] = bead2[0] = slice;
    
    /* No cutoff */
    for (bead1[1] = 0; bead1[1] < numParticles; bead1[1]++) {
        /* The loop over all other particles, to find the total wavefunction */
        for (bead2[1] = bead1[1]+1; bead2[1] < numParticles; bead2[1]++) {
            sep = path.getSeparation(bead2,bead1);
            r = sqrt(dot(sep,sep));
            psiT *= PsiTrial(r);
        } // bead2
    } // bead1
 
    return psiT;
}