LocalAction Class Reference

A base class to be inherited by actions that are local in imaginary time. More...

#include <action.h>

Inheritance diagram for LocalAction:

Collaboration diagram for LocalAction:

Public Member Functions
	LocalAction (const Path &, LookupTable &, PotentialBase *, PotentialBase *, WaveFunctionBase *, const TinyVector< double, 2 > &, const TinyVector< double, 2 > &, bool _local=true, string _name="Local", double _endFactor=1.0, int _period=1)
	Setup the path data members and local action factors.
virtual	~LocalAction ()
	Empty base constructor.
double	potentialAction ()
	Return the potential action for all time slices and all particles. More...
double	potentialAction (const beadLocator &)
	Return the potential action for a single bead indexed with beadIndex. More...
double	barePotentialAction (const beadLocator &)
	Return the bare potential action for a single bead indexed with beadIndex. More...
double	potentialActionCorrection (const beadLocator &)
	Return the potential action correction for a single bead. More...
double	potentialActionCorrection (const beadLocator &, const beadLocator &)
	Return the potential action correction for a path. More...
double	derivPotentialActionTau (int)
	The derivative of the potential action wrt tau on all links starting on slice. More...
double	derivPotentialActionLambda (int)
	The derivative of the potential action wrt lambda for all links starting on slice. More...
double	secondderivPotentialActionTau (int)
	The second derivative of the potential action wrt tau on all links starting on slice. More...
double	derivPotentialActionTau (int, double)
	The derivative of the potential action wrt tau on all links starting on slice within a cutoff radius. More...
double	derivPotentialActionLambda (int, double)
	The derivative of the potential action wrt lambda for all links starting on slice within a cutoff radius. More...
virtual dVec	gradPotentialAction (int slice)
double	rDOTgradUterm1 (int)
	Return the sum over particles at a given time slice of the gradient of the action potential for a single bead dotted into the bead's position. More...
double	rDOTgradUterm2 (int)
	Return the sum over particles at a given time slice of the gradient of the action potential for a single bead dotted into the bead's position. More...
virtual double	deltaDOTgradUterm1 (int slice)
virtual double	deltaDOTgradUterm2 (int slice)
virtual double	virKinCorr (int slice)
virtual TinyVector< double, 2 >	potential (int slice)
virtual double	potential (int slice, double maxR)
Public Member Functions inherited from ActionBase
	ActionBase (const Path &, LookupTable &, PotentialBase *, PotentialBase *, WaveFunctionBase *, bool _local=true, string _name="Base", double _endFactor=1.0, int _period=1)
	Setup the path data members and canonical re-weighting factors.
virtual	~ActionBase ()
	Empty base constructor.
string	getActionName ()
	Returns the action name.
double	kineticAction ()
	The full kinetic Action More...
double	kineticAction (const beadLocator &)
	The kinetic Action at a single slice More...
double	kineticAction (const beadLocator &, int wlLength)
	The kinetic Action for wlLength slices. More...
virtual double	potentialAction (const beadLocator &, const beadLocator &)
	Return the potential action for a path. More...
void	setShift (int _shift)
	The public method that sets the tau scaling factor.
int	getShift ()
	Get the tau scaling factor.
double	rho0 (const dVec &, const dVec &, int)
	The free-particle density matrix. More...
double	rho0 (const beadLocator &, const beadLocator &, int)
	The free-particle density matrix for two beadLocators with imaginary time separation M.
double	rho0 (const dVec &, const int)
	The free-particle density matrix for a given spatial and temporal separation. More...
double	ensembleWeight (const int)
	The ensemble particle number weighting factor. More...

Protected Member Functions
double	V (const beadLocator &)
	Returns the total value of the potential energy, including both the external potential, and that due to interactions for a single bead.
double	V (const int, const double)
	Returns the total value of the potential energy, including both the external potential, and that due to interactions for all particles in a single time slice within a certain cuttoff radius of the center of a cylinder. More...
TinyVector< double, 2 >	V (const int)
	Returns the external and interaction potential energy, for all particles on a single time slice. More...
double	gradVSquared (const beadLocator &)
	Return the gradient of the full potential squared for a single bead. More...
double	gradVSquared (const int)
	Return the gradient of the full potential squared for all beads at a single time slice. More...
double	gradVSquared (const int, const double)
	Return the gradient of the full potential squared for all beads at a single time slice restricted inside some cut-off radius. More...
double	Vnn (const beadLocator &)
	Returns the total value of the potential energy, including both the external potential, and that due to interactions for a single bead with all interactions confined to a single time slice using a nearest neighbor grid lookup table.
double	Vnn (const int)
	Returns the total value of the potential energy, including both the external potential, and that due to interactions for all particles in a single time slice using a lookup table which only considers particles within a sphere of some cutoff radius.
double	bareUnn (const beadLocator &, const beadLocator &)
double	gradVnnSquared (const beadLocator &)
	Return the gradient of the full potential squared for a single bead. More...
dVec	gradientV (const int)
	Return the gradient of the full potential for all beads at a single time slice. More...
dMat	tMatrix (const int)
	Returns the T-matrix needed to compute gradU. More...
double	RdotgradUterm1 (const int)
double	RdotgradUterm2 (const int)
double	deltadotgradUterm1 (const int)
	Return the sum over particles at a given time slice of the gradient of the action potential for a single bead dotted into the bead's position minus the center of mass of the WL that it belongs to. More...
double	deltadotgradUterm2 (const int)
	Return the sum over particles at a given time slice of the gradient of the action potential for a single bead dotted into the bead's position minus the center of mass of the WL that it belongs to. More...
double	virialKinCorrection (const int)
	Returns the value of the virial kinetic energy correction term. More...
