Commit b78c3e0a by Francois Gygi

compute dipole


git-svn-id: http://qboxcode.org/svn/qb/trunk@500 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 1d726db8
......@@ -3,7 +3,7 @@
// AtomSet.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: AtomSet.C,v 1.14 2007-01-27 23:44:59 fgygi Exp $
// $Id: AtomSet.C,v 1.15 2007-08-14 04:11:19 fgygi Exp $
#include "AtomSet.h"
#include "Species.h"
......@@ -383,6 +383,23 @@ void AtomSet::reset_vcm(void)
}
////////////////////////////////////////////////////////////////////////////////
D3vector AtomSet::dipole(void) const
{
D3vector sum;
for ( int is = 0; is < atom_list.size(); is++ )
{
double charge = species_list[is]->zval();
int i = 0;
for ( int ia = 0; ia < atom_list[is].size(); ia++ )
{
D3vector p = atom_list[is][ia]->position();
sum += charge * p;
}
}
return sum;
}
////////////////////////////////////////////////////////////////////////////////
void AtomSet::sync()
{
#if USE_MPI
......
......@@ -3,7 +3,7 @@
// AtomSet.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: AtomSet.h,v 1.15 2007-03-17 01:14:00 fgygi Exp $
// $Id: AtomSet.h,v 1.16 2007-08-14 04:11:19 fgygi Exp $
#ifndef ATOMSET_H
#define ATOMSET_H
......@@ -62,6 +62,7 @@ class AtomSet
void sync(void);
void reset_velocities(void);
D3vector vcm(void) const;
D3vector dipole(void) const;
void reset_vcm(void);
int size(void) const;
};
......
......@@ -3,7 +3,7 @@
// ComputeMLWFCmd.C:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: ComputeMLWFCmd.C,v 1.1 2007-08-13 21:26:27 fgygi Exp $
// $Id: ComputeMLWFCmd.C,v 1.2 2007-08-14 04:11:19 fgygi Exp $
#include "ComputeMLWFCmd.h"
#include<iostream>
......@@ -20,4 +20,33 @@ int ComputeMLWFCmd::action(int argc, char **argv)
mlwft->compute_transform();
mlwft->apply_transform(sd);
if ( ui->onpe0() )
{
cout << " <mlwf_set size=\"" << sd.nst() << "\">" << endl;
for ( int i = 0; i < sd.nst(); i++ )
{
D3vector ctr = mlwft->center(i);
double sp = mlwft->spread(i);
cout.setf(ios::fixed, ios::floatfield);
cout.setf(ios::right, ios::adjustfield);
cout << " <mlwf center=\"" << setprecision(6)
<< setw(12) << ctr.x
<< setw(12) << ctr.y
<< setw(12) << ctr.z
<< " \" spread=\" " << sp << " \"/>"
<< endl;
}
cout << " </mlwf_set>" << endl;
D3vector edipole = mlwft->dipole();
cout << " <electronic_dipole> " << edipole
<< " </electronic_dipole>" << endl;
D3vector idipole = s->atoms.dipole();
cout << " <ionic_dipole> " << idipole
<< " </ionic_dipole>" << endl;
cout << " <total_dipole> " << idipole + edipole
<< " </total_dipole>" << endl;
cout << " <total_dipole_length> " << length(idipole + edipole)
<< " </total_dipole_length>" << endl;
}
}
......@@ -3,7 +3,7 @@
// MLWFTransform.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: MLWFTransform.C,v 1.1 2007-08-13 21:26:27 fgygi Exp $
// $Id: MLWFTransform.C,v 1.2 2007-08-14 04:11:19 fgygi Exp $
#include <iostream>
#include <iomanip>
......@@ -174,37 +174,6 @@ void MLWFTransform::compute_transform(void)
a_[5]->ger(-1.0,cr,0,csz,0);
int nsweep = jade(maxsweep,tol,a_,*u_,adiag_);
// print Wannier center position
if ( ctxt_.onpe0() )
{
cout << " <mlwf_set size=\"" << adiag_[0].size() << "\">" << endl;
for ( int i = 0; i < adiag_[0].size(); i++ )
{
const double cx = adiag_[0][i];
const double sx = adiag_[1][i];
const double cy = adiag_[2][i];
const double sy = adiag_[3][i];
const double cz = adiag_[4][i];
const double sz = adiag_[5][i];
// Next lines: M_1_PI = 1.0/pi
const double t0 = 0.5 * M_1_PI * atan2(sx,cx);
const double t1 = 0.5 * M_1_PI * atan2(sy,cy);
const double t2 = 0.5 * M_1_PI * atan2(sz,cz);
const double x = t0*cell_.a(0).x + t1*cell_.a(1).x + t2*cell_.a(2).x;
const double y = t0*cell_.a(0).y + t1*cell_.a(1).y + t2*cell_.a(2).y;
const double z = t0*cell_.a(0).z + t1*cell_.a(1).z + t2*cell_.a(2).z;
D3vector ctr = center(i);
double sp = spread(i);
cout << " <mlwf center=\"" << setprecision(6)
<< setw(12) << ctr.x
<< setw(12) << ctr.y
<< setw(12) << ctr.z
<< "\" spread=\"" << sp << "\"/>"
<< endl;
}
cout << " </mlwf_set>" << endl;
}
}
////////////////////////////////////////////////////////////////////////////////
......@@ -299,11 +268,12 @@ double MLWFTransform::spread(void)
}
////////////////////////////////////////////////////////////////////////////////
D3vector MLWFTransform::center(void)
D3vector MLWFTransform::dipole(void)
{
// total electronic dipole
D3vector sum(0.0,0.0,0.0);
for ( int i = 0; i < sd_.nst(); i++ )
sum += center(i);
sum -= sd_.occ(i) * center(i);
return sum;
}
......
......@@ -3,7 +3,7 @@
// MLWFTransform.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: MLWFTransform.h,v 1.1 2007-08-13 21:26:27 fgygi Exp $
// $Id: MLWFTransform.h,v 1.2 2007-08-14 04:11:19 fgygi Exp $
#ifndef MLWFTRANSFORM_H
#define MLWFTRANSFORM_H
......@@ -42,7 +42,7 @@ class MLWFTransform
double spread(int i);
double spread(void);
D3vector center(int i);
D3vector center(void);
D3vector dipole(void);
MLWFTransform(const SlaterDet& sd);
~MLWFTransform(void);
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment