1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2008 The Regents of the University of California
4 //
5 // This file is part of Qbox
6 //
7 // Qbox is distributed under the terms of the GNU General Public License
8 // as published by the Free Software Foundation, either version 2 of
9 // the License, or (at your option) any later version.
10 // See the file COPYING in the root directory of this distribution
11 // or <http://www.gnu.org/licenses/>.
12 //
13 ////////////////////////////////////////////////////////////////////////////////
14 //
15 // PartialChargeCmd.cpp
16 //
17 ////////////////////////////////////////////////////////////////////////////////
18 
19 #include "PartialChargeCmd.h"
20 
21 #include <iostream>
22 #include <complex>
23 
24 #include "Context.h"
25 #include "Sample.h"
26 #include "Basis.h"
27 #include "Atom.h"
28 #include "ChargeDensity.h"
29 
30 using namespace std;
31 
32 ////////////////////////////////////////////////////////////////////////////////
action(int argc,char ** argv)33 int PartialChargeCmd::action(int argc, char **argv)
34 {
35   string usage("  Use: partial_charge [-spin {1|2}] name radius");
36 
37   // parse arguments
38   // plot [-spin {1|2}] name radius
39 
40   // ispin = 0: include both spins
41   // ispin = 1: include first spin
42   // ispin = 2: include second spin
43   int ispin = 0;
44   double radius = 0.0;
45   string atom_name;
46 
47   if ( !(argc==3 || argc==5) )
48   {
49     if ( ui->onpe0() )
50       cout << usage << endl;
51     return 1;
52   }
53 
54   const int nspin = s->wf.nspin();
55   // process arguments
56   int iarg = 1;
57   if ( !strcmp(argv[iarg],"-spin") )
58   {
59     if ( nspin != 2 )
60     {
61       if ( ui->onpe0() )
62         cout << "nspin = 1, cannot select spin" << endl;
63       return 1;
64     }
65     iarg++;
66     // process argument: ispin
67     if ( iarg==argc )
68     {
69       if ( ui->onpe0() )
70         cout << usage << endl;
71       return 1;
72     }
73     ispin = atoi(argv[iarg]);
74     if ( ispin < 1 || ispin > 2 )
75     {
76       if ( ui->onpe0() )
77         cout << " spin must be 1 or 2" <<  endl;
78       return 1;
79     }
80     iarg++;
81   }
82 
83   // argument must be the atom name followed by the radius
84   if ( iarg==argc )
85   {
86     if ( ui->onpe0() )
87       cout << usage << endl;
88     return 1;
89   }
90   atom_name = argv[iarg];
91   iarg++;
92   if ( iarg==argc )
93   {
94     if ( ui->onpe0() )
95       cout << usage << endl;
96     return 1;
97   }
98   radius = atof(argv[iarg]);
99 
100   // find atom and check validity of radius argument
101   Atom* a = s->atoms.findAtom(atom_name);
102   if ( a == 0 )
103   {
104     if ( ui->onpe0() )
105       cout << " PartialChargeCmd: atom " << atom_name << " not defined" << endl;
106     return 1;
107   }
108   if ( radius <= 0.0 )
109   {
110     if ( ui->onpe0() )
111       cout << " PartialChargeCmd: radius must be positive" << endl;
112     return 1;
113   }
114 
115   D3vector pos = a->position();
116   if ( ui->onpe0() )
117   {
118     cout << setprecision(5) << " Atom " << atom_name << " at " << pos << endl;
119     cout << " radius = " << radius << endl;
120   }
121 
122   const Context& ctxt = s->wf.sd_context();
123   ChargeDensity cd(s->wf);
124   Basis *vbasis = cd.vbasis();
125   cd.update_density();
126   const double omega = vbasis->cell().volume();
127 
128   double sum = 0.0;
129   const double fac = 4.0 * M_PI * radius * radius * radius / omega;
130 
131   // G=0 term
132   if ( nspin == 1 )
133   {
134     sum = fac * cd.rhog[0][0].real() / 3.0;
135   }
136   else
137   {
138     // nspin == 2
139     if ( ispin == 0 )
140       sum = fac * ( cd.rhog[0][0].real() + cd.rhog[1][0].real() ) / 3.0;
141     else if ( ispin == 1 )
142       sum = fac * cd.rhog[0][0].real() / 3.0;
143     else
144       sum = fac * cd.rhog[1][0].real() / 3.0;
145   }
146 
147   // Start sum after G=0 coeff if in first process row
148   int igstart = 0;
149   if ( ctxt.myrow() == 0 )
150   {
151     // skip G=0 coefficient
152     igstart = 1;
153   }
154 
155   const int ngloc = cd.rhog[0].size();
156   const double *gx = vbasis->gx_ptr(0);
157   const double *gy = vbasis->gx_ptr(1);
158   const double *gz = vbasis->gx_ptr(2);
159   const double *g = vbasis->g_ptr();
160   for ( int ig = igstart; ig < ngloc; ig++ )
161   {
162     // translate origin: compute exp(i G*pos )
163     double arg = pos.x*gx[ig] + pos.y*gy[ig] + pos.z*gz[ig];
164     double c = cos(arg);
165     double s = sin(arg);
166     // Re ( c(G) * (c + i*s )) = Re(c(G))*c - Im(c(G)*s
167 
168     complex<double> cg;
169     if ( nspin == 1 )
170     {
171       cg = cd.rhog[0][ig];
172     }
173     else
174     {
175       // nspin == 2
176       if ( ispin == 0 )
177         cg = cd.rhog[0][ig] + cd.rhog[1][ig];
178       else if ( ispin == 1 )
179         cg = cd.rhog[0][ig];
180       else
181         cg = cd.rhog[1][ig];
182     }
183 
184     // real part of coefficient of translated function
185     double ctrans_re = cg.real() * c - cg.imag() * s;
186     // product of norms: g * radius
187     double gr = g[ig]*radius;
188     // Bessel function j1(z) = sin(z)/z^2 - cos(z)/z
189     double j1gr = sin(gr)/(gr*gr) - cos(gr)/gr;
190     // factor 2 in next line: G and -G
191     sum += 2.0 * fac * ctrans_re * j1gr / gr;
192   }
193 
194   // accumulate sum across tasks
195   ctxt.dsum('c',1,1,&sum,1);
196 
197   if ( ui->onpe0() )
198   {
199     cout << " <partial_charge atom=\"" << atom_name
200          << "\" radius=\"" << radius;
201     if ( ispin > 0 )
202       cout << "\" spin=\"" << ispin;
203     cout << "\"> " << setprecision(6) << sum << " </partial_charge>" << endl;
204   }
205 
206   return 0;
207 }
208