dVec	gradU (const int)
	Returns the value of the gradient of the potential action for all beads at a given time slice. More...
Protected Member Functions inherited from ActionBase
double	tau ()
	The local shifted value of tau.
void	updateSepHist (const dVec &)
	Update the separation histogram. More...

Protected Attributes
int	eo
	Is a slice even or odd?
TinyVector< double, 2 >	VFactor
	The even/odd slice potential factor.
TinyVector< double, 2 >	gradVFactor
	The even/odd slice correction factor.
Protected Attributes inherited from ActionBase
string	name
	The name of the action.
LookupTable &	lookup
	We need a non-constant reference for updates.
const Path &	path
	A reference to the paths.
WaveFunctionBase *	waveFunctionPtr
	A pointer to a trial wave function object.
double	endFactor
	Mutiplictive factor of the potential action on ends.
int	shift
	The scaling factor for tau.
bool	canonical
	Are we in the canonical ensemble?
int	numBeads0
	The target number of beads.
bool	window
int	windowWidth
bool	gaussianEnsemble
double	gaussianEnsembleSD
beadLocator	bead2
beadLocator	bead3
dVec	sep
dVec	sep2
double	dSep

Additional Inherited Members
Data Fields inherited from ActionBase
const bool	local
	Is the action local in imaginary time?
const int	period
PotentialBase *	externalPtr
	The external potential.
PotentialBase *	interactionPtr
	The interaction potential.
Array< int, 1 >	sepHist
	A histogram of separations.
Array< int, 1 >	cylSepHist
	A histogram of separations for a cylinder.

Detailed Description

A base class to be inherited by actions that are local in imaginary time.

Definition at line 150 of file action.h.

Member Function Documentation

barePotentialAction()

double LocalAction::barePotentialAction ( const beadLocator &  beadIndex )

virtual

Return the bare potential action for a single bead indexed with beadIndex.

This action corresponds to the primitive approximation and may be used when attempting updates that use single slice rejection.

Reimplemented from ActionBase.

Definition at line 352 of file action.cpp.

                                                                     {
 
    double bareU = VFactor[beadIndex[0]%2]*tau()*Vnn(beadIndex);
 
#if PIGS
     /* We tack on a trial wave function and boundary piece if necessary */  
     if ( (beadIndex[0] == 0) || (beadIndex[0]== (constants()->numTimeSlices()-1)) ) {
         bareU *= 0.5*endFactor;
         bareU -= log(waveFunctionPtr->PsiTrial(beadIndex[0]));
     }
#endif
 
    return bareU;
}

Here is the call graph for this function:

deltadotgradUterm1()

double LocalAction::deltadotgradUterm1 ( const int  slice )

protected

Return the sum over particles at a given time slice of the gradient of the action potential for a single bead dotted into the bead's position minus the center of mass of the WL that it belongs to.

NOTE: This is the first term. These were split up because it is beneficial to be able to return them separately when computing the specific heat via the centroid virial estimator.

This includes both the external and interaction potentials.

Definition at line 1228 of file action.cpp.

                                                      {
 
    eo = slice % 2;
    int numParticles = path.numBeadsAtSlice(slice);
    int virialWindow = constants()->virialWindow();
 
    /* The two interacting particles */
    beadLocator bead1, beadNext, beadPrev, beadNextOld, beadPrevOld;
    bead1[0] = bead2[0] = slice;
 
    dVec gVi, gVe, gV, delta;
    double rDotgV = 0.0;
 
    /* We loop over the first bead */
    for (bead1[1] = 0; bead1[1] < numParticles; bead1[1]++) {
        gVi = 0.0;
        gVe = 0.0;
        gV = 0.0;
        delta = 0.0;
        /* Sum potential of bead1 interacting with all other beads at
         * a given time slice.*/
        for (bead2[1] = 0; bead2[1] < numParticles; bead2[1]++) {
 
            /* Avoid self interactions */
            if (!all(bead1==bead2)) {
 
                /* The interaction component of the force */
                gVi += interactionPtr->gradV(path.getSeparation(bead1,bead2));
            } 
        } // end bead2
        
        /* Now add the external component */
        gVe += externalPtr->gradV(path(bead1));
        
        /* Compute deviation of bead from COM of worldline, 
         * WITHOUT mirror image conv. */
        dVec runTotMore = 0.0;
        dVec runTotLess = 0.0;
        dVec COM = 0.0;
        dVec pos1 = path(bead1);
        beadNextOld = bead1;
        beadPrevOld = bead1;
        for (int gamma = 0; gamma < virialWindow; gamma++) {
            // move along worldline to next (previous) bead
            beadNext = path.next(bead1, gamma);
            beadPrev = path.prev(bead1, gamma);
            // keep running total of distance from first bead
            // to the current beads of interest 
            runTotMore += path.getSeparation(beadNext, beadNextOld);
            runTotLess += path.getSeparation(beadPrev, beadPrevOld);
            // update center of mass of WL
            COM += (pos1 + runTotMore) + (pos1 + runTotLess);
            // store current bead locations
            beadNextOld = beadNext;
            beadPrevOld = beadPrev;
        }
        COM /= (2.0*virialWindow);
 
        delta = pos1 - COM; // end delta computation.
        path.boxPtr->putInBC(delta);
        
        gV += (gVe + gVi);
 
        rDotgV += dot(gV, delta);
 
    } // end bead1
 
    return (VFactor[eo]*constants()->tau()*rDotgV);
}

Here is the call graph for this function:

deltadotgradUterm2()

double LocalAction::deltadotgradUterm2 ( const int  slice )

protected

