Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F102274640
wavepacket.h
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Feb 19, 00:12
Size
7 KB
Mime Type
text/x-c++
Expires
Fri, Feb 21, 00:12 (2 d)
Engine
blob
Format
Raw Data
Handle
24320777
Attached To
rLAMMPS lammps
wavepacket.h
View Options
# ifndef WAVEPACKET_H
# define WAVEPACKET_H
# ifndef _USE_MATH_DEFINES
# define _USE_MATH_DEFINES
# endif
# include <cmath>
# include <complex>
# include <functional>
# include "cvector_3.h"
using namespace std;
/** @file wpmd.h
@brief Classes to handle Gaussian Wave Packets. */
// Constants
const double MIN_EXP_ARG = -15.; // Minimum noticeable argument for exp function
class WavePacket;
///\en Template for v=der operation in \ref Wavepacket::int2phys_der()
template<class Type>
struct eq_second : public binary_function <Type, Type, Type> {
Type operator()(const Type& _Left, const Type& _Right) const{
return _Right;
}
};
///\en Template for v=-der operation in \ref Wavepacket::int2phys_der()
template<class Type>
struct eq_minus_second : public binary_function <Type, Type, Type> {
Type operator()(const Type& _Left, const Type& _Right) const{
return -_Right;
}
};
///\en Compares complex numbers on a per component basis.
/// \return \retval 0 if all component differences are 0 within tolerance \a tol (EQUAL),
/// \retval -1 for LESS
/// \retval 2 for GREATER
template< class CT >
int compare_compl(const CT &a, const CT& b, double tol=0.){
double dd=real(a)-real(b);
if(dd<-tol)
return -1;
if(dd>tol)
return 1;
dd=imag(a)-imag(b);
if(dd<-tol)
return -1;
if(dd>tol)
return 1;
return 0;
}
///\en Compares vectors on a per component basis.
/// \return \retval 0 if all component differences are 0 within tolerance \a tol (EQUAL),
/// \retval -1 for LESS
/// \retval 2 for GREATER
inline int compare_vec(const Vector_3 &a, const Vector_3& b, double tol=0.){
for(int i=0;i<3;i++){
double dd=a[i]-b[i];
if(dd<-tol)
return -1;
if(dd>tol)
return 1;
}
return 0;
}
/// wavepacket is w(x)=exp(-a*x^2+b*x+lz)
class WavePacket{
/// constructs a conjugate packet
friend WavePacket conj(const WavePacket &wp);
public:
cdouble a;
cVector_3 b;
cdouble lz;
WavePacket(const cdouble &sa=cdouble(1.,0.),const cVector_3 &sb=cVector_3(0.,0.), const cdouble &slz=0.): a(sa), b(sb), lz(slz){
}
WavePacket operator*(const WavePacket& other) const {
return WavePacket(a+other.a,b+other.b,lz+other.lz);
}
/// returns the integral of w(x) over 3D space
cdouble integral() const {
cdouble z = lz + b.norm2()/(4.*a);
if(real(z) < MIN_EXP_ARG) return 0.;
return pow(M_PI/a,3./2.)*exp(z);
}
/// init normalized packet with physical parameters: r0, p0, width, pw
/// w(x)=(3/2pi width^(3/4)exp[-(3/(4 width^2)-i pw/(2*width) )(x-r)^2+i p (x-r)]
void init(const double width=1., const Vector_3 &r=0., const Vector_3 &p=0., const double pw=0.){
a = (3./(2.*width) - cdouble(0.,1.)*pw)/(2.*width);
b = (2.*a)*r + cdouble(0.,1.)*p;
lz = log(3./(2.*M_PI*width*width))*(3./4.) + (-a*r.norm2() - cdouble(0.,1.)*(p*r));
}
/// init normalized packet with complex parameters a and b
/// w(x)=(3/2pi width^(3/4)exp[-(3/(4 width^2)-i pw/(2*width) )(x-r)^2+i p (x-r)]
void init(const cdouble &a_, const cVector_3 &b_){
a = a_;
b = b_;
Vector_3 r = get_r();
double r2 = r.norm2();
lz = cdouble( log( real(a)*(2./M_PI) )*(3./4.) - r2*real(a), r2*imag(a) - r*imag(b) );
}
cdouble operator()(const Vector_3 &x) const {
return exp(lz - a*x.norm2() + b*x);
}
/// ajusts lz so that Integral[w*(conj(w))] is 1 after this operation
WavePacket& normalize(){
//lz=0.;
//lz=overlap(conj(*this));
//lz=(-1./2)*log(lz);
Vector_3 r = get_r();
double r2 = r.norm2();
lz = cdouble( log( real(a)*(2./M_PI) )*(3./4.) - r2*real(a), r2*imag(a) - r*imag(b) );
return *this;
}
/// computes 3D overlap of wavepackets Inegral[w*other]
cdouble overlap(const WavePacket &other) const{
WavePacket wp=(*this)*other;
return wp.integral();
}
/// returns translated packet to the position of r'=r+dr
WavePacket translate(const Vector_3 &dr) const {
WavePacket wp(a,b,lz);
wp.b+=2*dr*a;
Vector_3 r=get_r();
double dr2=(r+dr).norm2()-r.norm2();
wp.lz+=-a*dr2-cdouble(0.,(get_p()*dr));
return wp;
}
/// width
double get_width() const{
return sqrt(3./4./real(a));
}
/// width momentum
double get_pwidth() const{
return -imag(a)*sqrt(3./real(a));
}
/// both width and width momentum
pair<double,double> get_width_pars() const{
double c=sqrt(3./2./real(a));
return make_pair(c/2,-imag(a)*c);
}
/// position
Vector_3 get_r() const {
return real(b)/(2*real(a));
}
/// momentum
Vector_3 get_p() const {
return imag(b) - real(b)*(imag(a)/real(a));
}
///\en Transforms derivatives of a function whith respect to WP parameters
/// from internal into physical representation, i. e.:\n
/// from df/d{are,aim,b0re,b0im,b1re,b1im,b2re,b2im} (8 values accessed by input iterator d_it in the given order)\n
/// to df/d{x0,x1,x2}, df/d{p0,p1,p2}, df/dw, df/dpw
/// The supplied inputs (val) are modified by op: val=op(val,phys_der).
/// Use operation=eq_second for the supplied inputs to be replaced by new physical derivative values.
/// The inpput and output locations may coinside, an internal buffer is used for transformation.
template<template<class A> class operation, class d_it, class dfdx_it, class dfdp_it, class dfdw_it, class dfdpw_it>
void int2phys_der(d_it dfdi_,dfdx_it dfdx, dfdp_it dfdp, dfdw_it dfdw, dfdpw_it dfdpw, double h_p=h_plank) const {
operation<double> op;
double dfdi[8], dfn[8];// internal buffer
for(int i=0;i<8;i++)
dfdi[i]=*dfdi_++;
double w=get_width();
Vector_3 r=get_r();
double t=3/(2*w*w*w);
dfn[6]= -t*dfdi[0]-imag(a)*dfdi[1]/w; //dw
dfn[7]=-dfdi[1]/(2*w*h_p); // dpw
for(int i=0;i<3;i++){
dfn[i]= 2*real(a)*dfdi[2+2*i]+2*imag(a)*dfdi[2+2*i+1];
dfn[3+i]= dfdi[2+2*i+1]*(/*m_electron*/1./h_p) ; //*(h_plank/m_electron);
dfn[7]+=-(r[i]*dfdi[2+2*i+1]/w)/h_p;
dfn[6]+=-2*r[i]*(t*dfdi[2+2*i]+imag(a)*dfdi[2+2*i+1]/w);
}
int i=0;
for(int k=0;k<3;k++)
*dfdx++=op(*dfdx,dfn[i++]);
for(int k=0;k<3;k++)
*dfdp++=op(*dfdp,dfn[i++]);
*dfdw=op(*dfdw,dfn[i++]);
*dfdpw=op(*dfdpw,dfn[i++]);
}
///\en Compares the wave packet to another on a per component basis.
/// \return \retval 0 if all component differences are 0 within tolerance \a tol (EQUAL),
/// \retval -1 for LESS
/// \retval 2 for GREATER
int compare(const WavePacket &other, double tol=0.) const {
int res=compare_compl(a,other.a, tol);
if(res)
return res;
for(int i=0;i<3;i++){
res=compare_compl(b[i],other.b[i]);
if(res)
return res;
}
return 0;
}
};
#if 0
double w=wk.get_width();
Vector_3 r=wk.get_r();
double t=3/(2*w*w*w);
fe_w[ic1+k1]+= t*E_der[s1][indw1+8*k1]+imag(wk.a)*E_der[s1][indw1+8*k1+1]/w;
fe_pw[ic1+k1]+=E_der[s1][indw1+8*k1+1]/(2*w*h_plank);
for(int i=0;i<3;i++){
fe_x[ic1+k1][i]+= -2*real(wk.a)*E_der[s1][indw1+8*k1+2+2*i]-2*imag(wk.a)*E_der[s1][indw1+8*k1+2+2*i+1];
fe_p[ic1+k1][i]+= (-E_der[s1][indw1+8*k1+2+2*i+1])*(m_electron/h_plank); //*(h_plank/m_electron);
fe_pw[ic1+k1]+=(r[i]*E_der[s1][indw1+8*k1+2+2*i+1]/w)/h_plank;
fe_w[ic1+k1]+=2*r[i]*(t*E_der[s1][indw1+8*k1+2+2*i]+imag(wk.a)*E_der[s1][indw1+8*k1+2+2*i+1]/w);
}
#endif
/// constructs a conjugate packet
inline WavePacket conj(const WavePacket &wp){
return WavePacket(conj(wp.a),conj(wp.b),conj(wp.lz));
}
# endif // WAVEPACKET_H
Event Timeline
Log In to Comment