245 lines
7.4 KiB
C++
245 lines
7.4 KiB
C++
# 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 {
|
|
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 {
|
|
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 with 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 input 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
|