Return the sum over particles at a given time slice of the gradient of the action potential for a single bead dotted into the bead's position minus the center of mass of the WL that it belongs to.

NOTE: This is the second term

This includes both the external and interaction potentials.

Definition at line 1307 of file action.cpp.

                                                      {
 
    eo = slice % 2;
    int numParticles = path.numBeadsAtSlice(slice);
    int virialWindow = constants()->virialWindow();
 
    /* The two interacting particles */
    beadLocator bead1, beadNext, beadPrev, beadNextOld, beadPrevOld;
    bead1[0] = bead2[0] = slice;
 
    dVec gVi, gVe, gV, g2V, delta;
 
    double term2 = 0.0;
    if (gradVFactor[eo] > EPS){
 
        /* constants for tMatrix */
        dVec rDiff = 0.0;
        double rmag = 0.0;
        double d2V = 0.0;
        double dV = 0.0;
        double dVe = 0.0;
        double dVi = 0.0;
        double g2Vi = 0.0;
        double g2Ve = 0.0;
 
        /* We loop over the first bead */
        for (bead1[1] = 0; bead1[1] < numParticles; bead1[1]++) {
            gV = 0.0;
            g2V = 0.0;
            delta = 0.0;
            dMat tMat = 0.0; // tMat(row, col)
            dVec gVdotT = 0.0;
 
            /* compute external potential derivatives */
            gVe = externalPtr->gradV(path(bead1));
            dVe = sqrt(dot(gVe,gVe));
            g2Ve = externalPtr->grad2V(path(bead1));
 
            gV += gVe;  // update full slice gradient at bead1
 
            /* Sum potential of bead1 interacting with all other beads at
             * a given time slice.*/
            for (bead2[1] = 0; bead2[1] < numParticles; bead2[1]++) {
 
                /* separation vector */
                rDiff = path.getSeparation(bead1, bead2);
                rmag = sqrt(dot(rDiff,rDiff));
 
                /* Avoid self interactions */
                if (!all(bead1==bead2)) {
 
                    /* Compute interaction potential derivatives */
                    gVi = interactionPtr->gradV(rDiff);
                    dVi = sqrt(dot(gVi,gVi));
                    g2Vi = interactionPtr->grad2V(rDiff);
                    
                    /* total derivatives between bead1 and bead2 at bead1 */
                    dV = dVi + dVe;
                    d2V = g2Vi + g2Ve;
 
                    /* compute the T-matrix for bead1 interacting with bead2 */
                    for (int a=0; a<NDIM; a++){
                        for (int b=0; b<NDIM; b++){
                            tMat(a,b) += rDiff(a)*rDiff(b)*d2V/(rmag*rmag)
                                - rDiff(a)*rDiff(b)*dV/pow(rmag,3);
                            if (a == b)
                                tMat(a,b) += dV/rmag;
                        }
                    }   // end T-matrix 
                    
                    gV += gVi;  // update full slice gradient at bead1
 
                }   
            }   // end bead2
 
            /* blitz++ product function seems broken in my current 
             * version, so this performs matrix-vector mult. 
             * Checked --MTG */
            for(int j=0; j<NDIM; j++){
                for(int i=0; i<NDIM; i++){
                    gVdotT(j) += gV(i)*tMat(j,i);
                }
            }  
            
            /* Compute deviation of bead from COM of worldline, 
             * WITHOUT mirror image conv.*/
            dVec runTotMore = 0.0;
            dVec runTotLess = 0.0;
            dVec COM = 0.0;
            dVec pos1 = path(bead1);
            beadNextOld = bead1;
            beadPrevOld = bead1;
            for (int gamma = 0; gamma < virialWindow; gamma++) {
                // move along worldline to next (previous) bead
                beadNext = path.next(bead1, gamma);
                beadPrev = path.prev(bead1, gamma);
                // keep running total of distance from first bead
                // to the current beads of interest
                runTotMore += path.getSeparation(beadNext, beadNextOld);
                runTotLess += path.getSeparation(beadPrev, beadPrevOld);
                // update center of mass of WL
                COM += (pos1 + runTotMore) + (pos1 + runTotLess);
                // store current bead locations
                beadNextOld = beadNext;
                beadPrevOld = beadPrev;
            }
            COM /= (2.0*virialWindow);
            delta = pos1 - COM; // end delta computation.
           
            path.boxPtr->putInBC(delta);
 
            term2 += dot(gVdotT, delta);
 
        }   // end bead1
 
        term2 *= (2.0 * gradVFactor[eo] * pow(tau(),3) * constants()->lambda());
    }
 
    return (term2);
}

Here is the call graph for this function:

derivPotentialActionLambda() [1/2]

double LocalAction::derivPotentialActionLambda ( int  slice )

virtual

The derivative of the potential action wrt lambda for all links starting on slice.

It is essential to have these slice overloaded values as we need to be careful of the factor of 1/2 or i<j in the full potential action.

Parameters

slice

the imaginary timeslice of the first link

Reimplemented from ActionBase.

Definition at line 469 of file action.cpp.

                                                         {
    eo = (slice % 2);
 
    /* As long as thers is a finite correction, we include it */
    if ( gradVFactor[eo] > EPS )
        return gradVFactor[eo] * tau() * tau() * tau() * gradVSquared(slice);
    else
        return 0.0;
}

Here is the call graph for this function:

derivPotentialActionLambda() [2/2]

double LocalAction::derivPotentialActionLambda ( int  slice, double  maxR  )

virtual

The derivative of the potential action wrt lambda for all links starting on slice within a cutoff radius.

It is essential to have these slice overloaded values as we need to be careful of the factor of 1/2 or i<j in the full potential action.

Parameters

