Commit 165a3690 by Francois Gygi

Added partial_charge command


git-svn-id: http://qboxcode.org/svn/qb/trunk@1763 cba15fb0-1239-40c8-b417-11db7ca47a34
parent bb48df30
......@@ -49,7 +49,7 @@ OBJECTS=qb.o AtomSet.o Atom.o Species.o \
ExtForceSet.o ExtForce.o AtomicExtForce.o PairExtForce.o \
GlobalExtForce.o \
uuid_str.o sampling.o CGOptimizer.o LineMinimizer.o \
ElectricEnthalpy.o \
ElectricEnthalpy.o PartialChargeCmd.o \
$(PLTOBJECTS)
CXXFLAGS += -DTARGET='"$(TARGET)"'
$(EXEC): $(OBJECTS)
......
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// PartialChargeCmd.C:
//
////////////////////////////////////////////////////////////////////////////////
#include "PartialChargeCmd.h"
#include <iostream>
#include <complex>
#include "Context.h"
#include "Sample.h"
#include "Basis.h"
#include "Atom.h"
#include "ChargeDensity.h"
using namespace std;
////////////////////////////////////////////////////////////////////////////////
int PartialChargeCmd::action(int argc, char **argv)
{
string usage(" Use: partial_charge [-spin {1|2}] name radius");
// parse arguments
// plot [-spin {1|2}] name radius
// ispin = 0: include both spins
// ispin = 1: include first spin
// ispin = 2: include second spin
int ispin = 0;
double radius = 0.0;
string atom_name;
if ( !(argc==3 || argc==5) )
{
if ( ui->onpe0() )
cout << usage << endl;
return 1;
}
const int nspin = s->wf.nspin();
// process arguments
int iarg = 1;
if ( !strcmp(argv[iarg],"-spin") )
{
if ( nspin != 2 )
{
if ( ui->onpe0() )
cout << "nspin = 1, cannot select spin" << endl;
return 1;
}
// process argument: ispin
iarg++;
if ( iarg==argc )
{
if ( ui->onpe0() )
cout << usage << endl;
return 1;
}
ispin = atoi(argv[iarg]);
if ( ispin < 1 || ispin > 2 )
{
if ( ui->onpe0() )
cout << " spin must be 1 or 2" << endl;
return 1;
}
}
else
{
// argument must be the atom name followed by the radius
atom_name = argv[iarg];
iarg++;
if ( iarg==argc )
{
if ( ui->onpe0() )
cout << usage << endl;
return 1;
}
radius = atof(argv[iarg]);
}
// find atom and check validity of radius argument
Atom* a = s->atoms.findAtom(atom_name);
if ( a == 0 )
{
if ( ui->onpe0() )
cout << " PartialChargeCmd: atom " << atom_name << " not defined" << endl;
return 1;
}
if ( radius <= 0.0 )
{
if ( ui->onpe0() )
cout << " PartialChargeCmd: radius must be positive" << endl;
return 1;
}
D3vector pos = a->position();
if ( ui->onpe0() )
{
cout << setprecision(5) << " Atom " << atom_name << " at " << pos << endl;
cout << " radius = " << radius << endl;
}
const Context& ctxt = *s->wf.spincontext();
ChargeDensity cd(s->wf);
Basis *vbasis = cd.vbasis();
cd.update_density();
const double omega = vbasis->cell().volume();
double sum = 0.0;
// non-spin-polarized calculation
const double fac = 4.0 * M_PI * radius * radius * radius / omega;
// G=0 term
sum = fac * cd.rhog[0][0].real() / 3.0;
// Start sum after G=0 coeff if in first process row
int igstart = 0;
if ( ctxt.myrow() == 0 )
{
// skip G=0 coefficient
igstart = 1;
}
const int ngloc = cd.rhog[0].size();
const double *gx = vbasis->gx_ptr(0);
const double *gy = vbasis->gx_ptr(1);
const double *gz = vbasis->gx_ptr(2);
const double *g = vbasis->g_ptr();
for ( int ig = igstart; ig < ngloc; ig++ )
{
// translate origin: compute exp(i G*pos )
double arg = pos.x*gx[ig] + pos.y*gy[ig] + pos.z*gz[ig];
double c = cos(arg);
double s = sin(arg);
// Re ( c(G) * (c + i*s )) = Re(c(G))*c - Im(c(G)*s
complex<double> cg;
if ( nspin == 1 )
{
cg = cd.rhog[0][ig];
}
else
{
// nspin == 2
if ( ispin == 0 )
cg = cd.rhog[0][ig] + cd.rhog[1][ig];
else if ( ispin == 1 )
cg = cd.rhog[0][ig];
else
cg = cd.rhog[1][ig];
}
// real part of coefficient of translated function
double ctrans_re = cg.real() * c - cg.imag() * s;
// product of norms: g * radius
double gr = g[ig]*radius;
// Bessel function j1(z) = sin(z)/z^2 - cos(z)/z
double j1gr = sin(gr)/(gr*gr) - cos(gr)/gr;
// factor 2 in next line: G and -G
sum += 2.0 * fac * ctrans_re * j1gr / gr;
}
// accumulate sum across tasks
ctxt.dsum('c',1,1,&sum,1);
if ( ui->onpe0() )
cout << " <partial_charge atom=\"" << atom_name
<< "\" radius=\"" << radius
<< "\"> " << setprecision(6) << sum << " </partial_charge>" << endl;
return 0;
}
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2015 The Regents of the University of California
//
// This file is part of Qbox
//
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// PartialChargeCmd.h:
//
////////////////////////////////////////////////////////////////////////////////
#ifndef PARTIALCHARGECMD_H
#define PARTIALCHARGECMD_H
#include <iostream>
#include <cstdlib>
#include <string>
using namespace std;
#include "UserInterface.h"
#include "Sample.h"
class PartialChargeCmd : public Cmd
{
public:
Sample *s;
PartialChargeCmd(Sample *sample) : s(sample) {};
const char *name(void) const { return "partial_charge"; }
const char *help_msg(void) const
{
return
"\n partial_charge\\n\n"
" syntax: partial_charge [-spin {1|2}] name radius\n"
" The partial_charge command computes the amount of charge\n"
" density contained in a sphere centered on an atom.\n"
" When using the -spin option, the charge of the given spin\n"
" is computed.\n\n";
}
int action(int argc, char **argv);
};
#endif
......@@ -59,6 +59,7 @@ using namespace std;
#include "ListSpeciesCmd.h"
#include "LoadCmd.h"
#include "MoveCmd.h"
#include "PartialChargeCmd.h"
#include "PlotCmd.h"
#include "PrintCmd.h"
#include "QuitCmd.h"
......@@ -275,6 +276,7 @@ int main(int argc, char **argv, char **envp)
ui.addCmd(new ListSpeciesCmd(s));
ui.addCmd(new LoadCmd(s));
ui.addCmd(new MoveCmd(s));
ui.addCmd(new PartialChargeCmd(s));
ui.addCmd(new PlotCmd(s));
ui.addCmd(new PrintCmd(s));
ui.addCmd(new QuitCmd(s));
......
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