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