slice	the imaginary timeslice of the first link
maxR	the cutoff radius

Reimplemented from ActionBase.

Definition at line 515 of file action.cpp.

                                                                      {
    eo = (slice % 2);
 
    /* As long as there is a finite correction, we include it */
    if ( gradVFactor[eo] > EPS )
        return gradVFactor[eo] * tau() * tau() * tau() * gradVSquared(slice,maxR);
    else
        return 0.0;
}

Here is the call graph for this function:

derivPotentialActionTau() [1/2]

double LocalAction::derivPotentialActionTau ( int  slice )

virtual

The derivative of the potential action wrt tau on all links starting on slice.

It is essential to have these slice overloaded values as we need to be careful of the factor of 1/2 or i<j in the full potential action.

Parameters

slice

the imaginary timeslice of the first link

Reimplemented from ActionBase.

Definition at line 422 of file action.cpp.

                                                      {
 
    double dU = 0.0;
    eo = (slice % 2);
 
    /* We first compute the base potential action piece */
    dU = VFactor[eo]*sum(V(slice));
 
    /* As long as there is a finite correction, we include it */
    if ( gradVFactor[eo] > EPS )
        dU += 3.0 * gradVFactor[eo] * tau() * tau() * constants()->lambda() 
            * gradVSquared(slice);
    return (dU);
 
}

Here is the call graph for this function:

derivPotentialActionTau() [2/2]

double LocalAction::derivPotentialActionTau ( int  slice, double  maxR  )

virtual

The derivative of the potential action wrt tau on all links starting on slice within a cutoff radius.

It is essential to have these slice overloaded values as we need to be careful of the factor of 1/2 or i<j in the full potential action.

Parameters

slice	the imaginary timeslice of the first link
maxR	the cutoff radius

Reimplemented from ActionBase.

Definition at line 489 of file action.cpp.

                                                                   {
 
    double dU = 0.0;
    eo = (slice % 2);
 
    /* We first compute the base potential action piece */
    dU = VFactor[eo]*V(slice,maxR);
 
    /* As long as there is a finite correction, we include it */
    if ( gradVFactor[eo] > EPS )
        dU += 3.0 * gradVFactor[eo] * tau() * tau() * constants()->lambda() 
            * gradVSquared(slice,maxR);
    return (dU);
 
}

Here is the call graph for this function:

gradientV()

dVec LocalAction::gradientV ( const int  slice )

protected

Return the gradient of the full potential for all beads at a single time slice.

This includes both the external and interaction potentials.

Definition at line 978 of file action.cpp.

                                           {
 
    int numParticles = path.numBeadsAtSlice(slice);
 
    /* The two interacting particles */
    beadLocator bead1;
    bead1[0] = bead2[0] = slice;
 
    dVec gV;                    // The 'force'
 
    /* We loop over the first bead */
    for (bead1[1] = 0; bead1[1] < numParticles; bead1[1]++) {
        gV = 0.0;
        /* Sum up potential for all other active beads in the system */
        for (bead2[1] = 0; bead2[1] < numParticles; bead2[1]++) {
 
            /* Avoid self interactions */
            if (!all(bead1==bead2)) {
 
                /* The interaction component of the force */
                gV += interactionPtr->gradV(path.getSeparation(bead1,bead2));
            } 
        } // end bead2
 
        /* Now add the external component */
        gV += externalPtr->gradV(path(bead1));
    } // end bead1
 
    return gV;
}

Here is the call graph for this function:

Here is the caller graph for this function:

gradU()

dVec LocalAction::gradU ( const int  slice )

protected

Returns the value of the gradient of the potential action for all beads at a given time slice.

This should be used for the pressure.

Definition at line 1458 of file action.cpp.

                                       {
 
    dVec gU2 = 0.0;
    dMat tM = tMatrix(slice);
    dVec gV = gradientV(slice);
 
    /* Combine the potential and a possible correction */
    eo = (slice % 2);
    dVec gU1 = VFactor[eo]*tau()*gradientV(slice);
 
    /* We only add the correction if it is finite */
    if ( gradVFactor[eo] > EPS ) {
        /* blitz++ product function seems broken in my current 
         * version, so this performs matrix-vector mult. */
        for(int j=0; j<NDIM; j++){
            for(int i=0; i<NDIM; i++){
                gU2(i) += gV(j)*tM(i,j);
            }
        }
        /* now scale by constants */
        gU2 *= 2.0 * gradVFactor[eo] * pow(tau(),3) * constants()->lambda();
    }
    
    return ( gU1+gU2 );
}

Here is the call graph for this function:

gradVnnSquared()

double LocalAction::gradVnnSquared ( const beadLocator &  bead1 )

protected

Return the gradient of the full potential squared for a single bead.

This includes both the external and interaction potentials using the nearest neighbor lookup table.

Definition at line 915 of file action.cpp.

                                                           {
 
    /* This should always be called after Vnn, such that the lookup table 
     * interaction list has been updated!!! */
 
    double totF2 = 0.0;     // The total force squared
 
    if (path.worm.beadOn(bead1)) {
 
        /* The 'forces' and particle separation */
        dVec Fext1,Fext2;
        dVec Fint1,Fint2,Fint3;
 
        /* Get the gradient squared part for the external potential*/
        Fext1 = externalPtr->gradV(path(bead1));
 
        /* We first loop over bead2's interacting with bead1 via the nn lookup table */
        Fint1 = 0.0;
        for (int n = 0; n < lookup.numBeads; n++) {
            bead2 = lookup.beadList(n);
 
            /* Eliminate self and null interactions */
            if ( !all(bead1 == bead2) ) {
 
                /* Get the separation between beads 1 and 2 and compute the terms in the
                 * gradient squared */
                sep = lookup.beadSep(n); 
                Fint2 = interactionPtr->gradV(sep);
                Fint1 -= Fint2;
                Fext2 = externalPtr->gradV(path(bead2));
 
                /* We now loop over bead3, this is the time-intensive part of the calculation */
                Fint3 = 0.0;
                for (int m = 0; m < lookup.numBeads; m++) {
                    bead3 = lookup.beadList(m);
 
                    /* Eliminate self-interactions */
                    if ( !all(bead3==bead2) && !all(bead3==bead1) ) {
                        sep = path.getSeparation(bead2,bead3);
                        Fint3 += interactionPtr->gradV(sep);
                    }
 
                } // end bead3
 
                totF2 += dot(Fint2,Fint2) + 2.0*dot(Fext2,Fint2) + 2.0*dot(Fint2,Fint3);
 
            } // bead2 != bead1
 
        } // end bead2
 
        totF2 += dot(Fext1,Fext1) + 2.0 * dot(Fext1,Fint1) + dot(Fint1,Fint1);
 
    } // bead1 on
 
    return totF2;
}

Here is the call graph for this function:

Here is the caller graph for this function:

gradVSquared() [1/3]

double LocalAction::gradVSquared ( const beadLocator &  bead1 )

protected

Return the gradient of the full potential squared for a single bead.

This includes both the external and interaction potentials.

Definition at line 762 of file action.cpp.

                                                         {
 
    double totF2 = 0.0;     // The total force squared
 
    if (path.worm.beadOn(bead1)) {
 
        /* All interacting beads are on the same slice. */
        bead2[0] = bead3[0] = bead1[0];
 
        /* The 'forces' and particle separation */
        dVec Fext1,Fext2;
        dVec Fint1,Fint2,Fint3;
 
        int numParticles = path.numBeadsAtSlice(bead1[0]);
    
        /* Get the gradient squared part for the external potential*/
        Fext1 = externalPtr->gradV(path(bead1));
 
        /* We loop through all beads and compute the forces between beads
         * 1 and 2 */
        Fint1 = 0.0;
        for (bead2[1] = 0; bead2[1] < numParticles; bead2[1]++) {
 
            if (!all(bead1==bead2)) {
 
                sep = path.getSeparation(bead2,bead1);
                Fint2 = interactionPtr->gradV(sep);
                Fint1 -= Fint2;
                Fext2 = externalPtr->gradV(path(bead2));
 
                /* There is a single term that depends on this additional interaction
                 * between beads 2 and 3.  This is where all the time is taken */
                Fint3 = 0.0;
                for (bead3[1] = 0; bead3[1] < numParticles; bead3[1]++) {
                    if ( !all(bead3==bead2) && !all(bead3==bead1) ) {
                        sep = path.getSeparation(bead2,bead3);
                        Fint3 += interactionPtr->gradV(sep);
                    }
                } // for bead3
 
                totF2 += dot(Fint2,Fint2) + 2.0*dot(Fext2,Fint2) + 2.0*dot(Fint2,Fint3);
 
            } //bead1 != bead2
 
        } // for bead2
 
        totF2 += dot(Fext1,Fext1) + 2.0 * dot(Fext1,Fint1) + dot(Fint1,Fint1);
 
    } // bead1 On
 
    return totF2;
}

Here is the call graph for this function:

Here is the caller graph for this function:

gradVSquared() [2/3]

double LocalAction::gradVSquared ( const int  slice )

protected

Return the gradient of the full potential squared for all beads at a single time slice.

This includes both the external and interaction potentials.

Definition at line 821 of file action.cpp.

                                                {
 
    double totF2 = 0.0;
 
    int numParticles = path.numBeadsAtSlice(slice);
 
    /* The two interacting particles */
    beadLocator bead1;
    bead1[0] = bead2[0] = slice;
 
    dVec F;                 // The 'force'
 
    /* We loop over the first bead */
    for (bead1[1] = 0; bead1[1] < numParticles; bead1[1]++) {
 
        F = 0.0;
        /* Sum up potential for all other active beads in the system */
        for (bead2[1] = 0; bead2[1] < numParticles; bead2[1]++) {
 
            /* Avoid self interactions */
            if (!all(bead1==bead2)) {
 
                /* The interaction component of the force */
                F += interactionPtr->gradV(path.getSeparation(bead1,bead2));
            } 
 
        } // end bead2
 
        /* Now add the external component */
        F += externalPtr->gradV(path(bead1));
 
        totF2 += dot(F,F);
    } // end bead1
 
    return totF2;
}

Here is the call graph for this function:

gradVSquared() [3/3]

double LocalAction::gradVSquared ( const int  slice, const double  maxR  )

protected

Return the gradient of the full potential squared for all beads at a single time slice restricted inside some cut-off radius.

This includes both the external and interaction potentials.

Definition at line 864 of file action.cpp.

                                                                   {
 
    double totF2 = 0.0;
 
    int numParticles = path.numBeadsAtSlice(slice);
 
    /* The two interacting particles */
    beadLocator bead1;
    bead1[0] = bead2[0] = slice;
    dVec r1;
    double r1sq;
 
    dVec F;                 // The 'force'
 
    /* We loop over the first bead */
    for (bead1[1] = 0; bead1[1] < numParticles; bead1[1]++) {
 
        r1 = path(bead1);
        r1sq = 0.0;
        for (int i = 0; i < NDIM-1; i++)
            r1sq += r1[i]*r1[i];
 
        if (r1sq < maxR*maxR) {
 
            F = 0.0;
            /* Sum up potential for all other active beads in the system */
            for (bead2[1] = 0; bead2[1] < numParticles; bead2[1]++) {
                /* Avoid self interactions */
                if (!all(bead1==bead2)) {
 
                    /* The interaction component of the force */
                    F += interactionPtr->gradV(path.getSeparation(bead1,bead2));
                } 
            } // end bead2
 
            /* Now add the external component */
            F += externalPtr->gradV(path(bead1));
 
            totF2 += dot(F,F);
        } // maxR
    } // end bead1
 
    return totF2;
}

Here is the call graph for this function:

potentialAction() [1/2]

double LocalAction::potentialAction ( )

virtual

Return the potential action for all time slices and all particles.

Computes the total potential energy by summing over all particles and time slices. We must be very careful with the factor of 1/2 or the fact that the sum is over i<j for particles to avoid double counting.

!!NB!! It is important to use V and gradVSquared here, as this is used for debugging purposes with check move and must match up with the single bead calculation.

Reimplemented from ActionBase.

Definition at line 289 of file action.cpp.

                                     {
 
    double totU = 0.0;
 
    /* Combine the potential and a possible correction */
    for (int slice = 0; slice < path.numTimeSlices; slice++) {
        eo = (slice % 2);
        totU += VFactor[eo]*tau()*sum(V(slice));
        
        /* We only add the correction if it is finite */
         if ( gradVFactor[eo] > EPS ) 
             totU += gradVFactor[eo] * tau() * tau() * tau() * constants()->lambda() 
                 * gradVSquared(slice);
    }
 
    return ( totU );
}

Here is the call graph for this function:

potentialAction() [2/2]

double LocalAction::potentialAction ( const beadLocator &  beadIndex )

virtual

Return the potential action for a single bead indexed with beadIndex.

We have to use the beadIndex form of gradVSquared here to ensure that bead indexed action changes are equal to total potential action changes. We use a pre-processor directive to determine whether or not to use the nearest neighbor lookup table.

!!NB!! If we are using checkMove to debug moves you MUST always include the correction term.

Parameters

beadIndex

the bead that we are computing the potential action for

Reimplemented from ActionBase.

Definition at line 320 of file action.cpp.

                                                                 {
 
    double bareU = 0.0;
    double corU = 0.0;
    eo = (beadIndex[0] % 2);
 
    /* We first compute the base potential action piece */
    bareU = VFactor[eo]*tau()*Vnn(beadIndex);
 
    /* If we have a finite correction and are not at the base level
     * of bisection (shift=1) we compute the full correction */
     if ( (shift == 1) && (gradVFactor[eo] > EPS) )
        corU = gradVFactor[eo] * tau() * tau() * tau() * constants()->lambda() 
            * gradVnnSquared(beadIndex);
     
#if PIGS
     /* We tack on a trial wave function and boundary piece if necessary */  
     if ( (beadIndex[0] == 0) || (beadIndex[0] == (constants()->numTimeSlices()-1)) ) {
         bareU *= 0.5*endFactor;
         bareU -= log(waveFunctionPtr->PsiTrial(beadIndex[0]));
     }
#endif
 
     return (bareU + corU);
}

Here is the call graph for this function:

potentialActionCorrection() [1/2]

double LocalAction::potentialActionCorrection ( const beadLocator &  beadIndex )

virtual

Return the potential action correction for a single bead.

Provided that the correction is small, we return its value.

Reimplemented from ActionBase.

Definition at line 372 of file action.cpp.

                                                                           {
 
    eo = (beadIndex[0] % 2);
    /* If we have a finite correction */
    if (gradVFactor[eo] > EPS) {
 
        /* We need to update the lookup table here */
        lookup.updateInteractionList(path,beadIndex);
 
        /* Compute the correction */
       return (gradVFactor[eo] * tau() * tau() * tau() * constants()->lambda() 
               * gradVnnSquared(beadIndex));
    }
    else 
        return 0.0;
 
}

Here is the call graph for this function:

Here is the caller graph for this function:

potentialActionCorrection() [2/2]

double LocalAction::potentialActionCorrection ( const beadLocator &  startBead, const beadLocator &  endBead  )

virtual

Return the potential action correction for a path.

Given a starting and ending bead, compute the potential action correction for the path.

Parameters

startBead	the starting bead
endBead	the end bead on the path

Reimplemented from ActionBase.

Definition at line 399 of file action.cpp.

                                    {
 
    double corU = 0.0;
    beadLocator beadIndex;
    beadIndex = startBead;
    do {
        corU += potentialActionCorrection(beadIndex);
        beadIndex = path.next(beadIndex);
    } while (!all(beadIndex==path.next(endBead)));
 
    return corU;
}

Here is the call graph for this function:

rDOTgradUterm1()

double LocalAction::rDOTgradUterm1 ( int  slice )

virtual

Return the sum over particles at a given time slice of the gradient of the action potential for a single bead dotted into the bead's position.

Used for regular virial energy and thermodynamic pressure.

NOTE: This is the first term. These were split up because it is beneficial to be able to return them separately when computing the specific heat via the centroid virial estimator.

This includes both the external and interaction potentials.

Reimplemented from ActionBase.

Definition at line 1079 of file action.cpp.

                                                  {
 
    eo = slice % 2;
    int numParticles = path.numBeadsAtSlice(slice);
 
    /* The two interacting particles */
    beadLocator bead1;
    bead1[0] = bead2[0] = slice;
 
    dVec gVi, gV;
    double rDotgV = 0.0;
 
    /* We loop over the first bead */
    for (bead1[1] = 0; bead1[1] < numParticles; bead1[1]++) {
        gVi = 0.0;
        /* Sum potential of bead1 interacting with all other beads at
         * a given time slice.*/
        for (bead2[1] = 0; bead2[1] < numParticles; bead2[1]++) {
 
            /* Avoid self interactions */
            if (!all(bead1==bead2)) {
 
                /* The interaction component of the force */
                gVi += interactionPtr->gradV(path.getSeparation(bead1,bead2));
            } 
        } // end bead2
        
        /* Now add the external component */
        gV = (externalPtr->gradV(path(bead1)) + gVi);
 
        rDotgV += dot(gV, path(bead1));
 
    } // end bead1
 
    return (VFactor[eo]*constants()->tau()*rDotgV);
}

Here is the call graph for this function:

rDOTgradUterm2()

double LocalAction::rDOTgradUterm2 ( int  slice )

virtual

Return the sum over particles at a given time slice of the gradient of the action potential for a single bead dotted into the bead's position.

Used for regular virial energy and thermodynamic pressure.

NOTE: This is the second term

This includes both the external and interaction potentials.

Reimplemented from ActionBase.

Definition at line 1126 of file action.cpp.

                                                  {
 
    eo = slice % 2;
    int numParticles = path.numBeadsAtSlice(slice);
 
    /* The two interacting particles */
    beadLocator bead1;
    bead1[0] = bead2[0] = slice;
 
    dVec gVi, gVe, gV, g2V;
 
    double term2 = 0.0;
    if (gradVFactor[eo] > EPS){
 
        /* constants for tMatrix */
        dVec rDiff = 0.0;
        double rmag = 0.0;
        double d2V = 0.0;
        double dV = 0.0;
        double dVe = 0.0;
        double dVi = 0.0;
        double g2Vi = 0.0;
        double g2Ve = 0.0;
 
        /* We loop over the first bead */
        for (bead1[1] = 0; bead1[1] < numParticles; bead1[1]++) {
            gV = 0.0;
            g2V = 0.0;
            dMat tMat = 0.0; // tMat(row, col)
            dVec gVdotT = 0.0;
 
            /* compute external potential derivatives */
            gVe = externalPtr->gradV(path(bead1));
            dVe = sqrt(dot(gVe,gVe));
            g2Ve = externalPtr->grad2V(path(bead1));
 
            gV += gVe;  // update full slice gradient at bead1
 
            /* Sum potential of bead1 interacting with all other beads at
             * a given time slice.*/
            for (bead2[1] = 0; bead2[1] < numParticles; bead2[1]++) {
 
                /* separation vector */
                rDiff = path.getSeparation(bead1, bead2);
                rmag = sqrt(dot(rDiff,rDiff));
 
                /* Avoid self interactions */
                if (!all(bead1==bead2)) {
 
                    /* Compute interaction potential derivatives */
                    gVi = interactionPtr->gradV(rDiff);
                    dVi = sqrt(dot(gVi,gVi));
                    g2Vi = interactionPtr->grad2V(rDiff);
                    
                    /* total derivatives between bead1 and bead2 at bead1 */
                    dV = dVi + dVe;
                    d2V = g2Vi + g2Ve;
 
                    /* compute the T-matrix for bead1 interacting with bead2 */
                    for (int a=0; a<NDIM; a++){
                        for (int b=0; b<NDIM; b++){
                            tMat(a,b) += rDiff(a)*rDiff(b)*d2V/(rmag*rmag)
                                - rDiff(a)*rDiff(b)*dV/pow(rmag,3);
                            if (a == b)
                                tMat(a,b) += dV/rmag;
                        }
                    }   // end T-matrix 
                    
                    gV += gVi;  // update full slice gradient at bead1
 
                }   
            }   // end bead2
 
            /* blitz++ product function seems broken in my current 
             * version, so this performs matrix-vector mult. 
             * Checked --MTG */
            for(int j=0; j<NDIM; j++){
                for(int i=0; i<NDIM; i++){
                    gVdotT(j) += gV(i)*tMat(j,i);
                }
            }  
            term2 += dot(gVdotT, path(bead1));
 
        }   // end bead1
 
        term2 *= (2.0 * gradVFactor[eo] * pow(tau(),3) * constants()->lambda());
    }
 
    return (term2);
}

Here is the call graph for this function:

secondderivPotentialActionTau()

double LocalAction::secondderivPotentialActionTau ( int  slice )

virtual

The second derivative of the potential action wrt tau on all links starting on slice.

It is essential to have these slice overloaded values as we need to be careful of the factor of 1/2 or i<j in the full potential action.

Parameters

slice

the imaginary timeslice of the first link

Reimplemented from ActionBase.

Definition at line 447 of file action.cpp.

                                                            {
 
    double d2U = 0.0;
    eo = (slice % 2);
 
    /* As long as there is a finite correction, we include it */
    if ( gradVFactor[eo] > EPS )
        d2U += 6.0 * gradVFactor[eo] * tau() * constants()->lambda() 
            * gradVSquared(slice);
    return (d2U);
 
}

Here is the call graph for this function:

tMatrix()

dMat LocalAction::tMatrix ( const int  slice )

protected

Returns the T-matrix needed to compute gradU.

This is used for the virial energy estimator in the TI or GSF action. CURRENTLY ONLY RETURNS VALUES FOR INTERACTIONS, NOT EXTERNAL.

Definition at line 1015 of file action.cpp.

                                         {
 
    int numParticles = path.numBeadsAtSlice(slice);
 
    dMat tMat = 0.0; // tMat(row, col)
 
    dVec rDiff = 0.0;
    double rmag = 0.0;
    double d2V = 0.0;
    double dV = 0.0;
    
    /* two interacting particles */
    beadLocator bead1;
    bead1[0] = bead2[0] = slice;
    
    for (int a=0; a<NDIM; a++){
        for (int b=0; b<NDIM; b++){
 
            /* loop over first bead */
            for (bead1[1]=0; bead1[1]<numParticles; bead1[1]++){
                
                /* loop over all other beads */
                for (bead2[1]=0; bead2[1]<numParticles; bead2[1]++){
 
                    /* avoid self-interactions */
                    if (!all(bead1==bead2)) {
                        
                        rDiff = path.getSeparation(bead1, bead2);
                        rmag = sqrt(dot(rDiff,rDiff));
                        d2V = interactionPtr->grad2V(rDiff);
                        d2V += externalPtr->grad2V(path(bead1));
                        dV = sqrt(dot(interactionPtr->gradV(rDiff)
                                    ,interactionPtr->gradV(rDiff)));
                        dV += sqrt(dot(externalPtr->gradV(path(bead1))
                                    ,externalPtr->gradV(path(bead1))));
 
                        tMat(a,b) += rDiff(a)*rDiff(b)*d2V/(rmag*rmag);
                        if (a != b)
                            tMat(a,b) -= (rDiff(a)*rDiff(b)/pow(rmag,3))*dV;
                        else
                            tMat(a,b) += (1.0/rmag - rDiff(a)*rDiff(b)/pow(rmag,3))*dV;
                    }
                } // end bead2
            } // end bead1
        }
    }
 
    tMat /= 2.0; // don't want to double count interactions
    
    return ( tMat ); 
}

Here is the call graph for this function:

Here is the caller graph for this function:

V() [1/2]

TinyVector< double, 2 > LocalAction::V ( const int  slice )

protected

Returns the external and interaction potential energy, for all particles on a single time slice.

This is really only used for either debugging or during the calculation of the potential energy. As such, we update the separation histogram here.

Definition at line 573 of file action.cpp.

                                                   {
 
    double totVint = 0.0;
    double totVext = 0.0;
 
    beadLocator bead1;
    bead1[0] = bead2[0] = slice;
 
    int numParticles = path.numBeadsAtSlice(slice);
 
    /* Initialize the separation histogram */
    sepHist = 0;
 
    /* Calculate the total potential, including external and interaction
     * effects*/
    for (bead1[1] = 0; bead1[1] < numParticles; bead1[1]++) {
 
            /* Get the state of bead 1 */
            beadState state1 = path.worm.getState(bead1);
 
            /* Evaluate the external potential */
            totVext += path.worm.factor(state1)*externalPtr->V(path(bead1));
 
            /* 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);
                updateSepHist(sep);
                totVint += path.worm.factor(state1,bead2) * interactionPtr->V(sep);
            } // bead2
 
    } // bead1
 
    /* Separate the external and interaction parts */ 
    return TinyVector<double,2>(totVext,totVint);
}

Here is the call graph for this function:

V() [2/2]

double LocalAction::V ( const int  slice, const double  maxR  )

protected

Returns the total value of the potential energy, including both the external potential, and that due to interactions for all particles in a single time slice within a certain cuttoff radius of the center of a cylinder.

This is really only used for either debugging or during the calculation of the potential energy. As such, we update the separation histogram here.

Definition at line 619 of file action.cpp.

                                                        {
 
    double totVint = 0.0;
    double totVext = 0.0;
    dVec r1,r2;
 
    double r1sq,r2sq;
 
    beadLocator bead1;
    bead1[0] = bead2[0] = slice;
 
    int numParticles = path.numBeadsAtSlice(slice);
 
    /* Initialize the separation histogram */
    cylSepHist = 0;
 
    /* Calculate the total potential, including external and interaction
     * effects*/
    for (bead1[1] = 0; bead1[1] < numParticles; bead1[1]++) {
 
        r1 = path(bead1);
 
        r1sq = 0.0;
        for (int i = 0; i < NDIM-1; i++)
            r1sq += r1[i]*r1[i];
 
        if (r1sq < maxR*maxR) {
 
            /* Evaluate the external potential */
            totVext += externalPtr->V(path(bead1));
 
            /* The loop over all other particles, to find the total interaction
             * potential */
            for (bead2[1] = bead1[1]+1; bead2[1] < numParticles; bead2[1]++) {
                r2 = path(bead2);
 
                r2sq = 0.0;
                for (int i = 0; i < NDIM-1; i++)
                    r2sq += r2[i]*r2[i];
 
                sep = path.getSeparation(bead2,bead1);
                if (r2sq < maxR*maxR) {
                    int nR = int(abs(sep[2])/dSep);
                    if (nR < NPCFSEP)
                        ++cylSepHist(nR);
                }
                totVint += interactionPtr->V(sep);
            } // bead2
 
        } // maxR
    } // bead1
 
    return ( totVext + totVint );
}

Here is the call graph for this function:

virialKinCorrection()

double LocalAction::virialKinCorrection ( const int  slice )

protected

Returns the value of the virial kinetic energy correction term.

See also

S. Jang, S. Jang and G.A. Voth, J. Chem. Phys. 115, 7832 (2001).

Definition at line 1433 of file action.cpp.

                                                       {
 
    double vKC = 0.0;
 
    /* Combine the potential and a possible correction */
    eo = (slice % 2);
 
    /* We only add the correction if it is finite */
    if ( gradVFactor[eo] > EPS ) {
        
        vKC += gradVSquared(slice);
        
        /* scale by constants */
        vKC *= gradVFactor[eo] * pow(tau(),3) * constants()->lambda();
    }
    
    return ( vKC );
}

Here is the call graph for this function:

---

The documentation for this class was generated from the following files:

• include/action.h
• src/action.cpp